Quadratic conservative scheme for relativistic Vlasov--Maxwell system

For more than half a century, most of the plasma scientists have encountered a violation of the conservation laws of charge, momentum, and energy whenever they have numerically solve the first-principle equations of kinetic plasmas, such as the relativistic Vlasov--Maxwell system. This fatal problem is brought by the fact that both the Vlasov and Maxwell equations are indirectly associated with the conservation laws by means of some mathematical manipulations. Here we propose a quadratic conservative scheme, which can strictly maintain the conservation laws by discretizing the relativistic Vlasov--Maxwell system. A discrete product rule and summation-by-parts are the key players in the construction of the quadratic conservative scheme. Numerical experiments of the relativistic two-stream instability and relativistic Weibel instability prove the validity of our computational theory, and the proposed strategy will open the doors to the first-principle studies of mesoscopic and macroscopic plasma physics.


Introduction
The relativistic Vlasov-Maxwell system has been regarded as the first-principle equations of weakly coupled plasmas. The relativistic Vlasov equation is where f = f (t, x, p) is the distribution function; t, x = [x, y, z] T , and p = [p x , p y , p z ] T are the time, space, and momentum, respectively; m and q are the particle mass and charge, respectively; E = [E x , E y , E z ] T and B = [B x , B y , B z ] T are the electric and magnetic field respectively; and c is the speed of light in vacuum. γ is the Lorentz factor described as This equation is coupled with the governing equations for the electromagnetic field: Maxwell's equations in Gaussian-cgs units: quadratic conservative schemes are based on some type of product rule in discrete form. Although the original Morinishi scheme was composed for incompressible fluid dynamics, this strategy has been extended to compressible fluid dynamics [20], and many other quadratic conservative schemes for compressible fluid dynamics have been proposed [21,22,23,24,25]. Further, the product rules in discrete form are useful in constructing conservative numerical methods for hyperbolic hydrodynamic equations in non-conservative formulation [26,27,28]. Such a strategy has been called a "structure-preserving" theory [29]. A conservative algorithm should be constructed for the relativistic Vlasov-Maxwell system using the structure-preserving strategy.
In this article, a quadratic conservative scheme is proposed for a relativistic Vlasov-Maxwell system, which is based on the finite-difference method and strictly maintains the conservation laws of charge, momentum, and energy. In Sec. 2, the quadratic conservative scheme for the relativistic Vlasov-Maxwell system is proposed. The theoretical proof of Gauss's law, solenoidal constraint of the magnetic field, and conservation laws of charge, momentum, and energy is given in Sec. 3. Some mathematical formulae used in this study are also introduced. Section 4 describes the experimental demonstration of the conservation property via the relativistic two-stream instability and relativistic Weibel instability. Section 5 gives the conclusions of this article.

Structure-preserving theory for relativistic Vlasov-Maxwell system
Before the quadratic conservative scheme for the relativistic Vlasov-Maxwell system is introduced, some important strategies for constructing the proposed scheme are described here. When proving the conservation laws of momentum and energy, the product rule and integration-by-parts are required both in differential and discrete forms. In addition, the commutative property of finite-difference operators is required to derive Gauss's law and solenoidal constraint of the magnetic field. Therefore, the finite-difference operators should have a linearity. Accordingly, we employed the Crank-Nicolson method for the temporal difference, and a 2nd-order central difference method for the spatial and momentum differences.
The quadratic conservative discretization of the relativistic Vlasov equation is where n, M y ], and j 3 ∈ [1, M z ] are the indices of t, x, y, z, p x , p y , and p z , respectively. ∆t, ∆x, ∆y, ∆z, ∆p x , ∆p y , and ∆p z are the grid intervals of t, x, y, z, p x , p y , and p z , respectively. The finite-difference operators and interpolation operators are defined as follows: where F is an arbitrary function. Moreover, the distribution function f must be maintained at the round-off level near the momentum boundaries: Therefore, the computational domain of momentum should be large enough to maintain Eqs. (8)-(10) Maxwell's equations (3) and (4) must be discretized as follows: When solving Eqs. (11)- (16), the current density J should be obtained from the distribution function f . Furthermore, Gauss's law Eq. (5) is required to derive the law of momentum conservation, which is associated with the charge density ρ. Therefore, these quantities are defined as follows: where ∆V = ∆p x ∆p y ∆p z . Note that the domain of summation covers the entire computational domain: These constitute the complete set of our quadratic conservative scheme.

Proof of conservation property
This section gives the proof of the exact conservation properties of charge, momentum, and energy for the proposed discretization method in the previous section. First, the discrete product rule, summationby-parts, and commutative laws of finite-differential operators are derived in Sec. 3.1. The law of charge conservation is derived in Sec. 3.2. In Sec. 3.3, Gauss's law and the solenoidal constraint of the magnetic field are obtained to prove the law of momentum conservation. The laws of momentum and energy conservation are derived in Sec. 3.4 and Sec. 3.5, respectively.

