Nonlinear oscillations in the Knudsen plasma diodes

Under analytical and numerical study are the kinetic phenomena inherent for nonlinear oscillations in 1D diodes in which charged particles are supplied from an emitter and move with no collision in the interelectrode spacing. An original numerical method (E,K code) is used to build up the ion and electron velocity distribution functions as well to study the structures forming in the course of the nonlinear oscillations both in the vacuum diode with monoenergetic electron beam and in the Knudsen diode with surface ionization. The oscillations are shown to occur owing to the Bursian-Pierce instability.

Time-dependent processes occurred in planar plasma diodes are considered. We suppose that the electrons and ions with a prescribed velocity distribution functions (DF) leave the emitter surface and travel collisionlessly between electrodes (the Knudsen regime) at a self-consistent electric field. Two type of devices are considered: 1. The vacuum diode with no ions and DF of the emitted electrons being near the monoenergetic one (the Bursian diode), 2. The Knudsen diode with surface ionization, where the emitted electrons and ions have the half-Maxwellian DF (the KDSI).
The devices exhibit the large amplitude oscillations of the electron current across the diode when a current of the entered electrons exceeds a certain threshold value. The paper presents a theory of such oscillations.

Numerical code
The highly nonlinear oscillations arisen in the diode are accompanied by an intense energy exchange between charged particles (electrons and ions) and the time-dependent electric field. Such process results in a strongly non-equilibrium DF and even its disruption. In order to describe correctly such processes the kinetic and Poisson's equations should be solved. Such a problem can be solved analytically only for very particular cases. Therefore, we solve this problem numerically.
The kinetic part of the problem is the most difficult one. The best known and used in nowadays plasma simulations are codes based on plasma model as a system of macroparticles (see, e.g., XPDP1 code [1]). These codes, however, have a disadvantage of being rather noisy, which is mainly due to smaller number of particles within each cell. They are presented as rather high-speed ones. However, as it is shown in Reference [2], only a strong decrease in size of space and time steps of computational grid and, consequently, an appreciable increase in the number of macroparticles makes it possible to obtain the reliable results with XPDP1 code, in particular, in situations involving the particle reflection.
We have developed a new numerical code, which is called the E, K code in the literature [3]- [6]. The code involves the fact that, in the collisionless case, the relationship takes place: f (τ, ζ(τ ), u(τ )) = f + 0 (τ 0 , u 0 ). (1) Here f and f + 0 are the DF at the point (ζ, τ ) and at the emitter, respectively, ζ(τ ) = ζ(τ ; τ 0 , u 0 ), u(τ ) = u(τ ; τ 0 , u 0 ), andζ(τ 0 ) = u 0 . In calculating the DF at the time τ = t, the electric field at all the points of the interelectrode gap is known at all preceding times, i.e., over entire time interval (0, t). The values of the field distribution ε(t, ζ) at the nodes of the spatio-temporal mesh are stored in computer memory. Calculating the DF at a mesh node is reduced to computing a number of the trajectories of test particles.
The main feature of the method is that the trajectory of each particle is calculated backward in time until it reaches the electrode surface. As a result, from a given arrival velocity u 1 , the velocity u 1 0 and time τ 1 0 at which a particle is injected from the emitter are determined, as well as the value of the DF at the velocity u 1 . Next, the code gives the trajectory for a new arrival velocity u 2 = u 1 + ∆u, and so on. The velocity step ∆u is chosen so that the difference between the values of the DF at the neighbouring trajectories does not exceed a given value: As a result, the step ∆u is chosen automatically from the gradient of the DF. This choice of ∆u provides the required accuracy in the calculation of the DF and its moments and, in particular, ensures a high accuracy in the determination of the points in which the DF has discontinuities. In calculating the trajectories, the electric field strength within each cell is approximated by a linear function of the coordinate and time, and the position of a particle and its velocity are calculated as the power series: For the coefficients a s , the simple recurrence relations were obtained: Here, Z = ζ(τ 1 ) and U =ζ(τ 1 ). This method for calculating the trajectories ensures that the electric field is continuous when the boundaries of the cells are crossed -a circumstance that is especially important for the trajectories of reflected particles -and simultaneously provides high accuracy in the calculation of the parameters of the trajectories.
In Reference [7], E, K code was checked for several electric fields, which allow to obtain analytical formulas for trajectories of particles, as well as for self-consistent time-dependent processes, which are described by analytical theory [8]. All cases under consideration demonstrated high accuracy of the code. Comparison of E,K code and XPDP1 code shown, that the former exceed considerably XPDP1 code in accuracy [7]. Besides, in order to obtain reliable results using XPDP1 code one has to make an average over a large number of time steps. However, its value is a free unknown parameter.

