Collisionless dissipation at the boundary sheath of magnetized low temperature plasmas

In partially magnetized technological plasmas like magnetron discharges or Hall-effect thrusters, the interaction of gyrating electrons with the plasma boundary sheath plays an important role. This study shows that the interaction can be described as a dissipative process. It is assumed that the Debye length λD is much smaller than the electron gyroradius rL , and rL in turn much smaller than the mean free path λ and the field gradient length l. Focusing on the scale of the gyroradius, the sheath is assumed as thin ( λD→0 ), collisions are neglected ( λ→∞ ), the magnetic field is taken as homogeneous, and electric fields (= potential gradients) in the bulk are neglected ( l→∞ ). The interaction of an electron with the electric field of the plasma boundary sheath is represented by specular reflection, v→v−2v⋅nn , where v is the electron velocity and n is the sheath normal. A previous study showed that this set-up leads to a scattering process where gyrophase invariant incoming velocity distributions are scattered into outgoing distributions with fractal dependence on the gyrophase (Krüger and Brinkmann 2017 Plasma Sources Sci. Technol. 26 115009). The present research focuses on the subsequent evolution. It is shown that the gyromotion results in fast phase mixing with respect to the gyrophase which quickly restores the gyrophase invariance. In combination, the two phases of the sheath interaction constitute a mapping within the space of gyrophase invariant distributions. This combined process is dissipative and leads to isotropization. As an explicit comparison shows, it acts similarly to elastic electron-neutral collisions and can even dominate the latter in low pressure plasmas. A simplified version of the operator that applies for weakly anisotropic distributions is suited for practical simulation work.


Introduction
Magnetized low temperature plasmas such as Hall effect thrusters, DC or RF magnetrons, and magnetically enhanced glow discharges are widely used in technology and industry [1][2][3][4][5][6]. They operate at low pressures, p = 0.1-1 Pa, and high plasma densities, n e = 10 18 -10 20 m −3 . The magnetic field, with flux densities of typically B = 10-100 mT, serves to confine the electrons by forcing them on narrow helices around the magnetic field lines (gyromotion) [7]. The gyroradius (Larmor radius) r L , for thermal electrons around 10 −4 m, is much smaller than the field gradient length l and the mean free path λ which both measure a few cm. The ion gyroradius, in contrast, is comparable to l and λ.
A suitable modeling framework is kinetic theory: Each species s (of mass m s and charge q s ) is represented by a distribution function f s . The evolution is governed by kinetic equations that take into account the influence of spatial gradients and the macroscopic fields E and B as well as the collisions with other particles and the walls. The electromagnetic fields, in turn, must be calculated by Maxwell's equations. Typically, the magnetic field is external and the electrostatic approximation E = −∇Φ applies; it suffices then to use Poisson's equation for the electrostatic potential Φ. Complemented by suitable initial and boundary conditions, kinetic models can describe virtually all magnetized technological plasmas [8,9].
Practical simulation, however, is a different matter: As sixdimensional systems of non-linear differential equations, kinetic descriptions are quite demanding, analytically and numerically. (The popular Particle-in-Cell/Monte Carlo collisions algorithm faces its own set of challenges for high-density discharges [10][11][12][13].) For fast and robust numerical models, physically justified simplifications are needed. With regard to the electron component, it would be an attractive strategy to use the characteristic scale relations discussed above. By considering the ratio of the gyroradius r L to the scales l and λ as a small number, asymptotic perturbation analysis may allow to separate the 'microscopic' gyromotion of the magnetized electrons from their 'macroscopic' drift motion. If it is possible to effectively integrate out the microscopic scales, all that remains to be considered is the macroscopic gyroinvariant dynamics on longer length and slower time scales, either in analytical models or in numerical simulations that are much less demanding. The drift-kinetic and gyrokinetic theories of high temperature plasma physics and plasma astrophysics demonstrate the viability of such a strategy and can serve as inspiration [14,15].
Unfortunately, magnetized technological plasmas differ considerably from their astrophysical or fusion cousins and a direct transfer of theories established in these areas is not possible. An important difference concerns the topology of the magnetic field: Fusion devices such as tokamaks and stellarators ensure long-time plasma confinement with magnetic field lines that wind around the apparatus many times and are essentially (or even literally) infinite [16]. Astrophysical plasmas typically have less regular magnetic field topologies, but also there the field lines often have no defined endpoints [17]. In technological plasmas, on the other hand, the field lines have well-defined endpoints at which they enter or leave the plasma domain. Electrons which approach these endpoints on their path interact with the boundary sheath, a thin zone with a strong electric field. The extension of the sheath boundary is governed by the Debye length λ D , typically λ D = 10 −6 -10 −5 m, much smaller than the gyroradius r L . The associated potential drop exceeds by far the typical electron energy (per electron charge). The magnetization condition fails, and the described simplifications may fail as well.
Our previous research has confirmed this suspicion. In two studies, we discussed the interaction of an incoming stream of electrons with a gyroinvariant distribution with a specularly reflecting boundary sheath [18] and a more physical, spatially extended Bohm sheath [19]. Both studies came to the same conclusion: An incoming gyroinvariant distribution is mapped onto an outgoing distribution that has pronounced structures in the gyrophase.
Does this mean that the envisioned program of complexity reduction cannot be carried out? In this study we will argue that this fear is not true. On the contrary, we will demonstrate the existence of a mechanism that restores the assumptions of gyrophase invariance dynamically. It will be shown that this mechanism-phase mixing [20]-has an important impact on the dynamics of magnetized technological plasmas, not only in terms of their modeling.
The rest of the work is organized as follows. Section II, a condensed review of [18], describes the set-up under study and casts the sheath-interaction of electrons in a one-particle picture. Section III switches to kinetic theory. First the sheath scattering operator S γ is recalled. We then describe the evolution of the outgoing electron distribution during the subsequent gyromotion by a transport operator T s , with s as the distance from the sheath along the field. For large values of s, T s essentially acts as gyroaverage · . The averaged operator S γ = S γ is proposed as an effective scattering operator. It is flux-conserving but dissipative and acts similarly to elastic collisions, as shown in section 5. Section 6 provides a practical approximation of the scattering operator, section 7 gives a short summary and outlook.

