Computer Simulation of Explosive Emission Processes in Strong Electromagnetic Fields

The paper discusses the problem of computer simulation of explosive emission processes in strong electromagnetic fields. In this work, a numerical method for calculating electron emission from the surface of a metal cathode for the axial symmetric geometry of a technical system is proposed. For modelling, the particle method and grid calculations of electromagnetic fields based on Maxwell’s equations are applied. The approach uses the representation of large smoothed Gaussian particles and implements calculations of the electromagnetic fields on Cartesian spatial grids. The software realization is directed towards on parallel computing. The calculation of the emission process in a coaxial cylindrical diode with magnetic self-insulation with the use of developed technique was carried out. The process is considered for two situations: with and without the longitudinal plasma layer. The calculations show that the presence of plasma leads to an increase in the emission current. This result is consistent with the results of experiments.


Introduction
Technical, environmental, economic and other systems that are studied by modern science most commonly use the methodology of mathematical and computer modeling.This is because a direct full scale experiment, due to many factors and many options, can be long, expensive, often either dangerous or impossible.The methodology of the mathematical and computer modeling allows us to understand the essence of a separate object or system, study their structure and internal connections, calculate the basic properties and laws of development, manage them with the help of external influences, predict the consequences of these influences [1,2].
For development of complex technical systems, the mathematical and computer modeling is used in a number of areas, including creation more efficient and promising designs.Based on extensive experimental information and developed theoretical methods, mathematical and computer modeling in some cases is a catalyst for the creation of new technologies and industries.Recently, the mathematical and computer modeling has been aimed not only at the quantitative description of processes in existing natural and artificial materials [3,4], but also at the creation of new materials with specified properties [5].One of the technologies for obtaining such materials is surface treatment with ion, electron and laser beams.For this purpose, electrical installations are created that generate one or another type of energy carriers.
A very promising direction at present is the creation of microwave radiation generators with a large controlled power in a given frequency range.One of the designs of such installations is a coaxial cylindrical diode with magnetic self-isolation [6].In such a device, the formation of a relativistic electron beam is produced by the explosive emission [7] of electrons from the surface of a cylindrical cathode.Under conditions of a strong electromagnetic field, the emitted electrons form a relativistic electron beam (REB).A part of the energy of this beam is transformed into radiation under interacting with a pre-created plasma layer.Explosive emission is accompanied by the formation of plasma on the cathode surface, consisting of atoms of the cathode material and gas molecules stored on its surface (for example, hydrogen).In some cases, a mechanism for the formation of plasma based on inert gases (for example, argon) is specially added to the system.Selection of design parameters and taking into account of plasma formation processes allows us to obtain the specified characteristics of the device.
The purpose of this work is to present a mathematical model, a numerical algorithm for its implementation and simulation results related to emission calculations in coaxial cylindrical diodes with magnetic self-isolation.The focus of attention in the work is the technique of reproducing argon plasma in a numerical algorithm and changing the emission parameters in the presence of a plasma layer.
The second section of the paper describes the computational geometry of the model problem and the basic equations.The third section discusses a numerical technique combining the grid approach and the particle method, as well as a parallel software implementation.The fourth section is devoted to verification of the numerical model of plasma.Moreover, the results of calculations of the emission process with and without a plasma layer are presented here.