The Bursian diode
About 30 years ago, the vacuum diode with directed electron beam attracted an attention in connection with the creation of high power microwave generators, such as vircator, reditron, and so on. Under certain conditions, a potential barrier, reflected a portion of electrons, is formed, named the virtual cathode (VC), and highly nonlinear oscillations arise in the diode.

Steady solutions
In plasma diodes with an electron beam it is convenient to use the kinetic energy of the electrons at the emitter W b and the beam Debye length λ D as energy and length units, respectively [9]: Here, the current density of the beam J b and the acceleration voltage V b = W b /e are expressed in amperes per square centimeter and volts, respectively. The steady solutions of the Bursian diode are determined only by two dimensionless values: an inter-electrode distance δ = d/λ D and the external voltage The PD with a single potential minimum is typical for the steady solutions of the diode. When the minimum potential |η m | turns into the initial electron energy, a portion of electrons is reflected and returns to the emitter. This is the solution equipped with a virtual cathode. It names a solution of electron reflection. If V C is fixed, the steady solutions lie on the dependence ε 0 (δ) on the {ε 0 , δ} plane ( Figure 1). The branches 1 and 2 correspond to solutions without electron reflection, and 3 -with reflection [9].
A stability of steady solutions without reflection with respect to the small perturbations has been studied by using a dispersion equation. The solutions lay on branch 1 are stable, whereas those on branch 2 are unstable. However, there is no dispersion equation for solutions with reflection. We have investigated the stability properties of such solutions with E, K code [6]. In calculations, the time scale is selected to be the time needed an electron to overcome λ D . The code is the exact one giving an opportunity to determine the eigenvalues (growth rate Γ and frequency Ω) of the main mode from simulation of a small perturbation evolution.
We have calculated such process for a set of gap values, and obtained the dependencies of eigenvalues on δ. One can see from Figure 2, that a threshold δ l exists below which all solutions on branch 3 are stable to small electron perturbations. For δ > δ l , the steady solutions are unstable to the small perturbations. As δ increases, the alternation of regions of stable and unstable steady solutions occurs. Thus, the steady-states with a partial reflection of electrons from the VC can be realized in the Bursian diode.