Configuration and one-particle picture
As stated, the magnetized plasmas under study are characterized by a Debye length λ D that is much smaller than the electron gyration radius r L , which in turn is much smaller than the gradient and curvature length l of the electric and magnetic fields and the mean free path λ. This research focuses on phenomena which take place at and in the immediate neighborhood of the boundary sheath of such plasmas. The size of the neighborhood is some ten gyroradii, typically a few millimeters. Assuming that this is still much smaller than the scales l and λ, electric fields (potential gradients), magnetic field gradients, and collisions can be neglected. The sheath is infinitesimally thin and located at the plane z = 0 of a suitably oriented system of Cartesian coordinates (x, y, z). The plasma occupies the half-space R 3 > = {(x, y, z)|z > 0}. The angle between the normal e z and the field direction b is γ; Figure 1. Geometry of the configuration under study. A homogeneous magnetic field is assumed, described by the gyrofrequency vector Ω. Its inclination with respect to the surface normal n ≡ ez is given by the angle γ.
we set b = sinγ e x + cosγ e z and define t 1 = cosγ e x − sinγ e z and t 2 = e y . Note that (t 1 , t 2 , b) is a right-handed system. For γ = 0 it coincides with (e x , e y , e z ). Sheath-parallel fields are excluded, |γ| < π/2 (figure 1). With Ω = eB/m e as the electron gyrofrequency, the equations of motion are: Under the stated assumptions, the interaction of an electron with the sheath at z = 0 can be described as a specular reflection that maps an incoming (−) to an outgoing state (+).
The one-particle picture which results from these equations of motion is shown in figure 2. An incoming electron approaches from infinity. Its trajectory is a helix around a field line with foot points X and Y. We introduce the guiding center R = r + b×v/Ω which moves on the field line and the gyration vector ρ = −b×v/Ω which points from there to the particle. Reflection at the boundary sheath scatters the electron onto a new helix, possibly repeatedly, until the electron reaches an outgoing helix and returns to infinity. Thus, a full trajectory H consists of a sequence of two or more helical sections The definition intervals of the sections are I 0 = (−∞, t 1 ] for the incoming electron, I i = [t i , t i+1 ], i = 1 . . . N − 1 for the momentarily trapped electron, and I N = [t N , ∞) for the electron that escapes to infinity. The dynamics is the hard-wall limit of a Hamiltonian system with invariance in t, x, and y. Consequently, three physical An example electron trajectory for the inclination γ = 3π/8 [18]. The electron arrives from infinity gyrating around the field line in the left with a pitch angle of γ in = 2.749 (sinγ in = 0.383). After three reflections at the sheath (indicated by the green disks), it emerges gyrating around the field line in the right (shift in y-direction), with the pitch angle set to γout = 0.684 (sinγout = 0.632). Note that the guiding center can go below the sheath. The virtual interaction described in section 2 moves the guiding center directly from its first formal encounter with the sheath (red disk left) to its last encounter (red disk right). Reproduced from [18]. © IOP Publishing Ltd. All rights reserved. quantities exist that are constants of both the free gyromotion and the sheath reflection. These are the Hamiltonian itself, by value the kinetic energy ϵ, and the canonical momenta in the invariant directions p x and p y , all per particle mass: The implicitly defined constants can be interpreted: v denotes the speed of the particle, while X =X and Y =Ŷ − v tanγ cosχ/Ω specify the foot point of the guiding field line. These quantities can be completed to an alternative coordinate system, the gyrocoordinates. The additional variables are the arc length s which locates the guiding center on its field line, and the pitch angle χ and the gyrophase φ which describe, together with the electron speed v, the electron velocity in field-aligned spherical coordinates: x =X + s sinγ + v cosγ sinχ sinφ /Ω, The Jacobian of the transformation is: In the new variables, the equations of motion are simplified. As already stated, v,X, andŶ are constants throughout. During the free gyromotion, the pitch angle χ is also constant, while the arc length s and the gyrophase φ have a constant time derivative: A helical section H i can thus be described by its domain I i and the constants (χ i ,φ i ,t i ), where the time functions s i (t) and φ(t) are given as: In a sheath reflection, the incoming pitch angle χ − and gyrophase φ − are mapped to the outgoing quantities χ + and φ + via: Assume an electron that approaches from infinity on its initial helical section H 0 ≡ H in witht 0 ≡t in ,χ 0 ≡χ in , andφ 0 ≡φ in . The angles (χ in ,φ in ) belong to the negative unit hemisphere. Before the sheath reflection event R i , the electron has the helical elements (χ i−1 ,φ i−1 ,t i−1 ). The reflection takes place when the z-coordinate is equal to zero, i.e. when: With χ − =χ i−1 and φ − =φ i−1 + Ω t i −t i−1 , the scattering law then provides the new helical elementsχ i = χ + ,φ i = φ + − sinφ + tanγ tanχ + , andt i = t i − sinφ + tanγ tanχ + /Ω. After a finite number N of sheath reflections, the electron returns to infinity witht N ≡t out ,χ N ≡χ out , andφ N ≡φ out . The angles (χ out ,φ out ) belong to the positive unit hemisphere S 2 + . The set of all electron trajectories defines a mapping, depending on the inclination angle γ, from the incoming helical elements (t in ,χ in ,φ in ) to the outgoing elements (t out ,χ out ,φ out ).
It can be decomposed into the angular scattering and the time shift: It was shown in [18] that the mappings are piecewise smooth but have, except for γ = 0, kinks when at least one of the reflection points z i = 0 is reached with vanishing velocity v z . In its differentiable points, the Jacobian of the mapping S γ is:

