S-PARMOS-a method for simulating single charged particle motion in external magnetic and electric fields

In this paper we present a new software package for computational simulation of single particle motion in the presence of static external electric and magnetic fields. This seemingly simple problem is in fact very complicated one. Namely, while analytic formulations of this problem are unambiguous, it is very difficult to predict even qualitatively a charged particle trajectory for an arbitrary combination of field parameters, except in simplified cases. Instead, one might perform a single particle motion simulation. Our software package represents a unique tool for this problem. A special new feature of our approach is constructing the instant Larmor Center Trajectory. For the case of slowly changing fields, the Larmor center trajectory reduces to the less general guiding center trajectory approximation. The possibilities for further investigating these two approaches by using our software might be of great interest for both educational and engineering purposes, especially in the areas of gaseous electronics and laboratory, fusion and space plasmas.


INTRODUCTION
Interactive computer simulations with the aid of graphics systems represent a powerful tool in modern education and research.Due to high nonlinearity some apparently simple systems from, e.g., Newtonian mechanics, are in fact, extremely complicated to characterize and understand.Their behavior is very difficult and frequently impossible to predict without empirical and/or numerical experiment.However, both empirical and numerical methods involve their own specific difficulties, e.g., uncertainties and sources of errors.Finally, when these problems are satisfactorily solved, transfer of the understanding for educational and engineering purposes remains a next difficult task.Visualization is one important aspect of this process.There are numerous examples of such attempts in various fields (e.g., regarding general methodology in learning and teaching purposes [1], simulations towards understanding the particular concepts of velocity and acceleration in projectile motion [2]).As argued in [2] computer simulations can significantly help students confronting their cognitive constraints and developing functional understanding of physics.Therefore visualization plays an important role in understanding physical phenomena.There are programs specialized in particular areas of physics, e.g.3D-interactive computer graphics system for waves and optics [3], and also a large scale virtual field laboratory for studies of environmental properties and processes using geostatic and scientific visualization techniques [4].Simulation of general Newtonian mechanics systems is supported by e.g., the physics simulator named PHYSIKI [5] and a physics simulation suite of various demonstration programs [6].There are simple user interface where all parameters of particles and environments can be set.In spite of the fact that both of them have fairly extensive user manuals with examples a great disadvantage of these packages is a rather poor graphic quality (old DOS graphics) which should be in our opinion one of the essential feature of simulations.Moreover, both these packages are designed for a wide number of various possible physical scenarios at an elementary level, rather than being suited for a detailed study of a particular class of problems.On contrary, for studying such a particular class of problems in electrodynamics i.e., single particle motion in external electric and magnetic fields, the Berkeley group developed a simple code ([7]) which is, at the moment, the only available open source program package of which we are aware.Unfortunately, this code was not upgraded for many years.Since this is an important subject of interest in many teaching and research activities and institutions, we find necessary to develop a new code for educational and research purposes with a high versatility, friendly user interface and capabilities for quick upgrading.A special new feature of our approach is constructing the instant Larmor Center Trajectory.For the case of slowly changing fields, the Larmor center trajectory reduces to the less general guiding center trajectory approximation.The concept of Larmor center trajectory is exact one.The question that arises from this concept is if a space element around the instant Larmor center can be used for describing the space-time averaged motion and interaction of a particle with its surroundings.
Our code is user friendly, open source, written in Matlab but easily translatable to other languages.It is intended to perpetual upgrading for better visualization of already solved problems, but primarily for opening and solving new questions with a capability for simplifying certain extremely complex and unpredictable charged particle motion in experimental and self-consistent electric and magnetic non-stationary fields.At the moment the code is still in a non-mature stage but nevertheless, it is already a very powerful tool for understanding particle motion under complex conditions and for attacking new very general physics questions.
In Section 2 we describe the theoretical backgrounds of single particle motion in external fields.The concept of Larmor center trajectory is explained in Section 3. The simulation code is described in Section 4. Section 5 is the presentation of the results and Section 6 contains conclusing remarks.