Nonlinear oscillations
In order to see the process termination when a small perturbation increases, it is necessary to study its nonlinear stage. We studied a regime with the VC at zero external voltage.
In numerical simulations, the time step was ∆τ = 0.05; the coordinate step was ∆ζ = δ/N with N = 200; and the electron DF at the emitter was chosen to be with a small velocity spread ∆. At each time moment τ p , firstly the DF of electrons is calculated with E, K code. Then having calculated the DF and electron density, we solve the Poisson's equation with the following boundary conditions for the potential:  The problem reduces to solving the set of linear algebraic equations Here, 1 ≤ k ≤ N and (∆n) p k is the electron charge in the layer (ζ p k−1 , ζ p k ), which is calculated with allowance for the electron fluxes ν through the cell boundaries, At each time step τ p , the self-consistent solution is achieved as follows. For the time τ p , the electric field distribution is known only at all preceding time moments τ ≤ τ p−1 , but not yet at τ = τ p . The field is extrapolated to the time interval (τ p−1 , τ p ), and the DF and electron density at the time τ p are then calculated. The density being known, the refined field distribution at the time τ p is calculated and the density is recalculated for the new field. As a rule, it is sufficient to carry out only one to two such iterations.
The nonlinear process ends, as a rule, with the periodic oscillations. The dependencies of the potential η m at the top of the VC on its position ζ m , when growth of a small perturbation finished, are closed curves (Figure 3,a-d) [10]. As the electrode gap increases, the shape of the curve changes from elliptic (near the threshold) to nearly triangular. Figure 3 shows such curves for the 1st instability region (δ l , δ r ).
The E, K code gives an opportunity to plot the electron velocity DF and density n e (ζ) at any time with high accuracy. Figure 4,a shows the DF projection on the velocity-coordinate plane at a fixed time. The DF has a very complicated structure. Integration of the DF over velocities yields the electron distribution within the gap. Figure 4,b shows that the electron density varies by several orders of magnitude. It has a sharp peak at the point corresponding to the left boundary of the region from which the electrons are reflected.
For each instant of time, we can obtain such fine parameters of the process as a position ζ M and potential η M of the point corresponding to the maximum of the density distribution and to plot the attractor in the {η M , ζ M } plane (Figure 3,e-h). The attractor is a closed curve of an oval shape. In analogy with the theory of finite-dimensional dynamic systems, we call such a curve the limiting cycle. In Figure 3, the heavy dot within the cycle corresponds to a steady solution. At the final stage of the process, the dependence η M on ζ M winds onto the limiting cycle from inside (Figure 3,g). Further, we should check: i) whether the relevant solutions are stable or not, and ii) whether they are unique for given external parameters.
In order to solve these problems we selected the initial conditions responsible for the points located beyond the cycles. In these processes, the dependence η M (ζ M ) wound like a helix onto the limiting cycle from outside (Figure 3,f). This cycle coincides with the limiting cycle obtained in simulating the development of a small perturbation from the steady solution. Thus, in the 1st instability region, the steady limiting cycles correspond to the regime of nonlinear periodic oscillations and each such cycle is unique. We calculated a set of the transient processes started from initial conditions located beyond the cycles, for the electrode distances beyond the instability threshold δ l . We found that, in these cases, the processes ended also with nonlinear periodic oscillations of the same type as those in the 1st instability region.
In the 2nd instability region (δ 2 l , δ 2 r ), we revealed a new class of oscillations of which amplitudes are substantially smaller than that of the oscillations described above. In Figure   25th  5, the new cycle is shown as the minor oval. However, for these gap values, large-amplitude oscillations exist also. We found that, at certain times, the convective current density increases in a jump like manner (Figure 6,a). The jump in the current develops on a time scale much shorter than the electron time-of-flight through the electrode gap. Such behaviour of the current density stems from the character of the evolution of the electron distribution in the electrode gap during the oscillations. The reason is that, at the times when the height of the VC has a minimum, a fairly large group of electrons overcomes the potential barrier. Further, these electrons form a sharp front in their spatial distribution (Figure 6,b). The front moves toward the collector, and when it reaches the collector, the convective current density increases sharply.

Long-lived electrons
Earlier in numerical simulations an unusual phenomenon, being characterized for oscillations, was revealed: existence of so-called long-lived electrons [11]. These are the electrons that reach the VC, which serves as a potential barrier for them, and oscillate with it during several periods. We have found a reason for the onset of long-lived electrons, and analysed their properties [10]. In order to study the features of the long-lived electrons we have considered a model of the electric field, which describes properly trajectory of the VC oscillation: Here, α 2 = −2η 0 m /(ζ 0 m ) 2 , η 0 m and ζ 0 m are the height and the position of the VC top in the steadystate case (κ = 0), and the parameters Ω and κ characterize the frequency and amplitude of the VC oscillations. The quantity α is determined by the curvature of the PD in the vicinity of the VC top. For an electron that is injected from the emitter with the velocity v 0 at the time τ 0 , we obtain the following expressions: Here, B = (1 + ϕ) −1 , ϕ = Ω/α, A ± = ζ 0 m ∓ v 0 /α + ζ 0 m √ Bκ cos [Ωτ 0 ∓ arctan(ϕ)] /2. These formulas show that, at α(t − τ 0 ) 1, the electrons can either overcome the potential barrier and reach the opposite electrode (A + > 0) or they can be reflected from the barrier and return to the emitter (A + < 0) that depends on the relationship between v 0 and τ 0 . It is also seen that a group of electrons exists that can oscillate in the vicinity of the VC for a fairly long time until they fall into the 1st or the 2nd group. These are just This relationship is deduced from the condition A + = 0, and does not depend on parameters of the arrival point (ζ, τ ). Trajectories of long-lived electrons are determined by last term in formulas (11). Hence, long-lived electrons are those beam electrons, the kinetic energy of which is close to zero when they reach the vicinity of the potential minimum. As a result, such electrons are reflected from the potential barrier, but they move slower than the VC (B < 1), which thus overtakes them, so they again catch up with the barrier and are again reflected from it, and so on. It is the reason why long-lived electrons exist. More over, any electron "lives" only a finite time. Indeed, the formulas (11) contain an exponentially increasing term. Therefore electron trajectory is unstable to small variations in the initial conditions. On sufficiently long time scales τ , this term will become greater than the remaining terms and the electron will come to one of the electrodes, no matter how close the point (v 0 , τ 0 ) is to curve (12) or how small the coefficient A + is.

