On the justification of the Poisson–Boltzmann equation in the context of technological plasmas

The Poisson–Boltzmann (PB) equation is a nonlinear differential equation for the electric potential that describes equilibria of conducting fluids. Its standard justification is based on a variational principle which characterizes the thermodynamic equilibrium of a system in contact with a heat reservoir as a minimum of the Helmholtz free energy. The PB equation is also employed in the context of technological plasmas. There, however, the standard justification is inapplicable: technological plasmas are neither in thermodynamic equilibrium nor in contact with heat reservoirs. This study presents an alternative variational principle which is based on the functionals of entropy, particle number, and electromagnetic enthalpy. It allows to justify the PB equation for a wide class of technological plasmas under realistic assumptions.


Introduction
The processes of plasma surface engineering-for example plasma activation, plasma etching, or plasma-enhanced chemical vapor deposition-are physically and chemically complex [1]. They involve a hierarchy of coupled subsystems-the power source, the transmission network, the process gas supply, the reaction chamber, the plasma, and the surface of the workpiece-, and incorporate phenomena on various length and time scales. Correspondingly difficult and resourceconsuming is their mathematical modeling and simulation. * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. This applies even when only individual submodels are considered. The plasma model, for example, which is the focus of this research, consists in its most general form of a particle sector with six-dimensional kinetic equations for the phase space distribution functions, and a field sector with the Maxwell equations for the electric and magnetic fields. The relevant length scales span nearly six orders of magnitude; the relevant time scales even more [2].
Within the range of currently available computational resources, resolution of the full model is beyond reach. For practical simulation work, drastic simplifications have to be introduced. The arsenal is plentiful and diverse [3]. Of course, all simplifications must be justified with sound physical arguments, and their applicability must be carefully checked.
The subject of this research project is a particularly powerful mathematical simplification, the Poisson-Boltzmann (PB) assumption. It replaces Maxwell's equations and the electron model by a nonlinear second order differential equation for the electric potential Φ. In standard form it reads as follows, where ε 0 is the electric constant, e the elementary charge, T e the spatially constant electron temperature, n i (r) the net ion density with a prescribed, not necessarily homogeneous distribution, andn a density constant that represents the electron density at the reference point of the potential: The PB equation is of importance also in other research areas. In their recent topical review 'nonlinear electrostatics: the PB equation', Gray and Stiles listed over 20 branches of natural science and engineering, and gave more than 60 citations [4]. They called the PB equation applicable for 'a conducting medium in thermal equilibrium, such as an electrolyte solution or a plasma', and gave a detailed derivation based on the extended Helmholtz free energy functional F.
For an electrolyte solution of singly charged anions and cations, the Gray and Stiles version of the PB equation reads (in our notation) as follows, where ε r is the relative permittivity of the solvent (e.g. ε r ≈ 80 in the case of water),n the bulk concentration of both species (when the fluid is neutral and the potential zero), and ρ f a fixed charge density, for example that of a molecule whose solvation behavior may be studied: (2) PB equations (1) and (2) share obvious similarities but have also profound differences: • Both equations invoke Poisson's equation to relate the Laplacian of the potential Φ to the charge density ρ. (Equation (1) for a plasma assumes vacuum, equation (2) takes the possible presence of a solvent with permittivity ε r into account.) The underlying assumptions, however, are physically different: equation (2) considers a static situation, and Φ is a true electrostatic potential. Equation (1), in contrast, only invokes the so-called 'electrostatic approximation': the time derivative of the magnetic field is neglected in Faraday's equation, ∇×E ≈ 0, allowing to set E = −∇Φ. It is not necessarily assumed that the situation is static. • Both equations express the charge density in functional dependence of the potential. In equation (2) all mobile charge carriers-positive and negative ions-are assumed to follow a Boltzmann distribution. In equation (1), this applies only to electrons; the description of the ion dynamics is not part of the model. This means that the ion density n i (r) is treated as fixed, similarly as the fixed charge density ρ f (r) in (2). • Both equations are of thermodynamic appearance. However, only equation (2) assumes a true thermodynamic equilibrium, where T is the temperature of the environment. As stated, reference [4] justifies equation (2) via a variational principle that employs the Helmholtz free energy F = U − TS, where U is the internal energy and S the entropy. References [5][6][7] also employ the free energy F. Equation (1), in contrast, applies to a dissipative system that is, even when stationary, far from thermodynamic equilibrium. The value of the electron temperature T e , which is much higher than the temperature of the ions and the neutrals, does not arise from a contact with the environment but from the balance of competing energy gain and loss processes.
It is evident that the thermodynamic logic of equation (2) is not suited to justify equation (1), at least not without additional arguments. It is our intention to provide such an argument. We will argue that it is not necessary to assume that the system is in 'true equilibrium'; it suffices to ask that it is in 'quasiequilibrium'. This is a reference to a time scale argument: Under precisely defined conditions, the dynamics of a wide class of technological plasmas can be divided into two wellseparated temporal regimes. When analyzed within the fast regime, the system is thermally insulated and its stable states prove as thermodynamic equilibria. Only in the slow regime these equilibria must be interpreted as dissipative steady states that are maintained by a constant supply of energy and matter. (Compare [8,9].) The rest of this manuscript is organized as follows. In section 2, we discuss the physics of a concrete example, a cylindrically symmetric inductively coupled plasma of laboratory scale, motivated by a discharge experimentally studied by Godyak et al [10]. Using a global model and a kinetic description in two-term-expansion, we can define the fast and slow temporal regimes explicitly. Section 3 provides a more abstract definition of the fast regime which covers also other discharge configurations. Section 4 discusses the corresponding balance equations of particle number, electromagnetic enthalpy, and entropy. Section 5 constructs a thermodynamic variational principle and shows that it allows to justify the PB equation. Summary and outlook conclude the study.