Kinetic description and moment analysis
To study the evolution of an ensemble of different electrons, a kinetic formulation is needed. In the original Cartesian system, the subject is the distribution function The kinetic equation reads: Electrons enter the configuration from 'infinity' with negative velocity in field direction. (Physically, 'infinity' is a region several ten gyration radii away from the boundary sheath.) The specular reflection at the sheath at z = 0 is described by the condition: In principle, the kinetic equation could be explicitly solved by the method of characteristics. This method is, however, of limited use, as the trajectories must be found numerically [18].
In the system of gyrocoordinates, the kinetic equation (25) assumes a more transparent form. Its subject is the distribution function f(X,Ŷ, s, v, χ, φ, t). As the variablesX,Ŷ, and v are constants of motion, they can be treated as parameters and suppressed in the notation. We also assume that f evolves slowly compared to the microscopic timescales that govern the sheath neighborhood and drop the time dependence. Furthermore, we split the distribution into the incoming function f − and the outgoing function f + defined on the subspaces R + 0 × S 2 ± , where S 2 ± are the positive and negative unit hemispheres, v|cosχ| In this picture, the electrons, represented by their gyrocenters, enter the domain from s = ∞ with pitch angles χ ∈ [π/2, π]. This is described by imposing asymptotic conditions on f − . The idea of integrating out the microscopic degrees of freedom would be realized in this model by assuming gyroinvariance, i.e. independence of the distribution f − on the gyrophase φ.
As equation (27) shows, this implies also independence of s.
The gyrophase invariance of f − , however, does not survive the interaction with the sheath. (We exclude the exception γ = 0 where the magnetic field is aligned with the sheath normal.) View mapping (22) as a scattering at s = 0: For each pair (χ, φ) ∈ S 2 + of outgoing elements, the pair − reconstructs the incoming elements they originate from. For a gyroinvariant incoming distribution, the outgoing distribution is: Thus, even if the incoming distribution is gyroinvariant, the outgoing distribution is not. As illustrated in figure 3, it generally shows a pronounced structure in the gyrophase φ which is-at least partially-of a fractal nature. (See [18] for a more detailed discussion.) The subsequent evolution amplifies this even further. Solving the differential equation (28), we get an analytical expression for the distribution f + at the position s. It is a 'twisted' copy of the distribution at s = 0; the function is rotated in the gyrophase φ by an angle that depends on the distance s and the pitch angle χ. Figure 4 provides an illustration: The information contained in a distribution f can be understood in terms of its moments. For a weight function g(χ, φ), symmetric under χ → π − χ, we define the flux moments via the following flux-oriented scalar products over the unit hemispheres S 2 ± : The most prominent flux is the moment of order zero, the particle flux. It is conserved both in the scattering and the gyromotion as no particles are lost or gained in the dynamics: The moments of first order describe the center of weight of the flux: The central moment of second order measures the anisotropicity of the flux, Higher-order moments can be constructed systematically. A complete family of moments is provided by the Zernike polynomials Z m n (χ, ϕ), the natural system of orthogonal functions associated with the scalar product (31). They are products of polynomials of degree n in sinχ with harmonics of degree |m| in φ, with integers n ⩾ 0, −n ⩽ m ⩽ n, and n − m even [22]. See figure 5 for the first few examples, and appendix A for more details: The fluxes corresponding to the Zernike polynomials Z m n (χ, ϕ) will be termed ψ m n ± (s). They include the moments (32)-(35): ψ 0 0 ± is equal to the particle flux ψ (n) ± and ψ (s) ± , i.e. describe the center of weight of the distribution function with respect to the gyrophase, ψ 0 2 ± is proportional to the anisotropy measure ψ (a) ± : For the gyroinvariant incoming distribution f − , the fluxes ψ m n − (s) are nonzero only if m = 0. The ψ 0 n − are spatially constant, as is the distribution itself. They are, however, not conserved in the scattering at the sheath, except for the particle flux ψ 0 0 . (The latter statement simply reflects particle conservation: All electrons which enter the sheath must eventually leave it.) As a consequence, the outgoing distribution f + as given by (29) has fluxes corresponding to all angular moments (37). The spatial evolution of these fluxes after the scattering is governed by the differential equations: For s larger than some ten r L , the phase-mixed distribution Ts f+ is practically gyroinvariant.
Clearly, the fluxes for m = 0 are spatially constant. The moments built with m = 0, however, evolve in s and quickly tend to zero. The slowest moments to disappear are ψ 1 1 + and ψ −1 1 + , but also they essentially vanish after a few ten gyroradii. For the example distribution figure 4, this is shown in figure 6. Thus, the outgoing distribution function converges, in some sense, to its own gyroaverage which is spatially constant: However, this limit cannot be interpreted in the usual sense. This follows from the fact that the norm of the distribution is also a constant of both sheath interaction and gyromotion. Physically, the square of the norm represents the flux of free energy; for a discussion of its relation to the thermodynamic potential of the same name, see [21]: The free energy flux can be seen as a measure of the structural content of the distribution. For any given particle flux, the distribution with the smallest free energy flux ψ ( f ) ± is isotropic, i.e. independent of the pitch angle χ and the gyrophase φ. The conservation of the free energy flux is a direct result of the Hamiltonian character of the dynamics under study. The gyroaveraged distribution, however, has a smaller free energy flux. The equality sign in the following relation is only assumed when the distribution is gyroinvariant already: According to (39), the distribution f + converges towardf + , while (41) asserts that their difference has a constant, generally non-vanishing norm. However, there is no contradiction: The convergence (39) holds 'in the weak sense'. For s → ∞, it is not the norm of the difference f + −f + that vanishes but its scalar product with all fixed functions such as the Z m n . This has an intuitive physical interpretation: Call the determination of a flux moment of type (37) as an 'observation' of the distribution f + , somehow realized by a technical device. As advanced as this device may be, it will always have a finite resolution in phase space. The distribution f + , however, develops increasingly finer structures in a dynamic process which is known as 'filamentation' or 'phase mixing' [20]. (There is a certain similarity to the phenomenon of Landau damping, which, however, has a different physical origin [24].) Ultimately, these fine structures cannot be resolved by the device, which will then only report observations of the phase average. In what might be called a 'thought experiment', the calculation of the norm (40) would be an 'observation of the distribution by itself'. This 'observation' could resolve the filamentation and report the full structural content. However, it would not be possible to realize it by any feasible device.

