Concepts, Capabilities, and Limitations of Global Models: A Review

For researchers wishing to generate an understanding of complex plasma systems, global models often present an attractive first step, mainly due to their ease of development and use. These volume averaged models are able to give descriptions of plasmas with complex chemical kinetics, and without the computationally intensive numerical methods required for spatially resolved models. This paper gives a tutorial on global modeling, including development and techniques, and provides a discussion on the issues and pitfalls that researchers should be aware of. Further discussion is provided in the form of two reviews on methods of extending global modeling techniques to encompass variations in either time or space.


Introduction
Low temperature plasmas have a large number of scientific and commercial applications, so an understanding of their properties is important for optimization of plasma based processes and technologies. Modeling a plasma discharge and its behavior can help greatly in gaining this comprehension, but depending on the precise properties of the system of interest, the appropriate method to simulate the plasma may vary. Global modeling represents a numerical method of describing plasma discharges, based on fluid equations, that neglects spatial derivatives in order to enhance computational efficiency. They are able to quickly predict spatially averaged plasma parameters such as densities or temperatures for systems that would otherwise be difficult to simulate, and relationships between key parameters can be explored across a broad range of system properties.
Global models are based on two types of equation: particle balance equations, written for each included species, and power balance equations, which are primarily used for electrons, but can be included for other species. Finding solutions to the resulting set of equations requires a range of information about the system of interest, including physical characteristics, power coupling method, and a list of species included and the reactions between them. The equations used to build a global model can appear simple, but the resulting system that they describe can behave in ways that are unusual and often unintuitive. Global models have therefore been widely used for analyzing the chemistry and identifying the main reactions in low temperature plasmas, as they allow for complex chemical reaction schemes, with a large number of species and reactions to be studied, usually without the associated long computation times of models inclusive of spatial resolution.
Due to the widespread use and high value to researchers of global models, a number of well known codes are available that will solve global models for a user defined set of conditions. One such program is the GlobalKin code, [1,2] which solves a global model using given rate coefficients, and also contains a module for finding the electron energy distribution function (EEDF) using a common two term Boltzmann approximation. A commercial application, Quantemol-P, is built around GlobalKin to provide a graphical user interface for increased usability. [3] Another powerful tool is ZDPlasKin, which is able to solve for the time evolution of species densities in a system with a user defined reaction scheme. [4] As with GlobalKin, it includes a solver for the EEDF, and is able to calculate gas temperatures. The outputs of both ZDPlasKin and GlobalKin are compatible with the PumpKin pathway analysis software. [5] The plasma modelling toolkit PLASIMO also provides users with a framework to create a global model. [6] Despite these tools existing, and usually being easy to use, it is still sometimes desirable to create a global model of a particular system. This may be for reasons of furthering the understanding of the underlying physics, or because the system to be described is unusual in some way.
In the first part of this tutorial review, an introduction on global models is given, along with a tutorial on the development of a simple model. Emphasis is placed on the key issues that need to be considered when constructing such a model, and generality is maintained so that the framework presented can be easily built upon. Following a description of the necessary considerations, the tutorial documents the formulation of a simple argon model, from equation derivation through to numerical solution and example results. A discussion is then given on the importance of reaction schemes, including the need to source rate coefficients and other data for a potentially large number of reactions, depending on the species being considered. The reactions that are necessary may include atomic or molecular excited states, for which reliable data are often difficult to come by, despite their importance in certain plasma systems. Finally, a word of caution is given as to some of the limitations of global models, including the implications of poor source data for reaction rate coefficients. As there can be a relatively small number of reactions that dominate the behavior of a system, small changes in these coefficients for specific reactions can lead to drastic differences in the resulting behavior. Due in part to these issues, and despite their widespread use, there are a number of caveats that must be considered while interpreting the results that global models provide.
The second part of this work looks beyond the common assumptions of temporal and spatial homogeneity normally used in the basic global model approach, and presents reviews of two classes of significant extensions. A large number of low temperature plasmas are driven by oscillatory electric fields, in which the power absorption is inherently time dependent. To include such effects in a model typically also requires some degree of spatial resolution. The first review section explores some of the techniques that have been used to include time dependent power deposition, while working around the need for spatial resolution and the associated computational cost.
The second review section investigates some applications of global models to systems with known high amounts of spatial non-uniformity, and the different ways developed by researchers to tackle this problem.

