HIGH ORDER GAUSS-SEIDEL SCHEMES FOR CHARGED PARTICLE DYNAMICS

. Gauss-Seidel projection methods are designed for achieving desirable long-term computational eﬃciency and reliability in micromagnetics simulations. While conventional Gauss-Seidel schemes are explicit, easy to use and furnish a better stability as compared to Euler’s method, their order of accuracy is only one. This paper proposes an improved Gauss-Seidel method-ology for particle simulations of magnetized plasmas. A novel new class of high order schemes are implemented via composition strategies. The new algorithms acquired are not only explicit and symmetric, but also volume-preserving to- gether with their adjoint schemes. They are highly favorable for long-term computations. The new high order schemes are then utilized for simulating charged particle motions under the Lorentz force. Our experiments indicate a remarkable satisfaction of the energy preservation and angular momentum conservation of the numerical methods in multi-scale plasma dynamics compu-tations.


1.
Introduction. Plasma is a collection of charged particles interacting with electromagnetic fields, whose sources can be either external or internal [1]. The most fundamental physical process in collective dynamics of magnetized plasmas is the motion of charged particles under the influence of electromagnetic fields. Plasma characteristics can be understood and analyzed in terms of a single particle motion satisfying the Lorentz force equation. Efficient numerical methods for simulating such single particle models may significantly benefit the study of plasma dynamic behaviors. Discussions in the topic can be found in numerous recent publications (see [16,22,24,9] and references therein, for example).
In the study of single particle motions and guiding center dynamics [12], highly stable integrations are necessary. Further, features to preserve certain geometric 2. The motion equation under Lorentz force. The dynamic motion of a charged particle in an electromagnetic field is governed by the Newton-Lorentz equation where x is the position of the particle, m is the mass, and q denotes the electric charge. For the convenience, we may assume that B and E are static. Thus, B = ∇ × A and E = −∇ϕ with A and ϕ being the potentials. Eq. (1) is equivalent to an Euler-Lagrangian equation with Recall [17,15,22]. Let the conjugate momentum p = mẋ+qA(x). Then the system (1) is with a Hamiltonian which enables the use of canonical symplectic algorithms. Since (3) is not separable, that is, it cannot be split into the form H(x, p) = T (p)+U (x), a general symplectic method for solving (1) cannot be explicit [8,4,23]. However, note that Eq. (1) can be separated naturally to where the energy is H(x, v) = mv · v/2 + qϕ(x). Letting z = x v , the system (4) can be rewritten as where K(z) = −qB(x) −mI mI 0 is a skew-symmetric matrix witĥ There is a non-canonical symplectic structure (K-symplectic structure) preserved by the flow of (4) [9,14]. It has been a difficult task to construct an explicit method that preserves the non-canonical structure K(z), and in the same time, is more efficient than symplectic method [10]. Denote According to Liouville's theorem, the corresponding volume is invariable along the solution flow of (4). A new class of volume-preserving methods was derived recently via splitting strategies [9].
3. Symmetric and volume-preserving methods. We consider a first order Gauss-Seidel scheme and its adjoint for (4). Based on them, a new class of explicit, symmetric and volume-preserving Gauss-Seidel schemes can be accomplished via compositions.
To this end, we rewrite (4) as Based on the Gauss-Seidel scheme for Eq. (7), we have Its adjoint method Φ * h can be comprised to Apparently, schemes (8) and (9) can be computed explicitly. Let us first show that (9) is volume-preserving. For this, we denote (9) as z k+1 = Φ * h (z k ), z k = [x k ; v k ]. Theorem 3.1. The scheme (9) is volume-preserving.
Proof. To see this, we only need to prove According to properties of a determinant of a square matrix, we have where and c = hq/m and B i = B i (x k ), i = 1, 2, 3; k = 0, 1, 2, . . . Therefore, This completes our proof.
In fact, Gauss-Seidel methods of any order can be constructed through compositions of Φ h or Φ * h [8]. Targeting at Lorentz force systems, the following result offers a simple and straightforward way to construct high order symmetric Gauss-Seidel methods.
Theorem 3.2. Based on first order Gauss-Seidel schemes (8) and (9), we have the following volume-preserving and symmetric formulas.
1. Second order scheme: Proof. a) It is firstly proved that the adjoint method has the same order as the original method, and, with a possible sign change, also the same leading error term. For the autonomous differential equationṡ let ϕ t be the exact flow and let Ψ h be a one-step method of order p satisfying Let e = y 0 − Ψ −h ϕ h (y 0 ) . It follows (14) that The . Using the Taylor expansion, we have Furthermore, if the method is symmetric, its (maximal) order is even. Because Ψ h = Ψ * h implies C(y 0 ) = (−1) p C(y 0 ), and therefore C(y 0 ) can be different from zero only for even p.
b) The adjoint method satisfies the usual properties such as h/2 is symmetric. Next, we prove that the method GS 2 h is of order 2. Based on part (a), we can assume that the first-order Gauss-Seidel method Φ h and its Let (17)-(18) that Similar to the proof of part (a), we have Therefore, we obtain which implies that the method GS 2 h is of order 2. c) The argument of order of the composition method GS h is symmetric and of order 2l. Now we prove that GS 2(l+1) h is symmetric and of order 2(l + 1). Since GS 2l h is symmetric, i.e., GS 2l is also symmetric. Assume that the 2l-order method In addition, we have Notice that 2α l + β l = 1 and 2α 2l+1 l + β 2l+1 l = 0, we obtain which implies that the method GS 2(l+1) h is at least of order 2l + 1. Because the method GS 2(l+1) h is symmetric, so the order of GS 2(l+1) h is up to 2l + 2. Note that, since the methods Φ h and Φ * h is volume-preserving, any composition method based on Φ h and Φ * h is still volume-preserving. The aforementioned explicit methods can be effective and efficient for long-term simulations of the multi-scale dynamics of plasmas. 4. Numerical experiments. We consider two respective cases in this section. 4.1. 2D dynamics in a static electromagnetic field. Firstly, we consider the 2D dynamics of a charged particle in a static electromagnetic field where the potentials are chosen to be A = [−y/2, x/2, 0] T , ϕ = 10 −2 R with R = x 2 + y 2 . In this example, the physical quantities are normalized by the system size a, the characteristic magnetic field B 0 , and the gyro-frequency ω 0 ≡ qB 0 /m of the particle. The Lagrangian is L = (Ṙ 2 + R 2ξ2 +ż 2 )/2 + R 2ξ /2 − ϕ in cylindrical coordinates (R, ξ, z). Due to ∂L/∂ξ ≡ 0, the angular momentum is invariant along the solution trajectory of system (4). It is known that the energy Starting from the initial position x 0 = [0, −1, 0] T with the initial velocity v 0 = [0.1, 0.01, 0] T , the analytic orbit of the charged particle is a spiraling circle with constant radius. The large circle corresponds to the ∇ · B drift and the E × B drift of the guiding center, and the small circle is the fast-scale gyro-motion.
We first apply the Runge-Kutta method of order four denoted by RK4 and the numerical results are shown in Fig. 1. The step size is h = π/10 which is the 1/20 of the characteristic gyro-period. Though RK4 has a fourth order of accuracy, the numerical error accumulation after 2.7 × 10 5 steps gives rise to a complete wrong solution orbit in Fig. 1(b), where the gyro-motion is numerically dissipated. The normalized errors of the invariants as a function of integration time are performed in Fig. 1(c) and (d). These figures show that the RK4 cannot preserve the invariants inherited by this system.
Next, we test the symmetric and volume-preserving numerical methods. Fig. 2 shows the numerical orbits of the second order methods after 5×10 5 steps in (a), and that of the fourth order methods after 2.5 × 10 5 steps in (b). The volume-preserving methods provide a correct and stable gyro-motion over a very long integration time.
To test accuracy of the methods, the exact solution is computed by the RK4 with a small step size h = 0.0001. In Fig. 3, the solution errors ||[x n , v n ]−[x(t n ), v(t n )]|| 2 are computed at time t = 30. It is shown that the method GS 2 h and the Boris methodG 2 h (Ref. [9]) have the convergence rate of order two; the method GS 4 h and the splitting method G 4 h (Ref. [9]) have the convergence rate of order four. It is also noticed that the method GS 2 h (or G 4 h ) is more accurate than the Boris method (or GS 4 h ). In Fig. 4, we plot the relative errors of the invariants H and p ξ as functions of the step size h. We observe that the invariant errors are about a scale of h 2 for h is best in preserving energy, while the method GS 4 h is best in preserving angular momentum.
In Fig. 5, we illustrate the relative errors of the energy H and the angular momentum p ξ computed by the four methods GS 2 h ,G 2 h , GS 4 h and G 4 h . It is observed that the errors are bounded over a long integration time. The errors of invariants become smaller when the higher order methods are used.