Interpretation in the frame of functional analysis
To describe the situation more formally, we resort to the framework of functional analysis. The incoming and outgoing distributions are elements of the complex Hilbert spaces H ± defined on the unit hemispheres S 2 ± , respectively [18]: The underlying flux-oriented scalar products are the complex versions of (31), In this notation, the virtual interaction of a distribution with the boundary sheath at s = 0 can be described by the scattering operator (for details, see [18]), The operator has the following properties: • It is unitary, S † γ = S −1 γ , and thus preserves the scalar product, S γ f, S γ f ′ + = f, f ′ − , and the norm, S γ f = f . This implies the conservation of the free energy flux. • It conserves the particle number flux and, as the kinetic energy ϵ is a constant of motion, also the flux of kinetic energy.
• It maps isotropic incoming distributions (which are constant with respect to both the pitch angle and the gyrophase) onto their outgoing equivalents.
For a matrix representation of the scattering operator, [18] used the orthonormal real Zernike polynomials Ψ j (χ, φ) described in appendix B. The matrix elements were found by numerical integration for selected values of the angle γ: The following observations were made: • The matrix S γi j is orthogonal, i.e. its transposed is equal to its inverse, S T γij = S γji = S −1 γij . • The element S γ11 is unity; all other elements in the first row and column vanish. • The absolute values of the elements in the higher rows or columns are bounded by unity. • If the inclination angle γ is equal to zero, S γi j is identical to the unit matrix.
Now, turn to the subsequent transport. Solution (30) can be used to define the operator T s which takes a distribution from the formal scattering plane at s = 0 to an s > 0: This operator rotates the distribution function in the gyrophase φ by an angle that depends on the distance s and the pitch angle χ. Some of its properties are: • It is unitary, T † s = T −1 s , and thus preserves the scalar product, Physically this represents the conservation of the flux of free energy, • It conserves the particle number flux and, as the kinetic energy ϵ is a constant of motion, also the flux of the kinetic energy, • It maps gyroinvariant distributionsf onto themselves, Further insight can be won by writing the Hilbert space H + as the direct sum of the rotational subspaces to the Fourier indices m ∈ Z, H + = ∞ m=−∞ H (m) + . The transport operator is then just a multiplication with a phase factor: For m = 0, T Thus, for any m = 0, the operator T (m) s converges in the weak sense to 0, and the operator T s itself converges, also in the weak sense, to the projection on the gyroinvariant subspace H (0) + , i.e. to the operator · that performs gyrophase averaging: For a quantitative assessment, we calculate the matrix elements of the operator T (m) s in the basis of the complex Zernike polynomials ψ nm (see appendix). As the subspaces for each m propagate independently, the operator is diagonal with respect to the index pairs m, m ′ : In the gyroinvariant subspace, i.e. for m = 0, the result is simply the identity: For subspaces with m = 0, the elements T m n ′ n (s) are given as: T m n ′ n (s) = 2 (n ′ + 1)(n + 1)ˆπ These expression evaluate, for growing values of |m|, to exceedingly lengthy functions of s. Appendix C gives explicit examples and figures. For our argumentation, it is only important that the matrix elements converge quickly to zero, asymptotically as: Altogether, we now get the following interaction scenario: Assume an incoming gyroinvariant distribution function in H (0) − with certain fluxes of particle number and kinetic free energy. That function is mapped onto an outgoing function with the same fluxes of particle number, energy, and free energy, but with a complex fractal dependence on the gyrophase ( figure 3). This function is then twisted by the transport operator and converges in the weak sense to an gyroinvariant distribution ( figure 4). Practically, this takes place within a few ten gyroradii. An effective mappingS γ can be defined which skips the intermediate step and maps the incoming gyroinvariant function directly to the outgoing one, see figure 7: This operator has a number of interesting properties: • The operator is not unitary and does not conserve the norm, i.e. the free energy flux. It is thus dissipative; the outgoing free energy flux is smaller than the incoming one: • The fluxes of particle number and kinetic energy are conserved. • The operator returns a isotropic function identically. In fact, the only distribution whose free energy is not dissipated is the isotropic one.
In the basis of the gyroinvariant Zernike polynomials, where n is even and m is equal to zero, the matrix elements are defined as follows: The structure of the matrix is as follows, where the functions are displayed in figure 8: One can make the following general statements: • The matrix is symmetric, and its elements are even functions of γ. • Elements for n = n ′ = 0 are unity, other elements in the first row or column vanish. • The elements of the other rows and columns are by absolute value smaller than unity. • One eigenvalue of the matrix is unity, all others are smaller. • For γ = 0, and for γ → π/2, the matrix becomes the unit matrix.

