### Computation in Plasma Physics

#### I. Introduction

Though it is not generally realized, Plasma Physics is the only major branch of physics which has come into being largely in the past generation (NRC(b), 1986).

Since the 1830’s, when Faraday first studied electrical discharges, to the 1920’s – when Langmuir discovered plasma oscillations, and Breit and Tuve reflected radio waves from the ionosphere – discharge phenomena in electrical devices, in the earth’s atmosphere and outer space were considered as belonging to separate fields.

In the subsequent period up to about 1950, investigations by a galaxy of scientists – Alfven, Appleton, Chandrasekhar, Landau, Saha, Spitzer, among many others – revealed that the behavior of these apparently different systems was rooted in the same common structure of matter, which is manifest when a gas or vapor is ionized.

Their investigations brought into prominence the essential property of ** plasmas**, a term first used by Langmuir, which is absent in the other three states –

*collective behavior*. They laid the foundation of modern plasma physics by investigating how ionospheric plasmas influence the propagation of radio waves, how activity in the sun effects the auroras and creates magnetic storms on the earth, what is the role of plasmas in the formation of stars and galaxies, how waves and charged particles in plasmas are coupled.

Experiments on laboratory gas-discharges also yielded new and accurate methods of measurement of plasma parameters such as temperature and charge density. The space program and the fusion program has stimulated the growth of research in plasmas. Along with the realization that understanding natural processes on a cosmic scale requires the solution to problems in plasma physics, it is now amply clear that progress in nuclear fusion is also dependent upon the solution of problems in plasma physics rather than in nuclear physics.

The result has been the birth of a new subject area called Plasma Science. Although this includes plasma physics, in its scope also lies the description of the large class of ionized matter which participates in atomic and molecular transport processes, chemical reactions, as well as the excitation of neutral materials and their interactions with ionized media.

Though the essential motivation – and funding – to study plasmas came from these research programs, new application areas of plasma technology soon appeared. These have become full fledged research programs themselves. Some of them are free-electron and x-ray lasers, plasma isotope separation methods, neutron sources, and accelerators based on collective effects.

The principles of plasma behaviour at the microscopic level are well established. The challenge lies in the developement of techniques that will elucidate the behaviour of plasmas at the macroscopic level knowing the principles which govern the microscopic behaviour. However, it is becoming increasingly clear that though the dynamical equations governing plasmas in many configurations are structurally simple, they are impossible to solve analytically, since they are highly non-linear, and involve infinite number of degrees of freedom (Armstrong, 1970).

Progress in plasma physics would thus be considerably enhanced if computational methods were employed. In a survey carried out by a committee appointed by the National Research Council (NRC(b), 1986), it was recognized that computational techniques have considerably helped in understanding the processes occurring inside magnetically confined fusion plasmas.

In the area of space and astrophysical plasmas, the committee stressed that large scale numerical models are expected to play a central role in interpreting astronomical observations, as well as in motivating new and different kinds of observations. Observing the fact that access to major computational facilities was limited, the committee recommended a national computational program dedicated to basic plasma physics, space physics and astrophysics, with the aim of maintaining the technology and providing access to facilities suitable for carrying out large-scale calculations and simulations.

The use of such facilities was expected to lead to significant advances in understanding hot and strongly coupled plasmas, such as those occurring in inertially confined fusion targets, or in red giants and neutron stars. The importance given to computational techniques in the fusion program has given birth to the proposal of a “Numerical Tokamak” (Dawson, 1993, 1995). By harnessing the power of modern supercomputers and parallel computers, plasma simulation can model the behavior of experimental tokamaks. These simulations will significantly improve the understanding of energy and heat transport in such devices, and could save money otherwise spent in building the next generation of these machines. The High Performance Computing and Communications Initiative of the US Department of Energy has declared the “Numerical Tokamak” a Grand Challenge, and a consortium of eleven institutions in the USA has taken up this task.

