Steady state stability features of a diode with electrons and positrons entering from opposite boundaries and moving without collisions in plasma are numerically studied. The most complex regime when charged particles are reflected from potential barriers is considered. This problem arises, in particular, when modeling pulsar diodes. A small perturbation evolution is studied. It has been established that at the initial stage of the process the perturbation amplitude changes in time according to an exponential law. It is shown that stationary solutions with a potential barrier for electrons located near the electron-emitting electrode and a potential barrier for positrons located near the opposite electrode are stable when the inter-electrode distance is below a certain threshold. As the inter-electrode distance increases, the solutions become unstable. Solutions of another type when barriers reflecting particles are located in the opposite to the emitting electrode parts of the gap are also studied. However, these solutions turned out to be unstable.

## I. INTRODUCTION

Plasma consisting of electrons and positrons occurs in various high-energy astrophysical objects, such as neutron stars and their surroundings, black holes and active galactic nuclei, and relativistic jets that produce gamma-ray bursts. In an attempt to better understand its properties, not only much theoretical work is published, but laboratory experiments also have been carried out (see, e.g., the review^{1}). One of the most complex astrophysical objects containing electron–positron plasma is pulsars. These radio sources were discovered more than 50 years ago. Despite the fact that the main models describing the structure of the pulsar magnetosphere and the processes occurring in it appeared right after their discovery, there is still no clear idea of either this radiation mechanism or the reason for switching between modes (see, e.g., Ref. 2).

Reference 3 suggests a hypothesis in which the pulsar radiation is caused by electric field oscillations in the diode plasma, arising due to the instability of steady state. This kind of instability is typical for diodes with collisionless plasma, and the mechanism of instability development is the same as that of the well-known Pierce instability, where an electron beam moves in the diode space through an immovable ion background.^{4} It should be noted that collective effects in electron–positron plasma created by a laser beam were studied, for example, in Ref. 5.

In Ref. 3, a model of a vacuum diode with counter flows of electrons and positrons entering from opposite boundaries is proposed. Steady states of such diodes are studied in detail in Ref. 6. It has been established that they can be divided into two main types depending on the charged particle movement type: (1) all particles reach the opposite electrode and (2) a particle portion is reflected from the potential barrier within the diode and returns to the emitting electrode. In Ref. 3, an equation for the electric field perturbation for stationary solutions of the first type is derived. For a uniform field distribution, an analytical solution to this equation is found, a dispersion equation is obtained, and dispersion branches are constructed. It has been proven that homogeneous solutions are stable when inter-electrode gap values are less than $2\u2009\pi \u2009\lambda D$, where $\lambda D$ is the Debye–Hückel length. The stability features of the first type state with a nonuniform field distribution were studied in Ref. 7. It is proven that all such solutions are unstable with respect to small perturbations.

In this work, we study numerically the stability features of stationary solutions for the regime with particle reflection from potential barriers. A numerical research method is used due to the fact that to date no equation for the electric field perturbation for such a regime has been obtained. The study comes down to calculating the initial stage of a small disturbance development. From the time characteristics of the process, the perturbation eigenmodes are determined, and by the sign of the growth rate, i.e., the real part of the complex frequency, it is determined whether the stationary solution is stable. To validate the results obtained, calculations are carried out using two numerical codes: EK-code and PIC-code.

## II. STATIONARY SOLUTIONS

In Refs. 3 and 7, we have obtained an analytical solution of equation for the electric field perturbation. We could use monokinetic velocity distribution functions (VDF) to describe the behavior of charged particles in the cases under consideration. However, in this paper, the stability problem is solved numerically, and it is necessary to use a non-monokinetic VDF. We consider that the non-relativistic electron flow enters from the left electrode, and the positron one enters from the right electrode with VDF slightly different from the monoenergetic one. We assume that the average velocities of incoming particles are the same, but of the opposite sign: $v\xafp,0=\u2212v\xafe,0$. In this case, the average energies of incoming particles are equal to $W0\u2261mev\xafe,02/2=mpv\xafp,02/2$, where $me$ and $mp$ are the electron and positron masses, respectively ( $me=mp=m$). We also assume that the particle densities at the electrodes are equal: $ne=np=n0$. Particles move in the inter-electrode gap without collisions. In addition, we consider that a particle that reaches any electrode is absorbed by it. The potential difference *U* between the electrodes is assumed to be zero.

We pass on to dimensionless quantities, choosing the particle energy $W0$ and the Debye–Hückel length $\lambda D=[(2\epsilon \u03030W0)/(e2ne,0)]1/2$ (*e* is the electron charge, and $\epsilon \u03030$ is the vacuum permittivity) as units of energy and length. For dimensionless coordinate, potential and electric field strength, velocity, and time, we have $\zeta =z/\lambda D$, $\eta =e\Phi /(2W0)$, $\epsilon =eE\lambda D/(2W0)$, $u=v/v\xafe,0$, and $\tau =t/(\lambda D/v\xafe,0)$.

At first, we obtain stationary solutions for the VDFs (1) and (2). If $\Delta $ is small, they should not differ much from the solutions corresponding to the case of $\delta $-shaped VDF of particles, described in detail in Ref. 6.

