KINETIC SIMULATION OF VACUUM PLASMA EXPANSION BEYOND THE “PLASMA APPROXIMATION”

Methods: The given results are obtained by the methods of deterministic modeling based on the numerical solution of the system of Vlasov-Poisson equations. Results: The provided comparison of the analytical expressions for the solution of kinetic equations in the "plasma approximation" and the numerical solutions of the Vlasov-Poisson equations system convincingly show the limitations of the "plasma approximation" in some important cases of the considered problem of plasma formation decay. Conclusion: The theoretical results of this work are of great importance for understanding the shortcomings of the "plasma approximation", which can manifest themselves in practical applications of computational plasma physics.


Introduction
The "plasma approximation" is known to be applicable to a low-frequency and steady-state phenomenon, popular among plasma scientists in various fields (Chen, 1984). Sometimes more appropriable term is used here "the plasma condition" (Nishikawa & Wakatani, 1990), i.e. the number of electrons in a Debye sphere is large enough to effect charge shielding. But utilizing it leads to inconsistencies in the equation of motion and prevents a proper, field-theoretic treatment of a condensed matter in the plasma state. This circumstance takes place due to the fact that if plasma reaches a quasi-neutral state ne  ni, its space charge is approximately equal to zero   0. According to Poisson equation, this leads to E = 0. But the "plasma approximation" states that E  0 and the electric field can be found elsewise (Chen, 1984). Such a separation of the initially consistent solution of plasma and field equations in the most cases leads to an ambiguous multivalued interpretation of the electric field definition. That seems to be the main methodological drawback of the "plasma approximation".
As only formal definitions equate with formal mathematics, so the "quasineutral" term correlates to "the plasma condition" ne  ni. It refers to the profound tendency of plasma electrons to change their positions as a response to the electrostatic potential of ions to exponentially attenuate the Coulomb field, and is often taken as the definition of the "plasma approximation." We know that the the traditional term "neutral" already embraces the implications of quasineutrality, as no discrete medium remains neutral on characteristic scales sufficiently smaller to resolve isolated charges. In quasineutral media, the microscopic field fluctuates strongly, but on the particle scale it averages out as the differential volume element grows. The author of the monograph (Chen, 1984) claims that "the plasma approximation is almost the same as the condition of GLASNIK / MILITARY TECHNICAL COURIER, 2022, Vol. 70, Issue 3 quasineutrality discussed earlier but has a more exact meaning … is a mathematical shortcut that one can use even for wave motions… it is usually possible to assume ne = ni and div E  0 at the same time".
The physical kinetics of plasma provides a different point of view on a "quasineutrality" concept given from more fundamental positions. Indeed, if we consider simple two component plasma that consists of least of electrons and single-charged ions, then the ensemble of each type of particles in terms of physical kinetics is characterized by its distribution function, here fe and fi, respectively. Hence the number density of each particle type is a special case of the distribution function zero-moment where (r, p)phase-space coordinates, ttime variable. Such zero-moments, like e.g.
(1), of the particle distribution function, do not characterize the microscopic state of the ensemble of particles. They just represent particular macroscopic characteristics of certain plasma components. That is why the approximation of ne  ni is the equality in a "weak form", i.e. the identity fe = fi does not follow from the quasineutrality condition. Frequently, depending on a particular physical problem based, one has to introduce extra conditions for two or more additional (higher) moments of the distribution function in order to satisfy the "plasma approximation". If we assume that the "plasma approximation" is a convenient computational approach to a number of physical problems, then one has to determine the limits of its use. This paper is aimed to clarify the details of the two-component vacuum plasma bunch expansion into free space by numerically solving the system of Poisson-Vlasov equations. Its main purpose is to show the features of this process without using the "plasma approximation". For simplicity, but without loss of generality, we solve the problem of plasma expansion in a one-dimensional Cartesian spatial configuration. The calculation results are compared with the exact self-similar solutions of the Vlasov equations with the "plasma approximation" (Dorozhkina & Semenov, 1998) pointing out possible shortcomings of the "plasma approximation".