Nonlinear oscillations in the KDSI
The Knudsen diode with surface ionization (KDSI) is a collisionless, one-dimensional, singleemitter plasma diode. The emitter serves as a source of electrons by Richardson emission and ions by surface ionization of neutral alkali atoms. The emitted electron and ion DFs are the half-Maxwellian ones. The Knudsen thermionic energy converter [12] and Q machine [13] are examples of the KDSI.

Steady solutions
The self-consistent steady solutions of the KDSI are determined by three dimensionless values: emitter neutralization ratio γ = n 0,+ i /n 0,+ e , electrode distance δ = d/λ D , and external voltage V C = eΦ C /kT E . Typical PD exhibits a quasineutral plasma region with practically constant potential, embedded between two space-charge sheaths adjacent to the electrodes. We have analysed stability properties of such solutions in overneutralizated regime (γ > 1), and built up a boundary of the region of stable solutions (Figure 7) [14]. Fat points show data of experiment carried out to determine the boundary beyond which the oscillations initiate [15]. One can see the excellent agreement between our theory and this experiment.

Nonlinear oscillations
A number of experiments demonstrated the large amplitude oscillations of electron current in the KDSI (see, e.g., [12,13]). Their period is comparable with the ion transit time of the electrode gap d. We have shown that the oscillations are characterized by the periodic sequences of two types of stages: the slow (related to the slow motion of the ions) and the fast (compared with the time for an electron to overcome the electrode distance) ones [4]. In calculating the slow stage, we assume that the electrons overcome the gap d before the ions cross the distance of the Debye length λ D . Then, we consider that to the moment when ions moved over distance λ D , the electrons and electric field in the electrode gap have already a time to redistribute and conform with a given ion distribution. Therefore the time scale is selected to be the time needed an ion to overcome λ D . The ion DF is calculated by E, K code. At each time-step τ p , the electron and electric field distributions for known ion background n i (τ p , ζ) are found from the stationary problem. In this case, an explicit expression for electron density can be obtained for PD of any type. It depends on a potential value at the given point, and also at the points of PD local minimum, being the points of electron reflection: n e = n e (η, η 1 , . . . , η l ). Inserting the electron density in the Poisson equation, we reduce the problem to solution of nonlinear ordinary differential equation of a 2nd-order: with boundary conditions: As a rule, this problem has several solutions. The PDs that can be realized as well as their stability properties, are analysed by the η, ε-diagram technique [14]. During the slow stage, the ions are redistributed in the electrode gap. At certain instants of this stage, the conditions of the Pierce type instability onset appear [16]: a potential jump forms in the emitter sheath, which accelerates strongly the electrons, and electron beam with a small spread in velocities propagates through the plasma. At such instants KDSI plasma differs from the Pierce diode one only due to the ion background, through which the electron beam goes. The ion background is nonuniform in the KDSI.
During the fast stage, the ion distribution may be considered as the steady one. And the time scale is selected as the time needed an electron to overcome λ D . Self-consistent simulations of this stage are similar to those in the vacuum diode. Figure 8,a is an example of calculation of the electron transient process which finishes in a state of the VC (curve 4). The initial PD is a distribution of the monotonic growth. Relevant state belongs to the boundary of a stable region (curve 1). Figure 8,b shows the maximum deviation ∆Φ of a potential from the initial one in time [17].
The instability gives an interesting physical phenomenon -a cut-off of current. It occurs in the case when the electron transient process finishes in a state for which the PD with the VC realizes: a region with negative potential forms at a distance from the emitter (Figure 9,a). As a result a strong limitation of electron current, passing through the diode, occurs. can speak that the current changes instantaneously in diode. This structure is stable. During the slow stage it moves to the collector with a characteristic ion velocity ∼ (2kT E /m i ) 1/2 and drives an electron current, passing through the diode. In Figure 9,b, these structures correspond to the parts of the solid curves with the current density less than j E . The vertical dashed lines denote the fast electron transitions. Oscillation period is about a time-of-flight of the thermal ions over the diode gap, and modulation depth β = [j max (t) − j min (t)] /j max (t) ∼ 0.5. These facts are corroborated by experiments [12,13]. From the DP evolution during an oscillation period (Figure 9,a) one can see that, within the time interval (12.7 < τ < 16.4) to the right of a potential minimum, a double layer exists, within which the potential changes by an amount of about an external voltage Φ C , and an essential separation of charges occurs at the distance of several λ D . The electron current, of which magnitude depends on a potential barrier height for electrons, passes through the layer.
Development of electron instability gives rise to a number of interesting physical phenomena. For example, during the fast stage, a potential well is formed (Figure 9,a, τ > 16.4) and accumulation of potential energy occurs. Then, during a slow stage, the latter converts into directed kinetic energy of the ions and thus results in ion acceleration up to relatively high energy. More over, the ions as beams had a very small velocity spread (Figure 10,a), and the ion kinetic energy localization occurs (Figure 10,b).  Figure 10. Evolution of an ion DF at the collector (a) and of ion kinetic energy during an oscillation period; γ = 5, δ = 16, V C = 5.