Steady-state solutions are determined by three dimensionless parameters: the inter-electrode distance $\delta =d/\lambda D$, the potential difference between the electrodes $V=eU/(2W0)$, and the electric field strength at the left electrode $\epsilon 0$. However, they are also slightly influenced by the velocity spread $\Delta $.

For a fixed potential value *V*, it is convenient to represent steady-state solutions by points on the ( $\epsilon 0,\delta $) plane. These points form separate curves, which are called solution branches.^{6} They correspond to wave-like potential distributions (PD), which for solutions belonging to different branches differ in the number of extrema and their positions relative to each other. The relation between the potential minimum and maximum $\eta m$ and $\eta M$ is given by the formula $\eta M=V\u2212\eta m$.^{6}

In the case of our interest, when $V=0$ and particles enter the plasma from opposite electrodes with identical charges, masses, and kinetic energies, the total charge in the inter-electrode gap should be equal to zero, and the PDs should have odd symmetry about the gap center.^{3} Due to the symmetry, the number of solution branches turns out to be less than those for the general case $V\u22600$, described in Ref. 6. Taking into account the symmetry, it also makes an opportunity to correct calculations when modeling time-dependent processes in the diode.

In Fig. 1, solution branches for $\Delta =0$ and $\Delta =0.01$ are shown. The solution branches for the mode without particle reflection from potential extrema are designated as $nk$ (index *k* is the number of extrema on the PD) in Fig. 1. In this case, $\eta m>\u2212(1\u2212\Delta )2/2$. In the case when electron reflection occurs, we call the potential minimum the virtual electron emitter (*e*-VE). Similarly, we call the potential maximum the virtual positron emitter (*p*-VE). There are two types of PD with particle reflection. When the *e*-VE is located to the left of the *p*-VE, the corresponding branches are designated $dk$ (the index *k* is the number of extrema on the PD lying between *e*-VE and *p*-VE) in Fig. 1. In the opposite case, when the *p*-VE is located to the left of the *e*-VE, the corresponding branches are designated $di,j$ (the index *i* is the number of PD minima lying to the left of the *p*-VE, and index *j* is the number of PD maxima lying to the right of the *e*-VE) in Fig. 1. Due to the symmetry, in the case of $V=0$, solutions belonging to branches $nk$ and $dk$ with $k=0,\u20092,\u20094,\u2009\u2026$, as well as $ds,s$ with $s=0,\u20091,\u2009\u2026$ can only occur. Examples of PD corresponding to branches $n2$ and $d0$ are shown in Fig. 2(a), and for branches $d00$ and $d11$ in Fig. 2(b). The potential distributions for solutions without reflection were obtained in Ref. 7.

^{8}

### A. $d0$ branch

The characteristic PD for the $d0$ branch is shown in Fig. 2(a) with a red curve. Since the solution is symmetrical to the middle of the inter-electrode gap, it is sufficient to consider the PD on the interval $(0,\delta /2)$.

To find the PD, we divide the interval $(0,\delta /2)$ into three segments: (1) from the left boundary ( $\zeta =0$) to the point $\zeta \u2212$, where the reflection of electrons starts (potential $\eta =\eta \u2212=\u2212(1\u2212\Delta )2/2$); (2) from the point $\zeta \u2212$ to the position of the potential minimum $\zeta m$, where the reflection of electron finishes; and (3) from the point $\zeta m$ to $\delta /2$, where all the unreflected electrons penetrate. In these three segments, the electron densities are described by various formulas [the first, second, and third of formula (5), respectively], which leads to different types of functions $D(\eta )$ in relation (7). The form of these functions is given in Table I. The positron density is determined by the first of formula (6).

Segment . | $D(\eta )$
. |
---|---|

$[0,\zeta \u2212]$ | $D1e(\eta )=[(1+\Delta )2+2\u2009\eta ]3/2+(\u22122\u2009\eta m+2\u2009\eta )3/2\u22122\u2009[(1\u2212\Delta )2+2\u2009\eta ]3/2+[(1+\Delta )2\u22122\u2009\eta ]3/2\u2212(\u22122\u2009\eta m\u22122\u2009\eta )3/2$ |

$[\zeta \u2212,\zeta m]$ | $D2e(\eta )=[(1+\Delta )2+2\u2009\eta ]3/2+(\u22122\u2009\eta m+2\u2009\eta )3/2+[(1+\Delta )2\u22122\u2009\eta ]3/2\u2212(\u22122\u2009\eta m\u22122\u2009\eta )3/2$ |

$[\zeta m,\delta /2]$ | $D3e(\eta )=[(1+\Delta )2+2\u2009\eta ]3/2\u2212(\u22122\u2009\eta m+2\u2009\eta )3/2+[(1+\Delta )2\u22122\u2009\eta ]3/2\u2212(\u22122\u2009\eta m\u22122\u2009\eta )3/2$ |

Segment . | $D(\eta )$
. |
---|---|

$[0,\zeta \u2212]$ | $D1e(\eta )=[(1+\Delta )2+2\u2009\eta ]3/2+(\u22122\u2009\eta m+2\u2009\eta )3/2\u22122\u2009[(1\u2212\Delta )2+2\u2009\eta ]3/2+[(1+\Delta )2\u22122\u2009\eta ]3/2\u2212(\u22122\u2009\eta m\u22122\u2009\eta )3/2$ |