Framework of a Simple Global Model
In order to build a model describing a plasma, a number of decisions must be made about what considerations to include, and what to exclude, while keeping in mind the limitations that follow from any assumptions made. One of the first considerations to make is the physical properties of the system to be modeled, as supposedly simple parameters such as the geometry of the plasma can have impacts on the behavior and results of the model; the ratio of vessel internal surface area to plasma volume is particularly important for global models, as will be discussed. For investigations of general properties, a simple planar or cylindrical geometry is usually sufficient. [7][8][9] However, if a particular device is to be modeled for the purposes of comparison with experiment, then more complicated geometries are often required, as well as considerations of gas flow and mechanical pumping systems. [10][11][12][13] Alongside geometrical considerations, the method by which power is deposited in the plasma must also be determined. It is important to know if the plasma is powered through a DC or AC drive, and if this is radio frequency (RF), or another waveform, and whether this is capacitively or inductively coupled, [14][15][16] or both. [17] Alternatives such a helicon drive [18] or electron-cyclotronresonance [19,20] are also possible.
Considerations must also be made about the plasma itself. The gas or gases being used to create the plasma are often dictated by the investigation being performed, but the precise properties of each gas, as well as those of any species created by the plasma, need to be taken into account. In particular, the energy distribution functions (EDFs) of each species are important, and can have significant effects on the model. Neutral particles are often taken to have a Maxwellian EDF with a particular temperature, although it is possible to model this temperature self consistently. [1,12,[21][22][23] Ionic species are not necessarily in thermal equilibrium, and so may deviate from a Maxwellian EDF; despite this they are often modeled as such due to the complications of including such a deviation. [11,16,17,24] Electrons can be described by a Maxwellian EDF only under certain specific conditions, and so it is often more appropriate to make a different assumption.
One of the more difficult considerations is the set of reactions, interactions and chemical products, including excited or metastable states, that the model includes. The interaction of the plasma with any walls must also be considered. Even models of apparently simple gases can contain very large numbers of reactions. Creating a suitable reaction scheme that encompasses all of the major effects, without becoming intractable, is a complex task, made more difficult by the requirement to find data for collisional cross sections or rate coefficients. This topic is complex, and will be discussed in more detail in Section 2.3.

Analytical Derivation
Once the necessary components of the model have been decided upon, equations must be obtained to describe the system. This section and the next detail the steps required to build a global model of a simple argon plasma, consisting only of electrons and positive ions in a neutral background of constant density and temperature, from generic equations through to generating results. The system is kept as general as possible, so that the reader may more easily apply the model to more complex systems.
As global models are about macroscopic descriptions of the plasma, they can be considered as a subtype of fluid models. Briefly, the fluid equations detail the evolution of the macroscopic properties of a system in both space and time, and can be used to describe the conservation of mass, momentum, energy, energy flux, and so on. For real systems, which are non-conservative, these equations must also consider gains and losses of their respective quantities, be this through particle interactions, losses to the wall, or other external effects.
The fluid equations used for a typical global model stem from the moments of the Boltzmann equation, the derivation of which is covered in many plasma physics theory texts. [25,26] The zeroth moment deals with the conservation of mass, or particles, and can be written @n @t þr r rÁ nu ðÞ ¼ dn dt ð1Þ The first moment relates the conservation of momentum to the forces acting on species, but is not included here as the removal of spatial considerations means that momentum is not resolved. The second moment, given as treats the conservation of energy. Both thermal and kinetic energies are included, as they are linked at a microscopic level. The symbols used in these and the following expressions can be found in Table 1.
In Equation (1), the two terms on the left hand side represent the change in particle density in time, and the spatial change of the particle flux, respectively. The right  (2), details the changes in volumetric energy density from internal sources on the left hand side, and external sources on the right. The internal sources comprise of changes in temperature in both space and time in the first two terms. The third and fourth terms track the change in internal energy density due to velocity divergence, and the final term on the left hand side represents changes in the transport of heat. Similarly to Equation (1), the term on the right hand side of Equation (2) contains changes to the total energy due to external sources.
In a basic global model, the assumption is made that the plasma is homogeneous, and so all spatial derivatives are zero. This leaves Equation (1) and (2) consisting of only temporal derivatives and effects due to external sources. Some debate exists as to the precise form of the external source terms, but they need to consider effects due to elastic and inelastic collisions with other species within the plasma, including particle and energy changes, as well as interactions with any walls.
In order to develop these terms, one must start to consider that each species in the system to be modeled requires its own set of equations, and that the interactions between species is contained within the collision terms on the right hand side of the equations. With this in mind, equations are presented hereon with subscripts e, i,o rg denoting the species being considered as electrons, ions, or background neutrals, respectively. A subscript a denotes any species.
The term for external changes in particle density needs to include particles created and lost through reactions, as well as those flowing into or out of the volume of interest. Flows into or out of the model can be physical flows, such as from a gas inlet or to a pumping system, or through interactions with the wall of the system, such as surface recombination or secondary electron emission. In the example model to be created, it is assumed that the neutral gas background is constant, so no gas flows or pumping are included, and that any excited or charged particles incident on the wall of the system are lost. Therefore, the term describing external changes to density dn=dt takes the form of the right hand side of In this expression, the first two terms on the right hand side contain reactions and (de)excitations that lead to volume gains or losses, respectively, and are of the form n a n b K R . The third term describes losses to the wall, and considers a flux of particles of type a, to be determined, onto a wall which has total surface area A.
The term for external changes to species energy needs to consider a variety of effects. Energy is gained by charged species through the application of external fields, be these radio frequency or otherwise, and a term for this must be considered. The precise form of this term depends on the heating mechanism being employed. There is also a transfer of energy between species, as they may have different temperatures, and so energetic species will lose energy to cooler ones through elastic collisions. Energy is gained and lost through chemical processes, due to both the gain and loss of particles and from the energy required or released by the reactions. The loss of particles to the wall must also be considered, as each escaping particle takes energy with it.
In this example, it is assumed that neutrals and ions have a constant temperature of 300 K, and so Equation (2) must only be derived for electrons. The transfer of energy due to elastic collisions is assumed to happen only between electrons and neutral particles, as these are the dominant elastic collision partners at low ionization fractions. In plasmas with a high degree of ionization, electron-electron collisions are also important for energy transfer. [27][28][29] The amount of energy transfer in electron-neutral collisions can be approximated through a hard sphere model. [25,26] The energy changes due to inelastic collisions can be found through a simple summation of the reactions involving electrons, and the gains or losses of electron energy as a result of each.
In this simple example, the power deposition is approximated by a constant value representing the time averaged power absorbed by the electrons, S abs . This is in general not a valid assumption, but is done here for ease of understanding and to allow for a simple solution scheme, as presented in Section 2.2.
Particularly in capacitively coupled discharges, there may be significant modulations of the deposited power in both space and time. The presence of non-ohmic heating mechanisms can also exacerbate the discrepancies between the true power deposition and a constant time averaged value. [30,31] By failing to account for these changes in deposited power, it is possible to miss temporal modulations of electron temperature. This can result in errors even in equilibrium values obtained from global models, as reaction rate coefficients generally have nonlinear dependencies on the electron energy. Discrepancies can also arise in results through the misunderstanding of what S abs represents, and confusion around the relationship between it and the power displayed for example on an RF amplifier, which may occur due to significant power absorption by ions in the sheath under certain conditions. [32] The considerations for energy gains and losses due to changing numbers of particles is similar to that given in Equation (3), however, there are additional energy losses associated with particles crossing the sheath. With the assumption of Maxwellian energy distributions for all species, each electron lost through the sheath takes 2k B T e of energy with it. [33] If it is assumed that ions enter the sheath with the Bohm velocity, then each ion takes 1 2 k B T e of kinetic energy, as well as being accelerated by the sheath voltage before being removed from the system. [33] The Bohm velocity is an analytical result that specifies a minimum ion velocity for the formation of a sheath in a simple electropositive plasma as being k B T e =m i ðÞ 1=2 . Although it is derived using the assumptions of zero ion temperature, [33] electrons in Boltzmann equilibrium with the plasma potential, and a collisionless sheath, it sees widespread use across a large range of plasma parameters. Extensions exist for multiple ion species, [34] including negative ions, [35] but the classical result will be used in this model. Although the Bohm velocity is a widely known result in low temperature plasma physics, the applicability of it to any system that deviates from the strict conditions from which it is derived is a topic of ongoing contention, [36][37][38][39][40][41] and researchers should use caution when using a value for the Bohm velocity for anything but the most basic of plasmas.
As it is assumed that the plasma is driven in this case by RF, then the mean value for the sheath voltage, V s must be used, as discussed later. By combining each of the above considerations with Equation (2) and (3), one obtains an expression for the rate of change of electron temperature 3 2 k B n e @T e @t ¼ S abs À 3 m e m g K eg n e n g k B T e À T g ÀÁ which can then be rearranged to give @T e e @t ¼ 2 3 S abs en e À 2 m e m g K eg n g T e e À T e This expression also uses the knowledge that the plasma must conserve current, so the wall fluxes of electrons and ions must be equal. To make computation easier, the change has also been made to defining temperatures and energies in eV, as both electron temperature and the energies of collisions are on the order of 1-10 in these units.
In Equation (3) and (4), there are a number of quantities that are yet to be defined. Most notably are the reactions, their reaction rate coefficients, and the energies gained or lost due to each reaction. These are given in Table 2 for a simplified argon example. The excitation reaction is considered only as an electron energy loss mechanism; the argon metastable states are not explicitly included. It is also assumed that the model is to be used in a pressure range where the effective rate coefficients for three-body collisions are negligibly small.
Also as yet undefined are the mean sheath voltage V s and the electron/ion wall flux G wall . For a plasma driven by RF in capacitive mode, it is possible to estimate a mean sheath voltage as a function of the ohmic power deposition and the electron-neutral elastic collision frequency where the assumption is made that the power is only deposited to the electrons through the ohmic channel. [33] For the flux of species to the wall, the ions are again assumed to enter the sheath with the Bohm velocity. With the assumption that no particles are created or destroyed inside the sheath, the flux of ions, and electrons leaving the plasma can be found using the Bohm velocity and the density of ions at the sheath edge. Assuming a collisional, quasineutral plasma contained between two electrodes with a small sheath, this density can be approximated as [24] where the n e given is the central ion/electron density. From this, the wall flux can be equated to which provides the last unknown of the system. The ratio of ion densities at the center and sheath edge, n s =n 0 , is commonly referred to as h l . The value given in Equation (6) is a classical result for a simple, single ion, electropositive plasma. However, as is the case for the Bohm velocity, the value is subject to change depending on the plasma conditions, and care must be taken to ensure that the value used is appropriate for the plasma being studied. Multiple values exist for different plasmas, and attempt to account for effects including electronegativity, multiply charged ions, or changing collisionality regimes. [43][44][45][46][47] The concept of the h l factor for various systems is explored more in Section 3.2.1.
These definitions leave the reaction scheme, input power, system pressure, and reactor geometry as inputs to the model. In this example, the plasma is assumed to be contained between two infinite planar electrodes separated by a distance l, so that the surface area to volume ratio A=V % 2=l. The electron and ion densities are considered equal, as discussed above, and so only a single density equation is required. Combining all of these aspects leads to a set of two differential equations that must be solved simultaneously, written as @n e @t ¼ K I n e n g À 2p n e eT e e m i n g K ig l 2 ð8Þ @T e e @t ¼ 2 3 S abs en e À 2 m e m g K eg n g T e e À T e g À 2 3 15:76 Â K I þ 12:14 Â K X ðÞ n g À T e e K I n g

Numerical Solution
Exact analytical solutions to global model equations are, in general, not possible, and so they must be found through numerical methods. A number of options are available, depending on the precise form of the equations obtained. For those with time dependent properties, such as models including gas flow [22] or the time dependence of power deposition, [45,48] it is most appropriate to perform a direct numerical integration of the equations over a given time period. However, for models with no time dependent inputs, it is usual to use a numerical root finding method to find equilibrium values, where all time derivatives are zero. If one is intending to use root finding, then it is prudent to rework the equations into a form that simplifies the numerical work needed. In the case of Equation (8) and (9) where n e and T e e are the only unknowns, it is possible to split the problem into two steps. Dividing Equation (8) by n e and setting the time derivative to zero yields an equation independent of n e . Although it is not possible to rearrange to solve for T e e due to the form of K I , it is possible to express the pressure-length product as a function of T e e , pl This also shows how the pressure and the length scale of the discharge are intrinsically linked, as is seen in classic analytical results, [24,33] and that the pressure-length product can be used as system property, as opposed to the length and pressure independently.
Using this it is possible to use numerical minimization to find the value of T e e that gives pl from the inputs p, l, and the assumption of T g ¼ 300 K. Once the value of T e e has been found, it is possible to rearrange Equation (9) to find n e analytically using the expression n e ¼ 2 3 e À1 S abs 2 m e m g K eg n g T e e À T e g þ 2 3 15:76K I þ 12:14K X ðÞ n g þT e e K I n g þ Using this two stage process, it is then possible to find the electron temperature and density for any given combination of system length, pressure, and input power density. This takes a very short time on a regular desktop or laptop computer, and can give thousands of solutions per second. Although this model is highly simplified, this property of rapid computation is common to most global models. Using them it is trivial to calculate trends and behaviors over a wide range of parameters, as is demonstrated in Figure 1.
As global models are based on the fluid equations, different effects can be isolated and investigated, such as the relative importance of elastic and inelastic collisions. This can be useful to identify the causes of unexpected or unintuitive behaviors. For example, in Figure 1a the electron density can be seen to be non-monotonic with pressure-length product. Analysis of the different energy loss mechanisms given in Equation (9) shows that at low pressures, electron energy is lost primarily to the wall and to Table 2. List of reactions for a simplified argon chemistry. [24,33,42] Expressions for K I , K X , and K eg are fits to data assuming a Maxwellian electron energy distribution function, from Gudmundsson. [42]  hi - q inelastic collisions. At high pressures, gas heating through elastic collisions is the leading cause of electron energy loss. This can be seen in Figure 2, which shows that between the two regimes, there is a mode transition, leading to the peak in electron density at intermediate pressures. This is also seen in Figure 3 across a wide range of input powers. Additionally, in Figure 3 it can be seen that the position of the peak density changes with power input. This is due to the dependence of the sheath voltage, found in the last term of Equation (11), on the input power, as given in Equation (5).
The example results presented in this section give an indication of how quickly a model can be developed and used. However, they also show how the models can easily be unwittingly used outside of the regime in which they are valid. For example, when the mean free path of particles becomes significant compared to the system length, one can no longer assume that the dominant effects are those dependent on collisions, such as ohmic heating. In this model the electron thermal mean free path is roughly 100% of the discharge length at a pressure-length product of 0.4 Pa Á m, and increases rapidly as the pressure drops, and so the results in the region below the dashed line in Figure 3 are mostly likely inaccurate.
Another important reason to be critical of these results is the limited number of species and reactions that are accounted for. The example model presented includes a metastable state of argon only as an energy loss term;  [49] Input power density was held at 1 kW Á m À3 .  resolving other effects would require including it as a separate species, and so increase the model complexity. As presented, the metastable energy loss channel affects the electron density only by one part in 10 5 , and so a naïv e interpretation would suggest that the effect of metastables is negligible. In fact, even in ''simple'' mono-atomic gases such as argon or helium, the importance of excited metastable states is widely known. [50][51][52] This is succinctly shown by global models that self-consistently include excited states of argon, which can contain up to 20 reactions. [53,54] It was found that the multi-step ionization process is the dominant ion creation pathway for a wide range of plasma parameters, accounting for nearly 70% of ionization under certain conditions. This is because each step requires a smaller amount of energy than ionization from the ground state, and so is more likely to occur due to the higher populations of lower energy electrons.
The introduction of molecular species into such a system again increases the complexity many times over. Considerations must be made for chemical products of the supplied gas mixture, as well as negative ions, excited states, and molecular dissociation.

Chemical Reaction Schemes
Global models containing detailed chemical kinetics are used by several communities to describe a variety of systems. Typical examples where complex chemistries are important are in atmospheric, astrochemical, and combustion modeling. Many species of interest in these fields overlap with those that are important in plasma modeling as well. There exist a number of databases from these communities, [55][56][57][58][59] providing a base for building complex chemical models. Typical species where significant overlap of interest exists between the discussed communities are those containing nitrogen (N), oxygen (O), hydrogen (H), or carbon (C). These species are important in many plasma assisted processes, such as the production of reactive species for biomedical applications, plasma etching/ deposition processes, or plasma assisted combustion. However, in picking the reactions relevant for the system to be described, one has to be careful to cover the important reaction mechanisms for all species of interest. As well as needing to ensure that all important reactions are considered, one must be aware of the dangers of inaccurate reaction rate coefficients. It is often possible to find multiple sources for a reaction rate coefficient or cross section, and it can be difficult to trace primary sources of data. This can lead to issues relating to unknown uncertainties in the values obtained, and if the reaction is one that plays a significant role, then large differences in reported behavior can occur. It is obvious that the more complex the system, the more difficult this task can become.
A topical example of this is atmospheric pressure plasmas (APPs) operated in noble gases, with different molecular admixtures included specifically to generate chemically active species. These plasmas are interesting for various applications, such as surface modification, etching, and also biomedical applications. [60][61][62][63][64][65][66] One of the more studied cases is an Ar or He discharge with molecular oxygen (O 2 ) admixture, [67][68][69][70] which is generally modeled using not more than 25 species and 373 reactions to describe the plasma chemistry.
Very recently, reaction schemes for water containing APPs have been established, such as for a noble gas containing a small amount of humidity (H 2 O). The molecular nature of H 2 O results in a large reaction scheme of 46 species and 577 reactions, [71] and with addition of O 2 , [72] this increases to 55 species and 855 chemical reactions. In addition, the highly complex case of humidity containing air plasmas or plasmas that contain impurities have also been investigated, in some cases with over 1800 reactions, some from the addition of carbon containing species. [1,48,[73][74][75] It becomes clear that, intuitively, the more species are part of the initial gas mixture, the higher the number of species and chemical reactions that must be considered. This is illustrated for APPs in Figure 4, where the number of reactions is plotted as a function of the number of species included in the model for a selection of published models. In general, similar scalings can be expected for other systems as the number of species considered is increased.
In addition to selecting the species and reactions included, rate coefficients for each reaction need to be found. This is not always a simple task, as these coefficients can depend on the gas temperature for heavy particle collisions, electron temperature for the electron collisions, and often the gas pressure for certain types of reaction. Furthermore, any analytical expression for a reaction rate coefficient can only ever be an approximation of the complex interactions taking place on a fundamental level, which often require quantum calculations to be described accurately. As a result of this, it is important to consider whether or not reaction rate coefficients obtained in the literature are suitable for the parameter range to be investigated by the model.
Ideally, coefficients for reactions between heavy particles should include dependencies on gas temperature and pressure, where these behaviors are known, as some can be highly variable with these parameters. For example, the density of ozone in atmospheric pressure plasmas is known to be heavily dependent on the gas temperature, due to the variation of the different rates of creation and destruction. [76,80] The formation of H 2 O 2 , usually depicted as is an example of a reaction that depends on both gas temperature and pressure in different ways. Although usually referred to as a three-body reaction, fundamentally this process proceeds as a sequence of two successive twobody reactions. These are given as which show the formation of an intermediate excited molecule. This excited and unstable molecule can subsequently decompose back into the original particles if the excess energy is not removed through collision with a third body. Therefore, the overall process is often treated as a three body reaction, to account for the need for a third body to remove the excess energy. At low pressures, the effective rate coefficient for Equation (12) will increase with increasing pressure, due to the greater stabilization rate of the intermediate excited state. This can be seen in the behavior of the low pressure reaction rate coefficient, k 0 , in Figure 5.
As the pressure increases, and the probability of collisional de-excitation of the intermediate excited species becomes greater, the dependence of the effective rate coefficient on pressure will become less pronounced. In the high pressure limit, the background gas density is high enough that practically all excited intermediates will be stabilized before decaying. Therefore, in this situation, the reaction becomes effectively two-body in nature, and can be described by a pressure independent rate coefficient, given as k 1 in Figure 5.
Although k 0 and k 1 are experimentally well characterized for many reactions, [55] the effective rate coefficients in the transition region are typically approximated using the Lindemann-Hinshelwood model, shown as the solid lines in Figure 5. [81] This method uses the ratio of the rate coefficients for the two known limiting cases to calculate a ''falloff'' curve that is able to bridge the two regimes. A method to account for possible broadening of these falloff curves due to multiple possible excited intermediates was proposed by Troe. [82] This broadening factor, F, is strictly dependent on temperature and pressure; however, the approximation of a constant central value F c is appropriate for most atmospheric pressure systems. [82] A dependence on temperature or pressure is of course not only possible for neutral-neutral interactions, but also for many reactions involving ionic species too. Therefore, one has to be careful to make sure that chosen rate coefficients match to the system being described, especially when taking rate coefficients from secondary publications, which may have created effective rate coefficients particular to their specific conditions. Even though, a large number of rate coefficients are available in the literature, many rate coefficients are still not known. This is particularly true in the case of reactions involving ions and excited states. In those cases, it is often possible to calculate rate coefficients through a variety of approximations.

The Role of Excited States
In general, excited states play an important role in many systems investigated by global models. [11,[83][84][85][86][87][88][89][90][91][92] Generally, Figure 5. Dependence of the calculated rate coefficient for the reaction in Equation (12) on pressure and temperature using N 2 as a third body, calculated from data from the IUPAC database. [55] Low and high pressure limits k 0 and k 1 are indicated as straight lines. [55] A. Hurlbatt et al. three different types of excitation exist: electronic, rotational, and vibrational. All of these can have a significant influence on the plasma dynamics. For example, electronic excitation can lead to the formation of metastable states, which transport energy and can help sustain plasmas. In atmospheric pressure plasmas, Penning ionization via collisions with metastable species is an important ionization process. [67] Even in low pressure systems, the presence of metastables can play an important role, for example in oxygen process plasmas, collisions with the O 2 1 D g ÀÁ excited state can be the dominant destruction mechanism for the O À negative ion, under certain conditions. [93][94][95][96] Another form of excitation is into vibrational states in plasmas containing molecules. These states are important in many systems as a crucial loss channel for the electron energy, as well as an important intermediate step in some reactions. Vibrationally excited states can lead to the easier dissociation of molecules such as CO 2 and CO, [86,88,92] N 2 , [97] and O 2 . [88][89][90][91] Additionally, vibrationally excited molecules often have larger dissociative attachment cross sections, and as such their presence can enhance the production of negative ions. An important example is the formation of H À ions in plasmas containing H 2 , such as those used in negative ion sources for fusion applications. [11,[98][99][100][101][102] A comparison between the results of global models, with and without vibrational states, has shown that there is better agreement between simulations and experiments when vibrational states are included, even at relatively low pressures. [85] Therefore, an accurate representation of the vibrational kinetics is often required to accurately predict the behavior of plasmas containing molecular species and improve industrial processes that depend upon them. [99,100] However, electron impact excitation cross sections and heavy particle rate coefficients for excitation and deexcitation of these excited states are not always experimentally measured, or available from full quantum mechanical calculations. It is possible to calculate rate coefficients for reactions where no data exists, through the use of either quantum mechanical approximations, or from existing measured cross sections. If the cross section for excitation from the molecular ground state to the first vibrational level is known, then it is possible to find estimates for the higher excited states through the use of an approximation on the relationship between the levels, often termed the Fridman approximation. [86,87] Another possibility is to obtain electron impact cross sections through computationally intensive techniques such as the R-matrix method, which has recently been employed for the calculation of electron impact cross sections involving vibrational states for several different molecules. [88,90,103,104]

Reduction of Reaction Schemes
It has been discussed that some plasma systems require large chemical reaction schemes to be described. Although a number of databases for reaction rates are available to construct extensive schemes, many rates are not known, or only poorly known, as previously mentioned. When constructing a model, it is tempting to include all reactions possible to describe a physical system. However, it can be difficult to realistically assess the accuracy of each reaction rate coefficient in schemes comprised of 1000's of reactions, and as such the uncertainty in the solution of the model is hard to quantify. This has been emphasized in recent investigations by Turner [69,105] and is demonstrated concisely in Figure 6. This shows a number of possible solutions to the time evolution of the density of helium metastables in a helium APP containing 0.1% oxygen. The reaction rate coefficients used to obtain each solution were generated randomly, but lie within the uncertainties attributed to each.
The recent work of Turner has shown that significant reduction of reaction schemes is possible without greatly influencing the final solution of the model. [105] In that work it was found that a He-O 2 reaction scheme could have the number of reactions reduced by 85% and still provide the same solution, within the error bars, when compared to the ''full'' reaction scheme. It was also found that the remaining rate coefficients contributed most strongly to the uncertainty in the final solution of the model. The method presented by Turner [105] opens up the possibility for the reduction of other reaction schemes, and could lead to an improvement in efficiency for many plasma models.
In addition to the Monte-Carlo approach of Turner, [69] other methods for reaction scheme reduction have been proposed. An algorithm was developed by Lehmann [106] to Figure 6. Possible solutions to the trajectories of the helium metastables in an APP, calculated with different values for reaction rate coefficients taken from within the uncertainty of each. ß IOP Publishing. Reproduced with permission. All rights reserved. [69] identify the significant pathways in a reaction scheme by performing an analytical investigation of the possible reactions and their relative rate coefficients. This algorithm has since been used to create the software tool PumpKin, [5] which is able to perform this analysis for user defined reaction schemes, and can identify which processes are important over different time scales.
These tools and other approaches can also help to reduce full reactions schemes to a set of the most significant reactions by disregarding those reactions which fall below a certain threshold. [68,71,75,107] This can then subsequently be used in more complex models where an overly large number of reactions would be of detriment to the execution time or some other aspect of the model. [14,108] These reduced reaction schemes should however be used with caution, as the dominant pathways can depend on the precise conditions under which the plasma is operated.
Such studies are also potentially useful for the reduction of uncertainty in the solutions of such models, as the reactions identified as contributing most to the uncertainty can be in principle be measured again, with the aim of reducing the uncertainty in the corresponding reaction rate coefficient. Experimental studies to reduce the uncertainty of key reaction rate coefficients are arguably more worthwhile endeavors than those to measure many reaction rates which have little consequence on the results of models.

Limitations and Pitfalls
Although global models are used extensively in low temperature plasma research in a very successful manner, it is necessary to emphasize that they do have a range of inherent limitations. In the preceding discussion, several limitations regarding reaction rate coefficients have been highlighted. Similar points could be raised regarding electron impact cross sections. These issues surrounding fundamental input data are common to all plasma modeling approaches; however, they take on special significance in global models where the primary aim is often to understand systems comprising of complex chemistries. In this context, inaccurate or unknown rate coefficients or cross sections can have a very significant effect on the conclusions of works utilizing global models. In particular, it is important to have a critical view of the results of models where a significant number of reaction rates or cross sections are not experimentally measured, but are instead estimated.
In addition to fundamental data concerns, global models have a number of specific limitations based on the assumptions used in their development. Some of the most limiting assumptions are those based on analytical or semiempirical expressions such as that for the Bohm velocity or the h l factor, both discussed in Section 2.1. Such expressions generally have a fixed range of conditions under which they provide acceptable solutions, and can be limited by several factors, such as the degree of electronegativity in the case of the Bohm velocity, for example.
Further, limitations may be imposed based on assumptions regarding the EEDF, which cannot be obtained selfconsistently in global models, but can have a significant effect on the solution they produce. [109][110][111][112] A realistic representation of the EEDF can be obtained under highly collisional conditions where the electrons can be assumed to be in equilibrium with the local reduced electric field using a two-term-approximation Boltzmann equation solver such as Bolsigþ. [28] The pressure at which this assumption is valid is dependent upon the details of the electron impact cross sections for the gas to be modeled. [113] In this approach, the Boltzmann solver uses electron impact cross sections as input to calculate an EEDF, from which electron impact rate coefficients can be derived for use in the global model. [1,4,114] Conversely, under very low pressure conditions on the order of a few Pa, and in highly ionized plasmas, such as in inductively coupled systems, the assumption of a Maxwellian EEDF is generally well justified. However, in the pressure range between these two extremes an EEDF must be chosen which is not necessarily physically well justified. As such the use of global models in this regime requires additional caution.

Beyond Spatiotemporal Averaging
The assumptions of temporal and spatial homogeneity inherent in the basic global model approach can provide significant limitations under certain conditions. In general, these assumptions are borne from the desire to obtain fast solutions of the model. However, provided suitable relations can be found to describe the relevant forms of the temporal and spatial inhomogeneities involved in the system, then global models can, in principle, be extended to deal with such systems in a computationally efficient manner. In Sections 3.1 and 3.2, examples of global models where particular emphasis has been placed on the treatment of temporal and spatial inhomogeneities will be discussed. The discussed treatments represent extensions of the basic global model approach in order to overcome the assumptions of temporal and spatial homogeneity, while still allowing for rapid solution of the model.

Time Varying Power Deposition
In low temperature plasmas, both the shape and effective temperature of the EEDF are closely tied to the spatiotemporal structure of the electric field for any given geometry. In some cases, it is reasonable to assume that neither the shape nor the temperature of the EEDF varies significantly in time. Under these conditions the time independent approximations of conventional global models apply well enough that their results can reproduce experimental values and trends relatively well. [115,116] A common example is the assumption of a Maxwellian EEDF in global models of low pressure inductively coupled plasmas (ICPs) where the effective temperature is solved for using an electron energy balance equation such as Equation (9). In these cases, the EEDF is generally well approximated as Maxwellian, and to a large degree spatially and temporally invariant, while the temperature must still be solved for. However, in cases where significant temporal variations in either the effective temperature or the shape of the EEDF are expected it is necessary to represent these variations in the model in order to achieve an accurate representation of the overall discharge dynamics.
To model such variations fully self-consistently generally requires higher order models inclusive of the solution of Poisson's equation in time and space, and in certain cases the inclusion of kinetic approaches to describe collision-less or non-local heating. [30,117,118] To include time variations in the EEDF in global models without significantly increasing the computational time or the model complexity, it is possible to approximate such variations using analytical expressions, [14,15,108] simplified numerical schemes which do not explicitly deal with spatial variations, [9,16,17,[119][120][121] or inputs from higher order models. [48,73,74] Here, two examples will be discussed where this kind of approach is necessary in order to properly describe the discharge dynamics. The first of these relates to the modeling of temporal instabilities in low pressure electronegative inductively coupled plasmas, and the second deals with temporal variations in the EEDF in radio-frequency atmospheric pressure plasma jets.

Temporal Instabilities
It is known that low pressure ICPs produced in electronegative gases, such as those typically used in plasma processing applications, are susceptible to temporal instabilities, which depend on the operating pressure, discharge power, and feed gas. [9,16,17,119,120] These instabilities typically manifest themselves as temporal variations in the discharge light emission, charged particle densities and electron temperature which in turn can have an effect on other plasma properties. As such, an understanding of the fundamental processes behind these instabilities is important for the optimization of plasma processing applications.
In order to treat these phenomena properly in a global model, it is necessary to first understand their origin, and to then develop a simplified model which is capable of predicting the main features of the instability based on this understanding. Lieberman et al. were the first to develop such a global model of this phenomena. [9] Their model, which has since been built upon and used to investigate several different electronegative plasmas, [16,17,119,120] proposed that the temporal oscillations resulted from the repeated transition of the discharge between capacitive and inductive modes. The basic theory developed to understand this phenomenon relates to the discharge power as a function of the electron density. It was proposed [9] that the total volumetric discharge power is composed of two parts: that due to capacitive coupling and that due to inductive coupling Physically, Equation (17) gives the power transfer between the RF coil, which represents the primary of a transformer, and the plasma, which represents the secondary. The parameter n I is considered the electron density for which the inductive power deposition is maximum. At low electron densities, the plasma represents a weakly conducting loop and acts like an open circuit, leading to low power transfer, while at very high electron densities the highly conducting plasma acts as a short circuit, again leading to lower transferred power. The value of n I represents the optimum point between the two extremes with regard to power transfer. The equivalent parameter for capacitive coupling is n C , although in this case the efficiency of capacitive power transfer simply decreases when n e > n C . The parameters R ind and R cap represent the equivalent resistance of the discharge as a result of inductive and capacitive coupling, respectively. In the original model, the values of n C , R ind , and R cap were free parameters which were varied in order to replicate phenomena observed experimentally; however, more recent publications have refined the understanding of these parameters and related them explicitly to other plasma properties allowing these models to offer enhanced predictive capabilities. [16,17,119] Equation (16) and (17) can be used to specify the input power with appropriately chosen values of I RF . Along with electron temperature dependent rate coefficients for the system, the governing equations of such a model can be solved in a time dependent manner, allowing for the time variation of the discharge properties to be predicted. An example of the results of such a global model applied to a Cl 2 plasma, reproduced from the work of Despiau-Pujo and Chabert [16] is given in Figure 7. Here, the characteristic oscillations in the charged particle densities can be seen, as a result of transitions between capacitive mode, where the electron density is low, and inductive mode, where the electron density is high.

RF Power Deposition at High Pressure
Another scenario in which a time varying power deposition is required in order to properly treat the discharge dynamics is in the simulation of RF atmospheric pressure plasma jets. The high electron-neutral collision frequencies in these systems mean that the electron temperature, and the shape of the EEDF, is strongly modulated at time-scales on the order of nanoseconds, as electrons rapidly lose energy through collisions with the neutral gas. In order to selfconsistently and efficiently describe such rapid variations in the electron temperature in a global model, it is necessary to formulate expressions describing the variation of electric fields in the plasma within the RF cycle, without the explicit solution of Poisson's equation in time and space. Examples of such schemes are described in detail by Lazzaroni et al. [14] and Niemi et al. [121] The cited works rely on the so-called homogeneous RF discharge model, in which the bulk plasma is assumed to oscillate throughout one RF cycle. With this model as a base the time varying motion of the plasma bulk, and consequently the plasma sheath, can be derived allowing for a variation of the electric field, or electron temperature throughout one RF cycle. From the time varying electron temperature ''effective'' rate coefficients for electron impact processes can be derived which account for the varying electron temperature throughout the RF cycle. These effective rate coefficients account for the highly nonlinear variation of electron impact rate coefficients with electron temperature. The importance of this approach, as discussed by Lazzaroni et al., [14,15] is emphasized by the fact that while the instantaneous EEDFs are all assumed to be Maxwellian in shape, the resultant time-averaged EEDF is non-Maxwellian. This means that the assumption of a time averaged EEDF, of any particular shape and temperature, will not necessarily be representative of the effective timeaveraged EEDF calculated with the knowledge of temporal variations. As a result, the electron impact rate coefficients calculated from EEDFs derived from such time-varying global models and those calculated from global models assuming no time variation will differ, leading to differences in the plasma parameters calculated by the models. [15]

Spatial Variation
In low temperature plasmas, the non-uniform spatial distribution of species densities and energies can have a significant impact on the plasma properties. For example, in CCPs, reaction rates can be modulated strongly in both space and time, [122,123] leading to a significant loss of information if spatial averaging is performed. Even in supposedly simple plasmas, considerations must be made as to the spatial properties of the system, such as how the aspect ratio of a cylindrical vessel affects the spatial distribution of species. [124,125] For this reason, models with the ability to resolve at least one spatial dimension often provide a more accurate picture of plasma phenomena than those without. However, as mentioned previously, the inclusion of spatial resolution tends to increase the computation time of the model due to the increased numerical complexity required. Thus, a number of methods have been developed to incorporate some degree of spatial information, without the associated computational load.

The h l Factor
The primary method of introducing spatial considerations into global models is, as mentioned in Section 2.1, the use of a factor to describe the relationship between the central ion density and that at the sheath edge. A number of different expressions exist for this factor, often termed h l , that relate to a particular set of conditions. They are typically derived from analytic assumptions [33,47,126,127] or through the finding of empirical relations. [11,24,44,45,128] For the most simple cases, those without negative ions and having a simple geometry, ''classical'' analytic expressions exist for h l that are dependent only on the collisionality regime of the system. [24,47] For high pressure systems, ions can be assumed to collide with neutrals at a constant rate giving [126] For low pressures, ion motion can be assumed to be collisionless, [127] giving rise to In the intermediate range, ion motion is collisional, but the collision rate is not constant. [128] From this concept, one can derive the expression These ''classical'' expressions for the h l factor in these three regimes apply to the case of an electropositive plasma, with a single positive ion species and a small sheath, contained between two parallel plate electrodes.
These expressions perform well within their defined ranges, but if a plasma transitions from one of these ranges to another, then the expression being used will be rendered invalid. They are also not applicable to more complex systems, such as those containing negative ions. To counteract this first issue, a heuristic fit linking the h l factor in the three pressure ranges discussed above has been proposed by Chabert and Braithwaite, [24] and is written h l; heu % 0:86 3 þ 1 2 This expression agrees well with the results from a 1D analytical model developed recently. [129] However, comparisons with a 1D particle-in-cell (PIC) simulation showed that discrepancies arose at higher pressures, greater than around 30 Pa in a 2 cm wide argon CCP due to the increasing non-uniformity of the ionization rate. [123] Although this and related expressions perform well for the plasmas they are designed to describe, different expressions are naturally required for more complicated systems.
The environment of electronegative plasmas poses a particular challenge for empirically describing the edge to center density ratio, as they can exhibit a number of different structures depending on the operating conditions. [127][128][129][130][131][132] A variety of expressions, based mainly on heuristic descriptions, have been developed over the years, but are often only valid for a very small subset of plasmas. [8,[133][134][135] General expressions have been developed more recently that attempt to provide a more broad description, applicable to a wider range of systems. [44][45][46] These typically consist of an ansatz of limited expressions, empirically coupled to provided a single expression that is able to describe a large range of systems. Similarly to the electropositive case, the constituents of these expressions are descriptions of low, high, and intermediate pressures, often taken from simple 1D models. [44,45] The result of Monahan and Turner [44] is a multicomponent expression, with the constituents written as h a % 0:86 for the case of equal ion temperatures. A subscript n denotes a property of the negative ions, and a subscript 0 refers to a central value.
The results from a global model using this expression was shown to perform reasonably well when compared to a significantly more complicated PIC simulation. [44] More recently, Chabert derived an expression that applies to the low pressure, low density regime of electronegative plasmas. [43] In that work, a 1D model was developed assuming cold positive ions and electrons and negative ions in Boltzmann equilibrium with the plasma potential. This model was then used to heuristically derive an expression for h l for this low pressure regime. The result, is similar to Equation (21), with additional terms to account for the effect of negative ions on both the overall positive ion density as well as the positive ion velocity. The author was careful to state that the resulting expression has a valid parameter range that is more limited than the results of Monahan and Turner [44] and Kim et al. [45] However, within the range of low pressure and density, and where the assumption of Boltzmann equilibrium for the negative ions is valid, the results of Equation (27) are more accurate than those provided by Equation (22). For systems with a particularly high degree of spatial non-uniformity, or those with unusual properties, the behavior may not be able to be captured by a simple density ratio. In particular, cases where there is a significant modulation of reaction rates in space, such as the aforementioned case of high pressure CCPs, [122,123] require some consideration of this. The effect of neutral dynamics is also of concern for high pressure systems, where gas temperatures may be elevated, [136,137] and high density plasmas where the high degrees of ionization lead to neutral gas depletion [138][139][140] through a variety of effects.

Analytical and Semi-Analytical Alternatives
Although expressions have been derived for h l that consider the impact of neutral gas depletion, [47] it is analytically difficult to account for all of the effects that contribute; the role of non-uniform gas temperature is particularly difficult. Due to the coupling of particle and power balance expressions that arises when neutral dynamics are considered, some degree of spatial resolution is required to satisfactorily describe such a system, as the effects of reduced gas density will not be evenly distributed across the discharge.
In order to capture this behavior, while maintaining the low computation times of volume averaged models, a number of authors have developed analytic [141,142] or semianalytic [143,144] models to describe the equilibrium state of a plasma subject to neutral gas depletion. Fructman et al. [144] developed an analytical model to describe an electropositive, quasineutral plasma in the high collisionality regime, with neutral dynamics. That work showed how the aforementioned coupling between particle and energy transport may lead to a counter-intuitive decrease in plasma density for an increase in deposited power under certain conditions. This analytic model was extended by Raimbault et al. [142] to provide solutions for the low and intermediate pressure regimes. The authors found that, in the low pressure cases, it was possible that there is an increase of the neutral density in the center, contrary to what one might expect, but admit that this could only occur in collisionless plasmas with a high degree of ionization, greater than 1%.
These analytic models are, by their nature, quick to solve, but are also restricted to describing general behaviors, due to the number of assumptions that must be made to reach a fully analytic solution. To provide more detail, semianalytical models for plasmas with neutral dynamics, including gas heating, have been developed for electropositive, [143] and later electronegative [144] plasmas in the high pressure limit. The conclusions of these works were that, as expected, the increase of gas temperature in the center of the discharge heightens the degree of neutral gas depletion that is observed.
Such investigation into the neutral properties would potentially not be possible without the spatial considerations of the models described. However, by implementing either analytical approximations, or a simple numerical integration in space, instead of an integration in time, these works were able to retain the short calculation time that makes global models so appealing.

Other Spatial Arrangements
For some systems, where spatial considerations are impractical or impossible to investigate using a 1D semianalytical technique, use of a global model in an unusual fashion may be the best option. A simple yet effective example is the study of gridded ion thrusters, for use as space propulsion devices. These systems are relatively simple in concept, and comprise of a plasma source bound on one side by a set of DC biased grids. Positive ions are accelerated by the grids, and leave the system with high energy, providing a high specific impulse, but low thrust, compared to combustion based propellants.
In order to provide a first analysis of such a system, and to estimate important parameters such as efficiencies, global models have been used. [12,13] The first example considered a xenon plasma driven as a cylindrical ICP capped at one end by a grid. The global model developed by Chabert et al. [12] has similarities with that developed in Section 2.1. However, the neutral species were also considered, with a source term defined by an input gas flow. The model used a different reaction scheme to describe xenon, and a selfconsistent implementation of the power deposition. The thruster aspect was implemented by altering the effective surface area to incorporate a ''semi-transparent'' grid, through which both ions and neutrals can escape. To account for the different interactions that ions and neutrals have with the grid, the transparency coefficients were different for each species, and the ions were accelerated to a speed depending on the grid voltage once they escape. This model was extended by Grondein et al. [13] to include iodine as a feed gas, in order to evaluate its properties as a novel propellant. This extension resulted in the consideration of six species, including negative ions, atomic iodine neutrals, and molecular positive ions, and considerations were made for the effect of the grid on the effective area seen by each species.
In both models, solutions were found by evolving the conservation equations in time until a steady state was reached. From these results, the thruster performance was evaluated. It was found in both cases that the efficiencies of power use and mass use behave in opposing ways. At low input powers and high gas input flows, the thrust power efficiency was high, as a large fraction of the total plasma power leaves as thrust. However, the mass use was inefficient, as the ion extraction rate is low compared to the gas input rate. At high powers, or low flow rates, the situation was reversed and mass use efficiency was favorable. This means that there is a compromise that must be made between efficiency in energy and efficiency in mass À two quantities particularly important for space flight applications. In the extended investigation of iodine, it was found that performance was similar between the two gases, with iodine performing better at low flow rates than xenon. This can be seen in Figure 8.

Coupled Global Models
There are systems where modifications to the global model framework do not provide a suitable description of the spatial arrangement. Again using the example of plasma based thrusters, the physical coupling of a power source, a plasma chamber, and a thruster outlet can be modeled through the numerical coupling of different sorts of models for each part. This was performed in the investigation of a proposed microthruster by Takao and Ono. [145] Their design consisted of a small dielectric plasma chamber, 1 mm in radius, contained in the end of a coaxial cable. This cable carried microwaves from an upstream source. On the end of the plasma chamber was a nozzle for generating increased thrust from the expanding plasma.
In order to model this system, the authors used a 2D electromagnetic model to calculate the interaction of the showing a comparison of overall thruster efficiencies using iodine or xenon as a propellant as a function of mass flow rate. Reprinted [13] with the permission of AIP Publishing. Figure 9. Diagram of the system used by Takao and Ono for the modeling of a microplasma thruster. The electromagnetic, global, and fluid models are depicted as EM, GM, and FM, respectively. Arrows show the transfer of data between the different parts of the system. ß IOP Publishing. Reproduced with permission. All rights reserved. [145] incident microwaves on the plasma, the properties of which were found using a global model. These two models were solved iteratively until the combined system reaches an equilibrium. The plasma properties were then given to a 2D fluid mechanics model to find the properties of the thruster exhaust, as depicted in Figure 9. This arrangement allowed the comparatively simple electromagnetic effects to be captured in 2D, while the more complex plasma system is solved by the global model, thus saving on computation time.
This system of coupled models was used to investigate the impact of adjusting the incident microwave power and frequency and on the power absorbed by the plasma, and the subsequent effects this had on the thrust efficiency of the device. It was found that performance increased greatly if the microwave frequency was tuned to produce standing waves within the device.
A conceptually similar analysis of a helicon plasma thruster was performed by Takahashi et al. [146] In this case, a global model with fixed input power was coupled to a 1D model of a magnetic nozzle, with the aim of investigating how the radius of such a device affects the performance. It was found that for a given power, length, and gas flow rate, the thrust from the device increased for increasing source radius. This agrees in principle with experimental data from the same publication, which investigated the relative performance of two devices differing only in their radii.
Thruster design is not the only situation in which a coupled global model is a sensible investigation method. Any system where there is a large but well defined difference between two regions can be well described by one or more coupled global models. For example, in the description of a high power impulse magnetron sputtering (HiPIMS) discharge, there are two distinct regions that can be described. Close to the surface, there is the interaction region of the plasma with the metal surface, often termed the ''racetrack.'' In the rest of the discharge vessel is an expanding plasma interacting with a magnetic field. These two regions can be well described by an ionization region model and a bulk plasma model, respectively. [147,148] These two models are independent global models, but coupled together describe the whole plasma of a HiPIMS system. This method of coupled global models compares well with experiment, [149] and has been used extensively in investigation of HiPIMS phenomena. [150][151][152][153]

Conclusions
This review has presented an overview of the use of global modeling in low temperature plasmas. Specific focus has been on the elements required to create a simple volume averaged model, and the extensions that can be applied to this basic framework in order to allow to variations in either time or space.
It has been shown that in order for a global model to be created successfully, considerations must be made not just about the components of the plasma, but also about the system in which it is contained. The interaction of the plasma with the surrounding environment, through effects such as power coupling or surface processes, plays a significant role in both the construction and outcome of the model. The choice of plasma components and the reactions is also not a trivial one. It has been shown that researchers can occasionally face great difficulty in deciding which species to consider, which reactions are important for the system of interest, and how to describe them with reaction rate coefficients. This last issue is particularly troublesome, and obtaining well described, accurate rate coefficients for reactions is unfortunately not always guaranteed.
In addition to problems obtaining reliable rate coefficients, it has been discussed how volume averaged models have a number of drawbacks due to the assumptions made in their development. In particular, the pressure range at which assumptions on the EEDF are valid places restrictions on the scope of where global models are accurate. This is on top of the obvious limitations arising from volume averaging, and the removal of important phenomena that this entails.
Despite these deficiencies, global models can be exceptionally simple and quick to develop and implement, as was shown in the first part of this review. The high quality and large quantity of results that can be obtained from these models, coupled with the rapid time of development and execution, means that they continue to see widespread use within the low temperature plasma physics community.
In the second part of this review, an overview of various extension schemes was given. These have been developed by the community to allow the rapid computation times of global models to be applied to systems where the assumptions of spatial uniformity or constant power deposition would give highly inaccurate results. It has been shown that it is possible to construct a model that is able to include time resolved power deposition, either through the dependence of the deposited power on plasma properties, or explicitly resolved over the RF cycle of the discharge.
It has been further shown that it is possible to perform fast investigations of systems that require some spatial information. This can be done by either altering the global model framework, treating the equations in space as opposed to time, or by coupling a global model with other models that incorporate the necessary spatial considerations.
Using these extended methods, global models are able to describe a wide variety of systems, include those that intuitively would require more complex analyses. Their use in low temperature plasma physics is likely to continue in growth, particularly with the increase in application oriented atmospheric pressure systems, where complex chemistries are the intention. There is the possibility though that these tools be used without a complete understanding of their limitations. In particular, the reliance of the chemical kinetics on uncertain rate coefficients means that care must be exercised if a global model is to be used in a predictive manner.