Simulating single-particle dynamics in magnetized plasmas: The RMF code

At the basic level, the motion of charged particles in plasmas must be understood by single-particle motion driven by the fields present at the particle at each instant. Particle-in-cell codes do this self-consistently, but are computationally expensive and time-consuming, evaluating the time evolution of billions of interacting particles. Test-particle codes, single-particle codes, are far quicker. They are highly accurate when the simulation time scale is far shorter than MHD or resistive time scales conditions often found in high temperature plasmas.

To aid the design of plasma diagnostics, analyze their results, and understand the physics of particle motion in complex fields, the Rotating Magnetic Field (RMF) code was written to calculate the trajectory of a charged particle in a specified electromagnetic field. The RMF code integrates Hamilton’s equations in cylindrical coordinates using the adaptive predictor-corrector Double-precision Variable-coefficient Ordinary Differential Equation (DVODE)¹ solver module for speed and accuracy. RMF has multiple capabilities for the field: a Hill’s vortex model² and a numerical Grad–Shafranov solution^3,4 for a Field Reversed Configuration (FRC); a rotating magnetic field;⁵ mirror coils; wires along the z-direction; static electrostatic fields normal to flux surfaces; and oscillatory magnetic-field-perpendicular grids to explore Fermi-acceleration physics.

Analyses include the computation of time-evolving energy histograms, Poincaré puncture plots, fast Fourier transforms of the energy, and Lyapunov exponents.⁶ Recent features include a synthetic diagnostic for simulating the observations of charge exchange neutrals⁷ and an array of grids to enable the study of a new form of Fermi acceleration.⁸ 

Particle motion is initialized by specifying the position (r, ϕ, z) and velocity using energy in eV and spherical angles about the z-direction. The six-dimensional state vector and many derived quantities are saved as functions of time. Time is measured in units of cyclotron periods for a selected magnetic field. A post-processing graphics code, XDRAW, is then used to plot up to 12 windows of any two quantities using different colors to denote successive time intervals. XDRAW provides extensive interactive features to zoom, digitize, compute slopes and ratios, and save postscript files. A single RMF simulation takes a few seconds. Multiple cases of RMF may be run in parallel and post-processing data mining run on the results.

The central feature of the RMF is numerical integration of Hamilton’s equations for the motion of a charged particle in a given electromagnetic field. The Hamiltonian is given in cylindrical coordinates by

where p_i are the conjugate momenta, q and m are the charge and mass of the charged particle, A is the magnetic vector potential, and φ is the electrostatic potential. The corresponding equations of motion are given by

They constitute a sixth-order coupled system of nonlinear ordinary differential equations. These equations are integrated numerically with DVODE set to a specified tolerance, typically 10⁻¹² in the H per step. (A long simulation, τ ∼ 10⁴, would have an accumulated error in H below 10⁻⁶.) The results are visualized with the XDRAW graphics code, which draws contour plots, including of the field, and plots of multiple curves of any two variables.

In the absence of a rotating magnetic field, the fields are independent of t and ϕ and the equations conserve the Hamiltonian H and conjugate angular momentum p_ϕ. We monitor H to verify that it is conserved to high accuracy. In the presence of a rotating magnetic field, these are no longer conserved separately, but there is a canonically transformed Hamiltonian K = H − ωp_ϕ, which we monitor for high accuracy.

Elastic or inelastic collisions are an option in RMF for which the Hamiltonian is re-initialized after each collision. A simulation is terminated when a particle’s trajectory extends beyond a specified region or when the specified simulation's duration is exceeded.

In RMF, a Field-Reversed Configuration (FRC) is represented by the Hill vortex model whose flux function, ψ, inside the separatrix is

where r_s is the separatrix radius at the midplane z = 0. An O-point (line) is at $ro=rs/2$⁠, ±z_s are the x-point locations, and κ ≡ z_s/r_s is the elongation. This satisfies a Grad–Shafranov equation,

with $ψ0=B0rs2/2$ being the confined magnetic flux, B₀ being the magnetic field at r = z = 0, maximum pressure $p0=2B02(1+1/4κ2)$⁠, and pressure profile p(ψ) = p₀(4ψ/ψ₀). The solution of Zakharov and Shafranov⁹ is used outside the separatrix to enforce non-negative pressure there. An alternative FRC representation may be implemented in RMF: the Grad–Shafranov equation numerically solved with specified boundary conditions.

