On Symmetrical Methods for Charged Particle Dynamics

In this paper, we use the scalar auxiliary variable (SAV) approach to rewrite the charged particle dynamics as a new family of ODE systems. The systems own a conserved energy. It is shown that a family of symmetrical methods is energy-conserving for a new ODE system but may not be for the original systems. Moreover, the methods have high-order accuracy. Numerical results are given to confirm the theoretical findings.


Introduction
In this paper, we are interested in investigating high-order numerical methods for solving the following charged particle dynamics [1][2][3][4]: q = p,ṗ = p × L(q) − ∇U(q), q(0) = q 0 , p(0) = p 0 , (q 0 , p 0 ) ∈ Ω ⊆ R 3 × R 3 , (1) where L(q) and −∇U(q) are, respectively, the magnetic and electric fields. In recent years, many energy-conserving numerical methods have been investigated within the framework of line integral methods [1,5] or Runge-Kutta methods [6][7][8][9]. The main results indicate that the method preserves polynomial Hamiltonians. For non-polynomial Hamiltonians, the energy errors appear in the long time simulation. As a result, the numerical methods may not exactly be energy-conserving. In this short paper, we propose a novel approach to developing some energy-conserving numerical methods for the charged particle dynamics. The main idea can be viewed as an extension of the recent scalar auxiliary variable approach [10][11][12][13][14] or the invariant energy quadratization approach [15]. Firstly, we introduce an auxiliary variable and obtain an auxiliary ordinary differential equation. Together with the original equations (1), we obtain a new family for the ODE system. The new system has the same solution as the original problem (1), and the energy of the new system can be rewritten as the quadratic invariants. Then, we studied the energy-conserving methods in the framework of Runge-Kutta methods. It can be seen that certain Runge-Kutta methods are energy-conserving for the ODE system, while the same Runge-Kutta methods cannot generally preserve the energy of the original system (1). Within the framework of Runge-Kutta methods, it is easy to find some other energy-conserving methods, e.g., Gauss methods and Hamiltonian boundary value methods [6][7][8][9].
The rest of the paper is arranged as follows. In Section 2, we propose the main approach and energy-conserving methods. In Section 3, we present several numerical results to confirm the effectiveness of the numerical methods. In Section 4, we summarize our results.
Together with (1), we have a new ODE system: where the initial condition of r is given by r(0) = U(q 0 ) + C 0 . The system has the same solutions for p and q, and the same conserved and continuous energy as (1), meaning that: This is true because, by the chain rule, it holds that: Let: Then, system (2) can be viewed as the standard ordinary differential equations, having the form:ẏ = f (y). Let: One can directly check that: and : where we have noted (3). We now consider an s-stage Runge-Kutta method to solve problem (4). Let (A, b, c) be a given Runge-Kutta method with vectors b = (b 1 , b 2 , · · · , b s ) T , c = (c 1 , c 2 , · · · , c s ) T and matrix A = (a ij ) s×s . Let h be the stepsize and define grid points t n by t n = nh, n = 0, 1, 2, · · · . Applying the s-stage Runge-Kutta method to solve problem (4) yields: where f j = f (Y n j ), y n and Y n i are the numerical approximations to y(t n ) and y(t n + c i h), respectively.
Proof. By (9) and definition (5), it holds that: where in the second equality, we have noted relation (8).

Remark 1.
The key in the proposed approach is to add an auxiliary ordinary differential equation. Then, the energy can be rewritten as a quadratic invariant. If the coefficients of a Runge-Kutta method satisfy (10), the method conserves quadratic invariants. For the quadratic-invariants conserving Runge-Kutta methods, we also refer readers to [6,7] or [8] for more details.

Remark 2.
The new system has the same solution as the original problem (1). However, we show that the proposed Runge-Kutta methods are energy-conserving for the ODE system, while the same Runge-Kutta methods cannot generally preserve the energy of the original system (1).

Remark 3.
The theoretical results can be similarly extended to the general Hamilton system of the form: where the Hamiltonian H(p, q) represents the total energy of the system. It is easy to confirm that H(p, q) = Const, along the solution curves of (13). By the proposed approach, a scalar variable r can always be introduced to modify the energy to a quadratic form. In this way, we have a new modified system by letting y = [p T , q T , r] T as follows: with the energy H(y) = y, y S , whereS is a symmetric matrix. Then, some symmetrical methods preserve the energy of (14).