Mathematical basis
In this article, the following finite-difference operator is defined to prove the conservation laws of momentum and energy: where F and G are the arbitrary functions. The equivalent operator has been used to construct kineticenergy-preserving schemes [19]. A product rule for the momentum dimensions is defined as follows: Obviously, Eq. (23) is also applicable to the spatial dimensions. A formula of summation-by-parts is obtained from Eq. (23): (24) can be written in a simpler form: Consequently, the constraint Eqs. (8)-(10) are enforced. Furthermore, another type of product rule is used for the time derivative: This formula is used to obtain time derivatives of the momentum of the electromagnetic field (E × B)/4πc and the energy of the electromagnetic field (E 2 + B 2 )/8π from Eqs. (11)- (16). To prove Gauss's law and the solenoidal constraint for the magnetic field, the commutative laws of finite-difference operators are derived as follows: δ δt Corresponding to the proof of energy conservation, some formulae related to the Lorentz factor γ are derived. From the definition of γ;

The law of charge conservation
The law of charge conservation is the 0th-order moment of the relativistic Vlasov equation. In differential form, this is described as follows: f dp x dp y dp z = 0, The corresponding equation in discrete form is derived from Eq. (7) as follows: By substituting Eqs. (17)- (20) into Eq. (33), the following expression can be obtained: Therefore, the law of charge conservation is strictly maintained, even in discrete form. The conservation laws of particle number and mass are also derived similarly.

Gauss's law and solenoidal constraint
Here we review the derivation of Gauss's law in differential form. The divergence of Eq. (3) is Substituting Eq. (32) into Eq. (35); Thus, Gauss's law is maintained if the following condition is satisfied at the initial state: To reproduce these operations in discrete form, Eqs. (11)-(13) are transformed into the following form using the commutative laws of finite-difference operators Eqs. (27) and (28): Therefore, Eq. (35) is automatically maintained by the above discretization. By substituting the law of charge conservation Eq. (34) into Eq. (38), a recurrence formula can be obtained: Hence, Gauss's law is strictly maintained even in discrete form if the law is satisfied at the initial state: Likewise, the solenoidal constraint of the magnetic field is strictly maintained, even in discrete form:

The law of momentum conservation
In this section, only the law of momentum conservation in the x-direction is discussed. Here we review the derivation of momentum conservation in differential form. The momentum of particles is described by the 1st-order moment of the relativistic Vlasov equation: f dp x dp y dp z The Maxwell's equations are transformed into the momentum of the electromagnetic field with the product rule: Finally, the law of momentum conservation is obtained because the terms on right-hand-side of Eqs. (42) and (43) cancel out: To obtain this relationship in discrete form, the following equations are derived by the summation-byparts of Eq. (25) To obtain the momentum of the electromagnetic field in discrete form, the temporal product rule Eq. (26), and the spatial product rule Eq. (23) are applied to Eqs. (12), (13), (15), and (16): where the Gauss's law Eq. (40) and the solenoidal constraint Eq. (41) are used to derive the right-hand-side. Therefore, the total momentum of charged particles and electromagnetic field is strictly conserved, even in discrete form because the terms on the right-hand-side of Eqs. (48) and (49) completely cancel out. Although the proof is omitted, the law of momentum conservation is also derived for the y-direction and z-direction, even in discrete form.

The law of energy conservation
Here we review the derivation of energy conservation in differential form. The energy of particles is described by the 2nd-order moment of the relativistic Vlasov equation: f dp x dp y dp z f dp x dp y dp z f dp x dp y dp z The Ampère-Maxwell and Faraday-Maxwell equations are transformed into the energy of the electromagnetic field by the product rule: Finally, the law of energy conservation is obtained because the terms on the right-hand-side of Eqs. (50) and (51) cancel out: To obtain this relationship in discrete form, the following relationships are derived by the summationby-parts of Eq. (25): Thus, the energy of particles is described using Eqs. (53)-(55) as follows, and it is ensured that the energy is not affected by the magnetic field: To obtain the energy of the electromagnetic field in discrete form, the temporal product rule Eq. where ω is the wave frequency, k is the wavenumber, and ω pe = (4πe 2 n e /m e ) −1/2 is the plasma frequency. The imaginary part of ω corresponds to the growth rate Γ of the instability. Solving Eq. (59) numerically, the most unstable mode and corresponding growth rate are obtained as follows: To calculate the most unstable mode, the length of a periodic domain L is set to be Lω pe /c = 18, and the upper/lower limits of the momentum domain are p/m i c = ±0.01 for protons and p/m e c = ±10 for electrons, respectively. The number of computational cells is 1024 × 1024. The temporal interval is given as c∆t/∆x = 1. The implicit method is implemented with the predictor-corrector method, and the number of iterations is 100 per time-step. A perturbation of wavelength L and amplitude 10 −5 n e is given to the electron density. The initial electric field is set to satisfy Gauss's law, i.e., Eq. (36). Figure 1 shows the time development of the electric field energy. The time is normalized by the plasma frequency (ω pe t). The energy of the electric field is amplified exponentially at the linear growth phase (30 ≤ ω pe t ≤ 70), and the numerical growth rate agrees well with the linear theory Eq. (60). Subsequently, the amplification of the electric field energy saturates and the instability enters the nonlinear regime. Figure  2 indicates the errors of global conservation. As shown in the theoretical proof in Sec. 3, all errors are strictly maintained at the round-off level, even if the instability has entered the nonlinear regime. Thus, the conservation property of the proposed strategy has been demonstrated experimentally. Note that the proposed scheme cannot maintain the conservation of the L1-norm of f due to the central difference. We should use many computational cells to mitigate the contamination of the numerical solutions by numerical dispersion. In Fig. 3, the distribution function of electrons becomes negative at ω pe t ∼ 100; this is a clear evidence that the proposed scheme cannot maintain the conservation of L1-norm. Moreover, the central difference does not include any numerical dissipation. If we employ the upwind difference, for example, the conservation of L1-norm might be maintained even in discrete form. However, the upwind difference will break up the quadratic conservative scheme since the mathematical formulae derived in Sec. 3.1 are no longer applicable. The best strategy to overcome this issue is to extend the proposed scheme to the Vlasov-Fokker-Planck-Maxwell system, and to introduce physical/artificial collision terms.