General terms
Let us consider a one-dimensional planar plasma bunch, consisting only of electrons and singly charged ions, located around the point x0, on the x-axis. The bunch has a localization region of the order of xc (spatial distribution half-width). We assume that both ions and electrons inside the bunch have Maxwellian velocity distributions with slightly different temperatures Te for electrons and Ti for ions. Assuming that the intial plasma is quasi-neutral, the particle distribution functions can be written in the following form ( ) where N0full number of particles in the bunch, me and mielectron and ion rest mass, respectively, (x, p)one-dimensional phase-space coordinates, Sbunch transversal cross section, xccharacteristic spatial scale of a plasma bunch, and Te,iinitial ion and electron temperatures, respectively. If the initial electron distribution is assumed to be a non-Maxwellian, then it evolves to a Maxwellian one due to electron-electron elastic collisions. The Maxwellian distribution conserves, since collisions no longer have any influence on the electron distribution function, because the relevant term for elastic collisions is zero for a Maxwellian distribution.
As the effects of collisions of the other kind are much more negligible with regard to electron-electrons, so the restriction to the Maxwellian initial distribution form is justified enough. The electron and ion distribution function comply with the collisionless Boltzmann (Vlasov) equations without a magnetic field (Vlasov, 1968) 0, 0, where E -the electric field vector component along the x-axis, qthe electron charge. During the expansion process, the assumption of a collisionless plasma is valid if the electron-electron collision time ee >> * = L/Cs, where L is the characteristic size and Cs = (qTe/m) 1/2 is the ion sound velocity. As where Coulomb log is   10 and A is the atomic number, we can estimate the condition for a number density ne << n* = Te 2 /(510 -12 A 1/2 L). For an antimony plasma (A = 51) at characteristic lengths of L = 1 cm and Te = 5 eV, the value of n* ~ 710 10 cm -3 (Baitin & Kuzanyan, 1998).
In most cases where the domination of space-charge effects is significant and global plasma quasi-neutrality condition is not met, the system of equations (3) is complemented by the Poisson's equation in order to account the electric field in a self-consistent way where  -electrostatic (electric) potential, 0vacuum dielectric permittivity, ne and ni are electron and ion number densities, respectively, that have to be found from (1).