Comparison with the action of elastic collisions
The dissipation mechanism discussed in this study is collisionless. It is instructive to compare its effect with that of collisional dissipation. For this purpose, we consider a configuration in which the magnetic field is oriented perpendicular to a boundary sheath located at s = 0. The collisionless dissipation effect is thus absent. Instead, elastic collisions are assumed to take place in the interval [0, l] on the s-axis-specifically, elastic electron-neutral collisions with constant mean free path λ described in the limit of massive scatter centers me/mn → 0. The corresponding kinetic model assumes gyroinvariant distributions f(s, χ) and reads: It is useful to separate the distribution into the incoming f − , defined in the interval [π/2, π], and the outgoing f + , defined in [0, π/2]. Measuring s in units of l and introducing the collision number ν = l/λ, the kinetic equations read: The isotropic source function S represents the scattered electrons: As the boundary condition at s = 0, specular reflection is assumed, From the right, particles enter with a gyroinvariant distribution The goal of the exercise is to determine the outgoing distribution In a first step, equations (63) and (67) can be solved for the function f − (s, χ): Next, the function f + (s, χ) can be found from equations (64) and (66), Evaluated at the point s = 1, this yields the outgoing distribution: We have introduced the weight function w for a compact notation, By evaluating relation (65) we can characterize the source function S(s) by a symmetric Fredholm integral equation of the second type: with the integral kernel given in terms of the exponential integral E 1 , Equations (71)-(74) constitute a mapping of the incoming gyroinvariant distribution f in to the outgoing distribution fout. This can be seen as a collision operator: The operator Cν has the following properties, which are indeed similar to those ofSγ: • The operator is not unitary and does not conserve the norm, i.e. the free energy flux. It is thus dissipative; the outgoing free energy is less than the incoming one: • The fluxes of particle number and energy are conserved.
• For isotropic incoming distributions, the operator returns the outgoing counterpart. Like the scattering operatorSγ, the operator Cν can be expanded into Zernike polynomials. The structure of the corresponding infinite matrix is as follows: A comparison of the properties of the operatorsSγ and Cν reveals the structural similarity, as we can make the following general observations: • The matrix is symmetric, and its elements are even functions of γ. • Elements for n = n ′ = 0 are unity, other elements in the first row or column vanish. • The elements of the other rows and columns are smaller than unity. • One eigenvalue of the matrix is unity, all others are smaller. • For ν = 0, the matrix becomes the unit matrix.
The matrix elements ofSγ and Cν are of comparable size for a collision number of ν ≈ 0.5 (see figure 9). For a typical laboratory magnetron, this translates into a gas pressure of about p = 1 Pa. Many magnetized technological plasma devices (such as magnetrons) work in that regime. It can be therefore assumed that collisional and collisionless dissipation will be of comparable importance and may even act in synergy.

