Study of adaptive symplectic methods for simulating charged particle dynamics

In plasma simulations, numerical methods with high computational efficiency and long-term stability are needed. In this paper, symplectic methods with adaptive time steps are constructed for simulating the dynamics of charged particles under the electromagnetic field. With specifically designed step size functions, the motion of charged particles confined in a Penning trap under three different magnetic fields is studied, and also the dynamics of runaway electrons in tokamaks is investigated. The numerical experiments are performed to show the efficiency of the new derived adaptive symplectic methods.


1.
Introduction. In studying the collective dynamics of magnetized plasmas, it is noticed that in the simulation of plasmas, the motion of charged particles under the influence of electromagnetic field is the most fundamental physical process. Many important phenomena in plasmas can be understood and analyzed in terms of single particle motion which satisfies the Lorentz force equations. It is also noticed that for plasma applications there are situations in which very longtime simulation is needed. For example, the simulation for dynamics of runaway electrons involves timescales spanning 11 orders of magnitude. This brings difficulties to both analytical and numerical studies. Conventional integrators, such as the fourth order Runge-Kutta method, often generate numerical results with significantly growing errors because of the lack of long-term stability and accuracy. A recent development is to implement geometric numerical methods [6] in computing the trajectory of particles. In contrary to conventional integrators, geometric numerical integrators have superior long-term stability. Not only can they guarantee the accurate simulation over long-time, but also preserve the important physical quantities such as energy, momentum, charge, etc. These geometric numerical methods have been successfully applied in solving the (relativistic) charged particle dynamics such as (variational) symplectic methods [19,27,29,28], volume-preserving algorithms [8,9], Poissonstructure preserving methods [13,10], etc.
Adaptive time integrations are a class of numerical methods whose time steps can be adapted during simulation. They have been well established and perform extremely well for many applications. As for Hamiltonian systems, the numerical results in [4] and [2] state that using the symplectic methods combined blindly with automatic step size control can deteriorate the correct qualitative behavior of numerical solutions. This poor performance can be explained by backward error analysis and the associated modified Hamiltonian. There have been a great of effort to overcome this problem. Among them a novel variable step size technique realized by time transformation is introduced in [5,21]. It is proved in [5] that the symplectic methods derived in this way can be used with adaptive time steps without losing qualitatively correct longtime behavior. This technique has also been analyzed in [22], and a model Hamiltonian system (the cubic oscillator) is investigated in the numerical experiments. The authors in [24] have extended the adaptive time stepping technique into the domain of variational integrators. In the application of plasma physics, there already have been some results about the adaptive time-stepping integrators [26,25]. In [26], the adaptive Boris-Buneman method is proposed for simulating the particle accelerations. Moreover, the adaptive symplectic methods are applied to simulate the dynamics of runaway electrons confined in tokamaks in [25].
It is clear that the key point in implementing the adaptive symplectic methods is to find an appropriate and cheap-to-compute time transformation which is called the step size function in this paper. Generally, for a given Hamiltonian system we choose arclength parameterization σ(z) = J∇H(z) −1 as the step size function which can be employed to any Hamiltonian system without having a priori knowledge of the solution to the system [6]. However, there are some shortcomings. For instance, for some systems the derivatives of step size function σ may become very complicated to calculate and program. And, the use of step size function may not improve the computation efficiency apparently. In this paper, we mainly study the application of adaptive symplectic methods to the charged particle dynamics. By analyzing the physical nature of the motion, the step size functions related to the particle's gyro-period are presented. In the numerical experiments, we study the motion of a single particle in a Penning trap. In this example, three kinds of magnetic fields: uniform, non-uniform and asymmetric are considered. The longtime simulation for the collective behavior of runaway electrons in tokamaks is also observed by means of parallel computing. The numerical results reveal the high efficiency and good accuracy of the new derived numerical methods.
The rest of this paper is organized as follows. In section 2, we formulate the Lagrangian, K-symplectic and Hamiltonian forms for the motion equations of charged particles under a given electromagnetic field. The idea of adaptive symplectic methods is introduced in section 3, and also the new step size function is proposed in this section. In section 4, we present and analyze the numerical results. Section 5 concludes this paper.
2. Equations of motion under Lorentz force. For a charged particle with rest mass m 0 and charge e, under the influence of electromagnetic field its dynamics is governed by the Newton-Lorentz equation where x ∈ R 3 denotes the position of the charged particle, and m = m 0 γ is the mass with γ = 1/ 1 − ẋ 2 /c 2 denoting the Lorentz factor. Here c is the speed of light in vacuum. Lorentz force F experienced by a particle is induced by a given electromagnetic field. The electric field E and magnetic field B can be described via the scalar potential ϕ(x, t) and vector potential A(x, t) as E = −∇ϕ(x, t) − ∂A(x,t) ∂t and B = ∇ × A(x, t). When the velocity of the particle is much smaller than the speed of light c, i.e., ẋ /c → 0, then Eq. (1) is reduced to the system in nonrelativistic case. In this paper, we assume the electromagnetic field to be static which leads to E = −∇ϕ(x) and B = ∇ × A(x). Lagrangian formulation.
The Newton-Lorentz equation (1) is the Euler-Lagrange equation obtained from the variational principle with Lagrangian function Denote the Lagrangian two-form on space (x,ẋ) by which is preserved by the flow of (2) and is the pullback of canonical symplectic form by Legendre transformation [17]. K-symplectic formulation. Eq. (1) can be written as the following first order ODE by introducing the momentum variable p where γ = 1 + p 2 /(m 2 0 c 2 ). We rewrite this system in a compact forṁ It is clear that dΩ = 0 for ∇ · B(x) = 0. This indicates that the skew-symmetric matrix K(z) defines a symplectic structure in (x, p). It can be proved that the phase flow φ t of the system (4) satisfies Hamiltonian formulation. Denote P the conjugate momentum defined as P = ∂L/∂ẋ, the charged particle dynamics (1) can also be described as a canonical Hamiltonian system in the coordinates (x, P) where γ = 1 + (P − eA(x)) 2 /(m 2 0 c 2 ). Let z = [x , P ] and we write the above system in the compact formż = J∇H(z), where J = 0 I −I 0 and H(x, P) = m 0 c 2 γ + eϕ(x). The canonical symplectic structure ω = dx i ∧ dP i is preserved by the flow of system (5).
Remark 1. The magnetic moment defined by µ(x,ẋ) = mẋ 2 ⊥ /2B is an important physical quantity in plasma physics. Here ⊥ represents the component which is orthogonal to the magnetic field line, and B = B(x) denotes the strength of the magnetic field. Let = ρ ∇B /B, where ρ = mẋ ⊥ /|e|B is the gyro-radius. In [25], it has been proved when E(x) ≡ 0 and the magnetic field varies slowly in space, i.e., 1, the magnetic moment µ is an adiabatic invariant which satisfies where < , > denotes the average along a gyro-period. For the case in which the magnetic field varies slowly in time w.r.t. gyro-period, the adiabatic invariance of the magnetic moment µ is also studied, and the result has been presented in [18]. Moreover, if the given magnetic field is strong and nonuniform, the authors in [1] and [7] proved that the magnetic moment µ for system (1) without considering the relativistic effects is nearly conserved over long-time. Specifically, for arbitrary N ≥ 1, there exists C N such that µ satisfies where C N is a constant which is independent of B.
As mentioned above the charged particle system (1) has some different formulations and geometric structures, thus the geometric integrators are good choices. Most geometric time integrators use constant time step size in numerical simulation while adaptive geometric integrators can be very effective from a computational point of view. We will present adaptive symplectic numerical methods for solving system (1) in the following section and the normalized physical quantities used in this paper has been listed in table 1.
3. Symplectic integrators with adaptive time steps. We start this section with introducing the basic idea of variable step size methods. Consider a differential equation Here, m 0 is the rest mass of a particle, e is the elementary charge, c is the speed of light, and B 0 is the given reference magnetic field. (In this paper m 0 = 9.1 × 10 −31 , e = 1.6 × 10 −19 , c = 3.0 × 10 8 .)