#### II. The early years

The roots of plasma simulation can be found in the methods used in designing vacuum tubes, oscillators and other first generation electronic devices. Numerical solutions to Poisson’s equation were commonly used to the map the electric field distribution in a triode; trajectories of electron streaming in a magnetron were obtained by solving the equations of motion on a desk calculator. They were, however, essentially solutions representing a static situation.

In a dynamic state represented by a collection of freely moving positive and negative charges subjected to external fields, it was realized that similar techniques could be used provided some modifications were made. These techniques were pioneered by scientists such as Buneman at Cambridge (Buneman, 1959), Dawson at Princeton (Dawson, 1960) in the late ’50s and early ’60s, and later extended by Birdsall, Langdon, Hockney, Eastwood and many others.

Initially attempts were made to model plasmas in one dimension by tracing the interaction of a few hundred charges constrained to move in a line. In these early calculations, the ions were kept static in the background while only the electrons were allowed to move. More elaborate programs were limited by the available computer resources. However, as computers grew in speed and storage capacity, the models also evolved. The simulation of warm plasmas became possible when the charges were assigned a distribution of velocities. Multi-component plasmas were simulated when ions were also allowed to move. With the extension of these models by adding two more velocity components, the basic plasma simulation model known as the 1d3v plasma (one dimension in length and three dimensions in velocity) as it is called, is still very widely used to day.

#### III. Extensions to the early models

A big step was taken when two dimensional models were introduced. The key concept which enabled this extension was the introduction of the mesh – or grid. This concept envisaged the division of the computational space into regular cells, which was superimposed over the plasma, so that charges moved from one cell to another.

By solving the field equations in the finite difference form, the calculation of the electric field became extremely economical. The charge density used to find the solution was taken from the positions of the charges, and assigned to the mesh. The theory of plasma simulation grew with understanding of how methods used for the assignment of charge to the mesh influenced the physics of the situation. The noise introduced by the spatial grid and limited number of particles were also recognized as potential error sources. The effect of the grid was to increase particle sizes, and reduce forces at very short distances, thus introducing a smoothening effect. It was demonstrated by Langdon (1970) that plasma theory also holds for finite sized particles, thus putting plasma simulation on a firm theoretical ground.

The early models allowed simulations only over limited space and time scales. Simulation over longer times became possible with the introduction of implicit differencing schemes (Langdon, Cohen, Friedman, 1983). Methods to increase the space scales in plasma simulation is an area of which is being addressed widely (Gibbons and Hewett, 1995).

#### IV. Modern plasma simulation codes

The state of the art in plasma simulation is represented by three-dimensional fully electromagnetic codes used to investigate fusion plasmas. These codes typically are used on supercomputers to model the plasmas existing inside inertial and magnetic confinement fusion devices. They integrate all the Maxwell equations, and take into account currents and magnetic fields generated by the charges, as well incorporating interactions between charges and radiation whose wavelength spans the full electromagnetic spectrum. Non-electromagnetic interactions such as atom-atom collisions and atom-wall interactions are also taken into account.

Computation in plasma physics has reached a mature and fully developed stage, in a manner which is not even hinted at in the standard text books on plasma physics. It is now an independent and recognized subfield of research called *Computational Plasma Physics*. Not only is it used provide insight into the essence of plasma phenomena, but also to model the behavior of whole plasma machines. It is one of the success stories of modern technology.

#### V. The analytical approach in Plasma Physics

Plasma Physics concerns itself with the behavior of charged particles subjected to their own and external electromagnetic fields. The particles may be electrons or ions, most often being a mixture of both in approximately equal number. They may be confined in a small volume or free to move over unimaginably large distances. Phenomena occurring in a plasma has much in common with other states of matter, notably interactions involving wave propagation, energy transfer, turbulence, instabilities, etc.