Numerical Experiments
In this section, we present some numerical results to confirm the theoretical findings.
We solve problems (1) and (2) by using the s-stage Gauss method with the stepsize h = 0.1. The discrepancies of the discrete energy are shown in Figure 1, where s-Gauss stands for the s-stage Gauss method here and below. Clearly, by solving (2), the discrepancies of the discrete energy are of the order of machine precision. However, by solving (1) with the same stepsize, the discrepancies of the discrete energy depend on the stepsizes and the orders of the methods. Then, we take the stepsize h = 1/8000 and consider the obtained numerical solutions as the reference solutions. We test the convergence of the numerical methods and list the maximum errors at t = 1 in Table 1. The results imply that the s-stage Gauss method is of the order of 2 s. These numerical results confirm the theoretical results and effectiveness of the methods.
The model has an important application in the investigation of the single particle motion in a non-uniform and static electromagnetic field [1,2]. We numerically solve problem (2) by using the s-stage Gauss methods. We show the maximum errors at t = 1 and convergence orders in Table 2, where the reference solutions are obtained by using the stepsize h = 1/8000. Then, we set the stepsize h = π/10 and solve problem (2) in the time interval [0, 10 3 π]. Figure 2 shows that the motion occurs in the (q 1 , q 2 )-plane. Moreover, Figure 3 shows that the discrepancies of discrete energy are of the order of machine precision. We also solved the original problem (1) by using the same s-stage Gauss method. From Figure 3, one can see that the discrepancies of discrete energy depend on the stepsizes and the orders of the methods. The numerical results further confirm the theoretical results of the present paper.
We numerically solve problem (2) by using the s-stage Gauss methods. We show the maximum errors at t = 1 and convergence orders in Table 3, where the reference solutions are derived by using the stepsize h = 1/8000. Then, we set stepsize h = 0.1 and solve problem (2). Figure 4 further shows that the discrepancies of the discrete energy are again of the order of machine precision. However, by solving the original problem (1) with the same method and stepsize, one can see from the figure that the discrepancies of the discrete energy depend on the stepsizes and the orders of the methods. The numerical results again confirm the theoretical results.  Example 4. The Hénon-Heiles Model is a well-known and frequently used model in certain fields of physics. For example, Hénon-Heiles equations arise in the description of a moving star in the axisymmetric potential of the galaxy. Consider the model mentioned in [7]: with the initial values: where p i , q i , i = 1, 2, are the i-th coordinates of p, q, respectively. This model has the conserved Hamiltonian: Using the approach, we let: We have the new system: with the energy: Firstly, we numerically solve model (17) by using s-stage Gauss methods. The maximum errors at t = 1 and convergence orders are shown in Table 4, where we obtain the reference solutions by Gauss methods with h = 1/8000. Secondly, we set stepsize h = 0.1 and solve the model (15) and (17) in the time interval [0, 1000]. Figure 5 shows the discrepancies of the discrete energy. We can see that the discrepancies are of machine precision by applying Gauss methods to system (17), while the same methods cannot preserve the energy of the original system (15). One can also see that the discrepancies depend on the orders of methods for the original system.
This system owns the Hamiltonian: Letting: which leads to the new system: with the energy: We use different Gauss methods to solve system (20) to test the accuracy of our schemes. We present the maximum errors at t = 1 and the convergence orders with the reference solutions calculated by stepsize h = 1/8000. The corresponding results are shown in Table 5, from which we can see that the accuracy property is held. Then, we set h = 0.1 and solve (19) and (20) by different Gauss methods. We show the discrepancies of the discrete energy in Figure 6, in which we can see that our schemes are energyconserving while the same methods are not for the original system.

Summary
In this paper, some energy-conserving numerical schemes were developed by introducing an auxiliary variable and an auxiliary ordinary differential equation. This approach, as an extension of the scalar auxiliary variable approach, was developed and applied for the charged particle dynamics. The new system has the same solution as the original system. It was shown that some symmetrical methods are energy-conserving for the new system, while the same methods cannot generally preserve the energy of the original system. Our numerical experiments confirmed the theoretical findings. Furthermore, the approach can be extended to other ODE systems or time-dependent partial differential equations.