Thermodynamic quasi-equilibria in high power magnetron discharges: a generalized Poisson–Boltzmann relation

The Poisson–Boltzmann (PB) equation is a nonlinear partial differential equation that describes the equilibria of conducting fluids. Using a thermodynamic variational principle based on the balances of particle number, entropy, and electromagnetic enthalpy, it can also be justified for a wide class of unmagnetized technological plasmas (Köhn et al 2021 Plasma Sources Sci. Technol. 30 105014). This study extends the variational principle and the resulting PB equation to high power magnetron discharges as used in planar high power pulsed magnetron sputtering. The example in focus is that of a circular high power magnetron. The discharge chamber and the magnetic field are assumed to be axisymmetric. The plasma dynamics need not share the symmetry. The domain is split into the ionization region close to the cathode where electrons are confined, i.e. can escape from their magnetic field lines only by slow processes such as drift and diffusion, and the outer region , where the electrons are largely free and the plasma is cold. With regard to the dynamics of the electrons and the electric field, a distinction is made between a fast thermodynamic and a slow dissipative temporal regime. The variational principle established for the thermodynamic regime is similar to its counterpart for unmagnetized plasmas but takes magnetic confinement explicitly into account by treating the infinitesimal flux tubes of as individual thermodynamic units. The obtained solutions satisfy a generalized PB relation and represent thermodynamic equilibria in the fast regime. However, in the slow regime, they must be interpreted as dissipative structures. The theoretical characterization of the dynamics is corroborated by experimental results on high power magnetrons published in the literature. These results are briefly discussed to provide additional support.