4.2.
2D dynamics in an axisymmetric tokamak geometry. In this subsection, we consider the motion of a charged particle in a 2-dimensional axisymmetric tokamak geometry without inductive electric field. The magnetic field in the toroidal coordinates (r, θ, ξ) is expressed as where B 0 = 1, R 0 = 1, and q = 2 are constant with their usual meanings. The corresponding vector potential A is chosen to be To apply the proposed methods in Section 3, we transform B (29) into the Cartesian coordinates (x, y, z) which is In this example, the physical quantities are normalized by the system size a, characteristic magnetic field B 0 , and the gyro-frequency ω 0 ≡ qB 0 /m. Here the Lagrangian is L = (Ṙ 2 + R 2ξ2 +ż 2 )/2 +Ṙz/(2R) + (R − 1) 2 + z 2 ξ /4 −ż ln R/2 in cylindrical coordinates (R, ξ, z). Due to ∂L/∂ξ ≡ 0, the angular momentum p ξ = R 2ξ + (R − 1) 2 + z 2 /4 is invariant along the solution trajectory of system (4). In the absence of electric field, the energy of the system is H = v · v/2, which   is a constant of motion. With the conserved quantities, the solution orbit projected on (R, z) space forms a closed orbit. Starting with the initial position x 0 = [1.05, 0, 0] T and the initial velocity v 0 = [0, 4.816e − 4, −2.059e − 3] T , the orbit projected on (R, z) space is a banana orbit, and it will turn into a transit orbit when the initial velocity is changed to v 0 = [0, 2 × 4.816e − 4, −2.059e − 3] T . Setting the step size h = π/10 which is the 1/20 of the gyro-period, we adopt the RK4 and the volume-preserving methods developed above. The simulation over 5 × 10 5 steps is shown in Fig. 6. For the results generated by the RK4 in Fig. 6(a) and (b), due to the numerical damping of energy the numerical orbits are not closed, and the gyro-motion is dissipated. The banana orbit gradually transformed into a circulating orbit and the transit orbit deviates to the right side. By comparison, in Fig. 6(c)-(d) the volume-preserving methods give correct and qualitatively better orbits over a very long integration time. We display the relative errors of the invariants H and p ξ in Fig. 7, which show that the proposed methods conserve the invariants well.  5. Concluding remarks. In this paper, we prove that both the Gauss-Seidel scheme and its adjoint scheme for the Lorenz force system are volume-preserving schemes. Actually they can also be derived by splitting the original system into several source-free solvable subsystems [7]. The Gauss-Seidel scheme is explicit and volume-preserving, but it is only of order 1. To improve order of accuracy of the Gauss-Seidel scheme, a class of high-order Gauss-Seidel schemes are constructed by the theory of composition methods. These methods are explicit, symmetric and volume-preserving, and thus are very effective for the long-term simulation of the multi-scale dynamics of plasmas. By numerical experiments, we verify the order of accuracy of the proposed schemes, supporting Theorem 3.2. Numerical experiments also confirm that the high-order Gauss-Seidel methods could conserve conservative properties well and have excellent performance in a long time computing. It is noting that the proposed high-order methods are explicit when they are applied for the separable source-free system.
Another new second-order Gauss-Seidel scheme can be constructed by changing the combinational order of the Gauss-Seidel scheme and its adjoint method, i.e.,

GS
Based on the second-order Gauss-Seidel scheme GS 2 h , the fourth-order method can be obtained in terms of the Theorem 3.2 where α 1 = 1/(2 − 2 1/3 ) and β 1 = 1 − 2α 1 < 0. By numerical calculation, the schemes (32) and (33) perform the same good long-term simulation behaviour as the proposed methods in Section 3. Thus this paper provides a series of efficient schemes for long time simulation of charged particle motion under the Lorentz force.
In addition, we can also apply other symmetric composition methods to develop new high-order methods (see II.4 and V.3.2 of [8]). In order to compute efficiently in a long time, Wang et al. [21] developed a Gauss-Seidel projection method for micromagnetics simulations. Based on the first-order splitting model, they carried out a combination of a Gauss-Seidel implementation of a fractional step implicit solver for the gyromagnetic term and the projection method for the heat flow of harmonic maps. However, their method is only of order 1. In fact, we believe that a second-order Gauss-Seidel projection method can be constructed by the Strang splitting and combining the second-order Gauss-Seidel scheme and the second-order projection method of Ref. [5]. It will be discussed carefully in our future work.