Low density plasmas are usually analyzed by studying the kinetics of the interaction of electromagnetic fields with a single charged particle. This is because collective effects are negligible. The equation of motion of the single particle is solved to obtain an expression describing its trajectory in space. Since the external fields may be spatially and temporally varying, the trajectory that is obtained varies in complexity accordingly. This technique finds an important extension in one method of plasma simulation known as the particle method. The single particle approach is not applicable to plasmas of high density.

High density implies that there is appreciable interaction between charges over larger distances. This interaction is what differentiates a plasma from other physical systems. The approach used in solving complex problems involving high density plasmas is completely different. The standard analytical model is called, somewhat inappropriately, the kinetic theory.

For the purpose of analysis, the finite number of discrete particles that constitute the plasma is considered to be infinitely subdivided into a continuous distribution of charge and mass. This distribution is smeared throughout the volume occupied by the plasma. The behavior of the plasma is described by the physical parameters which are defined at all points in space and time, and which can vary in a continuous fashion. The motion of the charges in any given small volume is replaced by a variation in time of the charge density and mass density in that volume. In other words, the plasma is thought to be a fluid.

The function which describes the distribution of charges in velocity space, position space and time is called the velocity distribution function, and it changes with time as the plasma evolves. The equation which governs the time rate of change of the distribution function is called the transport equation. This rate of change depends on many factors, notable ones being interparticle collisions and forces due to electrical and magnetic fields. Essentially, the study of the plasma involves finding solutions to the transport equation, under various physical conditions.

Since, in an experiment, one can observe only the cumulative action of a very large number of charges, the model must have ways to convert a microscopic distribution into macroscopic parameters. This is done by taking moments of the distribution function – essentially multiplying the distribution function by the required physical parameter and integrating over the velocity space. This procedure results in a partial differential equation specifying the behavior of that physical variable. The solution of this equation gives the required physical parameter. Unfortunately, the equation for each moment also includes another parameter which is dependent upon a higher moment. At no stage therefore is the set of equations complete. A solution is found by approximating the highest order moment and closing the series. The set of equations so formed are called the hydrodynamical equations, and different approximations used in the highest moment result in sets of equations having different characteristics (Chen, 1974; Seshadri, 1973).

Thus, the entity that represents the plasma is a mathematical model consistng of a set of equations. Unfortunately no such set of equations exists which can represent the plasma inside a discharge tube and also tackle fusion plasmas . The hydrodynamical equations describe the plasma only to some degree of accuracy and are applicable only over a limited range of parameters. The range of validity is dependent on the simplifying assumptions which have been made in order to arrive at the equations.

#### VI. Approaches To Plasma Simulation

Corresponding to the two different methods of analyzing plasma behavior, namely, the single particle approach and the kinetic approach, there are two basic methods of plasma simulation, namely, the particle method (Birdsall and Langdon, 1985; Hockney and Eastwood, 1981) and the fluid method (Armstrong, 1970).

In the fluid or continuum method of plasma simulation the time evolution of the velocity distribution function is calculated by coupling the Boltzmann or Vlasov transport equation with the Maxwell equations.

By “sampling” the plasma at a point in the phase space the program can give the value of the distribution function at that point, using which other physical parameters of interest can be calculated. Results so obtained are difficult to compare with experiment, but are directly comparable to theoretical results, and this is often the first method of simulation used by theoreticians.

Simple fluid models leave out much of the physics that is important in understanding the behavior of plasmas. By leaving out collisional kinetic effects, information about damping and nonlinear saturation of unstable modes is lost, and instabilities may be improperly simulated. In order to include these effects, therefore, the hybrid approach, which suitably combines the particle and fluid models, has been developed.

The fluid method has been used very successfully in understanding analytically the behavior of simple plasma systems. However, where complicated systems are concerned, the computational approach based on the particle method has had more success. The particle method is an extension of the technique of solving the equation of motion of a single particle. This technique is applied to all the charges that make up the plasma. The motion of each of the individual particles is separately calculated. The particles are tracked in time as they move in their own as well as externally applied electromagnetic fields, which are determined by the Maxwell equations.

