High order explicit Lorentz invariant volume-preserving algorithms for relativistic dynamics of charged particles

Abstract

Lorentz invariant structure-preserving algorithms possess reference-independent secular stability, which is vital for simulating relativistic multi-scale dynamical processes. The splitting method has been widely used to construct structure-preserving algorithms, but without exquisite considerations, it can easily break the Lorentz invariance of continuous systems. In this paper, we introduce a Lorentz invariant splitting technique to construct high order explicit Lorentz invariant volume-preserving algorithms (LIVPAs) for simulating charged particle dynamics. Using this method, long-term stable explicit LIVPAs with different orders are developed and their performances of Lorentz invariance and long-term stability are analyzed in detail.

Introduction

Structure-preserving algorithms possessing long-term stabilities play key roles in simulations of various fields [1], [2], [3], [4], [5], [6], [7]. In plasma studies, these advanced schemes have been applied in subjects including the secular dynamical simulations of charged particles [8], [9], [10], the long-term analysis of plasma kinetic processes [10], [11], [12], [13], [14], [15], the analysis of magnetohydrodynamical phenomena [16], [17], [18], and the simulations of nonlinear processes [19], [20]. The basic idea of constructing structure-preserving algorithms is to keep the fundamental geometry structures of ordinary differential equations (ODEs) or partial differential equations (PDEs) during discretization [3], [21], [22]. In other words, the one-step maps, approximating the continuous exact solutions, should inherit mathematical natures of original continuous systems, such as the volume-preserving property of source-free systems and the symplectic-preserving property of Hamiltonian systems [21]. These properties work as constrains for numerical systems, which can limit the global accumulations of numerical errors during long-term iterations [21], [22].

Explicit schemes, compared with implicit schemes, are more convenient for implementations and more efficient especially in cases of complex vector fields. Because the symplectic Runge-Kutta methods are often implicit [21], the Hamiltonian splitting technique has been applied widely to build explicit structure-preserving algorithms for charged particle dynamics. Written in canonical coordinates, the Hamiltonian equation of charged particles can be expressed as $J\dot{\mathbf{Z}}={\mathrm{\nabla }}_{\mathbf{Z}}H$, where Z is the canonical coordinate, J is the symplectic structure, and H is the Hamiltonian. Through the sum-split procedure of H, one can obtain canonical symplectic subsystems which can be discretized via standard symplectic methods, like the generating function method [23], [24]. Then the composites of subsystem algorithms are the canonical symplectic algorithm for the original system [25]. On the other hand, it has been find that the Lorentz force equation of charged particles possesses non-canonical symplectic structure, which can be written as $K\dot{\mathbf{z}}={\mathrm{\nabla }}_{\mathbf{z}}H$, where z is phase space coordinate and K denotes the K-symplectic structure. Similarly, the sum-split procedure of Hamiltonian produces the explicit K-symplectic algorithms [22], [26], [27]. Compared with the symplectic methods, volume-preserving algorithms (VPAs) impose looser constrains on discrete systems. However, VPAs still possess significant long-term stability and have been used widely in Particle-in-Cell simulations [28], [29], [30], [31], [32]. The volume-preserving systems can be expressed as $\dot{\mathbf{z}}=\mathbf{V}$, where the source-free vector field V satisfies ${\mathrm{\nabla }}_{\mathbf{z}}\cdot \mathbf{V}=0$. Through splitting V into several source-free sub-vectors, explicit VPAs can also be conveniently constructed [31], [32], [33], [34], [35], [36].

For relativistic dynamical systems, the Lorentz invariance is another fundamental property, which should also be kept after discretization. When the effects of the special relativity are considered, simulating a physical process in different Lorentz inertial frames can minimize the range of time and space scales and thus reduce the cost of calculation [37], [38]. In such cases, the Lorentz invariance of algorithms becomes very important. In 2008, Vay has pointed out the importance of Lorentz invariance for relativistic particles and constructed a new scheme that preserves the $\mathbf{E}\times \mathbf{B}$ velocity in different frames [38]. In 2017, Higuera and Cary improved Vay's method and built a VPA that also preserves $\mathbf{E}\times \mathbf{B}$ velocity [31]. The reference-independency of the numerical $\mathbf{E}\times \mathbf{B}$ velocity can be treated as one key outcome of Lorentz invariance of algorithms. Strictly speaking, the Lorentz invariant algorithms (LIAs) should produce reference-independent discrete numerical results when solving the same process in arbitrary Lorentz inertial frames. Equivalently, the difference equations of LIAs keep unchanged after all observables are transformed to that in another Lorentz inertial frame. Therefore, the definition of LIAs can be given as follows [39]. For a Lorentz invariant system F, an algorithm, $\mathcal{A}$, is Lorentz invariant if it satisfies that