Time scale analysis of an inductively coupled plasma
In order to motivate the assumptions of our abstract approach, we start with an explicit example and consider an inductively coupled argon plasma of laboratory size (see figure 1). Cylindrical symmetry is assumed. The plasma chamber V has a radius of R = 0.099 m and a height of H = 0.105 m; its volume is V = 3.2 × 10 −3 m 3 and its surface area A = 0.13 m 2 . The mantle has an area of 0.065 m 2 , the base (top and bottom together) an area of 0.065 m 2 . (All data are taken from the experimental study of Godyak et al [10].) The coil chamber V c , separated from the plasma by a thin dielectric window, is of similar size. V and V c together form the electromagnetic chamber with ideally conducting grounded walls. The discharge is driven by an RF power of P RF = 200 W at a frequency of f RF = 6.78 MHz. The gas pressure is p = 100 mTorr = 13.3 Pa. Assuming a gas temperature of T g = 650 K, the gas density is n g = 1.5 × 10 21 m −3 .
To calculate the plasma parameters, we turn to the global model provided in reference [1]. In its core is an ambipolar  [10]. It consists of an axially symmetric plasma chamber V of radius R = 9.9 cm and height H = 10.5 cm, topped by a coil chamber V c of similar size. The domains are separated by a thin dielectric window. Together, they form the electromagnetic chamber V ∪ V c ; it is assumed to have walls that are ideally metallic and grounded.
diffusion equation for the plasma (=electron and ion) density n e . It is assumed that the plasma is collisional, H/λ i ≈ R/λ i ≈ 150 1, with λ i = 6.7 × 10 −4 m as the ion mean free path, and its transport linear, with μ i = 2.6 m 2 V −1 s −1 as the mobility [11]. The ambipolar diffusion constant is then D a = μ i T e /e, where T e is the electron temperature. Expressing the plasma generation via the ionization frequency ν iz , the equation reads The eigenvalue problem posed by the stationary version of (3) is solved by a diffusion profile, where J 0 is the Bessel function of order zero and j 0,1 ≈ 2.4 its first zero: The constantn e signifies the density in the discharge center (r = 0, z = H/2) which differs from the average densityn e . The average-to-center ratio is n e /n e = 4J 1 j 0,1 /π j 0,1 ≈ 0.275. The corresponding eigenvalue ν iz is given as follows, where the implicitly defined quantity Λ is known as the effective diffusion length [1]: Further equations (not displayed here) describe the evolution of the volume-averaged densities of the 4s metastable and resonance levels and the 4p levels, accounting for excitation, radiative and collisional de-excitation, ionization, and particle losses to the chamber walls. An escape factor is introduced to account for the possible retrapping of radiation [12,13]. Lastly, assuming a spatially constant electron temperature, a balance for the mean electron energy equates the absorbed RF power to the energy losses from ionization and excitation and for the kinetic particle energies carried off to the walls.
The discharge is in a stationary state when all particle balances and the energy balance are fulfilled with vanishing time derivatives. Adjusting the escape factor so that the relative ratio of stepwise excitation is 60% (based on [14]), we obtain T e = 1.8 eV andn e = 7.3 × 10 17 m −3 . This is in good agreement with the experiments of Godyak, Piejack, and Alexandrovich who obtained an electron temperature of T e = 1.5 eV and a center density ofn e = 4.0 × 10 18 m −3 , the latter value corresponding to an average density ofn e = 1.1 × 10 18 m −3 [10].
It is interesting to explore the consequences of these results. About 99.5% of the electrons are in the elastic range with energies less than the excitation threshold of * = 11.55 eV. The Debye length is λ D = 1.2 × 10 −5 m, the electron elastic mean free path λ e = 6.4 × 10 −3 m, and the skin depth λ s = 6.2 × 10 −3 m. The electron plasma frequency is ω pe = 4.8 × 10 10 s −1 , the ion plasma frequency ω pi = 1.8 × 10 8 s −1 , and the collision frequency ν en = 8.7 × 10 7 s −1 . The energy relaxation length is λ = m i /2m e λ e = 1.2 m, which exceeds the discharge size. The Coulomb logarithm is ln Λ = 10.7, and the electron-electron Coulomb collision frequency for thermal electrons is ν ee = 3.6 × 10 7 s −1 . This exceeds by far the residual cooling rate due to elastic collisions 2ν en m e /m i = 2.4 × 10 3 s −1 . In the elastic range, the plasma electrons thus form a single thermodynamic body. Interestingly, this holds also in the boundary sheath, as the floating voltage is V sh = ln m i /2πm e + 1 T e /2e = 9.2V.
The assumption of a thin, collision-poor sheath, implicit in this formula, is justified because the ion mean free path λ i sufficiently exceeds the sheath thickness s = 3 × 10 −4 m (calculated with the Bohm model). The model can also be used to estimate time scales. In the vicinity of the stationary state, the equations can be linearized to describe the dynamics of superimposed perturbations. All eigenmodes have a negative real part; the equilibrium is therefore asymptotically stable. The electron temperature reacts with a time constant of τ = 1.0 μs, substantial heating thus requires several ten RF cycles, ω RF τ = 43. The reaction on plasma density fluctuations, governed by ν −1 iz = Λ 2 /D a = 150 μs, is a factor of 100 slower.
The global model implicitly assumes that thermodynamic equilibrium is always established. To analyze this in detail, a kinetic model is required. We adopt the two-term approxima- [1][2][3]15].) A scalar equation governs the isotropic part f (0) (r, z, v, t) of the electron velocity distribution. The collision term on the right accounts for electron-electron Coulomb interaction and for inelastic collisions with the residual cooling resulting from elastic electron-neutral collisions and electron-ion interaction (finite mass ratio effects) included. In order to properly account for the heating process, terms up to second order must be kept: The vector equation for the anisotropic part f (1) is split into its three components. The direct magnetic force is neglected. The friction coefficient ν m = ν en + ν ei = 1.2 × 10 8 s −1 accounts for the momentum loss due to elastic and Coulomb collisions: For the electromagnetic field, we start from Maxwell's equations under cylindrical symmetry. We introduce the potentials Φ and A, both defined in the entire electromagnetic chamber and assumed to be continuous through the separating window: The vector potential A is assumed to follow the Coulomb gauge, 1 r Poisson's equation then relates the scalar potential Φ to the net charge density ρ via The three components of the vector potential are governed by wave equations with the sum of the electron current density and the longitudinal displacement current acting as sources. Let us focus on the r and z components first: The inhomogeneities on the right are the total currents in r and z-direction, respectively. There cannot be any large scale contributions, due to the geometry of the set-up ( figure 1). Small scale current fluctuations, on the other hand, result in negligible values for A r and A z , due to the mathematical properties of the inverse of the d'Alembert (or wave) operator. Altogether, we can write the electric field E as follows, which is equivalent to adopting the electrostatic approximation in the r and z directions: The component A θ , however, is excited by the coil current and thus important: The next step of the argument serves to carefully separate the fast and the slow time scales. As stated above, the model includes the fast scales of plasma oscillation ω pe = 4.8 × 10 10 s −1 , collisional friction, ν m = 1.2 × 10 8 s −1 , and electron-electron thermalization ν ee = 3.6 × 10 7 s −1 . We now assume that this dynamics results in a relaxation of the system into a stable state, and that any further evolution is governed by the slow time constants of the global model. (A rigorous justification of these assumptions is subject of the later sections of this study.) It is important to understand that the presence of RF modulation does not conflict with this: although the RF modulation is fast, as ω RF = 2π f RF = 4.26 × 10 7 s −1 , its amplitude is weak, so that heating takes many RF cycles. As estimated above, ω RF τ = 43 1. Its worth-while to discuss this argument in more detail. The assumption of weak modulation can be carried over to the scalar potential Φ and the first order distributions f (1) r and f (1) z . For the former, this follows from Poisson's equation, and for the latter then from (7) and (9). Quasistatic behavior can be assumed, as the friction coefficient ν m exceeds the reaction speed of the slow dynamics by two orders of magnitude: We now turn to the dynamics in θ-direction which is directly excited by the induced RF. Equations (8) and (17) are linear in the quantities f (1) θ and A θ , while the distribution f (0) and the friction coefficient ν m are only weakly time-varying. It is therefore possible to make a time-harmonic ansatz in the RF frequency, (8) can then be solved to The ratio between the complex current density j θ and the complex electric field E θ is the complex conductivity of the plasma; it can be evaluated to Extending the definition of σ to the coil chamber (with a value of zero), a closed equation for the vector potential A θ , continuous across the domain boundary, can be given: We now return to equation (6). All terms that influence the evolution of f (0) vary only slowly, except the term which contains E θ f (1) θ . Here, however, we can split this term into its phase average and an average-free fluctuation in the frequency 2ω RF (with α being a constant): This relation motivates the definition of the heating operator T RF , a second order differential operator in the velocity v that acts on f (0) in a quasilinear way: Taking (18) and (19) as abbreviations, equation (6) turns into a parabolic equation for f (0) . The influence of the fluctuating heating term can be neglected; it merely creates a ripplef (0) that scales asf (0) / f (0) ∼ 1/(ω RF τ ) 1. This justifies the assumption of weak modulation. Together with the field equations (13) and (22), we now have a model of the dynamics: Suitable boundary condition must be posed. The simplest possibility is that of a perfectly absorbing wall at ∂V: all electrons which reach the wall of the plasma chamber are absorbed. Due to the presence of a boundary sheath, this applies only to a small number of electrons; most electrons are reflected by the sheath potential: Coupled to the electron dynamics are the equations for the potentials Φ and A θ which also must be equipped with boundary conditions. As the walls of the electromagnetic chamber are ideally metallic and grounded, one can demand Φ| ∂(V∪V c ) = 0.
As this point, the modeling procedure has achieved its goal, and a numerical solution of the equations could be sought. This is indeed the path that is chosen by many groups [15,16]. We, however, pursue a different approach: the analysis above has shown that the model still contains a number of different time scales. Let us group them in two regimes: • The slow dissipative regime consists of all phenomena related to the transport of ions, the generation and recombination of particles, the acquisition and dissipation of energy, and the charging of surfaces. This regime was already addressed by the global model. It acts on a time scale of τ = 1 μs or slower. • The fast thermodynamic time scale describes the transport of the electrons caused by drift and diffusion forces, and the electron thermalization via Coulomb interaction. It acts on the time scale of ν −1 ee = 30 ns or faster. We emphasize that is was only assumed that the discharge relaxes to a stable equilibrium on the fast time scale, it was not explicitly shown. We have to make up for what we have missed. It is straight-forward to define a model that only acts on the fast scale: the ion density must be assumed as stationary, and dissipative terms in the electron dynamics must be dropped. The resulting model reads, with f (1) r and f (1) z still given by (18) and (19): The boundary conditions must now guarantee that no electrons are absorbed at the walls. Note that the following assumptions imply Boltzmann equilibrium close to the walls: The Poisson equation and its boundary conditions continue to hold. The resulting model is compatible with the assumption of quasi-neutrality in the bulk and the emergence of a non-neutral sheath at the boundaries.