Forces on the particles are calculated using the Lorentz equation. These forces determine the change in position and velocity of the particles within a small time interval. As they move, their redistribution in space modifies the fields, which are recalculated periodically. The new fields are used in updating the particle positions and velocities. This is a self consistent technique, and its repeated application results in a simulation of the plasma behavior.

Physical parameters like currents and number densities, various energies and field quantities are “measured” at different points and times as the output of the simulation program. Results obtained using this method are directly comparable to experimental results, and it has been very successful in simulating a variety of non-linear plasma behavior.

In this simple approach, however, interparticle collisions which are responsible for energy transfer between charges and neutral particles cannot be accounted for. In order to include these effects, methods have been formulated by using Monte Carlo techniques (Birdsall, 1991). This makes it possible to investigate collisional processes within the particle approach. The appeal of the particle method lies in its preservation of the use of first principles in its calculations. These are beyond doubt. The use of artificially generated probabilities in what is essentially a deterministic model is a matter of concern to some purists.

#### VII. Basic References

*Armstrong, 1970*. **“Solutions of Vlasov’s equation by transform methods”**, `by T. P. Armstrong, R. C. Harding, G. Knorr, D. Montgomery in Methods in Computational Physics, volume 9`

`(Plasma Physics), edited by B. Alder, S. Fernbach, and M. Rotenberg, Academic Press, New York, 1970.`

*Buneman, 1959*. **“Dissipation of current in ionized media”**, `by O. Buneman, Phys. Rev., 115, 503, 1959.`

*Birdsall and Langdon, 1985*. **“Plasma Physics via Computer simulation”, **`by C. K. Birdsall and A. B. Langdon, McGraw-Hill Book Company, 1985.`

*Birdsall, 1991.* **“Particle-in-Cell Charged Particle Simulations, plus Monte Carlo Collisions with Neutral atoms, PIC-MCC”, ** `by C. K. Birdsall, IEEE Transactions on Plasma Sciences, 19, 65,`

`April, 1991.`

*Chen, 1974.* **“Introduction to Plasma Physics”**` by F.F.Chen, Plenum Press, New York, 1974.`

*Dawson, 1993.* **“High performance computing and Plasma Physics”, **`by Dawson, J. M., V. K. Decyk, R. Sydora, P. Liewer, in Physics Today, 46, 64, March 1993.`

*Dawson, 1960.* **“Plasma oscillations of a large number of electron beams”,** `by Dawson,J. M., Phys. Rev., 118, 381, 1960.`

*Gibbons and Hewett, 1995.*** “The Darwin direct implicit particle-in-cell (DADPIC) method for simulation of low-frequency plasma phenomena”**,` by M. R. Gibbons and D. W. Hewett, J. Comp. Phys, 120, 231, 1995.`

*Langdon, A. B., 1970.* **“Effects of the spatial grid in simulation Plasmas”,** `by A. B. Langdon, J. Comp. Phys, 6, 247, 1970.`

*Langdon, Cohen, Friedman, 1983*. **“Direct implicit large time step particle simulation of plasmas”,**` by A. B. Langdon, B. I. Cohen, A. Friedman, J Comp Phys , 51, 107, 1983`

*Hockney, Eastwood, 1981*. **“Computer simulation using particles”,**`by R. W. Hockney and J. W. Eastwood, McGraw-Hill Book Company, 1981.`

*NRC(a), 1986.* **“Scientific Interfaces and Technological applications”**, `a volume in the series Physics through the 1990s, National Academy Press, 1986.`

*NRC(b), 1986.*** “Plasmas and Fluids”**,` a volume in the series Physics through the 1990s,`

`National Academy Press, 1986.`

*Seshadri, 1973.*** “Introduction to Plasma Physics”**`by S. R. Seshadri,`

`American Elsevier Publishing Co.`