$[\zeta \u2212,\zeta m]$ | $D2e(\eta )=[(1+\Delta )2+2\u2009\eta ]3/2+(\u22122\u2009\eta m+2\u2009\eta )3/2+[(1+\Delta )2\u22122\u2009\eta ]3/2\u2212(\u22122\u2009\eta m\u22122\u2009\eta )3/2$ |

$[\zeta m,\delta /2]$ | $D3e(\eta )=[(1+\Delta )2+2\u2009\eta ]3/2\u2212(\u22122\u2009\eta m+2\u2009\eta )3/2+[(1+\Delta )2\u22122\u2009\eta ]3/2\u2212(\u22122\u2009\eta m\u22122\u2009\eta )3/2$ |

The PD on each of the segments is obtained by integrating relation (7) with the corresponding function *D* from the boundary of the segment up to a variable upper limit $\zeta $, where the potential is equal to $\eta (\zeta )$. The resulting dependencies $\zeta (\eta )$ are given in Table II. Note that in segments 2 and 3, one of the boundaries is the point of the potential minimum, at which the derivative vanishes and its value is not included in the integrand function.

Segment . | $\zeta (\eta )$
. |
---|---|

$[0,\zeta \u2212]$ | $\zeta =\u2212\u222b0\eta {\epsilon 02+[D1e(x)\u2212D1e(0)]/(3\u2009\Delta )}\u22121/2dx=\zeta 1e(\eta )$ |

$[\zeta \u2212,\zeta m]$ | $\zeta =\zeta m\u22123\u2009\Delta \u2009\u222b\eta m\eta [D2e(x)\u2212D2e(\eta m)]\u22121/2dx=\zeta m\u2212\zeta 2e(\eta )$ |

$[\zeta m,\delta /2]$ | $\zeta =\zeta m+3\u2009\Delta \u2009\u222b\eta m\eta [D3e(x)\u2212D3e(\eta m)]\u22121/2dx=\zeta m+\zeta 3e(\eta )$ |

Segment . | $\zeta (\eta )$
. |
---|---|

$[0,\zeta \u2212]$ | $\zeta =\u2212\u222b0\eta {\epsilon 02+[D1e(x)\u2212D1e(0)]/(3\u2009\Delta )}\u22121/2dx=\zeta 1e(\eta )$ |

$[\zeta \u2212,\zeta m]$ | $\zeta =\zeta m\u22123\u2009\Delta \u2009\u222b\eta m\eta [D2e(x)\u2212D2e(\eta m)]\u22121/2dx=\zeta m\u2212\zeta 2e(\eta )$ |

$[\zeta m,\delta /2]$ | $\zeta =\zeta m+3\u2009\Delta \u2009\u222b\eta m\eta [D3e(x)\u2212D3e(\eta m)]\u22121/2dx=\zeta m+\zeta 3e(\eta )$ |

An example of solution branches for $\Delta =0.01$ is shown in Fig. 1 in dashed lines. It can be seen that the $n2$ branch smoothly transforms into the $d0$ one as $\epsilon 0$ increases. The latter, with an increase in $\epsilon 0$, first goes in the direction of decreasing $\delta $, but then turns and goes in the direction of increasing $\delta $. Thus, there is a left bifurcation point $\delta BF$, which determines the left boundary of branch $d0$ along $\delta $. To find the dependence of the point $\delta BF$ on the spread value $\Delta $, we have calculated this branch for a number of $\Delta $ values. The results are shown in Fig. 3. It can be seen that the bifurcation point shifts toward increasing $\delta $ with rising in spread value.

In Table III, $\delta BF$ point coordinates for three spread values are tabulated. It turns out that this shift quite well satisfies the law: $\delta BF(a)(\Delta )=\delta BF(0)+(10/3)\u2009\Delta $ (cf. third and sixth columns of Table III).

n . | $\Delta $ . | $\delta BF$ . | $\epsilon 0,BF$ . | $\eta min,BF$ . | $\delta BF(a)$ . |
---|---|---|---|---|---|

1 | 0.00 | 2.1444 | 1.2412 | −0.5000 | 2.1444 |

2 | 0.01 | 2.4835 | 1.2579 | −0.494 25 | 2.4777 |

3 | 0.02 | 2.6116 | 1.2283 | −0.4886 | 2.6158 |

n . | $\Delta $ . | $\delta BF$ . | $\epsilon 0,BF$ . | $\eta min,BF$ . | $\delta BF(a)$ . |
---|---|---|---|---|---|

1 | 0.00 | 2.1444 | 1.2412 | −0.5000 | 2.1444 |

2 | 0.01 | 2.4835 | 1.2579 | −0.494 25 | 2.4777 |

3 | 0.02 | 2.6116 | 1.2283 | −0.4886 | 2.6158 |

### B. $d00$ branch