THEORETICAL BACKGROUND
The non-relativistic motion of a test particle with charge e and mass m in the presence of a combination of charged particle generated and external electric (E), magnetic (B) and other f (e.g., gravitational) fields, is given by Newton's equation of motion (see e.g., [8]): (1) where the fields are given by Maxwell's equations (2) with δ -the Kronecker delta and the dots denote differentiation with respect to time -∂ /∂t.Acceleration and radiative forces in this equation are ignored (this approach is known as the electrostatic approximation).The particle position r(t) is the unknown variable to be found for obtaining a complete solution.Since the analytic solution is not tractable in general (with exception of several simple cases) computer simulation is the only meaningful way to attack this problem.Although our code is capable for many body problem with mutual interaction here we deal only with a single particle assuming that the fields are known.Moreover here we deal only with the particular case of time independent electric and magnetic fields, ignoring the external fields of non-electromagnetic origin f(r, • r,t).In this case it is very important to understand the physical background of the particle motion in the electrostatic approximation.Namely, it can be readily shown that in the absence of electrostatic field the scalar product of • r with Eq.1 yields: ( where W = m • r • r/2 is the kinetic energy of the particle.This means that in the electrostatic approximation, the particle can not gain energy from magnetic field.Consequently, any possible particle energy gain or loss is due to electrostatic force. Although the conclusion about the energy conservation is a very useful one the motion of the particle in combined electric and magnetic field is still analytically unsolvable and qualitatively unpredictable in general.However, under certain conditions the particle averaged in short time scale motion can be well predicted.An important group of such a particle motion belongs to so called guiding center approximation.The guiding center approximation is a very useful concept, which works in those cases where we are able to find a moving (not necessarily inertial) frame of reference in which the particle trajectories are nearly circular.Consequently, the particle motion can be approximated by a time averaged center of rotation over the particle cycle.The prerequisite to apply this approximation is fullfillement of the so-called adiabatic conditions: (4) where ω B = eB/m is the particle orbital frequency around the instantaneous center of particle rotation, is the gyroradius of a particle with velocity and perpendicular to the direction of the magnetic field with gyrofrequency ω B = eB/m.These conditions simply mean that during one cycle of particle motion (in the coordinate system moving along guiding center) the magnetic field must be approximately constant both in space and time.
The first step in approximating the particle motion by using the concept of guiding center is separating the particle velocity into components parallel and perpendicular to the magnetic field.For an arbitrary vector A it can be written using elementary vector rules where b = B /B is the unit vector in the direction of the magnetic field, and by a b, we denote the standard dyadic product and 1 is the unit diagonal matrix.The concept of guiding center drift in the perpendicular direction has been rigorously described by employing the perturbation method ([9.10.11.12.13]).The general equation for the drift velocity v D involving the hierarchy up to the second order is: where the meaning of perpendicular and longitudinal terms is with respect to B. This equation may be considered as a result of averaging the particle motion over the cyclotron period.The first order term in the parenthesis on the right hand side eE + f includes all external forces (static and induced), the second order term (containing ) is due to the static magnetic field gradient and curvature and the last term (so called inertial due to the acceleration • v D ) is of second order.Ordering hierarchy means that more and more rigorous adiabatic criteria should be satisfied with respect to the order of particular terms for using the full approximation.Solution of Eq. 6 is so complicated that it masks the physics.For present purposes, it is enough to point out that the hierarchy of the drifts at least up to second order could be established in a more elementary manner (see e.g., [9]) as follows.
In homogenous magnetic and external fields, the guiding center drift in the direction normal to both external force and magnetic field lines is given by: (7) for general (e.g., gravitational) and electric field respectively.This is a striking result which shows that the resultant motion of the particle has no any component in the direction of the fields but the particle drifts in perpendicular direction.In the case of electric field the drift is also independent of the particle charge sign.There is no limitation regarding the drift velocity (note that this model neglects relativistic effects; if those are included the speed of light becomes an upper limit).In fact, there is no approximation made in deriving this classical drift.That is the reason for this kind of the drift to be considered as zero order.
In a combination of homogenous magnetic and non-homogenous electric fields, a correction term to previous formula arises yielding a total drift velocity: In the case of a non-homogenous magnetic field the drift velocity is (9) If the magnetic field lines cannot be approximated as straight lines, the curvature drift (10) appears.In fact, the curvature and non-homogenity of the magnetic lines in practice always appear together, so the resultant (gradient plus curvature) drift is: (11) These are just several most prominent examples of the particle drifts, i.e., the motion of the particle guiding center that can be expressed in analytic manner.These examples do not exhaust the variety of possible known or still unknown combinations of the external fields in which the particle motion can be time averaged.Moreover, the cases where approximation of a guiding center in non-inertial systems are solved are practically ignored in literature.For possible extension of the guiding center motion to more complicated cases some "empirical" experience is needed.Our idea is to open such a possibility by using the simulation codes with a more general concept of the instant particle motion center, i.e., the Larmor Center Trajectory.

LARMOR CENTER TRAJECTORY
In fact applicability of guiding center trajectory is rather limited to the space and time intervals in which the gyro-radius and gyro-frequency are constant.For a rigorous application of this approximation these quantities in general should be recalculated as the particle moves.Instead of doing this for tracing the particle motion in general nonhomogenous fields we introduce the definition and enable the software for calculating the instantaneous Larmor center trajectory, which is applicable even when the adiabatic conditions are not fulfilled.The position vector of the instant Larmor center can be found by using the Larmor formula .We calculate the Larmor radius with the center defined in the Frenet-Seret local coordinate system.The line r L (t) along which the Larmor radius vector origin moves is then given by: (15 where ds is the differential length of the particle trajectory (see [9]) and r(t) is the real particle trajectory.Our formula is very general one and applicable to both inertial and non-inertial systems.However, it is meaningful only if applied together with a simulation code.Then the general properties of the system of interest could be deduced from an exact solution in order to extract relevant information for possible approximate characterization of the system.

CODE DESCRIPTION
The simulation code S-PARMOS is an open source code of Ljubljana and Berkeley.It written in Matlab and can be easily upgraded by adding new combinations of stationary and non-stationary electric and magnetic fields.At the present stage these fields are analytic external ones.A self-consistent field as a result of many particle motion is possible to simulate as well, but this is not of special interest in this manuscript.There is a subroutine developed for employing the external field of arbitrary current carrying conductors and macroscopic charge carriers as well.However, a strategy for importing the externally generated fields by using other field generators, independent of our program, is adopted for future.At the moment the program is fairly documented by a user manual [15].The program has a graphical user interface (GUI) as illustrated in Fig. 1.The GUI window is currently divided into five parts: Magnetic field, Electric field, Particle information, Simulation and Field lines adjustment.In the "Magnetic field" and "Electric field" part, there is a pull-down menu, where the user can choose between different analytically defined fields, for example: field in z axis or x axis direction, exponentially decreasing field, or custom field.Once the field is chosen, the user is asked to set parameters which define the field.For each field a comment is written in order to aid the user in defining parameters.
In the "Field lines adjustment" the user can manipulate each field line length.With the "choose field" popup menu, the user chooses magnetic or electric field and with "choose line" selects the field line.In the current version, a maximum of 18 field lines can be drawn.When a field line is chosen, user can adjust its length by "set length" edit window.Length must be confirmed by pressing the return key.Inspection of lengths of individual field lines is simply done by selecting the appropriate field line number and the length will be shown in the "set length" edit window.In the "Simulation" part of GUI, there are two radio buttons for showing the trajectory and the guiding center.If switched on, additional traces of trajectory and Larmor center are drawn during the simulation.
During the simulation, two vectors are shown with their origin at the particle position.A vector for particle velocity and a vector for particle acceleration (or force).If the length of the vectors are inconvenient for viewing (due to scaling of the window versus magnitude of the vectors) the user can alter the scaling with two sliders in the "Simulation" part of GUI (the "adjust force vector visualization" slider and "adjust velocity vector visualization" slider).When all parameters are set, the button "Start simulation" should be pressed and a window showing the spatial location of the particle and trajectories is shown (Fig. 2).
An animation of this simulation is given in [16].User can see the field lines of magnetic (grey) and electric (blue) fields, the coordinate system and the particle with a force vector (magenta), velocity vector (yellow), trajectory (pink) and Larmor radius center (magenta).
During the simulation the user can switch on and off the particle trajectory, Larmor center and the vectors.The simulation can also be paused with the "Pause simulation" button and then restared with "Restart simulation" button.The "Start simulation", "Pause simulation" and "Restart simulation" buttons are one and the same.The viewpoint in this window can be manipulated with mouse or manipulation buttons.
The "Stop simulation" button stops the simulation completely.When the "start" button is pressed the fields are recalculated, the trajectory and the Larmour center are reset to zero length vectors.
For further data manipulation of a simulation, the user can save all the data in an ASCII file by selecting Save Data in File pull-down menu or by key combination Ctrl+D.The data is saved in three columns, each represents one coordinate.The data is saved as follows: particle trajectory, Larmor center, first magnetic field line, second magnetic field line, etc., Figure 1 Simulation GUI window is divided in to five parts: "Magnetic field", "Electric field", "Particle information", "Simulation" and "Field lines adjustment".first electric field line, second electric field line, etc.A description is written before each section of data (e.g."B line number 7").
At present moment this simulation is not yet a stand alone executable application, therefore it must be run under Matlab R14.As the simulation uses external vrml file the Virtual Toolbox should also be installed.
In Fig. 3 (animation [17]) we show an example of the Larmor center trajectory of a charged particle for a combination of non-homogenous magnetic fields which in a toroidal geometry which is at the same time characterized by a non-homogenity and the curvature of the magnetic field lines.
In tokamak devices, however, besides the toroidal magnetic field there is a poloidal field as well.In Fig. 4 (animation [18]) we show simulation in such a configuration with the so called safety factor defined as where a is the poloidal (local) radius and R is the toroidal (main) radius of the torus (The term "safety" is used because according Kruskal and Shafranov theory of the tokamak instabilities its value should be enough high i.e., always greater than unity for avoiding current-driven instabilities in tokamak devices.In geometrical sense safety factor is is the number of times a magnetic field line goes around a torus, so it is a local quantity which in typical tokamak devices ranges from near unity in the center of the plasma to 2-8 at the edge.).simulation in Fig 4 has been performed with q=3 at the main radius in the center of the plasma.Without further explanations of the tokamak plasmas we emphasize that the particle motion in Fig. 4, is so called "banana orbit" shaped.Namely, particles in certain range of initial energies become trapped in a region of tokamak due to magnetic field ingomogenity.Circular projection of this orbit, i.e., the guiding or Larmor center at a r-z plane in cylindrical system is banana-shaped.Particles with higher initial energies are untraped and rotate "indefinitely" along the magnetic field lines (animation [18]).In Figs. 5 (animation [19]) and 6 (animation [20]) there are two cases with the parameters chosen so as to illustrate the particle motion when the guiding center approximation a-priory fails.Example in Fig. 5 corresponds to so called magnetic mirror whereas example in Fig. 6 is addressed to magnetic cusp configuration.Both examples are created by using superposition of two magnetic monopoles.In spite of this theoretical idealization, these configurations are easily realized in real devices and are of high importance in space, laboratory, technological and fusion plasmas.We use these examples here to illustrate that the concept of Larmor radius center trajectory is still a well defined curve in spite of the fact that the adiabatic conditions are completely lost and makes no any sense to try any construction of the guiding center approximation.Of course, our result is just an illustration indicating that it is would be possible to reduce even adiabatic motion to an approximate path.How to do that for particular analytically given external fields cases is a future task that we could not solve within the framework of this work.

DISCUSSION
One fundamental problem of numerical simulation is, in general, the accumulation of numerical error, such as roundoff error and trunctation error, which are the consequence of discrete approximations made to the continuum equations.Numerical error can lead to nonphysical motions which vary depending on the specifics of the parameter regime, so this question should be investigated in more detail during the further code development.Important parameters for this future study are the time step, initial conditions, discretization method, numerical precision, etc.Based on preliminary investigations we are at this stage of development convinced that none of these parameters is crucial for reproducing the gobal particle motion behavior, and indeed can be suppressed or compensated by applying more advanced algorithms or filters.The principle development of this tool includes an educational module for insertion in graduate and advanced undergraduate courses in plasmas or electromagnetics.This will provide substantial synergy with current analytic approaches by allowing students to solve for and visualize fields in many types of plasma devices, from fusion tokamaks and stellarators to low temperature plasmas to accelerators to astrophysical plasmas.We anticipate that the educational benefits of the tool will extend well beyond the classroom to produce educational materials for the general public to improve understanding of such grand challenges as ITER as well as for use in describing the operational principles of commercial plasma devices.We intend to further upgrade the program by implementating a considerable number of new features.One of the most important ones, to be done urgently, is the possibility of importing the externally generated from other codes or from experimental data-basis magnetic and electric fields.This will enable potential users to work with arbitrary fields generated by independent computational tools, using our package without changing the source file and without knowing any programming.The second task is further implementation of time varying magnetic and electric fields as well as external electromagnetic waves.Finally, we intend to employ random field generators superimposed with the static fields for simulating and better understanding of simple particle motion in turbulent fields.Last but not least we are developing the virtual reality interfaces for a full 3D presentation of simulation for Virtual Realilty machines (e.g., Virtual Reality Passive System installed in LECAD laboratory).This system is used to create real 3D images and animations for better 3D object presentations.As vrml file contains 3D information an appropriate software viewer is able to generate stereo images that can be used with this system.These images are then fed separately into two projectors with common screen.Light from each image is predefined by a polarizer for a later separation.Both pictures are then reflected from the screen or passed through the screen (rear projection) without light polarization loss.A user wears special glasses with  polarizing filters which finally separate both pictures for each eye.In this way two computer generated pictures are combined by the user brain into a 3D image.Executable files and a version compatibile with Virtual Reality System will be available soon on the LECAD laboratory website.The current version of the source code is already available on request.

CONCLUSION
A new single particle motion simulation code for a wide range of applications has been developed, based on the general concept of Larmor center trajectory, which is of general importance for many applications.Under adiabatic particle motion the Larmor center trajectory reduces to the approximate guiding center trajectory.Our package is at the moment available with capabilities for simulating particle motion in stationary non-homogenous magnetic and electric fields, but extension to time varying external fields is straightforward and this work is in progress.At the moment the package deals with analytical external fields.The extension to the fields given in discrete points as a data obtained from other programs or experimental measurements is in progress as well.

Figure 2
Figure 2 wrl file viewer: the virtual world can be manipulated by a mouse or by navigation pannel at the bottom of the window.

Figure 3 Figure 4
Figure 3 Non-homogenic magnetic field in a torus.

ZFigure 5
Figure 5 Particle is trapped in the magnetic mirror, created by two magnetic monopoles aligned paralel with z axis.Monopoles have opposite sign.

Figure 6
Figure 6 Particle escapes from the magnetic cusp near the middle of both magnetic monopoles.Monopoles have equal sign.