Collisionless trap of electrons in the potential well
The PD with a hump is a potential well for electrons (Figure 10,a). Self-consistent collisionless electron trapping in the well, may result due to the energy exchange between electrons and the time-dependent electric field arising during the well formation. Really, when the well depth increases, an electron, moving in the well, loses its energy, to the electric field. With increasing gap length, the energy losses growth, and a threshold gap d * exists, in exceeding which these losses become larger than the initial energy of electron at the emitter. As a result, the electron reflects from a potential barrier in the collector region, and turns out to be trapped in the potential well [18]. The trapped electrons do not spread over the well as a time-independent structure, but they are a dynamic structure: forming a localized spatial bunch with sharp fronts, which bounces between the electrodes, being broaden, when passing the center region of the gap (Figure 11,a, curve 2), and shrank in the vicinity of electrodes (curves 1 and 3).
During the process, a charge of the bunch electrons decreases due to some electrons get away the electrodes. Such a time-dependent behaviour of the trapped electrons leads to oscillations of the PD (Figure 11,b), so that the well depth oscillations also (solid curve), being less than this one if the trapped electrons were absent (dashed curve). The electron trapping plays the crucial role in self-consistent calculations of nonlinear KDSI dynamics.

Conclusion
Thus, as in the diode with the beam-like velocity distribution function of emitted electrons, as in the Knudsen diode with surface ionization the large amplitude oscillations of the electron current can exist. They arise when a current of the incoming electrons exceeds certain threshold value which is of the order of the plasma beam Debye length. As a result, the Bursian-Pierce instability develops resulting in a current cut off. During the oscillation process a set of nonlinear structures form. For example, dynamical VC and DL exist, the long-lived and trapped electrons  Figure 11. Spatial distribution of trapped electrons at three time moments (a) and evolution of a well depth (b); δ = 1.97π.
appear, bunching charge particles occurs, and so on. We proposed to create fast electron switches on the basis of theory of nonlinear oscillations in the Bursian diode [8]. Results of this theory can be useful for microwave generators which convert oscillatory energy of the self-consistent electric field into microwave radiation. Used the theory of nonlinear oscillations in the KDSI we proposed new type of the themionic energy converter which converts heat energy directly to alternate current power [19].