The $d00$ branch corresponds to PD in which, unlike the $d0$ branch, the *p*-VE is located to the left of the *e*-VE [cf. Fig. 2(b), black curve, and Fig. 2(a), red curve]. This branch lies entirely within the negative $\epsilon 0$ region. It starts at the point where $n2$ branch ends. Here, the potential value at the maximum $\eta M=\eta +=(1\u2212\Delta )2/2$ is the smallest for the entire branch and corresponds to the boundary of the mode with reflection. Further, as $\epsilon 0$ increases, $\eta M$ value grows, while remaining less than $(1+\Delta )2/2$, up to the point where $\epsilon 0$ vanishes, and the $d00$ branch goes into $d11$ branch.

Considering the $d0$ branch is also sufficient to find the PD within the interval $(0,\delta /2)$. Let us divide this interval into three segments: (1) $[0,\zeta M]$, here the potential varies within $[0,\eta M]$, in this segment only positrons that have overcome the potential barrier occur; (2) $[\zeta M,\zeta +]$, here the potential varies within area $[\eta M,(1\u2212\Delta )2/2]$, a flow of reflected positrons is formed in this segment; and (3) $[\zeta +,\delta /2]$, here the potential varies within area $[(1\u2212\Delta )2/2,0]$, in this segment there are two positron flows: emitted from the right electrode and reflected. The positron density within these three regions is determined by the first, second, and third of formula (6), respectively. All three segments contain both direct and reflected electrons, whose densities are determined by the first of formula (5). Forms of the function $D(\eta )$ for the segments are given in Table IV.

Segment . | $D(\eta )$
. |
---|---|

$[0,\zeta M]$ | $D1p(\eta )=[(1+\Delta )2+2\u2009\eta ]3/2+(2\u2009\eta M+2\u2009\eta )3/2\u22122\u2009[(1\u2212\Delta )2+2\u2009\eta ]3/2+[(1+\Delta )2\u22122\u2009\eta ]3/2\u2212(2\u2009\eta M\u22122\u2009\eta )3/2$ |

$[\zeta M,\zeta +]$ | $D2p(\eta )=[(1+\Delta )2+2\u2009\eta ]3/2+[2\u2009(\eta M+\eta )]3/2\u22122[(1\u2212\Delta )2+2\u2009\eta ]3/2+[(1+\Delta )2\u22122\u2009\eta ]3/2+(2\u2009\eta M\u22122\u2009\eta )3/2$ |