Relativistic Weibel instability
In previous studies, the non-relativistic Weibel instability [32,33] was calculated as an electromagnetic or Vlasov-Maxwell test problem. Here we show the results of the relativistic Weibel instability calculated by SPUTNIK in the 1D3P mode. The initial distribution is given by the relativistic bi-Maxwellian distribution described as follows: where γ = 1 + p 2 /(mc) 2 , " " denotes the parallel direction, i.e., the x-direction, and "⊥" denotes the perpendicular direction, i.e., the y-and z-directions. Temperature anisotropy is given to electrons as α = 30, and α ⊥ = 5, while for isotropic background protons it is α = α ⊥ = 5. Again, α = mc 2 /k B T is a reciprocal of the temperature normalized by the rest mass energy. The electron distribution function at the initial state is shown in Fig. 4. In this situation, the dispersion relation for the relativistic Weibel instability [34] is given as    where K n is the modified Bessel function of the second kind of order n, ∆ = α /α ⊥ − 1, ξ ⊥ = α ⊥ + Γτ /ck, ξ = α + Γτ /ck, and ζ = ξ 2 + τ 2 . By solving Eq. (62) numerically, the most unstable mode and corresponding growth rate are obtained as follows: kc ω pe 0.92, Γ ω pe 0.144.
To calculate the most unstable mode, the length of a periodic domain L is set to be Lω pe /c = 7, and the upper/lower limits of the momentum domain are p/m i c = ±6 for protons and p/m e c = ±6 for electrons, respectively. The number of computational cells is 128 × 128 × 128 × 128. The temporal interval is given as c∆t/∆x = 1. The implicit method is implemented with the predictor-corrector method, and the number of iterations is 100 per time-step. A perturbation with a wavelength is L is given to B z . Figure 5 shows the time development of the magnetic field energy. The energy of the magnetic field is amplified exponentially at the linear growth phase (20 ≤ ω pe t ≤ 70), and the numerical growth rate agrees well with the linear theory Eq. (63). Subsequently, the amplification of the magnetic field energy saturates and the instability enters the nonlinear regime. Figure 6 indicates the errors of global conservation. As shown in the theoretical proof in Sec. 3, all errors are strictly maintained at the round-off level, even if the instability has entered the nonlinear regime. Thus, the conservation property of the proposed strategy has been demonstrated experimentally.

Conclusions
In this article, we have presented a quadratic conservative scheme for relativistic Vlasov-Maxwell system which is composed of a relativistic Vlasov equation, and Maxwell equations. The scheme is based on the finite-difference method, and thereby enables us to use non-periodic boundary conditions. We introduced  a Crank-Nicolson method for the temporal dimension, and a central difference for spatial and momentum dimensions. Before the discretization of the relativistic Vlasov equation, we have to convert it from a non-conservative formulation to a conservative formulation. The conservation property of the relativistic Vlasov-Maxwell system can be proven with a product rule, and integration-by-parts. Accordingly, we have also introduced some mathematical formulae for product rules in discrete form, and a summation-by-parts. All terms in the proposed scheme have been carefully designed so that every source term that is generated in the momentum and energy equations cancels out even in discrete form. We constructed a kinetic simulation code named SPUTNIK. Experimental demonstration of conservation property was performed via relativistic two-stream instability and relativistic Weibel instability. In both verification tests, SPUTNIK could maintain errors of global conservation of charge, momentum, and energy at a round-off level. It was confirmed that numerical growth rates agree well with linear stability theories of the instabilities. In the verification via the relativistic two-stream instability, an electron distribution function was found to become negative due to a numerical dispersion of the central difference, and the conservation property of L1-norm was clearly violated. By applying a quadratic conservative scheme for the relativistic Vlasov-Fokker-Planck-Maxwell system, this issue could be overcome. Recently, a mass-, momentum-, and energy-conserving scheme [35] for the Rosenbluth-Fokker-Planck equation [36] has been proposed. To construct the conservative Vlasov-Fokker-Planck-Maxwell scheme, the development of a fully conservative scheme for the relativistic Landau-Fokker-Planck equation based on Braams-Karney potential [37] may be required.