Technical system and basic equations
To create conditions for the occurrence of explosive emission, it is required to apply a strong electric field to the metal surface (cathode).In this case, large currents of the kilo ampere range are generated.One of the areas of application of such REBs is the development and creation of microwave generators based on coaxial diodes with magnetic self-isolation.In this paper, we consider one of the possible designs of such generator, shown in figure 1.It is assumed that the system is in a homogeneous magnetic field with moderate induction z B .Electron emission is initiated by a TEM-wave with a sufficiently large amplitude 0 E coming from the left boundary.The emerging relativistic electron beam interacts with the plasma layer generated at the preliminary stage by the anode-collector system.
A mathematical model describing the processes of REB formation and its further interaction with plasma may include various components depending on the characteristic time of the processes and geometric scales.In the case that we consider, the size range is measured in tens of centimeters, and the time range is units of nanoseconds.In this situation, the hydrodynamic approximation of the processes is not suitable due to the smallness of the time range.Therefore, we suggest on Maxwell's description of the evolution of electromagnetic fields, supplemented by the dynamics of charged particles under the action of the Lorentz force., .
Here D and B -electric and magnetic induction vectors, E and H -electric and magnetic field is the bulk charge density, negatively and positively charged particles, 0 q   and 0 q   -particle charges, n  and n volume concentrations of particles,     j j j is the total current density of the particles, To describe the dynamics of charged particles, the method of large smoothed particles applied in electrodynamic models (LSPE) is proposed.In this method, the equations of relativistic dynamics are used instead of the classical equations of continuity and momentum.These equations are written for separate charged particles that are clouds of ions and electrons.The equations have the following form: , .
q  and m  -radius-vector, momentum, velocity, charge and mass -the Dirac delta function describing the charge density of a particle of the sort (e.g., " "    -for plasma electrons, " "    -for positively charged plasma ions), N  is the number of particles of the kind  .
The systems of equations ( 2) and ( 3) are supplemented by the necessary boundary and initial conditions.In particular, the static parts of the electric and magnetic fields are set by stationary distributions and do not participate in the calculations.The dynamic parts of the electric and magnetic fields are absent at the initial moment.A TEM-wave is set on the left boundary: The following conditions are fulfilled at the boundaries of the conductors: On free boundaries, the conditions take the form: The initial conditions for particles are set to zero, excepting the situation when the generator initially contains a layer of plasma.In this case, some initial distribution of positively and negatively charged particles is given.We discuss this issue in more detail in section 4.
The boundary conditions on the emission surface for negatively charged particles that reappear as a result of emission have the form: Here ,k q  is the charge of a new negatively charged particle associated with an element of the emission surface (emission center) of an area S  , e is electron charge, e C is constant characterizing the emitter material, ,k  v and ,k  p are the initial velocity and momentum of a new large particle.