A practical approximation for use in simulations
The representation of the scattering operatorSγ in terms of Zernike polynomials is exact. Of course, the infinite series (77) must be truncated at some index, but the truncation error can be made as small as required. For use in simulations, however, an approximation may be desirable that is easier to handle. We can construct such an approximation under the reasonable assumption that the incoming distribution function f in is only slightly anisotropic. A consistent ansatz for such a nearly isotropic distribution function would be as follows, where the dominant isotropic part f (0) is amended by a small anisotropic correction f (1) . This is a two-term expansion into the gyroinvariant Zernike functions; note the similarity to the two-term expansion into Legendre polynomials [9]: Then, only the first two rows and columns of the scattering matrixS γnn ′ are of importance. All other diagonal elements are set equal toS γ22 , abbreviated asSγ, all off-diagonal elements are set equal to zero: This corresponds to the following integral operator, which returns the isotropic part of the distribution fully and all other contributions weakened by the factorSγ: With the constants a 0 = 2.688, a 1 = −8.544, a 2 = 2.514, a 3 = 29.075, a 4 = −18.757, a 5 = −38.783, a 6 = 21.224, the following polynomial is a reasonable fit (see figure 10): When employed in simulation work, the scattering should be realized as a random process, with the probabilities adjusted so that the expectation value corresponds to (80).