$${\mathcal{D}}_{\mathcal{A}}\circ {\mathcal{T}}_L\mathbf{F}={\mathcal{T}}_L\circ {\mathcal{D}}_{\mathcal{A}}\mathbf{F}$$

where, ${\mathcal{D}}_{\mathcal{A}}$ represents numerical discretization operation using $\mathcal{A}$, ${\mathcal{T}}_L$ denotes the transformation operation on variables based on the Lorentz matrix L, and “∘” is the composite operator. As an element of Lorentz group, L satisfies $Lg{L}^T=g$, where g is the Lorentz metric tensor [40].

The Lorentz invariance and the structure-preservation are two independent aspects for constructions of advanced algorithms. Symplectic-preserving or volume-preserving algorithms sometimes break the Lorentz invariance, while traditional algorithms like the Newton method and the Runge-Kutta method can produce Lorentz invariant schemes if they are applied to Lorentz invariant dynamical equations written in the 4-dimensional spacetime. Furthermore, if we combine the benefits of the structure-preservation and the Lorentz invariance, the resulting structure-preserving LIAs will possess reference-independent secular stability, which is very important for simulating relativistic multi-scale processes. To construct Lorentz invariant algorithms, one straightforward way is directly discretizing the Lorentz invariant dynamical equations in which all variables are written as geometric objects in 4-dimensional spacetime [40]. The integrity of these geometric objects, such as the spacetime coordinates and 4-dimensional tensors, should be kept during discretization [39]. For example, for the 4-dimensional Lorentz invariant Hamiltonian equation of charged particles [40], the symplectic-Euler method gives an explicit 1-order symplectic LIA [39]. Similarly, symplectic Runge-Kutta methods, such as the implicit mid-point symplectic scheme, also generate symplectic LIAs when applied to this system.

Although symplectic LIAs with different orders can be constructed using the symplectic Runge-Kutta methods, they are all implicit and their usage are not convenient. On the other hand, although Hamiltonian splitting technique can produce explicit high-order symplectic algorithms, it has to divide the Hamiltonian into several pieces to construct symplectic schemes for subsystems [23], [24], [25], [26]. The fine-grained splitting breaks the Lorentz invariance. Therefore, in this paper, we relax the constrain of symplectic-preservation and focus on the constructions of high order explicit Lorentz invariant volume-preserving algorithms (LIVPAs). Higher order schemes mean faster convergence rate, but, unfortunately, more CPU time consumed for one-step iteration. Generally speaking, the balance between the complexity and the convergence rate should be considered during practical simulations. Though it's not necessary to build algorithms with arbitrary high orders, schemes with order higher than 1 are still needed. From the view point of applications, the convergence rates of 1-order algorithms are sometimes too slow to be applied. Therefore, although a little complex than 1-order algorithms, 2-order and 4-order schemes, like the Crank-Nicolson method (2-order), the Boris method (2-order), and the 4-order Runge-Kutta method, have been widely applied in plasma studies. Especially, for simulations of processes with small scales, such as turbulence process and magnetic reconnection processes, high temporal and spatial resolutions are necessary. Though costing more CPU time for one-step iteration, high order schemes, with larger step length, can be cheaper than low order schemes.

Without exquisite considerations, the splitting technique can easily break the Lorentz invariance of the original system, which will be discussed in detail in Sec. 2. In this paper, however, we find a Lorentz invariant splitting procedure for constructing LIVPAs for charged particle dynamics. The splitting method is chosen for our kernel technique for two main reasons. First, through different composite of subsystem algorithms, high-order explicit schemes can be easily constructed [21]. Second, compared with the original systems, finding Lorentz invariant algorithms for subsystems is much easier, and it is also readily to understand that the composed scheme is Lorentz invariant if algorithms of all subsystems are Lorentz invariant. Although we split the original system into several pieces, the Lorentz invariance of algorithms still remains. Explicit LIVPAs with different orders are constructed and tested in detail. Compared with the explicit symplectic algorithms constructed via Hamiltonian splitting method, these LIVPAs can give reference-independent numerical results [24]. Compared with the Vay [38] and the Higuera-Cary [31] schemes that preserve the $\mathbf{E}\times \mathbf{B}$ velocity, LIVPAs show Lorentz invariance that is independent with the configurations of step length. Moreover, because of the preservation of phase space volume, LIVPAs show better secular stability than the Runge-Kutta method.

The rest part of this paper is organized as follows. In section 2, basic properties of Lorentz invariant algorithms and their relations with splitting method are discussed by use of a simple system. In section 3, the invariant form of the Lorentz force equation is analyzed. We introduce the Lorentz invariant splitting method and construct the LIVPAs in Sec. 4. In section 5, LIVPAs are applied in a typical field configuration and their numerical performances are studied. Finally, we conclude this paper in Sec. 6.