Numerical algorithm and parallel realization
Numerical implementation of the proposed model was carried out for the case of axial symmetry in coordinates (r, z).It is assumed that the condition of independence of the fields from the azimuthal variable  is fulfilled.The numerical implementation of equations ( 2) is based on the grid method of finite volume [8] on a curved Cartesian grid.He combined the explicit FDTD scheme [9,10] for solving Maxwell's dynamic equations and the well-known "cross" scheme for solving the Poisson equation.The Poisson equation is solved by the Jacobi iterative method with a diagonal transition matrix [11].Equations (3) are numerically implemented using the symmetric Adams scheme of the second order of accuracy [11][12][13].The description of the particle model requires special consideration.We assume that the charged particle in our model is a homogeneous cloud of ions (let's call it a "multi-ion") or electrons (let's call it a "multi-electron").This helps to describe REB electrons, ions and plasma electrons in a unified way.In the limit, when the charge tends to 1, these formations turn into ordinary ions and electrons.To describe the shape of the multi-ion and multi-electron, the following Gaussian representation for the probability density of particle location at the point ,k  r is used, approximating the Dirac function , ; , , A  depends implicitly on time and it is selected from the condition: The main motivation for our choice of particle model is the fact that the Gaussian representation under certain conditions is an exact solution of the equations of motion of the particles (3).
A parallel program was created based on the developed numerical methodology.The ANSI C and C++ programming languages, MPI and Open MP parallel computing standards were used for implementation of parallel code.The parallelization technology is based on the Domain Decomposition Method (DDM) [14,15] and Dynamic Load Balancing algorithms (DLB) of calculators [16].
DDM technology involves splitting a Cartesian cylindrical grid into compact domains of the same power.During calculations, each domain is attached to a specific MPI process, which sequentially performs all calculations according to the four blocks of the dynamic algorithm.At each fixed moment of time, not all domains can be occupied by particles.As a result, the distribution of grid domains is carried out only among MPI processes.When a certain number of particles enter a particular domain, they are processed by a group of parallel CPU threads organized with using the Open MP standard.If there are few particles, one thread is used.Otherwise, such a number of threads is added that provides approximate load equality among MPI processes.
Since particles dynamically move from one domain to another (that is, they move to another MPI process), it becomes necessary to correct the computational load.This correction is implemented by the DLB algorithm.With its help, at each step of the calculations, the running times of each MPI process are analyzed, the number of particles changes, the trajectories and shapes of which will be calculated on the next step.The algorithm assumes a smooth change of these characteristics (that is, it belongs to the class of diffusion methods of dynamic load balancing of calculators) and has confirmed its effectiveness.

Results
In this section, we consider two issues.The first one is associated with the clarification of the parameters of the numerical model of plasma.The second matter is devoted directly to the results of simulation of emission processes.
The analysis of the numerical model has already been carried out in [17].However, only the main properties of the numerical algorithm related to the quality of the implementation of the FDTD scheme and the particle model were considered there (8).The result of the consideration is the fact that the basic physics of the emission process in the framework of the proposed numerical implementation is transferred correctly at a qualitative level.By choosing the emission parameter e C in formula (7), numerical algorithm for a specific emitter material can be customized.As a result, it becomes possible to obtain not only a qualitative, but also a quantitative coincidence of the emission current with its values obtained by a full-scale experiment.
A similar consideration must be performed for calculating the emission with plasma layer.However, in this case it is necessary to clarify the parameters of the numerical model of plasma.These parameters are a grid resolution (spatial steps r h , z h and time step  ) and plasma parameters in a grid cell (the number of positively charged plasma centers -"multi-ions" inside the cell i M , and the number of negatively charged particles in the environment of the multi-ion -"multi-electrons" e M , the distance among a multi-ion and surrounding multi-electrons is ie R ).
To analyze the situation, a simplified model of design is shown in figure 2. It contains a plasma clot in the form of a square (see figure 2 on the left).In addition, in the figure 2, the initial arrangement of multi-ions and multi-electrons is given for  Let us now consider the choice of a criterion that will be used to calibrate the numerical model of plasma.The main problem of our model is to put data on the concentration and temperature of the plasma at equilibrium into it.For example, if it is known that a plasma layer is created in the generator system under development at the first stage, then it is possible to determine the density and temperature of the plasma from the parameters of its spectrum.In particular, the typical characteristics for argon plasma are the concentration of the order Here cell V is the three-dimensional volume of a grid cell, depending on its location along the radius.
The determination of the charges of specific multi-ions and multi-electrons for a fixed grid cell is made by the formulas: , , / , / ( ).
The plasma temperature can be put into our model by setting the arrangement of multi-ions and multi-electrons.Each specific configuration gives the system of a certain potential energy stored initially in an electric field.This potential energy must be related to the kinetic energy of plasma multi-ions and multi-electrons.The kinetic energy will appear as a result of plasma evolution in the absence of external influences after some time.
We will assume that in the system under development, the plasma created at the first stage is ideal and nonrelativistic.However, the equilibrium in it is established in some finite time.Then, near the equilibrium point, we can assume that the plasma temperature obeys the usual laws of thermodynamics [18]: v -masses and speeds of particles, p m and p v -average masses and particle velocities.Within the framework of the one- temperature plasma model, formulas (12) relate the velocities of multi-electrons and their total temperature.Let us consider the potential energy of the electric field stored in the initial arrangement of multiions and multi-electrons.It can be defined by considering the interaction energy P E of one multi-ion with one multi-electron [19]: e l e l e l e e q q E r r z z l Here i q and e q -multi-ion and multi-electron charges, 0  is the vacuum permittivity, , e l r -two extreme positions of a multi-electron when moving around a multi-ion,   0 0 , r z is the multi-ion position.The magnitude  when moving along the axis z can be estimated as follows: Here  actually coincides with the radius ie R , 0 G -geometric parameter.Similarly, we can consider the motion of a multi-electron along the transverse coordinate.As a result, it is possible to relate the temperature of the average multi-electron in quasi-equilibrium with the initial potential energy using the approximate expression: Let's move on to a numerical analysis of the calculation results of plasma equilibrium for various starting parameters.To do this, we consider the model geometry is shown in figure 3.There is no voltage on the electrodes (marked in red and blue in the figure).The plasma concentration is  We will vary the parameters of the space-time grid and the parameters of the plasma model in order to obtain the temperature of electrons close to ps (this ratio is chosen from the condition of calculating electromagnetic fields with amplitudes about 500kV/cm and 1Tl in the nanosecond time range).The results of calculating the temperature dynamics of multi-ions are shown in figure 4. They show that consistent mesh refinement on space and time leads to an increase in the stability range of the calculation, but only up to a certain limit.Further refinement is inexpedient, since it does not lead to the desired result.In addition, mesh refinement leads to a non-linear increase in computational costs.In particular, if the meshes are refined by 2 times, an 8-fold increase in computational costs is expected.These costs increase by more than 10 times, due to the presence of a calculation block associated with plasma particles.The second series of calculations is related to the analysis of the effect of plasma parameters on the stability of calculations in case of the fixation of grid parameters , at which the equilibrium temperature of the ions does not exceed 300 K during a time of at least 20 ns.
The third series of calculations is devoted modeling the development of REB with and without plasma.For this, the parameters are set: the amplitude of the TEM wave incoming from the left was 511 kV/cm, the pulse shape is trapezoidal (see below), the pulse duration is 2ns, and the constant magnetic field has a component equal to 0.025 T. The thickness of the plasma layer is 8 mm, length 25 cm.
Consider the evolution of the emission process.The distributions of electrons at characteristic moments of time are shown in figures 6 and 7.As can be seen from the figures, electron emission starts sequentially from the left end of the emitter and moves along its entire surface.Moving away from the emitter, the electrons are accelerated by the field and either reach the anode or at first, break through the plasma layer and go to the anode.As can be seen, in the second case, the beam electrons move the plasma electrons towards the anode.It is also observed in experiments.It should also be noted that there is the collector on the right closes the chain and stimulates secondary emission at the right end of the cathode.The end of the pulse leads to a redistribution of electron flows towards the collector and a general fading of the emission processes.
The shape of the TEM-wave pulse and the evolution of the total current of emission in the two considered situations are shown in figure 8.By range of values, they correspond to theoretical estimates and calculated data of the work [6].Herewith some increase in the emission current is observed in the second case, stimulated by the presence of the plasma layer.Quantitative differences may be due to the absence of accurate data on plasma parameters, materials of emitter and collector.Nevertheless, the main characteristics of the emission process are reproduced in the numerical model presented by us quite realistically.

Conclusion
The problem of modeling explosive emission processes in strong electromagnetic fields was considered.To solve it, a new computer model and a numerical algorithm were developed.The model combines the methods of grids and particles.A parallel software implementation of the developed numerical technique was carried out.As a verification of the computer model, the problem of the evolution of argon plasma in the absence of external influences was chosen.This allowed us to clarify the mechanism for numerically specifying the plasma and to establish the true accuracy of its reproduction by the code.Using the developed software tools, preliminary calculations of the emission process in a coaxial diode with magnetic self-insulation were carried out.Comparison of the calculation results with the works of other researchers confirmed the effectiveness of the proposed numerical technique.

4 iM
 (multi-ions are located at the nodes of a square Cartesian grid, 8 multi-electrons surround each multi-ion) and for two values 0

Figure 2 .
Figure 2. Model geometry (on the left) and two variants of the ion environment (in the center and on the right) about plasma concentration can be easily integrated into our model by setting the total charges of multi-ions and multi-electrons:

3 .
The starting velocities of multi-ions and multi-electrons are equal to 0.

Figure 3 .
Figure 3. Model geometry (on the left) and initial arrangement of multi-ions and multi-electrons (on the right).

8 eM
the same time, we will monitor the temperature of the multi-ions.This temperature starts from zero values and evolves slowly, but it is this temperature shows the time of the beginning of the development of the calculation instability.The first series of calculations is related to the analysis of the influence of grid resolution.For this, the plasma parameters are chosen 0 and the condition of matching the steps on time with the value of the steps on space is taken into account h

Figure 4 .
Figure 4.The evolution of temperature of the plasma ions for

9 Figure 5 .
Figure 5.The evolution of the temperature of the plasma ions depending on the number of multi ions in a grid cell 1, 2,4,8 i M  (curves 1-4).

Figure 6 .
Figure 6.Electron density distributions at times t=0.5, 1.0, 1.5 ns (from top to bottom) in the absence of a plasma layer.

Figure 7 .
Figure 7. Electron density distributions at times t=0.5, 1.0, 1.5 ns (from top to bottom) in the presence of plasma

Figure 8 .
Figure 8. Dependences of the total emission current on time (bottom) for a trapezoidal TEM-wave pulse (top).Numbers 1 and 2 correspond to the calculation results with and without the plasma layer.