Summary and outlook
Magnetized technological plasmas play an important role in numerous innovative fields from aerospace engineering to surface modification. Their modeling and simulation are of great theoretical and practical interest. Unfortunately, the most fundamental modeling approach, kinetic theory, involves a high level of mathematical complexity: In their most general form, kinetic models are systems of coupled nonlinear partial differential Equations situated in a six-dimensional phase space. This makes them analytically and numerically demanding and not directly suitable for practical simulation work. 'Drift kinetics' and 'gyrokinetics' -not individual theories but rather families of theories-pursue a specific strategy to reduce this complexity. Assuming strong magnetization, i.e. a regime where the gyroradius r L is the smallest of all length scales and the gyrofrequency Ω the largest of all frequencies, they eliminate ('integrate out') the microscopic gyromotion but keep the macroscopic drift. The focus on large length scales and slow timescales enables analytical models that are less complex and simulations that are less demanding.
In high-temperature plasma physics, this reduction strategy has proved highly successful. However, low-temperature plasmas are different. First, of course, only the light electrons can be magnetized, ions do not see the magnetic field at all. But the situation is not entirely clear for the electrons either. Namely, in technological plasmas, the magnetic field lines are not essentially infinite, as they are in typical fusion plasmas, but have defined endpoints at the chamber walls. Electrons approaching these points on their trajectories interact with the plasma boundary sheath. As the scale of the sheath is small and its electric field is strong, this interaction violates the magnetization condition. It is therefore not obvious a priori whether the described reduction strategy would work.
This manuscript, together with a previous publication [18], studies the situation in detail. The emphasis is put on the scale of the gyroradius r L . Assuming the Debye length λ D to be much smaller than r L , the sheath is considered as flat and the interaction represented by the law of specular reflection. Taking, in turn, the mean free path λ and the field gradient scale l as much larger than r L , collisions and the electric field are neglected and the magnetic field is treated as constant. The resulting model has only one free parameter, the inclination γ of the magnetic field direction b relative to the sheath normal n. (All other model parameters disappear in a dimensionless formulation.) Our earlier study describes the electron interaction with the sheath as a virtual scattering of the guiding centers at the plane s = 0. The scattering operator Sγ conserves the fluxes of particle number and energy, as well as the flux of free energy. However, it maps incoming gyroinvariant distributions onto outgoing distributions that have a pronounced, even fractal, structure in the gyrophase. At first glance, thus, the strategy of 'integrating the microscopic degrees of freedom' is indeed not feasible.
In this study, we widened our scope to include the subsequent free gyration of the electrons. Interpreted as a dynamic equation in the parameter s, the distance of the guiding center from the virtual scattering plane, its action could be represented by an evolution operator Ts. This operator also conserves the fluxes of particle number, kinetic energy, and free energy. However, it has a 'twisting' character: It rotates distribution functions in the gyrophase φ by an amount that depends on the pitch angle χ. With growing distances s from the sheath, the twisting gets stronger and the distribution function becomes increasingly filamentary. Asymptotically, any finite observer will be unable to resolve the structures in the gyrophase and will report only the gyro-invariant information contained in the scattered distribution. Practically, this process of phase mixing takes place within some ten gyroradii.
Combining the processes of sheath scattering and subsequent phase mixing, it was possible to define an effective scattering operatorSγ = Sγ which maps gyroinvariant incoming functions into outgoing functions that are also gyroinvariant. The effective operatorSγ conserves the particle and energy flux, but does in general not conserve the flux of free energy. It is, therefore, dissipative, analogous to the action of elastic electron-neutral collisions [21].
In summary, we were able to demonstrate that the envisioned program to simplify models of magnetized technological plasmas by 'integrating out the gyromotion' is not pointless. The additional problems posed by presence of the plasma boundary sheath can be overcome. This insight will help to establish predictive simulation models of magnetized technological plasmas that are less demanding than fully kinetic description. We believe that the effective collision operator described in section 6 may be a first step.
As a final remark, we emphasize that the significance of our results goes beyond modeling. Like the phenomenon of Landau damping-which has a different origin but also involves phase mixing-the effect of collisionless dissipation at the boundary layer of magnetized technical plasmas is of physical importance. Under certain conditions, especially at pressures below one Pa, it can even outweigh the influence of collisional interaction.