The RMF FRC work thus far reported using Hill’s vortex model. Landsman et al.¹⁰ examined the stability and periodicity of the three charged-particle orbit classes¹¹ in FRCs.¹² μ, the magnetic moment, herein defined as the perpendicular energy divided by the local magnetic field, was not conserved when particle orbits approached the X-points,¹³ a phenomenon studied earlier in mirror configurations^14,15 and geophysical plasmas.^16,17

Particles are accelerated by the induced electric field created the rotating magnetic field. Its vector potential is derived from vacuum conditions, ∇ × ∇ × A = 0. There are two types of parity,¹⁸ depending on the field’s symmetry about z = 0. In RMF, the odd parity’s vector potential is given by¹⁹ 

while the even parity’s vector potential is given by

where B_o and B_e are the odd and even magnetic amplitudes, I₀ and I₁ are modified Bessel functions, k = π/2z_s, ξ = kr, ψ = ϕ − ϕ₀ − ωt, and ω is the frequency of the rotating magnetic field.

Studies of the effects of rotating magnetic fields on charged particles in FRCs revealed a number of phenomena:

1. Ions were heated to energies near $mi(ωrs)2/2$ when 0.2 > ω/ω_ci < 2 and slightly lower when −0.2 > ω/ω_ci > −2 for B_o > 0.01B₀ while electrons were heated to energies up to $qπrs2ωBo$ for lower relative frequencies, ω/ω_ce ∼ |0.001|. There was a B_o/B₀-dependent threshold for heating.²⁰ 

2. Heated ions circulated around the FRC in betatron orbits nearly in synchronism with the RMF_o.

3. High energy ions and electrons circulated in opposite directions in betatron orbits, independent of the direction of the rotating magnetic field’s rotation.

RMF models the field as an ensemble of circular coils, each a filament of radius a carrying current I. A single coil, coaxial with the z axis and centered at z = 0, creates a dipole field with vector potential

with

and complete elliptic integrals

The elliptic integrals are represented by polynomials and logs of m.

By positioning two coaxial coils at z = ±z_c (with centers on the z axis) and currents flowing in the same direction, a mirror configuration is formed.

The purpose of the mirror field configuration was to study the lack of μ conservation in it. An example is shown in Fig. 1.

The particle was initially confined, as shown in Fig. 1 (top). However, after a few axial transits, μ fell to such a low value [Fig. 1 (bottom)] that the particle was lost through a mirror throat. Changes in μ occurred when the particle crossed the midplane. For most of the time, the ratio ρ_i/R_c was small, $<1%$⁠, which, according to the common theory, should result in μ conservation. Nevertheless, μ was not conserved. Theoretical^15,21,22 and RMF studies showed the importance of the field-parallel particle velocity.

Similar behavior, lack of μ conservation, was seen for ion motion in a dipole field configuration [see Fig. 2].

Guiding center theory shows that near a straight current-carrying wire a drift will develop parallel to the wire. Changes in μ do not occur^15,21,22 because a changing field curvature is required. That occurs in FRC and mirror devices but not for a single straight wire.

The vector potential of a single straight wire of length l is given by

where r is the distance from the wire. As Figs. 3 and 4 show, both high and low energy particles, those with ρ/r_c ∼ 0.5 and 0.05, show μ conservation, supporting the aforementioned theories and strongly in contrast with the orbits in mirror machines, FRCs, and dipoles.

The structure and several applications of the RMF code were described and note was made of current uses of the RMF code as a synthetic diagnostic for charge exchange measurements of ion energy distributions in high-β plasmas and to simulate field-parallel Fermi acceleration in mirror devices. Field options in the code allow studies of particle trajectories in the fields of static FRCs, axisymmetric mirrors, dipoles, cusps, and current-carrying wires, with or without static electric fields having a variety of geometries. Time-varying options include rotating magnetic fields and B-parallel electric fields created by transparent grids in the plasma.

This work was performed under DOE [Contract No. DE-AC02-09CH11466].

The authors have no conflicts to disclose.

The code output used to generate the graphs in this study is available from the corresponding author upon reasonable request. The figures and numerical results in this paper are openly available at http://arks.princeton.edu/ark:/88435/dsp01x920g025r, Ref. 23.