A more general model of the fast regime
Equipped with insight on a concrete example, we now proceed with our general program. For broader applicability, we define a more general model which avoids the assumption of cylindrical symmetry and the two-term-expansion. Moreover, we allow for the presence of electrodes with a non-vanishing constant potential (see figure 2): • The model describes only the domain V, i.e. the interior of the plasma chamber. Physical processes that take place in other regions of the device are not considered. • The model describes only the electrons and the field. Ions are assumed to be stationary. Neutral gas dynamics and surface processes are not addressed. • Electron-electron Coulomb collisions are included, as well as elastic electron-neutral and electron/ion collisions under the assumption of immovable scattering centers. • Dissipation, heating, and chemistry in the plasma or at the walls are not considered. • It is assumed that the electrostatic approximation holds.
The plasma chamber is a simply connected spatial domain V ⊂ R 3 with closed boundary ∂V. The electrons are described by a distribution function f = f (r, v, t) under a kinetic equation with time derivative and convective terms on the left and the terms of elastic scattering as well as electron-electron and electron-ion Coulomb interaction on the right: (32) Poisson's equation describes the electric potential: The boundary condition for the electrons assumes specular reflection: The boundary conditions for the potential assume that the domain boundary ∂V consists of grounded areas G, electrodes E k with k = 1 . . . N E , and insulating surfaces I: • Vanishing potential at the grounded surfaces G: Φ| G = 0, • Fixed potential values at the electrodes E k : • Constant surface charges at the insulating surfaces I:

The balance equations of the fast regime
Central for our further argumentation are the global balance equations of the fast regime. The first to consider is the electron number balance equation which makes a statement on the total number of electrons in the system, Integrating the kinetic equation (32) over velocity space with weight unity and making use of the fact that the collision terms on the right are particle conserving, we obtain the electron continuity equation with vanishing source term, Here, the electron density and the electron flux have their standard definition, Integrating over the domain of the discharge and utilizing the reflecting boundary conditions, we obtain the statement that the total electron number is conserved: The electron energy balance makes a statement on the kinetic electron energy: Integrating the kinetic equation (32) over velocity space with weight 1 2 m e v 2 and utilizing that the collision terms are energy conserving, we obtain the local energy balance The energy density and the energy flux are defined as Integrating over the domain of the discharge and utilizing the reflecting boundary conditions, we see the electron energy is only changed by the action of the field: We account for the energy of the electric field in electrostatic approximation, The corresponding local energy balance reads [17]: Integration over V yields the total power fed into the plasma.
Note that the normal current vanishes at ∂V due to the reflection conditions: The internal energy U is the sum of the electron energy W e and the field energy W el , It is not conserved because energy can be exchanged with the powered sources: (49) Here, the current I k is counted positive when it flows into the discharge, Due to the reflecting boundary conditions, the current accumulates as a surface charge: Charge and current are related as Under the assumption that the voltages V k are constant with respect to the fast regime, we can write the surface term as a total time derivative: (53) This motivates to define a quantity which, due to the analogy to classical thermodynamics, will be called the electromagnetic enthalpy. (Compare with [18].) The enthalpy is conserved under the action of the fast dynamics: Lastly, we study the negative entropy −S which is sometimes called the 'negentropy' [19]. It is a nonlinear functional of f , where Δμ measures the phase space elementary cell [20]. (Δμ is a constant of dimension m 6 s −3 to render the argument of the logarithm dimensionless. A quantum mechanics based choice would be Δμ = (h/m e ) 3 , where h is Planck's constant.) Boltzmann's constant k B is omitted, i.e. set equal to unity; it only serves to convert between units of energy (eV) and temperature (K): The corresponding local balance equation is Its right side is non-positive, which is a form of the Boltzmann H-theorem that holds for all physically reasonable collision terms [21]. (It is, actually, a manifestation of the second law of thermodynamics, a general physics principle with deep roots. See, for example, [24].) Owing to the presence of electron-electron Coulomb collisions, the equality sign is only assumed for distribution functions that are locally Maxwellian.
Integrating over the domain and taking the reflecting boundary conditions into account yields the statement that the negative entropy is a non-increasing quantity:

A variational principle for the fast regime
The results of section 4 allow to characterize the equilibria of the fast regime in terms of a variational principle. It was shown that the negative entropy is a non-increasing quantity, while the electron number and the enthalpy are conserved. Lyapunov's second argument then assures that the constrained minima of the negative entropy yield stable equilibria [22,23]. To employ that argument formally, all quantities must be expressed as explicit functionals of the distribution function f . For N this is already achieved, For the enthalpy H, we employ the formal solution of the Poisson equation (see appendix A). As defined there, G(r, r ) is the Greens function of the system, Φ (o) (r) the outer potential, and the charge density ρ stands for The enthalpy H can now be written as follows. The last term, the energy of the outer field, does not depend on the distribution f but only on the boundary conditions: The negative entropy is a functional of f as well: Using a standard technique, the unconstrained functional to be consulted is the following, where Λ N and Λ H are Lagrange multipliers and the prime indicates that the last term of H was dropped as being independent of f [24]: The first variation of L with respect to f is In an equilibrium, the first variation vanishes for all δ f . Reintroducing the potential Φ(r), see appendix A, the condition collapses to Solving for f yields the equilibrium condition (66) Renaming the constants, this gives a Maxwellian distribution in Boltzmann equilibrium: The corresponding electron density expression is: The second variation δ 2 L is positive definite. For the second term of (68), see the appendix A. All extrema of L are thus stable equilibria: Altogether, this demonstrates that under the stated assumptions, the PB equation holds for the fast regime: The boundary conditions are: • Vanishing potential at the grounded surfaces G: Φ| G = 0, • Fixed potential values at the electrodes E k : Φ| E k = V k , • Constant surface charges at the insulating surfaces I: n·∇Φ| I = σ(r)/ε 0 .

Summary and outlook
The subject of this study was the justification of the PB equation in the context of technological plasmas. Equation (1) is similar to the PB relations of other fields, as it describes an equilibrium of mobile charge carriers in terms of a nonlinear differential equation for the electric potential. There are, however, fundamental physical differences. The plasma PB equation covers only electrons, not all charge carriers, and does not address a system in thermal contact with the environment. Consequently, the standard justification of the PB equation via a variational principle in the free energy does not apply.
In search of an alternative justification, we invoked the concept of a 'quasi-equilibrium'. Under certain assumptions-sufficiently low pressure and sufficiently high electron density-a two-time-scale formalism could be set up. On the fast scale, a restricted electron dynamics holds that includes elastic and Coulomb interaction but excludes heating and collisional loss. The ion density and the electrode voltages are frozen. The regime conserves the electron number and the enthalpy that includes the electron and field energies and accounts for the energy exchange with the voltage sources. The associated entropy follows the H-theorem. A variational principle results that allows to characterize the equilibria of the fast regime. They are stable, have Maxwellian electron energy distributions, and follow the PB equation. As emphasized, these equilibria are 'true' thermodynamic equilibria only in the restricted dynamics of the fast regime. Physically, they are quasi-equilibria that are actually dissipative when analyzed within the slow regime. In particular, the steady state value of the electron temperature T e results from competing energy gain and loss processes.
Of course, the presented argumentation fails when the underlying assumptions are not met. For once, the scale separation is typically invalid for electrons which have enough energy to excite neutrals or to leave the plasma. However, these electrons form only a small minority; their charge density is negligible and does not undermine the validity of the PB equation. More important is the question whether the electrons in the elastic energy range below the first excitation level fulfill the assumptions. As stated, this requires the electron density to be sufficiently high so that the energy relaxation rate-measured by the Coulomb collision frequency ν ee -exceeds the cooling rate due to elastic interaction with neutral and ions. The investigations of section 2 and the measurements of reference [10] show that inductively coupled plasmas of moderate pressure and high power are good examples.
It may be possible to transfer the essence of our approach also to other discharge regimes. Two research directions come to mind: one could study discharges of lower plasma density where Coulomb collisions are less effective and the electron distribution is not Maxwellian. 'Thermodynamics' as such would not apply, but many aspects of it would continue to hold, provided the pressure is low enough so that the discharges are in the non-local regime [15]. In a different context-active plasma resonance spectroscopy (APRS)-members of our group analyzed such a situation and established a modified version of the PB equation [26]. A second research direction could address the regime of magnetized technological plasmas: magnetically confined high density plasmas as employed, for example, in Hall thrusters [27], magnetically enhanced hollow cathode arc discharges [28], or high power magnetrons [29,30]. The modeling of such discharges, however, bears additional difficulties which arise from the complex geometry of the magnetic field and the emergence of self-organized structures [31]. The need for making inroads in the field was emphasized in the 2017 Plasma Roadmap [32]. We are working on addressing the problems with our specific methods [33][34][35][36]. and its normal vanishes for boundaries with Neumann conditions, n ·∇ G(r, r )| r ∈I = 0.
Green's second theorem states: Formulating it for primed variables, taking φ(r ) = Φ(r ), ψ(r ) = G(r, r ), and making use of the boundary conditions, we obtain the formal solution of the Poisson problem: The first contribution is termed the inner potential Φ (i) (r); it represents the potential due to the volume charge density ρ, under homogeneous boundary conditions: • Vanishing potential at grounded surfaces G: Φ (i) G = 0, • Vanishing potential at electrodes E k : Φ (i) The last two terms form the outer potential Φ (o) (r). It fulfills the boundary conditions of the original problem and corresponds to a vanishing charge density in the domain: ∇ G(r, r )·d 2 r + 1 ε 0 I G(r, r )σ(r )d 2 r.
(A6) Several useful properties can be demonstrated. The first insight is that the energy of the inner field can be expressed as a quadratic form in the charge density: Now consider the scalar product of the inner and the outer field. The first term of the second identity rests on boundary conditions, the second term involves the Poisson equation: Thus, we have an identity that is used in the main text: (A9)