$[\zeta +,\delta /2]$ | $D3p(\eta )=[(1+\Delta )2+2\u2009\eta ]3/2+[2\u2009(\eta M+\eta )]3/2\u22122[(1\u2212\Delta )2+2\u2009\eta ]3/2+[(1+\Delta )2\u22122\u2009\eta ]3/2+[(2\u2009(\eta M\u2212\eta )]3/2\u22122[(1\u2212\Delta )2\u22122\u2009\eta ]3/2$ |

Segment . | $D(\eta )$
. |
---|---|

$[0,\zeta M]$ | $D1p(\eta )=[(1+\Delta )2+2\u2009\eta ]3/2+(2\u2009\eta M+2\u2009\eta )3/2\u22122\u2009[(1\u2212\Delta )2+2\u2009\eta ]3/2+[(1+\Delta )2\u22122\u2009\eta ]3/2\u2212(2\u2009\eta M\u22122\u2009\eta )3/2$ |

$[\zeta M,\zeta +]$ | $D2p(\eta )=[(1+\Delta )2+2\u2009\eta ]3/2+[2\u2009(\eta M+\eta )]3/2\u22122[(1\u2212\Delta )2+2\u2009\eta ]3/2+[(1+\Delta )2\u22122\u2009\eta ]3/2+(2\u2009\eta M\u22122\u2009\eta )3/2$ |

$[\zeta +,\delta /2]$ | $D3p(\eta )=[(1+\Delta )2+2\u2009\eta ]3/2+[2\u2009(\eta M+\eta )]3/2\u22122[(1\u2212\Delta )2+2\u2009\eta ]3/2+[(1+\Delta )2\u22122\u2009\eta ]3/2+[(2\u2009(\eta M\u2212\eta )]3/2\u22122[(1\u2212\Delta )2\u22122\u2009\eta ]3/2$ |

Segment . | $\zeta (\eta )$
. |
---|---|

$[0,\zeta M]$ | $\zeta =\u222b0\eta {\epsilon 02+[D1p(x)\u2212D1p(0)]/(3\u2009\Delta )}\u22121/2dx=\zeta 1p(\eta )$ |

$[\zeta M,\zeta +]$ | $\zeta =\zeta M+3\u2009\Delta \u222b\eta \eta M[D2p(x)\u2212D2p(\eta M)]\u22121/2dx=\zeta M+\zeta 2p(\eta )$ |

$[\zeta +,\delta /2]$ | $\zeta =\zeta ++3\u2009\Delta \u222b\eta (1\u2212\Delta )2/2{D3p(x)\u2212D3p[(1\u2212\Delta )2/2]+D2p[(1\u2212\Delta )2/2]\u2212D2p(\eta M)\u2009}\u22121/2dx=\zeta ++\zeta 3p(\eta )$ |

Segment . | $\zeta (\eta )$
. |
---|---|

$[0,\zeta M]$ | $\zeta =\u222b0\eta {\epsilon 02+[D1p(x)\u2212D1p(0)]/(3\u2009\Delta )}\u22121/2dx=\zeta 1p(\eta )$ |

$[\zeta M,\zeta +]$ | $\zeta =\zeta M+3\u2009\Delta \u222b\eta \eta M[D2p(x)\u2212D2p(\eta M)]\u22121/2dx=\zeta M+\zeta 2p(\eta )$ |

$[\zeta +,\delta /2]$ | $\zeta =\zeta ++3\u2009\Delta \u222b\eta (1\u2212\Delta )2/2{D3p(x)\u2212D3p[(1\u2212\Delta )2/2]+D2p[(1\u2212\Delta )2/2]\u2212D2p(\eta M)\u2009}\u22121/2dx=\zeta ++\zeta 3p(\eta )$ |

### C. $d11$ branch

The $d11$ branch continuously extends $d00$ one into the region of positive $\epsilon 0$ values. It includes PDs whose minimum $\eta m1$ is located near the left boundary. Its absolute value is less than the absolute minimum $|\eta m|$ [Fig. 2(b), red curve]. To the right of the minimum, the potential grows with increasing coordinate and reaches a maximum $\eta M$. The *p*-VE is located here. Then the potential decreases, reaches absolute minimum $\eta m$, and so on. As well as on the $d00$ branch, electron reflection occurs to the right of the middle of the gap, i.e., at $\zeta >\delta /2$.

Within the segment ( $0,\zeta M$), Eq. (12) with the function $D(\eta )\u2261D1p(\eta )$ is still valid. Moreover, in the region lying between the point behind the local minimum $\zeta 0$, where the potential vanishes, and $\delta /2$, the equations and formulas for the PD given above for the $d00$ branch remain valid. The $\epsilon 0$ value included in them now makes sense of the field strength at point $\zeta 0$. However, the PD in the segment $(0,\zeta 0)$ is symmetrical with respect to the local minimum point $\zeta m1$; therefore, the field strength at the point $\zeta 0$ is equal to $\u2212\epsilon 0$, and the value of $\epsilon 02$ in formulas describing the PD within the interval $(\zeta 0,\delta /2)$ does not change. Therefore, it is only sufficient to find PD within the region $(0,\zeta 0)$.

Varying $|\eta m1|$ from 0 (with $\epsilon 0=0$) to the maximum value determined by formula (16), we construct $d11$ branch. Figure 1 shows the solution branches for diodes with spreads $\Delta =0$ and $\Delta =0.01$ in the region of $\delta <6$. Comparison of solid and dashed curves shows that a slight spread of the VDF leads to a shift of these dependencies to the right and a slight stretch along the $\delta $ axis.

## III. NUMERICAL CALCULATIONS OF INSTABILITY DEVELOPMENT

When researching stationary solution stability features in electron–positron diode, we study the evolution of a small perturbation of the stationary electric field distribution. If the solution turns out to be unstable, the perturbation evolves to a solution different from the steady state one. In the opposite case, the process returns to the initial state. For simulation, we used two different numerical codes: PIC-code (implementing calculations by the particle method) and EK-code.

The PIC-code simulates the movement of a huge number of electrons emitted from the left boundary and positrons emitted from the right boundary, and moving under electric field set at the grid nodes. The grid consists of $N\zeta $ nodes with an equal distance $h\zeta $ between them. The charge density at the nodes is calculated in accordance with the “cloud-in-cell” model (linear particle contribution to the density at neighboring nodes),^{9} unit density corresponds to $N0$ particles in a cell. The electric field at the nodes is found as a solution of the Poisson equation, and between the nodes the field is approximated by a linear function.^{9} The position of a particle at each subsequent time moment is calculated by the “leapfrog” method,^{9} and the time step $h\tau $ is chosen to be constant. At each step, particles that reach any electrode are removed from the calculation, and particles that fly out from the electrodes with random velocities distributed in accordance with a given VDF are added. In the calculations described below, a uniform velocity distribution at the electrodes was specified in the interval $[1\u2212\Delta ,1+\Delta ]$ for electrons and in the interval $[\u2212(1+\Delta ),\u2212(1\u2212\Delta )]$ for positrons. The moment of particle departure from the electrode was considered random and uniformly distributed within the time step.

The numerical algorithm implemented in the EK-code is described in detail in the papers,^{8,10} where it is modified to calculate processes in a diode with counter flows of negatively and positively charged particles. At each time step, the charged particle VDF, their densities, and electric field distribution are sequentially determined. The particle VDF at a given point $(\tau i,\zeta k)$ is built by calculating a set of trajectories ending at this point. The computation is carried out backward in time from $\tau i$ up to the moment $\tau 0$ when the particle left the electrode. In this case, the electric field distributions at all previous times ( $\tau <\tau i$) are considered known. The field distribution at each time moment is found from the Poisson equation, into the right side of which the electron and positron densities calculated for this time moment are substituted. Since the distribution of the electric field is unknown at the last time interval ( $\tau i\u22121,\tau i$), its extrapolation to this interval is used when the densities are calculated. Then iterations are carried out with alternate recalculation of the field and particle densities.

When a small perturbation evolution was calculated by the EK-code, the electric field distributions were set equal to the stationary distribution over a sufficiently long period of time $\tau s$, so that particle density distribution equal to the stationary one had time to adjust in the inter-electrode gap. Then, at the time step immediately following the moment $\tau =\tau s$, the time-dependent self-consistent computation of the field and particle distributions starts. In this case, at the moment $\tau >\tau s$, a small random perturbation of the field distribution arose due to the finite accuracy of the calculations. However, sometimes another initial condition is used: during the time interval $(0,\tau s)$, an exponentially growing sinusoidal perturbation with an amplitude much less than unity was superimposed on the stationary field distribution.

When using the PIC-code, in some cases, the field distribution within the diode was also set equal to the stationary distribution during the time interval $(0,\tau s)$, and after that the time-dependent self-consistent calculation started. In other cases, the self-consistent calculation started at the time $\tau =0$, while at each point in the gap, a stationary field distribution (with or without a perturbation imposed on it) and the initial particle VDFs corresponding to the initial field distribution were specified. At each point, the initial VDF was of rectangular distribution form: the particles were uniformly distributed over velocities within the interval, the boundaries of which were determined from the stationary law of energy conservation.

First, we investigated the stability features of solutions belonging to $d00$ and $d11$ branches. In these calculations, the velocity spread was $\Delta =0.01$. The $d11$ branch continues the $d00$ branch as the inter-electrode gap length $\delta $ increases. Diodes with $\delta $ values equal to 2.7, 2.8, 3.0, and 3.5 were considered for the $d00$ branch, and $\delta =$ 5.0 and 5.6 were considered for the $d11$ branch. These values cover almost the entire range of possible diode lengths for which solutions of these types exist. Note that at $\delta =5.6$ there are two stationary solutions belonging to the branch $d11$, which differ in the value of $\epsilon 0$ (see Fig. 1). In all the cases studied, the solutions turned out to be unstable.

Figure 4 shows the dependencies of the potential maximum value $\eta M$ and the potential perturbation amplitude $\eta \u0303M$ on time at the initial stage of the process, obtained as a result of the EK-code modeling. The dependencies in Figs. 4(a) and 4(b) relate to the $d00$ branch. In all cases, $\eta M(\tau )$ oscillates, and the oscillation amplitude grows with time. During the oscillation process, the PD shape is rearranged several times, wherein the absolute maximum of the distribution is located either to the left or to the right of the diode middle. The growth of the perturbation amplitude $\eta \u0303M$ for values of the order of 0.01–0.1 is close to exponential: for $\delta =3.5$, the envelope $\eta \u0303M$ is a straight line on a logarithmic scale, and for smaller values, length dependence is less regular [Fig. 4(b)]. Figures 4(c) and 4(d) show the dependencies of $\eta M(\tau )$ and $\eta \u0303M(\tau )$ for solutions corresponding to $d11$ branch. It can be seen that at the initial stage oscillations $\eta M(\tau )$ occur, and the envelope of the dependence $\eta \u0303M(\tau )$ resembles an exponential.

These conclusions were confirmed by PIC-code calculations performed for the same diode lengths with varying shapes of potential perturbation at the initial time (sinusoidal or random due to errors in the numerical calculation). In all cases, the PD turned out to be unstable, and increasing oscillation amplitudes were observed at the initial stage of evolution; however, shape and amplitude of the oscillations strongly depended on the initial perturbation. Nevertheless, in all cases, the oscillation amplitude increased more or less exponentially with similar exponents. Figure 5 shows the PIC-code results for the initial stage of the PD evolution for the solutions belonging to the $d00$ branch for $\delta =3.5$ and $d11$ one for $\delta =5.0$.

The time dependence of the field strength at the left electrode $\epsilon (\tau ,0)$ does not coincide with the dependencies $\eta M(\tau )$, $\eta \u0303M(\tau )$, which is natural, since at this stage the distribution form is continuously restructuring. Note that different scenarios for the development of the initial perturbation are possible, so we can only state the instability of the solutions studied, but their change over time under various forms of perturbation may differ from the dependencies in Fig. 4.

In addition, the stability features of the solution belonging to the lower part of the $d0$ branch (below the bifurcation point) were studied at $\delta =2.55$, $\epsilon 0\u22481.027$ for $\Delta =0.01$. The solution turned out to be unstable; the process of the perturbation evolution ended with a uniform distribution of the potential.

Based on the calculation results, we determined the growth rates $\Gamma $ and the oscillation frequencies $\omega $ in Eq. (20). The $\Gamma (\delta )$ dependencies obtained using the EK-code are presented in Fig. 7(a). The oscillation frequency changes slightly in the studied range of diode lengths, since it is determined by the characteristic time of flight of charged particles through the diode.

The stability threshold $\delta thd$ (the value of $\delta $ at which the growth rate $\Gamma $ changes sign) shifts to the right as $\Delta $ increases. The dependence of $\delta thd$ on $\Delta $ for small $\Delta $ can be approximated by a linear function: $\delta thd=2.6+30\u2009\Delta $ [Fig. 7(b)]. Note that we determine the position of the stability threshold with an accuracy of 0.05.

The calculations have shown that in diodes with a length $\delta $ exceeding the stability threshold, the exponential growth of the oscillation amplitude continues to values of the order of several hundredths. After this, the increase in amplitude begins to slow down, and eventually the oscillations turn into periodic ones. At values of $\delta $ being near the stability threshold, the dependencies $\eta (\tau )$ and $\zeta (\tau )$ are close to sinusoidal; with increasing $\delta $, the shape of the oscillations is distorted.

We built the dependencies $\eta M(\zeta M)$ (the maximum value on the PD $\eta M$ and its position $\zeta M$ depend on time as a parameter). When the solution is unstable, this dependence is an unwinding spiral (with turns of irregular shape) at the stage of exponential growth of the oscillation amplitude. Then the distance between the turns gradually decreases, and, in the end, the curve reaches a closed loop that corresponds to steady-state nonlinear oscillations. Such loops are shown in Figs. 8(a) and 8(b) for $\delta =3$ and 3.5.

Note that a number of other time-dependent processes, which we also simulated, also end in nonlinear oscillations. In particular, at $\delta =3.0$ and 4.0, we considered an evolution of the perturbation of steady states corresponding to the $n2$ branch for various types of perturbation. In one of the possible scenarios for the development of the perturbation the process ended with nonlinear periodic oscillations in the vicinity of the $d0$ branch.^{7} The same oscillations were established in calculations with a zero initial condition, when it was assumed that before the start of the calculation, there were no charged particles in the diode. The cycle of steady-state oscillations of the PD maximum at $\delta =4.0$ in the case of the absence of particles in the diode at the initial time is shown in Fig. 8(c).

We plotted the dependence of the amplitude of steady-state oscillations on the diode length for two VDF half-width values: $\Delta =0.01$ and 0.015. In the vicinity of the threshold, the amplitude of oscillations grows with increasing $\delta $. Figure 9 shows how the value of the oscillation amplitude $\eta Mmax(\delta )\u2212\eta Mmin$ increases as $\delta $ moves away from the stability threshold $\delta \u2212\delta thd$. On a logarithmic scale, it is clear that this is a power-law dependence with an exponent close to 0.5. This indicates that the position of the stability threshold $\delta thd$ is the Hopf bifurcation point.

In the paper,^{11} where nonlinear oscillations in a diode with flows of electrons and ions were studied, it has been proven that long-lived ions appear in the plasma. These ions “live” for several periods of plasma oscillation, being, in a way, captured by an oscillating ionic potential barrier that periodically appears on their way. It is natural to assume that similar processes should occur in a diode with electron–positron plasma, and long-lived particles of both signs can also exist. To verify this assumption, we plotted the electron and positron distribution functions in terms of velocities and coordinates at different times and made sure that, as in the case of electron–ion plasma, the VDF at each point of the diode is a set of narrowing velocity intervals with a constant value of the VDF, outside of which this function is equal to zero. Over time, the number of intervals is getting larger and larger, which means the capture of charged particles into dynamic potential wells formed by moving PD extrema. Figure 10 shows the positron VDF at time moments corresponding to different points of the cycle shown in Fig. 8(b). On the scale used, one can distinguish only the widest intervals; examples of intervals structure at a fixed coordinate are shown in the insets. On panel (a), the interval *1a* splits to *4a, 5a* at $\zeta =0.86$, and the interval *6a* splits to *7a, 8a* at $\zeta =2.87$. On panel (b), the interval *1b* splits to *8b,9b* at $\zeta =2.47$. Remember that the same numbers on panels (a) and (b) designate different intervals, we did not try to set up a correspondence between the panels.

The presence of several narrow velocity intervals is illustrated in Tables VI and VII. The tables show the boundaries of the positron velocity intervals in the cross sections $\zeta =0.42$ and $\zeta =3.22$ at the bottom point of the cycle shown in Fig. 8(b). The interval numbers in the tables correspond to curve numbers in Fig. 10(b). At $\zeta =3.22$, the interval *5b* is too narrow to determine its boundaries correctly ( $\u223c10\u22128$), so we do not present their values in Table VII. The interval widths and the distance between intervals decrease exponentially with the velocity increase, except for the lowest intervals, where particles cannot yet be considered as long-lived. The same dependence had been observed in the case of long-lived ions.^{11}

n . | $umin$ . | $umax$ . | $umax\u2212umin$ . |
---|---|---|---|

1b | −1.653 819 8 | −1.649 713 1 | 0.004 106 711 9 |

2b | −1.600 249 4 | −1.597 954 6 | 0.002 294 804 7 |

3b | −1.597 048 1 | −1.596 887 6 | 0.000 160 494 1 |

4b | −1.596 670 0 | −1.596 649 5 | 0.000 020 542 3 |

5b | −1.596 641 5 | −1.596 640 4 | 0.000 001 120 4 |

n . | $umin$ . | $umax$ . | $umax\u2212umin$ . |
---|---|---|---|

1b | −1.653 819 8 | −1.649 713 1 | 0.004 106 711 9 |

2b | −1.600 249 4 | −1.597 954 6 | 0.002 294 804 7 |

3b | −1.597 048 1 | −1.596 887 6 | 0.000 160 494 1 |

4b | −1.596 670 0 | −1.596 649 5 | 0.000 020 542 3 |

5b | −1.596 641 5 | −1.596 640 4 | 0.000 001 120 4 |

n . | $umin$ . | $umax$ . | $umax\u2212umin$ . |
---|---|---|---|

6b | −0.462 520 86 | −0.407 972 07 | 0.054 548 798 |

7b | 0.163 014 80 | 0.221 614 12 | 0.058 599 325 |

8b | 0.256 868 25 | 0.266 720 51 | 0.009 852 254 |

9b | 0.270 251 44 | 0.270 606 68 | 0.000 355 231 |

2b | 0.271 069 52 | 0.271 114 51 | 0.000 044 993 |

3b | 0.271 133 33 | 0.271 135 76 | 0.000 002 426 |

4b | 0.271 139 06 | 0.271 139 54 | 0.000 000 474 |

n . | $umin$ . | $umax$ . | $umax\u2212umin$ . |
---|---|---|---|

6b | −0.462 520 86 | −0.407 972 07 | 0.054 548 798 |

7b | 0.163 014 80 | 0.221 614 12 | 0.058 599 325 |

8b | 0.256 868 25 | 0.266 720 51 | 0.009 852 254 |

9b | 0.270 251 44 | 0.270 606 68 | 0.000 355 231 |

2b | 0.271 069 52 | 0.271 114 51 | 0.000 044 993 |

3b | 0.271 133 33 | 0.271 135 76 | 0.000 002 426 |

4b | 0.271 139 06 | 0.271 139 54 | 0.000 000 474 |

## IV. CONCLUSION

We investigated the stability features of main steady state types in a diode with electrons and positrons counter flows.^{3,6,7} In the latter two papers, this problem was studied for the regime without reflection of charged particles from potential extrema, i.e., for $nk$ branches with $k=0,\u20092,\u20094,\u2009\u2026$. It was established analytically and confirmed numerically that for homogeneous solutions there is a threshold for the dimensionless inter-electrode gap $\delta $ (or, what is the same, for the current density), above which instability develops in the diode. The threshold $\delta th0$ turned out to be equal to $2\u2009\pi \u2009\lambda D$, which is $2$ higher than the known Pierce one.^{12} On the other hand, all inhomogeneous stationary solutions turn out to be unstable with respect to small perturbations. In this case, the time-dependent process either ends in a new state, characterized by nonlinear oscillations that occur around the stationary solution corresponding to the $ds$ branches, or reaches a homogeneous solution (if $d<2\u2009\pi \u2009\lambda D$).

In the paper, we studied the stability features of steady states in the regime with particle reflection from potential barriers, i.e., solutions corresponding to $d0$, $d00$, and $d11$ branches. Since the studies were carried out numerically, not a $\delta $-shaped VDF, but a distribution function with a small spread $\Delta $ in velocities was used. All stationary solutions have been described. They qualitatively coincided with that corresponding to the $\delta $-shaped VDF. Distinction is observed only in the vicinity of the left bifurcation point (see Fig. 3). A law has been found according to which the bifurcation point shifts along $\delta $ as the spread $\Delta $ increases.

The calculations of stationary solution small perturbation development were performed. It has been established that the solutions for the $d00$ and $d11$ branches are unstable, while solutions from the $d0$ branch turn out to be stable when the inter-electrode gap is less than a certain threshold value $\delta thd$. With increasing half-width $\Delta $ of the VDF, the position of the stability threshold $\delta thd$ shifts toward larger diode lengths, and the dependence $\delta thd(\Delta )$ is close to a linear function at a small $\Delta $. The nonlinear stage of instability development ends with nonlinear oscillations in the vicinity of the $d0$ branch. The process of instability development for solutions corresponding to branches without particle reflection can (under certain conditions) terminate by the same oscillations.

Studying the discovered nonlinear oscillations has shown that for the studied diode lengths oscillation amplitude grows with increasing in $\delta $, and in the vicinity of the threshold, the increasing law turns out to be close to the dependence $(\delta \u2212\delta thd)1/2$, i.e., the Hopf bifurcation takes place. The calculations of nonlinear oscillations have shown that the phase plane for charged particles of both types has a rather complex configuration and demonstrates the appearance of both long-lived electrons and positrons, i.e., particles that oscillate in the vicinity of their VEs during several oscillation periods until they fly out onto an electrode. In the literature on studying processes in electron–positron plasma, we have not found any mention of the long-lived particles.

To date, an analytical theory of diode stability has not been created for the regime with charged particle reflection from potential extrema. Therefore, during our research, we used two numerical codes: EK-code and PIC-code. Both codes made it possible to determine the main eigenmodes from the calculations, i.e., modes that have the largest growth rate in instability case and the smallest (in absolute value) dumping factor in dumping one. The growth rates obtained by both codes have completely coincided with each other. Thus, we have demonstrated the validation of the results obtained.

In further studies of electron–positron plasma diodes, we intend to consider the relativistic regime. It should more adequately describe, for example, the behavior of plasma in a pulsar diode.^{2} Although the main features of the processes in a diode are well described by the non-relativistic mode considered.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**L. A. Bakaleinikov:** Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). **V. I. Kuznetsov:** Conceptualization (lead); Formal analysis (lead); Investigation (equal); Methodology (lead); Project administration (lead); Supervision (lead); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). **E. Yu. Flegontova:** Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **D. P. Barsukov:** Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **I. K. Morozov:** Investigation (supporting); Software (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (supporting).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*Computer Simulation Using Particles*