Data availability statement
No new data were created or analysed in this study.  The table provides the Noll series for the 20 entries.   j  1  2  3  4  5  6  7  8  9  10  n , m  0,0  1,1  1,−1  2,0  2,−2  2,2  3,−1  3,1  3,−3  3,3  j  11  12  13  14  15  16  17  18  19 The basis functions Ψnm(χ, φ) (real) and Ψ nm (χ, φ) (complex) are defined on the incoming and outgoing hemispheres S 2 ± for integer index pairs (n, m) ∈ I, where: I = (n, m) ∈ Z 2 n ⩾ 0, |m| ⩽ n and n − m even . (A2) The real Zernike functions are normalized versions of the Z m n introduced above, while the complex Zernike functions are: Clearly, the relations are: and The corresponding orthogonality relations are and To enable a matrix representation of linear operators, is it useful to enumerate the orthonormal basis functions Ψnm and Ψ nm by a composite index j. The Noll series [23] provides a bijective mapping N → I of the positive integers into the index set, with j → (n j , m j ). For a given n, functions with a lower value of m are ordered first. Ψnm with cosine terms in φ and ψnm with m > 0 correspond to even j; Ψnm with sine terms in φ and Ψ nm with m < 0 are denoted by odd j. Gyroinvariant functions that are independent of φ have an even n, an m of zero, and the composite index j = (n 2 + n)/2 + 1 is odd. See table 1 for the first 20 entries of the Noll mapping. An algebraic representation is as follows, where the brackets denote the floor function, also called the Gauss-bracket: The renumbered basis functions are now defined as: Ψ j (χ, φ) = Ψ nj mj (χ, φ); (A12) they form a complete set on the Hilbert spaces and obey the orthogonality relations Ψ i |Ψ j σ = Ψ i |Ψ j σ = δ ij for i and j ∈ N. (A13)

Appendix B. Proof of statement (52)
Statement (52) asserts that, for any functions f(χ) and f ′ (χ), and for any integer m = 0, the scalar product f ′ , T (m) s f + converges to zero. Abbreviating h(χ) = 2π f * ′ (χ) f(χ) and measuring the length s in units of v/Ω, the integral in question is of the form: We will first demonstrate that, for functions f ′ and f that are continuously differentiable, the integral (52) converges to zero for both s → 0 and m → 0. Employing the identity: we obtain first, and then after some intermediate steps: Obviously, as the function h is continuously differentiable, the term in the bracket is finite. As it is independent of s and m, the integral converges to zero in both limits: The set of continuously differentiable functions is dense in the Hilbert spaces H ± under study so the statement holds in general.

Appendix C. Evaluation of the matrix element integrals
The matrix elements of the transport operator contain integrals in the Zernike polynomials. As above, s-measured in v/Ωis a non-negative real number, m is a non-negative integer, and n, n ′ are non-negative integers greater or equal to m with even difference. We define the real and the imaginary part of the matrix elements T m n ′ n (s) of equation (56) as: The quantities A 0 nn ′ (s) and B 0 nn ′ (s) are easily evaluated: For m > 0, the functions are nontrivial. Using the integral substitution ζ = ms sec(mχ) they can be written in terms of powers of s, trigonometric functions, and trigonometric integrals. The first few functions are, for |m| ⩽ 3 (see figures 11 and 12):