Introduction
Magnetron sputter deposition is a vacuum-based technique for depositing both metallic and dielectric (composite) thin films that is used in numerous industrial applications [1][2][3][4]. At the atomic level, sputtering is the process by which atoms, to be deposited on a substrate, are released from a negatively biased target (cathode) via bombardment by high-energy ions. Magnetron sputtering uses a magnetic field to trap the cathodeemitted secondary electrons and the plasma-born thermal electrons above the 'race track', thereby increasing the ionization rate of the gas, the flux of bombarding ions, and the sputtering and deposition rates. Ions, in contrast, are not trapped due to their high mass. A relatively recent development is high power impulse magnetron sputtering (HiPIMS) which utilizes power densities of the order of kW cm −2 in short pulses of tens to hundreds of µs at duty cycles of 0.5%-5% [5][6][7]. (Peak power and duty cycle are chosen so as to maintain an average cathode power similar to conventional direct current sputtering.) HiPIMS sources have plasma densities of up to 10 20 m −3 , leading to a high degree of ionization of the sputtered metal atoms and a high degree of dissociation of the process gas. The resulting ions are accelerated onto the substrate and hit there with relatively high energy, improving the quality of the deposited films [8][9][10].
Despite significant progress in recent years [11][12][13][14][15][16][17][18][19][20][21][22], comprehensive modeling and simulation of HiPIMS is still in its infancy. Given the complex mechanisms which are involved, this statement may not be very surprising. Multi-scale and multi-physics modeling is needed. In order to predict the properties of the deposited thin film, it is necessary to consider not only surface physics and chemistry, but also the processes in the discharge volume involving energetic and thermal electrons, neutral and ionized atoms and molecules (including the material sputtered off the target), the electric and magnetic fields, and possibly UV radiation. A number of individual models must be defined and coupled together. The resulting overall model is difficult to analyze and computationally highly demanding, as is evident from figure 1. This is all the more true as it has turned out that an in-depth analysis of the HiPIMS process requires a fully three-dimensional approach. In fact, it has been found that even when the discharge geometry has a certain symmetry, e.g. is axisymmetric (as in this investigation), the resulting plasma dynamics may not share it, contrary to what appears to the naked eye. For example, localized ionization zones (spokes) can form in these plasmas which propagate at several km s −1 and give rise to potential fluctuations in the MHz range [23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42].
In this situation, mathematical simplifications are welcome, provided they are justified by sound physical arguments. This study addresses the Poisson-Boltzmann (PB) assumption, a simplification that circumvents the difficult task of solving the fundamental description of the thermal electron component and the electromagnetic field-a six-dimensional kinetic equation nonlinearly coupled to Maxwell's equations-by assuming a Maxwellian energy distribution function for the electrons and setting up a three-dimensional elliptic differential equation for the electric potential Φ. Schematic of a simulation framework for the HiPIMS process: the objective is to predict the properties of the thin film deposited on the substrate in dependence of the characteristics of the target (which acts as cathode), the electrical power, and the supply of process gas and feed gas. The physical and chemical processes on the surfaces are mediated by the physical and chemical processes in discharge volume. Of importance are the species of the plasmaions, energetic and thermal electrons-the neutral gas, the electric and magnetic fields, and possibly the UV radiation. Multi-scale and multi-physics modeling is obviously required for a full description of the full process. The considerations of this study concern the self-consistent interaction of the thermal electrons, the magnetic field, and the electric potential.
In other branches of science, for instance in physical chemistry, the PB equation describes conductive fluids (electrolytes) that are in thermal equilibrium with their environment [43]. It can be justified by a variational principle based on the Helmholtz free energy F = U − TS, where U and S are the energy and entropy of the fluid in the considered domain, respectively, and T is the temperature of the environment. However, technological plasmas are far from thermodynamic equilibrium and their stationary states are instead dissipative structures [44]. They require a constant supply of electrical power. In particular, the value of the electron temperature does not result from equilibration with the environment. It is typically much higher than the temperature of the gas and the wall and arises from the balance of competing electron energy gain and loss processes.
In a recent study, we gave an alternative justification of the PB equation which applies to a large class of unmagnetized technological plasmas [45]. For a realistic example we focused on an inductively coupled argon discharge investigated experimentally by Godyak et al [46]. (Chamber radius R = 99 mm, chamber height H = 105 mm, argon gas pressure p = 13.3 Pa, radio frequency (RF) power P = 200 W at a frequency of f = 6.78 MHz, electron density n e = 4 × 10 18 m −3 .) In standard notation-ε 0 is the electric constant, e the electron charge, m e the electron mass, T e the spatially and temporally constant electron temperature, n i (r) the given ion density, and n a density constant-the PB equation reads The corresponding electron density and Maxwell-Boltzmann distribution function are f(r, v) = n e (r) The key idea of our approach was to replace the concept of thermodynamic equilibrium by the concept of thermodynamic quasi-equilibrium. Calculating the relevant time scales, we found two well-separated temporal regimes, a fast (30 ns or faster) 'thermodynamic' and a slow (1 µs or slower) 'dissipative' regime. The fast regime covers electromagnetic forces, the motion of electrons across the discharge domain, electron-electron Coulomb interaction, and elastic interactions with ions and neutrals in the limit of immobile scattering centers. The slow regime consists of all phenomena related to inelastic collisions, the generation and recombination of particles, particle losses to the wall, and the ion dynamics.
In the thermodynamic regime, the electron fluid is adiabatically isolated both from other particles and from the chamber walls; neither particle exchange nor heat exchange take place. This ensures the conservation of the particle number N and the electromagnetic enthalpy H. (The enthalpy, defined in analogy to the classical thermodynamic potential of the same name, accounts for the energy of the plasma and the energy exchange with external voltage sources. We note in passing that the conservation of H is related to the fact that the voltage sources were assumed as ideal. This assumption will be critically reexamined later in this study.) Furthermore, the second law of thermodynamics holds, which implies that the entropy S of the electron component is non-decreasing. Combining the statements, a variational principle was established to show that the system formed by the electrons and the electric potential relaxes to a thermodynamic equilibrium on the fast scale, described by the PB equation. The variational principle also determines the form of the electron distribution function as that of a Maxwellian in Boltzmann equilibrium with a constant electron temperature T e . For the chosen example, this is confirmed by Langmuir probe measurements [46]. Seen in the slow regime, however, the stationary states are no thermodynamic equilibria. Instead, they are dissipative structures. There is a steady flow of energy through the system, entering as electromagnetic energy and exiting as kinetic particle energy, heat, or radiation. Entropy is constantly being produced and transported out. Moreover, the free parameters of the PB equation and its solution, the density constantn and the electron temperature T e , are determined by competing gain and loss processes. In summary, we were able to classify the stationary states under consideration as thermodynamic quasiequilibria.
The concept of thermodynamic quasi-equilibrium provides insight and, at the same time, offers a mathematical simplification. Two conditions must be met for its physical validity. First, there must be enough interaction for the electron component to act as a single entity. This sets a size constraint: the system must be smaller than its own energy relaxation length. Second, the time scale separation must be sufficiently pronounced. The thermal relaxation within the electron component must act much faster than the dissipation processes which link the electron component to the environment. This yields a condition on the plasma density. In the inductively coupled plasma used as example, it must be greater than about 10 18 m −3 . (The time-scale separation applies only to thermal electrons. Coulomb interaction decreases with energy and inelastic collisions increase, so deviations from Maxwell's distribution may appear above a certain threshold energy. However, this does not affect the validity of the PB equation as long as the threshold energy is sufficiently high.) As stated, these results were derived for unmagnetized plasmas. Can they be transferred to magnetized plasmas, specifically to high power magnetron discharges? The answer to this question is not obvious. On the one hand, the thermalization time requirement is clearly met: at a (realistic) plasma density of n e = 5 × 10 19 m −3 and an electron temperature of T e = 3 eV, the electron Coulomb thermalization time is below 2 ns. This is very fast compared to the dissipative time scales which are in range of several µs. (See below for a detailed discussion.) There is also experimental and theoretical evidence that the electron energy distribution in high power magnetrons is indeed Maxwellian [16, 18-22, 25, 37, 47-51], excluding, of course, the cathode-born minority of energetic secondaries. On the other hand, the thermal electron component in a high power magnetron does not constitute a coherent thermodynamic entity. Transport across the magnetic field lines is suppressed; drift and diffusion processes are slow. The individual magnetic flux tubes are in fact adiabatically isolated from each other and form self-contained thermodynamic units that interact only through electrical forces.
The aim of this research is to examine the matter in more detail and to show that the quasi-thermodynamic approach can also be justified for high power magnetron discharges. We will use many of the ideas and results from our previous work, but now we will explicitly consider the magnetization of part of the discharge. Our manuscript is structured as follows: after this introduction, section 2 will examine the geometrical and physical properties of the discharge, with a specific focus on the temporal dynamics. Section 3 will introduce a thermodynamic model that covers both the magnetized and non-magnetized domains. In section 4, the balance equations will be discussed. In section 5, a new variational principle will be established, and a modified PB equation will be derived. Finally, section 6 will provide a brief summary and outlook to conclude the study.

Geometrical and physical considerations
We consider a circular planar magnetron with a static axisymmetric magnetic field B as shown schematically in Figure 2. Cross section of a circular planar magnetron and its magnetic field topology (to scale), mounted in a cylindrical plasma chamber (not to scale). The chamber wall and the anode (black) are at zero potential; their plasma-facing surfaces are G = E 0 . The cathode (red), the target (grey), and the target clamp (brown) are at negative potential; their plasma-facing surfaces are C ≡ E 1 . A further electrode E 2 (blue) is also on fixed potential. Insulators I (grey) separate the electrodes. The magnetic field lines are shown in black; the field strength is indicated by the background color and by the dashed lines (in mT). The domain V is split into an ummagnetized region Vum (green) and a magnetized region Vm. The part of the latter which consists of field lines from cathode to cathode forms the ionization region V ir (yellow). The complement is the outer region Vor = V/V ir . The boundary between the regions, i.e. the top of the ionization zone (yellow line), is called T ; the bottom of the ionization zone (part of the cathode) is B. figure 2 (and as used in many experimental publications cited above). We assume that the magnetron operates in a vacuum chamber-mathematically a simply connected spatial domain V-that shares its axisymmetry. This idealization is reasonable: the immediate neighborhood of an axisymmetric magnetron is usually axisymmetric as well; any infrastructure that is not (such as gas inlets and outlets) is distant and of little influence. The boundary ∂V consists of the grounded surface G = E 0 , the anode and the chamber wall, and the strongly biased cathode C = E 1 . Inserted into the chamber wall may be a number of axisymmetric electrodes E k , k = 2 . . . N E , to represent, for example, biased substrates [3]. All electrodes are separated by insulating segments I. The magnetic field emerges from an appropriate assembly of permanent dipole and ring magnets mounted below the cathode. It is axisymmetric and time independent. (A method to construct an analytical field model from Hall probe data is given by Krüger et al [52].) In natural cylindrical coordinates (r, θ, z) the magnetic field has the following form, where the last identity gives the representation by the flux function Ψ(r, z) and ensures that the field is divergence-free: The flux function equals, up to a factor of 2π, the magnetic flux through an axis-aligned circle of radius r and coordinate z, The vacuum character of the field imposes a constraint: It is useful to supplement the cylindrical coordinates (r, θ, z) by an alternative system of field-aligned flux coordinates (ψ, ϑ, s) [53]. They are created by amending the flux Ψ(r, z) with an arc length function S(r, z) for measuring distances along the magnetic field lines. S(r, z) can be constructed via integrals along the field lines, or alternatively characterized as a solution of the partial differential equation The azimuth θ remains unchanged but gets renamed: The back transformation enlists the functions R(ψ, s) and Z(ψ, s) to express the cylindrical coordinates in terms of the flux coordinates: The spatial volume element transforms as follows, which shows that the transformation is regular except at magnetic null points and on the center axis: An associated velocity transformation maps the Cartesian coordinates (v r , v θ , v z ) to field-aligned spherical coordinates (v, χ, φ), where χ is the pitch angle and φ the gyrophase: The back transformation is The velocity volume element transforms as Together, the variables (ψ, ϑ, s, v, χ, φ) are called the (position and velocity) gyrocoordinates. We will use the symbols r and v to denote particle positions and velocities in abstract form, to be translated into cylindrical or gyrocoordinates as appropriate.
The presence of a magnetic field B drastically affects the dynamics of the discharge: it prevents the electrons from moving freely through the plasma and forces them to gyrate in tight spirals around the magnetic field lines. The width of the spirals, less than 100 µm, is measured by the Larmor (gyro) radius r L . Motion across the field lines is thus restricted, more precisely impossible without higher order transport processes such as drift or diffusion. Motion along the field lines is free. In the course of this study, we will actually neglect all 'finite-Larmor-radius effects' and identify the position of an electron with its guiding center. The electron then moves along its field line, i.e. its coordinates ψ and ϑ are constant and its arc parameter s evolves with the parallel velocity, The equations of motion for the velocity gyrocoordinates are The exact validity conditions for these equations of motion will be discussed further below. As a consequence of these equations, the total energy ϵ tot and the magnetic moment µ of the electrons are constants of motion, Electrons are magnetized when their gyroradius is small compared to other relevant scales such as the gradient length L ∇B = B/|∇B| and the curvature length L κ = |κ| −1 of the magnetic field and the collisional mean free path λ. However, magnetization criteria based on a direct comparison of the length scales mentioned turn out to be rather unwieldy. For an algebraically simple criterion, we assume a minimum field value of, say, B min = 5 mT, and define the magnetized region as (see figure 2): Obviously, magnetization is a local property. For a global description of how the magnetic field restricts the movement of the electrons, we must take the field topology into account. Each magnetic field line segment which lies in V m is labeled by a flux ψ and an azimuth ϑ. Its endpoints have the arc coordinates s lo (ψ) and s hi (ψ); they either represent an intersection with the material chamber wall ∂V or a transition into the unmagnetized region V um = V \ V m . (In its further course, a field line may re-enter V m , which could be accounted for by numbering the segments by a discrete index n. For simplicity, we will omit that detail in the discussion.) End points lying on the cathode are referred to as active: high-energy secondary electrons can be injected there, and electrons that encounter them are reflected onto their field lines. End points lying on grounded areas, other electrodes, or insulating segments are inactive: they do not inject secondaries and reflect only thermal electrons, but not energetic ones. Finally, end points lying on a transition to V um are open: electrons can leave the domain V m and become free, perhaps returning later moving on a different field line. Only field line segments with two active end points can take part in the actual discharge. (Otherwise, at this low pressurefar to the left of the Paschen minimum-ignition is not possible at the typically applied cathode voltages.) All active segments together form the ionization region V ir ; its complement V or = V \ V ir is the outer region. The ionization region has the form of a half solid torus, see figure 2. Its boundary ∂V ir is formed by an annular section B of the cathode C and a flux surface T . Denoting the flux value of T as ψ lo and the maximal flux on B as ψ hi , the ionization region can be given in flux coordinates as Next, we analyze the temporal regimes of a typical high power magnetron discharge. The discharge is operated in argon at a pressure of p = 1 Pa and a temperature of T g = 800 K. The cathode radius is R C = 25 mm but a more typical scale is L = 10 mm which reflects the field line length in the ionization region and the size of the observed structures [23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42]. We select two sets of parameters, one for the ionization region, the other for the outer region. For a point in  Altogether, our analysis indicates that there are four wellseparated groups of time scales (or temporal regimes) present in the discharge dynamics: • Both regions exhibit a very fast regime-around 10 −11 sgiven by the inverse values of the electron plasma frequency ω pe . The corresponding dynamical phenomena are collective electron oscillations and waves. We assume that they are quiescent, i.e. that the electron fluid is in force equilibrium with the electric field. This is supported by experimental evidence on the spectrum of fluctuations in high power magnetron plasmas [23,33]. • The second, also very fast regime around 10 −10 s is only present in the magnetized region: the inverse of the gyrofrequency Ω marks the individual gyromotion of the electrons. It is much faster than the collision time which reflects the magnetization of the electrons. We assume that there is no collective motion on that scale, again following [23,33]. • The third regime ranges from 10 −9 s to 10 −8 s (in the zone V ir ) or 10 −7 s (in the zone V or ), identified as the thermodynamic time scale. It contains the inverse collision frequency ν −1 c and the internal energy relaxation time τ eq , both dominated by the Coulomb interaction, and also the time τ diff that an electron needs to diffusively traverse the respective region. In V ir this only applies to the motion along the field, outside the motion is not restricted. Also, for both regions, the inverse of the ion plasma frequency ω pi belongs to this regime. The identity ω pi = c s /λ D , however, implies that dynamics on the time scale of ω −1 pi would require spatial structures on the scale of the Debye length λ D . Such small scales are only present in the boundary sheath, while the length scales of the plasma are given by L. Since the structure of the sheath is controlled by the ion flux which passes through it, and since this flux can vary only slowly, we assume that all ion motion takes place on the time scale τ i = L/c s or slower, again supported by [23,33]. The capacitive time constant of the discharge chamber is τ cap = 2.6 ns, according to (A2) evaluated with a cathode radius of R C = 25 mm and an inner RF resistance of the voltage source R i = 10 Ω.
• The fourth and slowest time scale, lasting some 10 −6 s or longer, is the dissipative regime. This time frame encompasses all phenomena related to ion transport, particle generation and recombination, energy gain and loss processes, and surface charging. Central to this regime is the energy confinement time τ ec , which measures the duration for electrons to lose kinetic energy to ions or neutral gas. The cross Altogether, the time scale analysis reveals that a HiPIMS discharge is in many aspects similar to the unmagnetized inductively coupled plasma (ICP) studied earlier [45]. There is a pronounced separation between a fast thermodynamic and a slower dissipative time scale, which invites to generalize the quasi-thermodynamic approach of that work. However, there are also clear differences: the magnetized region contains an additional fast time scale, namely that of the gyromotion. While we may assume that it is not collectively excited, it surely will influence the individual electron motion and must be considered. Directly related is the fact that, in the active region, only the transport along the field lines is fast, while the cross field transport counts as slow. We will consider the implications of these facts in the next sections.

Setup of the model
To develop a model for the thermodynamic regime, we start from a kinetic equation that describes the evolution of the thermal electron distribution f = f(r, v, t) on the fast scales.
Adopting the electrostatic approximation E = −∇Φ, it reads in standard notation [54]: The collision operator on the right side contains only electronelectron Coulomb collisions, electron-ion Coulomb collisions in the limit of immobile ions, and elastic electron-neutral collisions in the limit of immobile neutrals. At the reactor wall, specular reflection is assumed; particle loss through the plasma boundary sheath is a slow dissipative effect and neglected.
Here, n is the normal of ∂V and U − denotes the set of reflected velocities with v · n < 0, The kinetic equation couples to the Poisson equation which calculates the electric potential Φ from the charge density ρ resulting from the net ion density n i and the electron density n e . (The charge density of the energetic electrons is viewed as a negligible correction term. Their number, controlled by the ions which impinge the cathode, is much smaller than that of the thermal electrons, and their distribution is not greatly affected by the plasma potential.) The net ion density is not necessarily axisymmetric, but it is assumed to be stationary on the fast time scales we are considering. We assume a set of mixed boundary conditions which account for ground G ≡ E 0 (the grounded section of the chamber wall and the anode cover), the negatively biased cathode C ≡ E 1 , and possibly other electrodes E k with k = 2 . . . N E . The electrodes are separated by insulating areas I which possibly carry surface charges σ. Altogether, the Poisson equation and its boundary conditions read It is important to note that our current approach differs from our previous one of [45] in that we no longer assume the voltage sources to be ideal. Instead, we take into account the internal resistance R k of each source, which causes a difference between the electrode voltage V k and the no-load voltageV k : The electrode currents I k are considered positive when they flow into the discharge chamber. Due to the assumption of specular reflection, only the displacement current contributes: To write the model in dimensionless form, we employ a normalization that reflects the status in the ionization region. We use L,T e ,n, andB as typical values for the length scale, the electron temperature, the plasma density, and the magnetic field strength, respectively. Useful derived scales are the thermal speedv = T e /m e 1/2 , the gyrofrequencyΩ = eB/m e , and the Debye lengthλ D = ε 0Te /e 2n 1/2 . We set . Furthermore, we define two dimensionless numbers that indicate the ratio of the frequency unit to the gyrofrequency and the ratio of the Debye length to the length unit: Dropping the tilde immediately, the kinetic equation reads and Poisson's equation assumes the form

Formal solution of Poisson's equation
The formal solution of Poisson's equation with mixed boundary conditions is given in [45]. Converted into dimensionless form, it reads It employs the Green's function G(r, r ′ ), defined as the solution of Poisson's equation for a unit point charge at r ′ ∈ V under homogeneous boundary conditions [55]: The Green's function is symmetric with respect to the exchange of it arguments, The electric potential Φ is usefully split into two contributions, the inner potential Φ (i) and the outer potential Φ (o) . The inner potential arises from the charge density distribution ρ in the domain V under homogeneous boundary conditions. As the charge density ρ itself, it can vary on the thermodynamic time scale: Its formal representation of the inner potential in terms of the Greens function is the first term of the sum given in (39): The outer potential arises from the voltages applied to the electrodes and from the surface charges on the dielectric surfaces under the assumption that the charge density vanishes: It is represented by the second and third term of (39). Because of the dependence of the V k on the electrode currents I k , it can also vary on the thermodynamic time scale: We also define the load-free outer potentialΦ (o) which depends on theV k and on σ and is thus constant on the thermodynamic time scale: For certain analyses, it is useful to further divide the outer potential. We introduce the influence function Ψ l for an electrode l, which is defined as the solution of and has the formal representation Finally, we define the surface-charge induced potential as the solution of: Condition (31) implies that the surface charges are constant on the thermodynamic timescale. Therefore, the corresponding potential is also constant, it is represented as Collecting the expressions, we can write the full formal solution of the potential as and the corresponding electric field as For later use, we collect a set of additional results. The inner electric field E (i) = −∇Φ (i) is orthogonal to the gradient of the influence functions, as the last integral vanishes:ˆV The same applies to the surface charge field E (s) = −∇Φ (s) : The product of the inner field E (i) and the surface charge field E (s) does not vanish, A corollary of these results is the identitŷ Finally, we define the capacitance coefficients of the electrodes by the following symmetric and positive definite matrix [56]:

Kinetic equation for the thermodynamic time scale
In the magnetized domain V m the electron dynamics contains not only the thermodynamic time scale but also the much faster scale of gyromotion. We will eliminate this scale with the help of a two-time-scale formalism. We first write the kinetic equation in gyrocoordinates, where the distribution function is f = f(ψ, ϑ, s, v, χ, φ, t) and the potential Φ = Φ(ψ, ϑ, s, t).
The magnetic field B, the cylinder coordinate functions R and Z, the field line curvature κ, and the gauge connection Γ all depend only on ψ and s: Using abbreviations of obvious meaning, the kinetic equation reads Clearly, the structure of this kinetic equation invites to apply a two-time-scale formalism. We define the fast gyro time scale t −1 = η −1 t and the slower thermodynamic time scale t 0 = t. (This is a different distinction than between the thermodynamic and dissipative regimes.) All dynamical functions depend on both scales; the time derivative is split accordingly: For the distribution function and the potential we specify a series in η: Sorting equation (59) by powers of η leads to two equations. The first of them is The general solution of (63) is any periodic function of φ − Bt −1 . However, we adopt the assumption that f (0) is gyrophase invariant; it is then also not a function of the fast time t −1 . Neither is, per Poisson's equation, the potential Φ (0) . The second equation then reads This equation can be solved for f (1) in the time scale t −1 . The explicit form of that solution is not of interest here, except for the fact that it possibly contains a term that is linear in t −1 . For the two-time-scale expansion to be consistent, such a linear divergence must be avoided. This poses a constraint on the terms that are constant with respect to the gyrophase: Restoring the full expressions and omitting the now unnecessary indices, we obtain a kinetic equation for the distribution function f = f(ψ, ϑ, s, v, χ, t) on the thermodynamic time scale: The kinetic equation requires boundary conditions at the endpoints s lo (ψ) and s hi (ψ). Physically, we continue to assume specular reflection, but it is not possible to harmonize this assumption with the condition of gyroinvariance. This insight was the result of recent research of our group [57][58][59]. Appendix B briefly reviews the material needed for this work. In a nutshell: specular reflection at the boundary destroys the gyrophase invariance of the incoming distribution function (unless it is isotropic, i.e. independent of the pitch angle). The outgoing distribution function shows complex, even fractal structures in the gyrophase. However, these structures are dissolved after a short distance by the process of phase mixing. Altogether, an effective boundary condition can be defined that maps incoming gyroinvariant distributions directly into outgoing gyroinvariant distributions:

Thermodynamically relevant local balance equations
For our later analysis, we also need to consider the balances of the thermodynamically relevant velocity moments of the distribution function f. These are the electron density n, the particle flux density Γ n , the energy density w, and the energy flux density Γ w , According to our assumptions above, the corresponding collision integrals are equal to zero. Consequently, the balances of particle number and energy read The energy balance of the electric field can be formulated as ∂ ∂t Adding the two energy balances results in ∂ ∂t A further important balance is that of the entropy S. In unitless notation, its phase space density is −f(ln(∆ µ f ) − 1), where ∆ µ is a numerical constant that is of no further influence. (Choosing ∆ µ = (h/m ev L) 3 would yield the semiclassical Sackur-Tetrode-equation for the entropy of a monoatomic gas [60].) We follow the convention to set ∆ µ equal to unity [54]. We thus define the local entropy density as and the corresponding flux density as The production density of entropy is a non-negative quantity. This is what the second law of thermodynamics states, but it can also be shown explicitly (as Boltzmann H-theorem [54]) for all collision terms that appear in the kinetic equation (37): This can be combined to the local entropy balance equation In the ionization region we work in gyrocoordinates. Integrating the distribution function over the velocity v, the pitch angle χ, and taking the factor 2π for the gyrophase into account, we obtain the electron density n, the flux density Γ n in the direction of the magnetic field, the energy density, and the energy flux density in the field direction, respectively: Again, the collisional production of the particle number and the energy vanishes identically, and we obtain the balances of particle number and energy as ∂ ∂t ∂ ∂t In a similar way, we obtain the entropy density and the entropy flux density along the magnetic field lines, respectively: Integrating the kinetic equation (65) with the weight ln f over velocity space yields the entropy balance equation for an infinitesimal flux tube, (87)

Balance equations of the thermodynamic regime
In order to further characterize the thermodynamic regime, we will derive a set of three balance equations that relate to the electron particle number N, the electron entropy S, and the electromagnetic enthalpy H. The first two balances will be considered globally in the region V or and field line (flux tube) resolved in the region V ir . The last balance will be given globally for the discharge chamber V. The particle numbers are conserved quantities, the entropies follow the second law of thermodynamics and are therefore non-decreasing, and the enthalpy is non-increasing.

Electron number balance
We integrate the continuity equation (71) over the region V or . Applying Gauß' theorem, the second term becomes a surface integral over ∂V or . This surface consists of segments of the chamber wall ∂V which are specularly reflecting, and of the flux surface T (the top of the ionization region) through which particle flux is prohibited due to the magnetization of V ir . The surface integral thus vanishes and we can state that, on the thermodynamic time scale, the particle number N or in the outer region, is a conserved quantity, so actually not a function of t: In the ionization region, we integrate the particle number balance (83) over the field line. The end points s lo (ψ) and s hi (ψ) lie on the cathode where boundary condition (66) holds. We find that the electron number per infinitesimal flux tube, is conserved on the thermodynamic time scale, so actually only a function of ψ and ϑ:

Entropy balance
By integrating the local entropy balance over V or we obtain the balance equation for S or , the entropy of the outer domain defined as The appearing surface term is non-negative. At the chamber wall, the entropy flux vanishes due to the reflecting boundary conditions. At the flux surface T this cannot be assumed. However, as the particle flux vanishes, we can simply postulate (on grounds of the second law) that the net entropy flux into the region V or is non-negative. Altogether, the total entropy of the outer region is a non-decreasing quantity: For the ionization region, we again use gyrocoordinates. The corresponding balance law is formulated for the entropy per flux tube, which is related to the total entropy in the ionization region via Integrating the entropy balance over the arc length parameter s gives The second term involves the negative entropy flux which has to be evaluated at the field line endpoints s lo (ψ) and s hi (ψ). Here, we again refer to our previous publications [57][58][59], where we have studied the interaction of a gyrating charge with a boundary sheath in detail. Appendix B gives a short overview. The bottom line is that an effective boundary condition can be defined which is dissipative in the sense that the outgoing entropy flux is larger than the incoming one. The second term in equation (97) is thus always non-positive. Altogether, it can be stated that the entropy of an electron distribution on a flux tube in the ionization region is an increasing quantity on the thermodynamic time scale:

Electromagnetic enthalpy balance
Finally, we discuss the concept of electromagnetic enthalpy H. This quantity measures the total energy of the discharge, including the plasma and electric field, adjusted for the energy exchange with external voltage sources. The concept was introduced in our earlier work [45]. However, as already stated, we have come to the opinion that the assumption of ideal voltage sources made there was an oversimplification; real voltage sources have internal resistance. To account for this we will modify the approach and now incorporate the effect of internal resistance in our analysis. We begin by integrating the total energy balance (74) over the discharge volume V and applying the Gauß theorem. The reflection boundary condition (31) implies that the matter flux to the wall vanishes, leaving only the displacement current: Since the displacement current cannot pass through the insulating sections I of the wall, the integral is limited to the electrodes E k , where the potential Φ is equal to V k : d dtˆV The surface integral terms in the summation correspond to the negative electrode currents. By expressing the electrode voltages in terms of the load-free voltagesV k and the resistive voltage drop in the source, V k =V k − R k I k , we obtain d dtˆV Employing the zero-load outer potentialΦ (o) as defined above, we get d dtˆV AsΦ (o) is time-invariant in the dynamics under study, the time derivative in the second term can be drawn in front of the integral. This defines a temporally non-increasing quantity, termed the preliminary (indicated by the prime) electromagnetic enthalpy, as The electrical field can be written E = E (i) − ∇Φ (o) . Utilizing equation (56), we can express the preliminary enthalpy as The first term of the integral, energy of the inner field, can be written aŝ Utilizing the capacitance coefficients introduced above, the second term becomeŝ The third term, finally, is dropped as temporally constant.
After some re-arranging of terms, we define the electromagnetic enthalpy in its final form as By construction, the electromagnetic enthalpy is a nonincreasing quantity, The decrease of the enthalpy is a phenomenon that takes place in the thermodynamic regime; it is controlled by capacitive time constant τ cap of the discharge.

Variational principle
In the previous section, we analyzed the balances of electron number, entropy, and electromagnetic enthalpy in the thermodynamic regime. We found that the infinitesimal magnetic flux tubes of the ionization region V ir , denoted by labels ψ and ϑ, behave as adiabatically insulated thermodynamic units. The outer region V or forms a thermodynamic unit by itself. In each of these adiabatically insulated units, the particle number is conserved and the second law of thermodynamics applies, meaning that the entropy is always non-decreasing. Additionally, the system enthalpy is non-increasing. This description bears resemblance to the fundamental scenario of classical thermodynamics. Applying Lyapunov's second method, we should be able to conclude that discharge states characterized by extremal values of the negative entropies and the enthalpy represent stable equilibria.
As a result of this reasoning, we encounter a non-standard optimization problem with an infinite number of objectives and an infinite number of constraints. Specifically, we seek to maximize the spatially distributed entropies and to minimize the enthalpy, while satisfying the requirement to keep the particle number quantities constant. However, the problem can be mapped onto a standard optimization problem which involves only one single objective. This is achieved by introducing a positive but otherwise arbitrary weight function T(r) which is constant in V or but depends on the magnetic field line labels in V ir , We then define a combined functional, denoted by the symbol G because of an obvious formal analogy to the Gibbs free energy of classical thermodynamics: This functional is also a temporally non-increasing quantity: The minima of the functional G with respect to a variation of f under the particle number constraints (89) and (91) represent stable equilibria of the thermodynamic regime. (The double-sum term over the I k is then automatically minimal, namely equal to zero.) These minima could easily be constructed directly. However, it is advantageous to carry out the minimization in a step-wise manner. In a first step, we thus minimize the functional G under the stronger condition that the electron density at each point is equal to a prescribed (and aptly named) positive quantity: The charge density ρ is then fixed, and it suffices to consider only certain elements of G. We introduce a positive Lagrangian multiplier function λ n (r) and write the corresponding unconstrained functional as The first variation is It is identically zero under the condition Solving for f and choosing the Lagrange parameter λ n (r) to meet the constraint (112) results in a Maxwellian with the prescribed temperature T(r), where the outer region is characterized by a global temperature T or , while in the ionization region each infinitesimal flux tube has its own temperature T(ψ, ϑ): The second variation is positive definite which shows that the Maxwellian f is indeed a unique minimum of the functional L: We are now left with a purely spatial problem, for which the functional G reads This target functional needs to be minimized under the number constraints (89) and (91). We note that certain parts of G are fixed, and define the streamlined version We define a Lagrangian multiplier function, constant in the outer region V or but field line dependent in the ionization zone V ir : The corresponding unconstrained functional is then +ˆV T(r) n(ln n − 1) d 3 r +ˆV λ(r) n d 3 r. (121) Using the symmetry of the Green's function, its first variation reads: It is identically zero under the following condition, where we have re-introduced the inner potential Φ (i) via (43) and set Φ = Solving for n and renaming the function λ(r), we obtain the Boltzmann relation The density prefactorsn can also be expressed in terms of the particle number quantities N. For the outer region, we obtain and for the ionization zonê Finally, we mention that the second variation is positive definite, the described solution is indeed a minimum of the functional L ′ : Combining the results, we see that the equilibrium distribution function is a Maxwellian that obeys the Boltzmann relation. The outer region V or is characterized by a constant temperature and a constant density prefactor, while in the ionization region V ir each infinitesimal flux tube has its own temperature and prefactor: The obtained distribution function and potential are indeed solutions of the kinetic model. The left side of the kinetic equation vanishes identically for any function of the total energy. In the ionization zone, derivatives are taken along the magnetic field, i.e. with respect to the arc parameter s, so that additional dependence on the field line labels ψ and ϑ is allowed. The collision term on the right vanishes for any Maxwellian. The boundary conditions are met by functions that are isotropic. Finally, also the Poisson equation is fulfilled.

Summary and outlook
In this study, we employed a quasi-thermodynamic approach to investigate the dynamics of high power magnetron discharges. Our goal was to demonstrate that these devices can be described by a modified PB equation, even when being far from thermodynamic equilibrium. Our research yielded several key findings, which we will review in this section.
To begin, we analyzed the discharge geometry and divided its domain V into two regions, the ionization region V ir and the outer region V or . The ionization region is characterized by the presence of field lines, identified by their flux ψ and azimuth ϑ, that are magnetized throughout and have both endpoints on the cathode. The outer region, as the complement of the ionization region, encompasses all points that do not lie on those active field lines. Our analysis of the dynamics showed that the time scales involved are more complicated than those observed in unmagnetized high density plasmas [45]. An additional fast time scale, namely the electron gyration time Ω −1 ∼ 0.1 ns, became present in the magnetized regions. While this could be addressed using a two-time-scale method, the magnetization also hinders the transport of electrons across active magnetic field lines, resulting in adiabatic insulation of the infinitesimal flux tubes of the ionization region. After some discussion we again found two well-separated temporal regimes, a fast (10 ns or faster) 'thermodynamic' and a slow (1 µs or slower) 'dissipative' regime. The fast regime covers electron waves and oscillations, electron-electron Coulomb interactions, and elastic interactions with ions and neutrals in the limit of immobile scattering centers. The slow regime represents inelastic collisions, generation and recombination of particles, crossfield drift/diffusion in the ionization region, electron losses to the chamber wall, and the ion dynamics. (The density of the minority of energetic electrons also changes slowly: they are produced through ion-induced secondary emission at the cathode with a rate that depends on the flux of impinging ions, and have lifetimes of a few microseconds. Moreover, as they account for only a fraction of the overall charge density and possess energies significantly higher than the energy fluctuations induced by the potential, their interaction with the potential is weak.) The result of our considerations was a consistent characterization of the thermodynamic regime of a high power magnetron. The electron number and entropy balances needed to be formulated separately for the flux tubes of the ionization region and for the outer region. The electromagnetic enthalpy balance was formulated for the entire chamber.
Combining the balance equations, we could set up a thermodynamic variational principle. A functional of the distribution function f was established, the minima of which fulfilled the kinetic equation and the boundary conditions. In physical notations they read where the electron density is in modified Boltzmann equilibrium: The density prefactorsn are related to the particle numbers. In the outer region we havê and in the ionization region, The electric potentialΦ fulfills the Poisson equation and the load-free boundary conditions: Together, (129)-(133) state that the electric potential in the ionization zone of a high power magnetron discharge obeys a modified PB equation and that the thermal electrons follow a Maxwell distribution with an electron temperature T e and a density prefactorn that can vary from field line to field line. In the outer region, a conventional PB equation applies and the temperature and the density prefactor are uniform.
The conclusions of this study were based on the analysis of a mathematical model using a number of length and time scale arguments. As always, it is important to verify the validity of the assumptions made and see whether the conclusions can be supported by other means. The following observations on high-power magnetrons have been reported: • The localized ionization zone (spokes) that develop in high power magnetron discharges are of macroscopic size [23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42]. They measure typically a centimeter, which is comparable to the length and curvature radius of a magnetic field line, as well as the mean free path. The spoke size is much larger than the gyroradius and Debye length. At the same time, their dynamics takes place on the scale of microseconds, not nanoseconds. This is evidence that the length and time scale assumptions we made in this study are correct. • Langmuir probe and laser Thomson scattering measurements demonstrate the Maxwellian distribution of the electrons in high power magnetron discharges [16, 18-22, 25, 37, 47-51]. Numerical simulations with particle-incell/Monte Carlo collisions codes [19] and with direct solutions of global kinetic models point in the same direction [16,20]. (Experimental detection of energetic electrons is challenging due to their low density. Numerical simulations face no such difficulties. They suggest that the assumption of a Maxwellian electron energy distribution function (EEDF) holds true for electrons with energies up to several tens of eV.) • Recent measurements by Dubois et al determined the electron drift velocity to be some percents of the thermal velocity [61]. This can be interpreted as an indication that the distribution function is essentially isotropic, and deviations only occur in a higher-order perturbation theory in η (which is not addressed here). • Assuming that the electron density n e is not drastically different from the ion density n i , Boltzmann equilibrium, equation (130), implies that the potential drop along the magnetic field lines is at most some T e /e. In contrast, the voltage drop across the field lines can be substantial. This was reported by many studies [62][63][64]. Importantly, it seems to hold not only for the time-averaged potential, but also for the time resolved potential which can be very dynamic in connection with the spoke phenomenon [34]. • The light emission from the observed spokes can show very strong gradients perpendicular to the field lines, but not along them [31]. This agrees with our result that the electron temperature is constant along the field lines. Of course the light emission can also be influenced by the energetic electrons, but then similar arguments apply.
Considering the experimental and theoretical evidence presented, we confidently conclude that our approach is solidly supported. Nevertheless, some issues require further attention. Furthermore, we intend to explore potential applications of our findings. The following list provides an overview of our planned future activities: • Our initial kinetic model, equations (37) and (38), has two smallness parameters, η and ε. The first number represents the ratio of the frequency unit L/v the to the gyrofrequency Ω, the second stands for the ratio of the Debye length λ D to the macroscopic scale length L. Our investigation has utilized the smallness of η, but the smallness of ε was not exploited. What could we conclude from postulating that λ D L? For unmagnetized discharges, the answer is clear: The discharge self-organizes into a quasi-neutral plasma on the scale L and an electron-depleted sheath on the scale λ D . The voltage differences in the plasma are at most a few T e /e, while the sheath voltages are large. For magnetized discharges such as high power magnetrons, the situation is less clear. There is an electron-depleted sheath, but is there a (uniformly) quasi-neutral plasma? The literature on this point is divided. Liebig and Bradley found quasineutrality experimentally confirmed to a high degree [62]. In contrast, Anders et al stated that 'Future work will require replacing the quasineutrality assumption with a solution of the Poisson equation' [31]. A study currently in preparation will offer a definitive answer to this open question. • Our modified Boltzmann-Poisson equation fixes the algebraic form of the equilibrium. However, specific values of the form functions T e andn are not provided by the theory. In fact, any positive value of the electron temperature and the density prefactor would define a suitable equilibrium. Which mechanism determines them in physical reality? To address this issue, we will have to go beyond the thermodynamic regime and analyze the dissipative regime by suitable methods. • Finally, we strive to link our research with the field of high-temperature plasma physics. In a recent publication, Negulescu studied fully ionized magnetized plasmas with mass disparate particles (electrons and ions) and rederived the electron Boltzmann relation [65]. How do our results match hers? At a first glance, there is some overlap, but differences concern the geometry and the boundary conditions. It will be one task of our future work to examine these connections more closely.

Data availability statement
No new data were created or analyzed in this study.
For a cathode with an area of A C = πR 2 C , the cathode sheath capacitance is C C = ε 0 A C /s. Assuming an inner source resistance R S , the corresponding capacitive time constant is The electrons of the discharge interact with each other and with the massive particles. With lnΛ = ln 12πn e λ 3 D as Coulomb logarithm, the Maxwellian-averaged electron-electron and electron-ion collision frequencies are For the electron-neutral interactions we rely on [68]. The rate constant for elastic scattering can be approximated as (A6) The motion of the electrons (if they are magnetized, along the magnetic field) follows ∂j ∂t = e 2 n i m e E − (K el n n + ν ei ) j.
To assess the thermal relaxation within the electron fluid, we assume an initial anisotropy with T e∥ = T e + ∆T e and T e⊥ = T e − ∆T e /2. The evolution of ∆T e is then given by [69] ∂∆T e ∂t = − 3 5 √ 2ν ee + 2ν ee ∆T e , Finally, we consider the energy balance of the electrons, with P as the heating power: ∂ ∂t 3 2 n e T e = P− K in E in + K iz E iz + 3 m e m n K el n n T e +3 m e m i ν ei n i T e n e .
Based on these considerations, we define the slow-down time τ sl , the equilibration time τ eq , and the energy confinement time τ ec as τ sl = (K el n n + ν ei ) −1 , (A10) The thermal electrons are heated via several channels. One channel results from the interaction with the energetic secondaries which start at an energy eV sh and are cooled down. If the energy ranges are well separated, the power is as follows, where the final expression approximates the energetic electron distribution function as a constant: (A14) In the ionization region V ir , the transport of the electrons along the field lines is diffusive; we calculate the diffusion constant and the diffusion time as follows; the same expressions also apply for the diffusion in the outer region V or : The diffusional transport across the field lines is highly suppressed: The drift velocity and the corresponding time constant are estimated via

Appendix B. Boundary conditions at closed field line end points
The boundary conditions (66) for the kinetic equation (65) at the end points s lo and s hi are based on insight into the interaction of magnetized electrons with a boundary sheath. Our group investigated this interaction in a series of studies [57][58][59]. We adopted the scaling λ D r L λ ≈ L as in this manuscript, but focused on the scale of the gyroradius r L . The Debye length λ D and thus the sheath thickness were set to zero. The mean free path λ and the characteristic system length L were set to infinity, i.e. the plasma was assumed to be collisionless and free of an electrical field and the magnetic field to be homogeneous. The inclination of the field to the sheath normal n was called γ. The trajectories of the electrons became regular helices, and their interaction with the sheath specular reflections. Below we summarize some results which are relevant for the current study. Without loss of generality we consider a left endpoint s lo , its coordinate shifted to s = 0. An electron is assumed incident from the right, having the velocity v, the pitch angle χ in , and the gyrophase φ in at the reference position s = 0. Being incoming means cos χ in < 0; the pair (χ in , φ in ) belongs to the negative unit hemisphere S − . After possibly several specular reflections at the sheath the electron becomes outgoing. The velocity is still v, but the pitch angle has changed to χ out with cos χ out > 0, and the gyrophase to φ out at the position s = 0. (Actually, the electron returns on a field line shifted by some r L , but this finite-Larmor-radius effect only needs to be taken into account when going beyond the expansion order used in this investigation.) The pair (χ out , φ out ) belongs to the positive unit hemisphere S + . An effective interaction S directly maps the incoming onto the outgoing angles: S : S − → S + , (χ in , φ in ) → (χ out , φ out ) = (χ γ (χ in , φ in ), φ γ (χ in , φ in )). (B1) The functions χ γ (χ, φ) and φ γ (χ, φ) must be constructed numerically. (This reflects the fact that the intersections of a helix with a plane cannot be given in a closed analytical form.) The construction was carried out in [58], in form of a parameter study for all values of γ. The functions turned out to be highly complex and even fractal in nature. However, they are differentiable except on a set of measure zero, and where they are, the Jacobian fulfils ∂χ γ ∂χ ∂φ γ ∂χ ∂χ γ ∂φ ∂φ γ ∂φ = sin χ cos χ sin χ γ cos χ γ = − sin χ cos χ sin χ γ cos χ γ . (B2) With these insights, the specular reflection condition (31) can be cast in gyrocoordinates. The outgoing distribution function f (out) (χ, φ), for (χ, φ) ∈ S + , results from the incoming