"Plasma approximation"
Regarding the problem of interest, the first paper where the "plasma approximation" and the associated theoretical approach have been introduced is Gurevich's paper (Gurevich et al, 1966). The case of a halfinfinite plasma with a sharp boundary has been considered by using the Boltzmann electron distribution in order to find the electric field under the quasi-neutrality assumption where n0is the constant characteristic number density of ions. This assumption is not appropriate, since the electrostatic potential is nonstationary and the electron distribution function is different from a Boltzmann for a collisionless plasma. So, the total electron energy changes in time and the electron thermal energy becomes a source of the ion component acceleration in the expanding plasma. In a bounded collisionless plasma, the electrons are trapped and oscillate in a potential well, which is formed in order to satisfy a quasi-neutrality condition. Since the parameters of the potential well change during plasma expansion, so the variation of the electron energy becomes significant. A more complicated approach to the implementation of the "plasma approximation" is based directly on the physical kinetics principles (Dorozhkina & Semenov, 1998).
pf dp p f dp qEn t m x pf dp p f dp qEn t m x In order to obtain the electric field strength from expression (6), the authors of the approach (Dorozhkina & Semenov, 1998) impose an additional approximation, namely, they assume that the plasma bunch is "currentless", so 0, −    e i q pf dp q pf dp (7) in this case by subtracting the second equation from the first one in (6), the following expression for the electric field is obtained As it could be estimated, the given formulas (8) and (9) represent significantly different electric field values for the same plasma parameters. Namely, formula (9) gives a stronger electric field, which is an order of magnitude higher than the similar value obtained by formula (8). These arguments demonstrate the inconsistency of the "plasma approximation" concept. The whole point is the ambiguity of the definition of the electric field, and there are other approaches to obtaining the electric field in "plasma approximation" leading to different electric field estimations which are all different (Baitin & Kuzanyan, 1998 The presented discussion highlights the main disadvantages of the so-called "plasma approximation". First, the variety of the electric field representations in the "plasma approximation" is determined by different theoretical approaches (liquid, kinetic or particle-in-cell) in use, obviously depriving the unambiguity of such approaches. Secondly, the calculation of the electric field requires some additional approximations that come far beyond the basic "quasi-neutrality" condition which is find to be insuffucient. In the scientific literature, the use of such approximations is given without sufficient justification (Dorozhkina & Semenov, 1998). Finally, there is no unambiguous way to identify physical situations where the electric field in the "plasma approximation" smoothly transforms into the field determined from the Poisson's equation or Maxwell's system of equations in the transition regions.

Numerical simulation
The obvious difficulties in choosing the correct formulation of the "plasma approximation" lead to the fact that the most accurate plasma dynamics has to be explained in terms of the complete Vlasov-Poisson system solution (3)-(4) where the electric field is determined in a selfconsistent way. Here we choose a one-dimensional formulation of the problem in the Cartesian coordinates. Its advantages are obvious: the obtained solution results can be directly compared to the analytical formulas in the "plasma approximation" from (Dorozhkina & Semenov, 1998).
The direct numerical integration of (4) by using the trapezoidal or Simpson methods leads to significant inaccuracies associated with the accumulation of errors. To accurately determine the electric field and potential in this work, we used the advanced fourth order method (Knorr et al, 1980). As the computational phase space is restricted to the finite spatial interval x [xmin, xmax], so we apply E = 0 boundary conditions at the both sides of it.
In this paper, the system of partial differential equations (3) was solved numerically on a rectangular uniform phase-space grid (x, p) having 5000 per 2001 grid points for electrons and ions. The Vlasov equations have been solved by using the high-order Cheng-Knorr semi-Lagrangian method similar to that previously used (Zubarev et al, 2020, Kozhevnikov et al, 2021 As an example, here we consider a two-component metallic plasma consisting of electrons and single-charged antimony ions Sb + . This plasma components are typical for vacuum discharge in diodes with antimony cathodes (Anders, 1997). For this plasma composition, a number of numerical calculation series have been carried out. The simulations were processed for a wide range of total number of particles N0 = 10 7 -10 13 . In each calculation, it was assumed that the plasma bunch had a characteristic scale of xc = 100 m, while the plasma was assumed to be nonequilibrium, i.e. Te = 1 eV, Ti = 5 eV that corresponds to real cathode plasma emission (Bugaev et al, 1975). For computational purposes, we restrict the spatial boundaries of the computational phase space to xmin = -2 cm, xmax = 2 cm. It is sufficient to simulate the ionized state behavior far enough from the computational borders. The obtained results were compared with the analytical solutions of the Vlasov equations with an electric field in the "plasma approximation" (8) (Dorozhkina & Semenov, 1998). Figure 1 shows the comparisons of the number densities distributions of electrons and ions at a time point of t = 500 ns for the cases of plasma expansion consisting of different initial number of particles. A quasi-neutral plasma distribution profile is shown with a black line in accordance with the analytical solution in the "plasma approximation" (Dorozhkina & Semenov, 1998). The results of the numerical calculations (without the "plasma approximation") -the number densities distribution profiles -are given for both electrons and ions to show the difference in their spatial distributions.
The first two plots in Figure 1 correspond to the decay of ionized states with a small number of particles -N0 = 10 7 and N0 = 10 9 . Such initial distributions of charged particles cannot be called "plasma" due to the fact that the Debye length (Chen, 1984) is much larger than the characteristic scale of the bunch, i.e. D >> xc. In the first case, the Debye length is much greater than xc, in the second case it has the same order of magnitude. In both cases, the numerical calculation shows a violation of the initially given quasi-neutrality conditions and a significant deviation from the "plasma approximation" profiles.
For a denser plasma (N0 = 10 11 ), a different situation is observed. This case corresponds to true plasma decay D << xc. The electron and ions number density distribution profiles in the Vlasov-Poisson model are close to the quasi-neutral profile obtained from the "plasma approximation". For the specified time point (t = 500 ns), quasi-neutrality is not significantly disturbed over the entire length of the plasma bunch. This situation is most fully described by the "plasma approximation": the plasma bunch expands with thermal velocities without significant loss of the initial quasi-neutrality. At the very beginning of the process, the more mobile and thermalized electronic component of the plasma is displaced with regard to the electrically neutral state, which leads to the appearance of a weak electric field close to the field in the "plasma approximation" in this case. Thus, plasma tends to remain quasi-neutral: if the ions move, then the electrons will follow them, and the electric field is adjusting to maintain the neutrality in accordance with the displacement of electrons and ions. Finally, the most important case considered here represents the decay of a dense plasma bunch. The corresponding result of these calculations is given at the fourth plot in Figure 1. For a plasma bunch with the total number of particles equal to N0 = 10 13 , other regularities are observed. First of all, one can find a similarly to the previous case N0 = 10 11 . For N0 = 10 13 the plasma bunch expands in time while preserving quasi-neutrality. But resulting from the calculations without the "plasma approximation", plasma decays faster. The profiles of a quasi-neutral plasma in the "plasma approximation" and that one obtained as a numerical solution of the Vlasov-Poisson equations (without the "plasma approximation") noticeably differ. In comparison with the decay of the ionized state (in N0 = 10 9 ), where the quasi-neutrality is noticeably violated, the decay of the dense plasma is more intense with regard to the reference calculation (black line in Figure 1).
Since the bunch quasi-neutrality is not violated (for N0 = 10 13 ), it can be assumed that in this case other factors affect the plasma decay process, which are not taken into account by the "plasma approximation" (Dorozhkina & Semenov, 1998). As we have already mentioned previously, the kinetic formulation of the "plasma approximation" requires an additional "currentless" approximation (7). The "currentless" approximation is introduced independently of the quasi-neutrality approximation. It leads to a linearization of kinetic equations (3), so that the plasma number densities (black lines in Figure 1) for different number of particles have the same characteristic width and differ only by the scale factor of the curve magnitude. In reality, the system of Vlasov-Poisson equations is essentially non-linear and its solutions for various parameters do not scale. The calculations show that condition (7) is violated for a dense plasma. It leads to faster plasma expansion due to the nonlinearly increasing influence of electron and ion currents without affecting the total bunch quasi-neutrality.

Conclusions
The results presented in this paper represent a counterexample convincingly showing the groundlessness of the "plasma approximation" for solving particular non-stationary plasma physics problems. The existing contradictions in the use of the "plasma approximation" can be formulated as the following theoretical statements: • In all of the existing "plasma approximation" formulations, the electric field at t = 0 is nonzero, while the initial plasma is quasineutral, so the electric field initially is E = 0; • Accurate numerical simulation (without the "plasma approximation") shows that for some cases of the quasi-neutral plasma bunch decay a local violation of plasma electrical neutrality appears. This leads to the electric field redistribution and affects further plasma dynamics, while in the "plasma approximation" plasma quasi-neutrality is preserved every time; • Finally, the known "plasma approximations" are ambiguous. They require additional physical approximations depending on the theoretical approach in use. Such approximations in some cases make the "plasma approximation" less accurate.
Following the comparisons given in this paper, it can be argued that the "plasma approximation" can be used to study plasma decay without an external electric field only for a restricted range of dense plasma parameters. In other cases, the initial quasi-neutrality conditions do not guarantee the preservation of electrical neutrality during the whole ionized state decay process making the use of the "plasma approximation" unacceptable.