Quantities Symbols
Non-relativistic Relativistic

Units Units
Time due to a time transformation t → τ, dt/dτ = σ(z). Applying a numerical method with constant time step ∆τ to the transformed system (8) leads to a numerical discretization of the original system (7) with variable step size defined by ∆t n = (n+1)∆τ n∆τ σ(z(λ))dλ ≈ σ(z n )∆τ . This strategy is illustrated in Fig. 1. For the Hamiltonian system (5), we define a new systeṁ where E(z) = σ(z) (H(z) − H(z 0 )) with a positive function σ(z) > 0 and a given initial value z 0 . Denote z(t) andz(t) be the solutions of systems (5) and (9) respectively. Starting from the same initial value z 0 , the solutionz of (9) at τ is corresponding to the solution z of (5) at t, that isz(τ ) = z(t) with t = τ 0 σ(z(λ))dλ. Moreover, the two systems share the same invariants. Also, the above analysis holds for the K-symplectic system and the corresponding transformed systeṁ z = K −1 ∇E(z) is still a K-symplectic system. At the discretization level, applying a symplectic integrator with a constant time step ∆τ > 0 to Eq. (9) gives which yields a numerical approximation of Eq. (5) on non-equidistant grids. This is what we called adaptive symplectic method. The nonlinear function σ involved in the time transformation is called the step size function. It is usually chosen as the arclength parameterization σ(z) = J∇H(z) −1 in [11] without a priori knowledge of the system. However, with t : τ : this choice the numerical computation may not be improved apparently because of the complicated calculation for the derivatives of σ. Thus, choosing an appropriate step size function for a specific problem is crucial.
Charged particle dynamics in plasma usually involves a wide range of temporal and spatial scales, especially for the relativistic situation. The motion of the particle in an external magnetic field is a superposition of motion along B and circular motion around the guiding center with gyro-frequency f c = eB/(2πm). It is known that the gyro-period T c = 1/f c determines the smallest time scale of the particle dynamics. In this paper, we choose the function σ(z) = 1/B for non-relativistic case and σ(z) = γ/B for relativistic case, i.e., the step size is adapted by the gyro-period. These two step size functions could avoid using too large or too small step size in comparing with gyro-period and thus, can balance the stability and computational cost well. Because with time step size larger than the gyro-period, the dynamical behavior shorter than time scale of the gyro-period will be lost and the numerical computation may be unstable. While the numerical computation will be very timeconsuming if a very small step size is used. In implementation, the Lorentz factor γ = (H − eϕ)/m 0 c 2 can be replaced with γ = (H 0 − eϕ)/m 0 c 2 by the fact that the Hamiltonian is constant along the solution. Then σ can be expressed as a P−independent function σ = (H 0 − eϕ)/Bm 0 c 2 , which simplifies the expression of Eq. (9).
In this paper, we mainly focus on the adaptive symplectic methods. The dimensionless transformed Hamiltonian systems can be expressed explicitly using the above step size function σ, Then we use constant step size symplectic methods to solve the above system to realize adaptive symplectic methods for original system (7). In fact any symplectic methods, such as symplectic Euler method, Störmer-Verlet method, symplecic Runge-Kutta methods, etc, can be used and the excellent longtime behavior of the adaptive symplectic methods can be explained using the backward error analysis in the following diagram.
Remark 2. For Hamiltonian system (5), denote a variable step symplectic integrator of order r by z n+1 = Ψ(∆t n , z n ). Then the numerical solution can be taken as the solution of modified Hamiltonian systeṁ Illustration of backward analysis for adaptive symplectic method.
with the modified Hamiltonian where ∆t n is connected to ∆τ by ∆t n = (n+1)∆τ n∆τσ (∆τ, z(λ))dλ ≈σ(∆τ, z n )∆τ . More details can be found in [5]. 4. Numerical experiments. In this section, we illustrate the accuracy and efficiency of adaptive symplectic numerical methods through simulating the motion of charged particles under a given electromagnetic field. In the first experiment, we investigate the motion of a single particle in a Penning trap under three different magnetic fields. This problem is modeled by a non-relativistic Newton-Lorentz equation. For the relativistic case, we study the motion of a single charged particle under an uniform electromagnetic field, and also the dynamics of runaway electrons in tokamaks.
Experiment I: Charged particles in a Penning trap. A Penning trap [14] is a device to confine charged particles via a homogeneous axial magnetic field and a inhomogeneous quadrupole electric field. If the trap is not ideal, the magnetic field in this trap might be slightly inhomogeneous or misaligned with the geometry of the trap [14]. This device works well to measure precisely the properties of ions and stable subatomic particles, and has been used in many laboratories.
In this experiment, we take the electric field to be a static quadrupole field which is The corresponding electric field potential can be expressed as ϕ(x) = −5(  Penning trap with uniform magnetic field. We first consider the ideal Penning trap with a uniform magnetic field B(x) = [0, 0, 100] whose vector potential is given by A(x) = [0, 100x, 0] . We simulate the particle trajectory in this kind of electromagnetic field and present the results in Fig. 3. As the magnetic field is uniform, the step size function σ is in fact constant with which the adaptive numerical method is reduced to a fixed time step method. Penning trap with magnetic bottle. We then consider the motion of charged particles under a nonuniform magnetic field which describes a magnetic bottle [23]. In this case, the magnetic field has a form of and the corresponding vector potential is We employ the implicit midpoint rule with fixed and adaptive time steps resp., and present the numerical trajectories in Fig. 4. From this figure, we can see that the adaptive symplectic method can provide the numerical trajectory as perfect as the fixed symplectic method does. The relative error of the Hamiltonian H computed by the adaptive symplectic method to T = 200 is shown in Fig. 5. It can be observed that the Hamiltonian error is bounded over longtime simulation. We also compute the error of magnetic moment µ which is presented in Fig. 6. As the magnetic moment for this case is an adiabatic invariant, the result in Fig. 6 demonstrates a good conservation. In order to illustrate the efficiency of adaptive symplectic methods, we compare the implicit midpoint rule with fixed and variable time steps by calculating the iteration steps up to the same time T = 100. Fig. 7 shows that using the adaptive time steps strategy permits fewer steps than using the fixed time step in the purpose of reaching the same error. This allows the improvement of computational efficiency by 180% in this example. Penning trap with asymmetric magnetic field. For the penning trap problem, we also consider the asymmetric magnetic field described as with vector potential The particle trajectory is shown in Fig. 8. In Figs. 9 and 10, we show the time evolution of the errors of Hamiltonian and magnetic moment computed by the adaptive and fixed time step midpoint rule. It is observed that both methods can bound the errors over long time, and our adaptive method can give more accurate results than the corresponding fixed step method. We also calculate the number of iteration steps for the symplectic method with adaptive and fixed time step up to T = 100 which is shown in Fig. 11. From this figure, a significant fewer iteration steps by using the adaptive symplectic method for reaching the same error can be observed. This implies the high efficiency and accuracy for using the adaptive time steps method. For this example, the computational efficiency with the adaptive symplectic method is improved about 400%. Furthermore, the adaptive symplectic For the adaptive time step method, we investigate numerically that the numerical solution can be simulated well until ∆τ = 1.98 is taken. In Fig. 12, for the adaptive implicit rule we plot the time step ∆t n ≈ ∆τ σ n with ∆τ = 1.98. The range of time step changing from 0.0118 to 0.0254 can be investigated. With the step size limit ∆τ = 1.98, the relative error of Hamiltonian and the absolute error of magnetic moment are displayed in Figs. 13 and 14, which indicates that the numerical errors computed by our adaptive method can still be bounded over longtime simulation even with larger time step. Experiment II: Relativistic charged particle dynamics. It is known that the charged particle dynamics is a multi-scale problem which usually spans a wide range of temporal and spatial scales. Thus, an efficient and stable numerical algorithm is demanding for simulating this physical problem. For studying the dynamics of relativistic system, we use the symplectic methods with step size function σ = γ/B. Here, γ is the Lorentz factor. Motion of relativistic particle in uniform electromagnetic field. To validate the adaptive symplectic methods for relativistic case, we first perform the simulation for a single particle in a uniform electromagnetic field with B 0 = 1, E 0 = 0.005. The corresponding potentials of the electromagnetic field can be expressed as For a particle initialized with velocity v x = 10 −2 , v y = v z = 0, and position y = 0.01, x = z = 0, two symplectic integrators including implicit midpoint rule and Gauss method of order 4 are used in this test. The coefficients of the 4th order Gauss method are shown in the following table. We apply the Gauss methods of fixed step and variable step by starting with the same time step ∆t = 2π/(100B 0 ) in this test. For comparison, we plot the motion of particle over the first gyroperiod and also the last gyro-period when the simulation time reaches T = 1000 and the orbit is exhibited in Fig. 15. As the gyro-period increases along the time, using the numerical method with fixed time step needs much more iteration steps to complete the simulation for the last gyro-period orbit than for the first gyro-period orbit. This results in the dramatic increasement of computational cost using the fixed time step methods. In comparison, the iteration steps of using adaptive time step method for every gyro-period are uniform, as the adaptive strategy permits time steps changing according to the gyro-period. In table 3, we illustrate the steps for simulating the motion over various gyroperiods. For numerical methods with fixed steps, it is observed that starting with 100 iteration steps for following the motion of first gyro-period, the number of iteration steps grows rapidly to reach 5740 when the particle moves to its 150th gyro-period. In Fig. 16, we present the relativistic error of Hamiltonian computed by adaptive Gauss method and adaptive implicit midpoint rule. The adaptive Gauss method shows the smaller error than the adaptive implicit midpoint rule. Collective behavior of runaway electrons in Tokamak. Tokamak is another important device used to confine charged particles. The runaway particles are some particles which experience an unlimited "runaway" acceleration during plasma disruption and fast shutdown [3]. As the large number of runaway particles possess high energy and large velocity, they may cause a great damage to the reactor. Thus, it is vital and necessary to simulate the dynamics of runaway particles over a long time.
In paper [25], we have applied the symplectic integrators with adaptive time step to the REs dynamics. The presented numerical results have verified the high efficiency and accuracy of the adaptive symplectic method. In this section, we further investigate the collective behavior of runaway particles with a new step size function σ = γ/B which has been slightly modified from the one presented in [25]. To describe the collective behaviors of a huge number of interacting particles, we start with the following Langevin equations which can date back to [15], where N denotes the number of particles, x i ∈ R 3 and v i ∈ R 3 are the positions and velocities of particles, resp., F(·) is external force which denotes the Lorentz force in our experiment, and W i are the independent Brownian motions. The diffusion parameters β and σ are the viscosity and the thermal diffusivity coefficients, respectively, which are related by σ = βκT /m, with κ being Boltzmann constant, T the temperature of the surrounding medium and m the mass of a particle. In this experiment, we assume the interaction between particles is so small that it can be ignored, and there is no diffusion effect, i.e., β = σ = 0. Then the Langevin equation is reduced to the system of N Newton-Lorentz equations. This permits us to use the numerical algorithms introduced in Section 3. The collective behaviors of REs is simulated by using the cluster LSSC-IV 4-th in LSEC 1 up to T = 5 · 10 10 .
In practical computation, we use 3456 cores with 30 particles for each. The magnetic and electric fields in tokamak can be expressed with toroidal coordinate (r, θ, ξ) as (see [25] for more details) where R = x 2 + y 2 , ξ, and z are radial distance, azimuth, and height of the cylindrical coordinate system respectively. R 0 is the major radius, and q is the safety factor. Using e ξ , e θ , e R and e z to denote the unit vectors along toroidal, poloidal, radial, and z directions. Constants B 0 and E 0 denote the magnitude of the toroidal magnetic field and the electric field at magnetic axis respectively. Here, we take the vector potential A and scalar potential ϕ as For a given electromagnetic field, it is known that the corresponding potential is not unique. To determine uniquely the potential, we can force the gauge fixing conditions. The usually used gauges include Coulomb gauge which requires ∇ · A = 0, Lorentz gauge which needs ∇ · A + ∂ϕ/∂t = 0 and Weyl gauge with ϕ = 0 [12]. It can be calculated that the above potential A satisfies ∇ · A = 0 which implies it is the Coulomb gauge. It can be seen that the potentials in Eq. (16) do not depend on time. We can also achieve the time-dependent potentials for the above electromagnetic field in tokamak by requiring the Weyl gauge condition [20]. In this case ϕ = 0, and the magnetic potential A is taken as A = A − Et with A and E defined in (16) and (15). This gives E = − ∂ A ∂t which implies that the electric field can be considered as an inductive field induced from loop voltage.
In this experiment, we use the initial values used in [16] which are sampled uniformly from 0 to 2π yielding a quadratic profile. On the poloidal plane, the distribution function of initial values is where r b = 2m is the outer boundary of the initial distribution. The initial kinetic energy is sampled as a normal distribution with expectation µ = 47.5Mev and standard deviation σ = 0.25Mev. The initial toroidal angle ξ 0 and poloidal angle θ 0 are sampled uniformly from 0 to 2π. The initial pitch-angle is uniformly sampled from 0 to 0.3, and the initial gyro-phase is distributed uniformly from 0 to 2π. For the adaptive symplectic integrators, the relative energy error keeps bounded in a small region. This fact is clearly demonstrated in Fig 17, where the energy error for a single runaway electron is plotted. To reveal the collective behaviors of runaway electrons, the evolutions of density distribution for the REs in the poloidal plane are displayed in Figs. 18 and 19. The difference between Fig. 18 and Fig. 19 is that the directions of two electromagnetic fields are opposite. It shows that the REs in an ideal configuration keep drifting outwards and finally strike the first wall of the device. This is consistent with experimental observations. It is also observed that changing the direction of electromagnetic field leads to the change of drift direction of the REs. In Fig. 20, we demonstrate the high efficiency of the adaptive symplectic method with step size function σ = γ/B. It shows that the computational efficiency of using the adaptive symplectic methods in our experiment can be improved about 10 times compared with the fixed time step method. 5. Concluding remarks. It is well-known that the geometric numerical integrators play an important role for simulating the Hamiltonian systems. In this paper, we combine the symplectic methods with an adaptive time step strategy for solving the charged particle dynamics. The adaptive time step strategy used in this paper is proposed in [6]. It can be realized by using a time transformation, i.e., a step size function. Here, we propose a new kind of step size function by following the dynamical nature of the charged particle system. This step size function depends on the gyro-period of particles. The new derived adaptive symplectic method has been applied to investigate the motion of one single charged particle in Penning trap, and the long-term dynamics of more than 10 5 runaway electrons in tokamaks. The numerical experiments show the good conservation for the invariants of system including energy and magnetic moment, and also illustrate the high efficiency of our numerical methods.