Explicit high-order non-canonical symplectic particle-in-cell algorithms for Vlasov-Maxwell systems

Xiao, Jianyuan; Qin, Hong; Liu, Jian; He, Yang; Zhang, Ruili; Sun, Yajuan

doi:10.1063/1.4935904

Explicit high-order non-canonical symplectic particle-in-cell algorithms for classical particle-field systems governed by the Vlasov-Maxwell equations are developed. The algorithms conserve a discrete non-canonical symplectic structure derived from the Lagrangian of the particle-field system, which is naturally discrete in particles. The electromagnetic field is spatially discretized using the method of discrete exterior calculus with high-order interpolating differential forms for a cubic grid. The resulting time-domain Lagrangian assumes a non-canonical symplectic structure. It is also gauge invariant and conserves charge. The system is then solved using a structure-preserving splitting method discovered by He et al. [preprint arXiv:1505.06076 (2015)], which produces five exactly soluble sub-systems, and high-order structure-preserving algorithms follow by combinations. The explicit, high-order, and conservative nature of the algorithms is especially suitable for long-term simulations of particle-field systems with extremely large number of degrees of freedom on massively parallel supercomputers. The algorithms have been tested and verified by the two physics problems, i.e., the nonlinear Landau damping and the electron Bernstein wave.

I. INTRODUCTION

The importance of numerical solutions for the Vlasov-Maxwell (VM) system cannot be overemphasized. In most cases, important and interesting characteristics of the VM system are the long-term behaviors and multi-scale structures, which demand long-term accuracy and fidelity of numerical calculations. Conventional algorithms for the VM systems used in general do not preserve the geometric structures of the physical systems, such as the local energy-momentum conservation law and the symplectic structure. For these algorithms, the truncation errors are small only for each time-step. For example, the truncation error of a fourth order Runge-Kutta method is of the fifth order of the step-size for each time-step. However, numerical errors from different time-steps accumulate coherently with time, and the long-term simulation results are not reliable. To overcome this difficulty, a series of geometric algorithms, which preserve the geometric structures of theoretical models in plasma physics have been developed recently.

At the single particle level, canonical Hamiltonian equation for charged particle dynamics can be integrated using the standard canonical symplectic integrators developed in the late 1980s.^1–9 Since the Hamiltonian expressed in terms of the canonical momentum is not separable, it is believed that symplectic algorithms applicable are in general implicit. Recent studies show that this is not the case, and high order explicit symplectic algorithms for charged particle dynamics have been discovered.^10–12

For the most-studied guiding center dynamics in magnetized plasmas, a non-canonical variational symplectic integrator has been developed and applied.^13–19 It is also recently discovered that the popular Boris algorithm is actually a volume-preserving algorithm.²⁰ This revelation stimulated new research activities.^21,22 For example, high-order volume-preserving methods²³ and relativistic volume-preserving methods²⁴ have been worked out systematically.

For collective dynamics of the particle-field system governed by the Vlasov-Maxwell equations,^25,26 Squire et al.^27,28 constructed the first geometric, structure-preserving algorithm by discretizing a geometric variational principle.²⁵ It has been applied in simulation studies of nonlinear radio-frequency waves in magnetized plasmas.^29,30 Similar methods apply to Vlasov-Poisson system as well.^31–33 We can also discretize directly the Poisson structures of the Vlasov-Maxwell system. A canonical symplectic particle-in-cell (PIC) algorithm is found by discretizing the canonical Poisson bracket,³⁴ and non-canonical symplectic methods are being developed using the powerful Hamiltonian splitting technique^11,35,36 that preserves the non-canonical Morrison-Marsden-Weinstein bracket^37–40 for the VM equations. Of course, geometric structure-preserving algorithms are expected for reduced systems as well. For example, an structure-preserving algorithm has been developed for ideal MHD equations,⁴¹ and applied to study current sheet formation in an ideal plasma without resistivity.⁴² The superiority of these geometric algorithms has been demonstrated. This should not be surprising because geometric algorithms are built on the more fundamental field-theoretical formalism and are directly linked to the perfect form, i.e., the variational principle of physics. The fact that the most elegant form of theory is also the most effective algorithm is philosophically satisfactory.

In this paper, we present an explicit, high-order, non-canonical symplectic PIC algorithm for the Vlasov-Maxwell system. The algorithm conserves a discrete non-canonical symplectic structure derived from the Lagrangian of the particle-field system.^25,26 The Lagrangian is naturally discrete in particles, and the electromagnetic field is discretized using the method of discrete exterior calculus (DEC). An important technique for interpolating differential forms over several grid cells are developed, which generalizes the construction of Whitney forms to higher orders. The resulting Lagrangian is continuous in time and assumes a non-canonical symplectic structure, the dimension of which is finite but large. Because the electromagnetic field is interpolated as differential forms, the time-domain Lagrangian is also gauge invariant and conserves charge. From this Lagrangian, we can readily derive the non-canonical symplectic structure for the dynamics, and the system is solved using a splitting method discovered by He et al.^11,12 The splitting produces five exactly soluble sub-systems, and high-order structure-preserving algorithms follow by combinations. We note that for previous symplectic PIC methods,^27–30,34 high-order algorithms are implicit in general, and only the first order method can be made explicit.³⁴ The explicit and high-order nature of the symplectic algorithms developed in the present study makes it especially suitable for long-term simulations of particle-field systems with extremely large number of degrees of freedom on massively parallel supercomputers.

The paper is organized as follows. The non-canonical symplectic PIC algorithm is derived in Sec. II with an Appendix on Whitney forms and their generalization to high orders. In Sec. III, the developed algorithm is tested and verified by two physics problems, i.e., the nonlinear Landau damping and the electron Bernstein wave.

II. NON-CANONICAL SYMPLECTIC PARTICLE-IN-CELL ALGORITHMS

We start from the Lagrangian of a collection of charged particles and electromagnetic field^25,26

\begin{matrix} L = ∭ d x (\frac{ϵ_{0}}{2} {(- \dot{A} (x) - \nabla ϕ (x))}^{2} - \frac{1}{2 μ_{0}} {(\nabla \times A (x))}^{2} \\ + \sum_{s} δ (x - x_{s}) (\frac{1}{2} m_{s} {\dot{x}}_{s}^{2} + q_{s} A (x) \cdot {\dot{x}}_{s} - q_{s} ϕ (x))), \end{matrix}

(1)

where $A (x)$ and $ϕ (x)$ are the vector and scalar potentials of the electromagnetic field, $x_{s}$ ⁠, m_s, and q_s denote the location, mass, and charge of the s-th particle, and ϵ₀ and μ₀ are the permittivity and permeability in vacuum. We let $ϵ_{0} = μ_{0} = 1$ to simplify the notation.

This Lagrangian is naturally discrete in particles, and we choose to discretize the electromagnetic field in a cubic mesh. To preserve the symplectic structure of the system, the method of DEC⁴³ is used. The DEC theory in cubic meshes can be found in Ref. 44. For field-particle interaction, the interpolation function is used to obtain continuous fields from discrete fields. The spatially discretized Lagrangian L_sd can be written as follows:

\begin{matrix} L_{s d} = \frac{1}{2} (\sum_{J} {(- {\dot{A}}_{J} - \sum_{I} \nabla_{d J I} ϕ_{I})}^{2} - \sum_{K} {(\sum_{J} {curl}_{d K J} A_{J})}^{2}) Δ V + \sum_{s} (\frac{1}{2} m_{s} {\dot{x}}_{s}^{2} + q_{s} ({\dot{x}}_{s} \cdot \sum_{J} W_{σ_{1 J}} (x_{s}) A_{J} - \sum_{I} W_{σ_{0 I}} (x_{s}) ϕ_{I})), \end{matrix}

(2)

where ΔV is the volume of each cell, integers I, J, and K are indices of grid points, and $\nabla_{d}$ and ${curl}_{d}$ are the discrete gradient and curl operators, which are linear operators on the discrete fields $A_{J}$ and $ϕ_{I}$ ⁠. Functions $W_{σ_{0 J}}$ and $W_{σ_{1 I}}$ are interpolation functions for 0-forms (e.g., scalar potential) and 1-forms (e.g., vector potential), respectively. They should be viewed as maps operating on the discrete 0-form $ϕ_{I}$ and 1-form $A_{J}$ to generate continuous forms. More precisely, $W_{σ_{1 J}} (x_{s}) A_{J}$ are the components of the continuous 1-form interpolated from the discrete 1-form $A_{J}$ ⁠. The idea of form interpolation maps is due to Whitney, and the interpolated forms are called Whitney forms. The original Whitney forms⁴⁵ are first order and only for forms in simplicial meshes (e.g., triangle and tetrahedron meshes). In the present study, we have developed high-order interpolation maps for a cubic mesh. The details of the construction of $\nabla_{d}$ ⁠, ${curl}_{d}, W_{σ_{0 I}}, W_{σ_{1 J}}$ ⁠, and the interpolating function for 2-forms (e.g., magnetic fields) $W_{σ_{2 K}}$ are presented in the Appendix. The major new feature of the form interpolation method adopted here is that the interpolation for electric field and magnetic field are different. Even for components in different directions of the same field, the interpolation functions are not the same. This is very different from traditional cubic interpolations used in conventional PIC methods,^46–48 where the same interpolation function is used for all components of electromagnetic fields. The advanced form interpolation method developed in the present study guarantees that the geometric properties of the continuous system are preserved by the discretized system.

The action integral is

S = \int d t L_{s d},

(3)

and the dynamic equations are obtained from Hamilton's principle

\frac{δ S}{δ A_{J}} = 0,

(4)

\frac{δ S}{δ ϕ_{I}} = 0,

(5)

\frac{δ S}{δ x_{s}} = 0 .

(6)

Equations (4) and (5) are Maxwell's equations, and Eq. (6) is Newton's equation with the Lorentz force for the s-th particle. For the dynamics to be gauge independent,¹⁶ it requires that the discrete differential operators and interpolation functions satisfy the following relations:

\nabla \sum_{I} W_{σ_{0 I}} (x) ϕ_{I} = \sum_{I, J} W_{σ_{1 J}} (x) \nabla_{d J I} ϕ_{I},

(7)

\nabla \times \sum_{J} W_{σ_{1 J}} (x) A_{J} = \sum_{J, K} W_{σ_{2 K}} (x) {curl}_{d K J} A_{J} .

(8)

The gauge independence of this spatially discretized system implies that the dynamics conserves charge automatically.

Since the dynamics are gauge independent, we can choose any gauge that is convenient. For simplicity, the temporal gauge, i.e., $ϕ_{I} = 0$ ⁠, is adopted in the present study. To obtain the non-canonical symplectic structure and Poisson bracket, we look at the Lagrangian 1-form γ for the spatially discretized system defined by $S = \int γ$ ⁠. Let $q = [A_{J}, x_{s}]$ ⁠, and the Lagrangian 1-form can be written as

γ = \frac{\partial L_{s d}}{\partial \dot{q}} d q - H d t,

(9)

where $d$ denotes the exterior derivative. In Eq. (9),

\frac{\partial L_{s d}}{\partial \dot{q}} = [{\dot{A}}_{J} Δ V, m_{s} {\dot{x}}_{s} + q_{s} \sum_{J} W_{σ_{1 J}} (x_{s}) A_{J}],

(10)

and

H = \frac{\partial L_{s d}}{\partial \dot{q}} {\dot{q}}^{T} - L_{s d},

(11)

= \frac{1}{2} Δ V (\sum_{J} {\dot{A}}_{J}^{2} + \sum_{K} {(\sum_{J} {curl}_{d K J} A_{J})}^{2}) + \sum_{s} \frac{1}{2} m_{s} {\dot{x}}_{s}^{2},

(12)

is the Hamiltonian. The dynamical equation of the system can be written as^25,49

i_{[\dot{q}, \ddot{q}, 1]} d γ = 0,

(13)

where $[\dot{q}, \ddot{q}, 1]$ represent vector field $\dot{q} \frac{\partial}{\partial q} + \ddot{q} \frac{\partial}{\partial \dot{q}} + \frac{\partial}{\partial t} .$ The non-canonical symplectic structure is

Ω = d (\frac{\partial L_{s d}}{\partial \dot{q}} d q),

(14)

and the dynamical equation (13) is equivalent to

\frac{d}{d t} [\begin{matrix} A_{J} \\ x_{s} \\ {\dot{A}}_{J} \\ {\dot{x}}_{s} \end{matrix}] = Ω^{- 1} [\begin{matrix} \frac{\partial}{\partial A_{J}} \\ \frac{\partial}{\partial x_{s}} \\ \frac{\partial}{\partial {\dot{A}}_{J}} \\ \frac{\partial}{\partial {\dot{x}}_{s}} \end{matrix}] H .

(15)

The corresponding non-canonical Poisson bracket is

{F, G} = [\frac{\partial F}{\partial A_{J}}, \frac{\partial F}{\partial x_{s}}, \frac{\partial F}{\partial {\dot{A}}_{J}}, \frac{\partial F}{\partial {\dot{x}}_{s}}] Ω^{- 1} {[\frac{\partial G}{\partial A_{J}}, \frac{\partial G}{\partial x_{s}}, \frac{\partial G}{\partial {\dot{A}}_{J}}, \frac{\partial G}{\partial {\dot{x}}_{s}}]}^{T},

(16)

or more specifically

\begin{matrix} {F, G} = \frac{1}{Δ V} \sum_{J} (\frac{\partial F}{\partial A_{J}} \cdot \frac{\partial G}{\partial {\dot{A}}_{J}} - \frac{\partial F}{\partial {\dot{A}}_{J}} \cdot \frac{\partial G}{\partial A_{J}}) + \sum_{s} \frac{1}{m_{s}} (\frac{\partial F}{\partial x_{s}} \cdot \frac{\partial G}{\partial {\dot{x}}_{s}} - \frac{\partial F}{\partial {\dot{x}}_{s}} \cdot \frac{\partial G}{\partial x_{s}}) + \sum_{s} \frac{q_{s}}{m_{s} Δ V} (\frac{\partial G}{\partial {\dot{x}}_{s}} \cdot \sum_{J} W_{σ_{1 J}} (x_{s}) \frac{\partial F}{\partial {\dot{A}}_{J}} - \frac{\partial F}{\partial {\dot{x}}_{s}} \cdot \sum_{J} W_{σ_{1 J}} (x_{s}) \frac{\partial G}{\partial {\dot{A}}_{J}}) \\ - \sum_{s} \sum_{J} \frac{q_{s}}{m_{s}^{2}} \frac{\partial F}{\partial {\dot{x}}_{s}} \cdot [\nabla \times W_{σ_{1 J}} (x_{s}) A_{J}] \times \frac{\partial G}{\partial {\dot{x}}_{s}} . \end{matrix}

(17)

Now, we introduce two new variables $E_{J}$ and $B_{K}$ ⁠, which are the discrete electric field and magnetic field

E_{J} = - {\dot{A}}_{J},

(18)

B_{K} = \sum_{J} {curl}_{d K J} A_{J} .

(19)

In terms of $E_{J}$ and $B_{K}$ ⁠, the partial derivatives with respect to $A_{J}$ and ${\dot{A}}_{J}$ are

\frac{\partial F}{\partial A_{J}} = \sum_{K} \frac{\partial F}{\partial B_{K}} {curl}_{d K J},

(20)

\frac{\partial F}{\partial {\dot{A}}_{J}} = - \frac{\partial F}{\partial E_{J}} .

(21)

Note that ${curl}_{d K J}$ in Eq. (20) is a matrix. Using Eq. (8), the Poisson bracket in terms of $E_{J}$ and $B_{K}$ can be written as

\begin{matrix} {F, G} = \frac{1}{Δ V} \sum_{J} (\frac{\partial F}{\partial E_{J}} \cdot \sum_{K} \frac{\partial G}{\partial B_{K}} {curl}_{d K J} - \sum_{K} \frac{\partial F}{\partial B_{K}} {curl}_{d K J} \cdot \frac{\partial G}{\partial E_{J}}) \\ + \sum_{s} \frac{1}{m_{s}} (\frac{\partial F}{\partial x_{s}} \cdot \frac{\partial G}{\partial {\dot{x}}_{s}} - \frac{\partial F}{\partial {\dot{x}}_{s}} \cdot \frac{\partial G}{\partial x_{s}}) \\ + \sum_{s} \frac{q_{s}}{m_{s} Δ V} (\frac{\partial F}{\partial {\dot{x}}_{s}} \cdot \sum_{J} W_{σ_{1 J}} (x_{s}) \frac{\partial G}{\partial E_{J}} - \frac{\partial G}{\partial {\dot{x}}_{s}} \cdot \sum_{J} W_{σ_{1 J}} (x_{s}) \frac{\partial F}{\partial E_{J}}) + \\ - \sum_{s} \frac{q_{s}}{m_{s}^{2}} \frac{\partial F}{\partial {\dot{x}}_{s}} \cdot [\sum_{K} W_{σ_{2 K}} (x_{s}) B_{K}] \times \frac{\partial G}{\partial {\dot{x}}_{s}}, \end{matrix}

(22)

and the Hamiltonian is

H = \frac{1}{2} (Δ V \sum_{J} E_{J}^{2} + Δ V \sum_{K} B_{K}^{2} + \sum_{s} m_{s} {\dot{x}}_{s}^{2}) .

(23)

The evolution equations is

\dot{F} = {F, H},

(24)

where

F = [E_{J}, B_{K}, x_{s}, {\dot{x}}_{s}] .

(25)

This is a Hamiltonian system with a non-canonical symplectic structure or Poisson bracket. In general, symplectic integrators for non-canonical systems are difficult to construct. However, using the splitting method discovered by He et al.,^11,12 we have found explicit high-order symplectic algorithms for this Hamiltonian system that preserve its non-canonical symplectic structure. We note that splitting method had been applied to the Vlasov equation previously without the context of symplectic structure.⁵⁰ We split the Hamiltonian in Eq. (23) into five parts

H = H_{E} + H_{B} + H_{x} + H_{y} + H_{z},

(26)

H_{E} = \frac{1}{2} Δ V \sum_{J} E_{J}^{2},

(27)

H_{B} = \frac{1}{2} Δ V \sum_{K} B_{K}^{2},

(28)

H_{r} = \frac{1}{2} \sum_{s} m_{s} {\dot{r}}_{s}^{2}, for r in x, y, z .

(29)

It turns out that the sub-system generated by each part can be solved exactly, and high-order symplectic algorithms follow by combination. The evolution equation for H_E is $\dot{F} = {F, H_{E}}$ ⁠, which can be written as

{\dot{E}}_{J} = 0,

(30)

{\dot{B}}_{K} = - \sum_{J} {curl}_{d K J} E_{J},

(31)

{\dot{x}}_{s} = 0,

(32)

{\ddot{x}}_{s} = \frac{q_{s}}{m_{s}} \sum_{J} W_{σ_{1 J}} (x_{s}) E_{J} .

(33)

The exact solution $Θ_{E} (Δ t)$ for any time step $Δ t$ is

E_{J} (t + Δ t) = E_{J} (t),

(34)

B_{K} (t + Δ t) = B_{K} (t) - Δ t \sum_{J} {curl}_{d K J} E_{J} (t),

(35)

x_{s} (t + Δ t) = x_{s} (t),

(36)

{\dot{x}}_{s} (t + Δ t) = {\dot{x}}_{s} (t) + \frac{q_{s}}{m_{s}} Δ t \sum_{J} W_{σ_{1 J}} (x_{s} (t)) E_{J} (t) .

(37)

The evolution equation for H_B is $\dot{F} = {F, H_{B}}$ ⁠, or

{\dot{E}}_{J} = \sum_{K} {curl}_{d K J} B_{K},

(38)

{\dot{B}}_{K} = 0,

(39)

{\dot{x}}_{s} = 0,

(40)

{\ddot{x}}_{s} = 0,

(41)

whose exact solution $Θ_{B} (Δ t)$ is

E_{J} (t + Δ t) = E_{J} (t) + Δ t \sum_{K} {curl}_{d K J} B_{K} (t),

(42)

B_{K} (t + Δ t) = B_{K} (t),

(43)

x_{s} (t + Δ t) = x_{s} (t),

(44)

{\dot{x}}_{s} (t + Δ t) = {\dot{x}}_{s} (t) .

(45)

The evolution equation for H_x is $\dot{F} = {F, H_{x}}$ ⁠, or

{\dot{E}}_{J} = - \sum_{s} \frac{q_{s}}{Δ V} {\dot{x}}_{s} e_{x} W_{σ_{1 J}} (x_{s}),

(46)

{\dot{B}}_{K} = 0,

(47)

{\dot{x}}_{s} = {\dot{x}}_{s} e_{x},

(48)

{\ddot{x}}_{s} = {\dot{x}}_{s} e_{x} \times \sum_{K} W_{σ_{2 K}} (x_{s}) B_{K} .

(49)

The exact solution $Θ_{x} (Δ t)$ of this sub-system can also be computed as

E_{J} (t + Δ t) = E_{J} (t) - \int_{0}^{Δ t} d t' \sum_{s} \frac{q_{s}}{Δ V} {\dot{x}}_{s} (t) e_{x} W_{σ_{1 J}} (x_{s} (t) + {\dot{x}}_{s} (t) t' e_{x}),

(50)

B_{K} (t + Δ t) = B_{K} (t),

(51)

x_{s} (t + Δ t) = x_{s} (t) + Δ t {\dot{x}}_{s} (t) e_{x},

(52)

{\dot{x}}_{s} (t + Δ t) = {\dot{x}}_{s} (t) + {\dot{x}}_{s} (t) e_{x} \times \int_{0}^{Δ t} d t' \sum_{K} W_{σ_{2 K}} (x_{s} (t) + {\dot{x}}_{s} (t) t' e_{x}) B_{K} (t) .

(53)

Exact solutions $Θ_{y} (Δ t)$ and $Θ_{z} (Δ t)$ for sub-systems corresponding to H_y and H_z are obtained in a similar manner. These exact solutions for sub-systems are then combined to construct symplectic integrators for the original non-canonical Hamiltonian system specified by Eqs. (22) and (23). For example, a first order scheme can be constructed as

Θ_{1} (Δ t) = Θ_{E} (Δ t) Θ_{B} (Δ t) Θ_{x} (Δ t) Θ_{y} (Δ t) Θ_{z} (Δ t),

(54)

and a second order symmetric scheme is

Θ_{2} (Δ t) = Θ_{x} (Δ t / 2) Θ_{y} (Δ t / 2) Θ_{z} (Δ t / 2) Θ_{B} (Δ t / 2) Θ_{E} (Δ t)

Θ_{B} (Δ t / 2) Θ_{z} (Δ t / 2) Θ_{y} (Δ t / 2) Θ_{x} (Δ t / 2) .

(55)

An algorithm with order $2 (l + 1)$ can be constructed in the following way:

Θ_{2 (l + 1)} (Δ t) = Θ_{2 l} (α_{l} Δ t) Θ_{2 l} (β_{l} Δ t) Θ_{2 l} (α_{l} Δ t),

(56)

α_{l} = 1 / (2 - 2^{1 / (2 l + 1)}),

(57)

β_{l} = 1 - 2 α_{l} .

(58)

III. NUMERICAL EXAMPLES

We have implemented the second-order non-canonical symplectic PIC algorithm described above using the C programming language. The code is parallelized using MPI and OpenMP. To test the algorithm, two physics problems are simulated. The first problem is the nonlinear Landau damping of an electrostatic wave in a hot plasma, which has been investigated theoretically^51–53 and numerically.^17,54,55 The density of electron n_e is $1.2116 \times 10^{16} m^{- 3}$ ⁠, and the electron velocity is Maxwellian distributed with thermal speed $v_{T} = 0.1 c$ ⁠, where c is the speed of light in vacuum. The computation is carried out in a $672 \times 1 \times 1$ cubic mesh, and the size of each grid cell is $Δ l = 2.4355 \times 10^{- 4} m$ ⁠. There is no external electromagnetic field, and there are 40 000 sample particles in each cell when unperturbed. The initial electric field is $E_{1} = E_{1} cos (k x) e_{x}$ ⁠, where $k = 2 π / 224 Δ l$ is the wave number, and the amplitude is $E_{1} = 36$ kV/m. The simulation is carried out for 15 000 time-steps, and the electric field is recorded during the simulation. We plot the evolution of electric field to observe the Landau damping phenomenon (Figs. 1 and 2). The theoretical damping rate of electric field is $ω_{i} = - 1.3223 \times 10^{9} rad / s$ ⁠, and it is evident that the simulation result agrees well with the theory.

FIG. 1.

View large Download slide

The time evolution of an electrostatic wave in a hot plasma.

FIG. 2.

View large Download slide

Logarithmic plot of the time evolution of the absolute value of the electric field. The slope of the solid green line is the theoretical damping rate.

Another test problem is the dispersion relation of electron Bernstein waves.⁵⁶ In this problem, an electromagnetic wave propagates perpendicularly to a uniform external magnetic field $B_{0} = B_{0} e_{z}$ with $B_{0} = 5.13$ T. Other system parameters are

n_{e} = 2.4 \times 10^{20} m^{- 3},

(59)

v_{T} = 0.07 c,

(60)

Δ l = 2.5 \times 10^{- 5} m,

(61)

Δ t = \frac{Δ l}{2 c} .

(62)

The computation domain is a $768 \times 1 \times 1$ cubic mesh, and the averaged number of sample points per grid is 4000. An initial electromagnetic perturbation is imposed, and after simulating 6000 time steps, the space-time spectrum of E_x is plotted in Fig. 3, which shows that the dispersion relation simulated matches the theoretical curve perfectly.

FIG. 3.

View large Download slide

The dispersion relation of electron Bernstein wave (contours) obtained by the non-canonical symplectic PIC method. The red dots are calculated from the analytical dispersion relation.

As a symplectic method, the non-canonical symplectic PIC algorithm is expected to have good long-term properties. To demonstrate that, we run another test, where the system parameters are same as the previous problem for the Bernstein wave, except that the number of sample particles is 40 per grid and the simulation domain is a $48 \times 1 \times 1$ cubic mesh. Both the second and the first order split methods are tested. The simulation is run for 2.5 × 10⁶ time-steps, and the evolution of total energy errors is plotted in Fig. 4. It is clear that the energy errors are bounded within a small value during the entire long-term simulation for both split methods, and the second order method is more accurate than the first order method.

FIG. 4.

View large Download slide

Total energy error as a function of time for a magnetized hot plasma obtained by the second (a) and the first (b) order non-canonical symplectic PIC method. Both energy errors are bounded within a small value.

IV. SUMMARY AND DISCUSSION

We have developed and tested a non-canonical symplectic PIC algorithm for the VM system. The non-canonical symplectic structure is obtained by discretizing the electromagnetic field of the particle-field Lagrangian using the method of discrete exterior calculus. A high-order interpolating method for differential forms is developed to render smooth interpolations of the electromagnetic field. The effectiveness and conservative nature of the algorithm has been verified by the physics problems of nonlinear Landau damping and electron Bernstein wave.

ACKNOWLEDGMENTS

This research was supported by ITER-China Program (2015GB111003, 2014GB124005, and 2013GB111000), JSPS-NRF-NSFC A3 Foresight Program in the field of Plasma Physics (NSFC-11261140328), the National Science Foundation of China (11575186, 11575185, 11505185, and 11505186), the CAS Program for Interdisciplinary Collaboration Team, the Geo-Algorithmic Plasma Simulator (GAPS) Project, and the U.S. Department of Energy (DE-AC02-09CH11466).

APPENDIX: HIGH ORDER INTERPOLATION FORMS FOR A CUBIC MESH

The interpolation forms for a cubic mesh is inspired by the Whitney forms,⁴⁵ which were originally developed as the first order interpolation forms over simplicial complex and has become an important tool in DEC.^43,45,57 In DEC theory, the discrete forms are defined on chains.⁴³ For example, the discrete 0-forms are defined on vertexes of the grid, the discrete 1-forms are defined on edges, and the discrete differential operators such as $\nabla_{d}$ and ${curl}_{d}$ are discrete exterior derivatives $d_{d}$ acted on these discrete forms. The Whitney map $ϕ_{W}$ is a map that allows us to define continuous forms based on these discrete forms.⁴³ With this map, the following relation holds for any discrete form α:

ϕ_{W} d_{d} α = d ϕ_{W} α,

(A1)

where $d$ is the continuous exterior derivative.

DEC solvers for Maxwell equations in cubic meshes are given by Stern et al.⁴⁴ For our purpose, we need to construct appropriate discrete differential operators $\nabla_{d}$ ⁠, ${curl}_{d}, {div}_{d}$ ⁠, as well as interpolation functions $W_{σ_{0 I}} (x)$ ⁠, $W_{σ_{1 J}} (x), W_{σ_{2 K}} (x)$ ⁠, and $W_{σ_{3 L}} (x)$ in a cubic mesh such that

\nabla \sum_{I} W_{σ_{0 I}} (x) ϕ_{I} = \sum_{I, J} W_{σ_{1 J}} (x) \nabla_{d J I} ϕ_{I},

(A2)

\nabla \times \sum_{J} W_{σ_{1 J}} (x) A_{J} = \sum_{J, K} W_{σ_{2 K}} (x) {curl}_{d K J} A_{J},

(A3)

\nabla \cdot \sum_{K} W_{σ_{2 K}} (x) B_{K} = \sum_{K, L} W_{σ_{3 L}} (x) {div}_{d L K} B_{K},

(A4)

hold for any $ϕ_{I}, A_{J}$ ⁠, and $B_{K}$ ⁠.

To accomplish this goal, we start from choosing an interpolation function $W_{σ_{0 I}} (x) = W_{σ_{0}} (x - x_{I})$ for 0-forms (e.g., scalar potential) as follows:

W_{σ_{0}} (x) = W_{1} (x) W_{1} (y) W_{1} (z),

(A5)

where $x_{I}$ is the coordinate of the I-th grid vertex and the cell size is chosen to be 1 for simplicity. It is required that $W_{1} (x)$ satisfies the following conditions:

W_{1} (x) = 0, if | x | > = 1,

(A6)

W_{1} (x) + W_{1} (x - 1) = 1, if 0 \leq x < 1 .

(A7)

For example, W₁ can be chosen to be piece-wise linear over one grid cell, i.e.,

W_{1} (x) = {\begin{matrix} 0, & | x | \geq 1, \\ 1 - | x |, & | x | < 1 . \end{matrix}

(A8)

For x, y, and z in $[i', i' + 1], [j', j' + 1], [k', k' + 1]$ ⁠, the x component of the left hand side of Eq. (A2)T_x is

\begin{matrix} T_{x} = \sum_{i \in {i', i' + 1}, j \in {j', j' + 1}, k \in {k', k' + 1}} ϕ_{i, j, k} W_{1}^{'} (x - i) W_{1} (y - j) W_{1} (z - k), \\ = \sum_{j \in {j', j' + 1}, k \in {k', k' + 1}} (ϕ_{i', j, k} W_{1}^{'} (x - i') + ϕ_{i' + 1, j, k} W_{1}^{'} (x - i' - 1)) W_{1} (y - j) W_{1} (z - k), \\ = \sum_{j \in {j', j' + 1}, k \in {k', k' + 1}} (ϕ_{i', j, k} - ϕ_{i' + 1, j, k}) W_{1}^{'} (x - i') W_{1} (y - j) W_{1} (z - k) . \end{matrix}

(A9)

The y and z components can be also deduced in the similar way

T_{y} = \sum_{i \in {i', i' + 1}, k \in {k', k' + 1}} (ϕ_{i, j', k} - ϕ_{i, j' + 1, k}) W_{1} (x - i) W_{1}^{'} (y - j') W_{1} (z - k),

(A10)

T_{z} = \sum_{i \in {i', i' + 1}, j \in {j', j' + 1}} (ϕ_{i, j, k'} - ϕ_{i, j, k' + 1}) W_{1} (x - i) W_{1} (y - j) W_{1}^{'} (z - k') .

(A11)

Equations (A9)–(A11) indicate that $W_{1}^{'} (x - i) W_{1} (y - j) W_{1} (z - k) d x$ ⁠, $W_{1} (x - i) W_{1}^{'} (y - j) W_{1} (z - k) d y$ ⁠, and $W_{1} (x - i) W_{1} (y - j) W_{1}^{'} (z - k) d z$ can be viewed as the bases for 1-form interpolation map, and that the discrete gradient operator $\nabla_{d}$ can be defined as linear operator on $ϕ_{I}$ as

{(\nabla_{d} ϕ)}_{i, j, k} = [ϕ_{i + 1, j, k} - ϕ_{i, j, k}, ϕ_{i, j + 1, k} - ϕ_{i, j, k}, ϕ_{i, j, k + 1} - ϕ_{i, j, k}] .

(A12)

For a given discrete 1-form field $A_{I}$ ⁠, the interpolated 1-form field is

\begin{matrix} \sum_{i, j, k} [A_{x i, j, k} W_{1}^{(1)} (x - i) W_{1} (y - j) W_{1} (z - k) d x \\ + A_{y i, j, k} W_{1} (x - i) W_{1}^{(1)} (y - j) W_{1} (z - k) d y \\ + A_{z i, j, k} W_{1} (x - i) W_{1} (y - j) W_{1}^{(1)} (z - k) d z], \end{matrix}

(A13)

where

W_{1}^{(1)} (x) = {\begin{array}{l} - W_{1}^{'} (x), & 0 \leq x < 1 \\ 0, & otherwise \end{array}

(A14)

is the one-cell interpolation function. The components of this interpolated 1-form field are written as

\sum_{i, j, k} W_{σ_{1, i, j, k}} (x) A_{i, j, k} \equiv \sum_{i, j, k} {[\begin{matrix} A_{x i, j, k} W_{1}^{(1)} (x - i) W_{1} (y - j) W_{1} (z - k) \\ A_{y i, j, k} W_{1} (x - i) W_{1}^{(1)} (y - j) W_{1} (z - k) \\ A_{z i, j, k} W_{1} (x - i) W_{1} (y - j) W_{1}^{(1)} (z - k) \end{matrix}]}^{T} .

(A15)

By the same procedure, we find that the discrete differential operators ${curl}_{d}$ and ${div}_{d}$ should be defined as

{({curl}_{d} A)}_{i, j, k} = {[\begin{matrix} (A_{z i, j + 1, k} - A_{z i, j, k}) - (A_{y i, j, k + 1} - A_{y i, j, k}) \\ (A_{x i, j, k + 1} - A_{x i, j, k}) - (A_{z i + 1, j, k} - A_{z i, j, k}) \\ (A_{y i + 1, j, k} - A_{y i, j, k}) - (A_{x i, j + 1, k} - A_{x i, j, k}) \end{matrix}]}^{T},

(A16)

{({div}_{d} B)}_{i, j, k} = [(B_{x i + 1, j, k} - B_{x i, j, k}) + (B_{y i, j + 1, k} - B_{y i, j, k}) + (B_{z i, j, k + 1} - B_{z i, j, k})] .

(A17)

For a given discrete 2-form field $B_{I}$ ⁠, the interpolated 2-form field is

\begin{matrix} \sum_{i, j, k} [B_{x i, j, k} W_{1} (x - i) W_{1}^{(1)} (y - j) W_{1}^{(1)} (z - k) d y \land d z \\ + B_{y i, j, k} W_{1}^{(1)} (x - i) W_{1} (y - j) W_{1}^{(1)} (z - k) d z \land d x \\ + B_{z i, j, k} W_{1}^{(1)} (x - i) W_{1}^{(1)} (y - j) W_{1} (z - k) d x \land d y], \end{matrix}

(A18)

whose components can be written as

\sum_{i, j, k} W_{σ_{2, i, j, k}} (x) B_{i, j, k} \equiv \sum_{i, j, k} {[\begin{matrix} B_{x i, j, k} W_{1} (x - i) W_{1}^{(1)} (y - j) W_{1}^{(1)} (z - k) \\ B_{y i, j, k} W_{1}^{(1)} (x - i) W_{1} (y - j) W_{1}^{(1)} (z - k) \\ B_{z i, j, k} W_{1}^{(1)} (x - i) W_{1}^{(1)} (y - j) W_{1} (z - k) \end{matrix}]}^{T} .

(A19)

For a discrete 3-form field ρ_I, the interpolated 3-form field is

\sum_{i, j, k} ρ_{i, j, k} W_{1}^{(1)} (x - i) W_{1}^{(1)} (y - j) W_{1}^{(1)} (z - k) d x \land d y \land d z,

(A20)

and the corresponding scalar is denoted as

\sum_{i, j, k} W_{σ_{3, i, j, k}} (x) ρ_{i, j, k} \equiv \sum_{i, j, k} ρ_{i, j, k} W_{1}^{(1)} (x - i) W_{1}^{(1)} (y - j) W_{1}^{(1)} (z - k) .

(A21)

We can verify that

{({curl}_{d} \nabla_{d} ϕ)}_{i, j, k} = 0, for any i, j, k and ϕ,

(A22)

{({div}_{d} {curl}_{d} A)}_{i, j, k} = 0, for any i, j, k and A .

(A23)

In Ref. 44, the discrete electromagnetic fields in cubic mesh seem different from ours on first look. But the difference is merely in the notation for indices. For example, for the discrete 1-form, we can alternatively use half-integer indices to rewrite Eq. (A12) as

\begin{matrix} [{(\nabla_{d} ϕ)}_{x i + 1 / 2, j, k}, {(\nabla_{d} ϕ)}_{y i, j + 1 / 2, k}, {(\nabla_{d} ϕ)}_{z i, j, k + 1 / 2}] \\ = [ϕ_{i + 1, j, k} - ϕ_{i, j, k}, ϕ_{i, j + 1, k} - ϕ_{i, j, k}, ϕ_{i, j, k + 1} - ϕ_{i, j, k}], \end{matrix}

(A24)

which is then identical with the notation in Ref. 44.

The above interpolation forms are defined over one grid cell. For the simulations reported here, in order to achieve higher accuracy, we have developed and deployed high-order interpolation forms over two grid cells. The interpolation 0-forms are

W_{σ_{0}} (x) = W_{1} (x) W_{1} (y) W_{1} (z),

(A25)

where $W_{1} (x)$ satisfies

W_{1} (x) = 0, if | x | > = 2,

(A26)

W_{1} (x + 1) + W_{1} (x) + W_{1} (x - 1) + W_{1} (x - 2) = 1, if 0 \leq x < 1 .

(A27)

The W₁ adopted in the algorithm is

W_{1} (x) = {\begin{array}{l} 0, & x > 2, \\ \frac{15}{1024} x^{8} - \frac{15}{128} x^{7} + \frac{49}{128} x^{6} - \frac{21}{35} x^{5} + \frac{35}{64} x^{4} - x + 1, & 1 < x \leq 2, \\ - \frac{15}{1024} x^{8} - \frac{15}{128} x^{7} + \frac{7}{16} x^{6} - \frac{21}{35} x^{5} + \frac{175}{256} x^{4} - \frac{105}{128} x^{2} + \frac{337}{512}, & 0 < x \leq 1, \\ - \frac{15}{1024} x^{8} + \frac{15}{128} x^{7} + \frac{7}{16} x^{6} + \frac{21}{35} x^{5} + \frac{175}{256} x^{4} - \frac{105}{128} x^{2} + \frac{337}{512}, & - 1 < x \leq 0, \\ \frac{15}{1024} x^{8} + \frac{15}{128} x^{7} + \frac{49}{128} x^{6} + \frac{21}{35} x^{5} + \frac{35}{64} x^{4} - x + 1, & - 2 < x \leq - 1, \\ 0, & x \leq - 2 . \end{array}

(A28)

It can be proved that this piece-wise polynomial function is 3rd order continuous in the whole space. For this two-cell interpolation scheme, the $\nabla_{d}$ defined in Eq. (A12), ${curl}_{d}$ defined in Eq. (A16), and ${div}_{d}$ defined in Eq. (A17) remain the same, but the function $W_{1}^{(1)}$ in Eqs. (A13), (A15), and (A18)–(A21) needs to be replaced by the $W_{1}^{(2)}$ function defined as

W_{1}^{(2)} (x) = - {\begin{array}{l} W_{1}^{'} (x) + W_{1}^{'} (x + 1) + W_{1}^{'} (x + 2), & - 1 \leq x < 2, \\ 0, & otherwise . \end{array}

(A29)

It can be proved that Eqs. (A2)–(A4), (A22), and (A23) hold for this two-cell interpolation scheme.

References

1.

R. D.

Ruth

,

IEEE Trans. Nucl. Sci.

30

,

2669

(

1983

).

https://doi.org/10.1109/TNS.1983.4332919

Google Scholar

Crossref

2.

K.

Feng

, in

Proceedings of the 1984 Beijing Symposium on Differential Geometry and Differential Equations

, edited by

K.

Feng

(

Science Press

,

1985

), pp.

42

–

58

.

Google Scholar

3.

K.

Feng

,

J. Comput. Math.

4

,

279

(

1986

).

Google Scholar

4.

K.

Feng

and

M.

Qin

,

Symplectic Geometric Algorithms for Hamiltonian Systems

(

Springer-Verlag

,

2010

).

Google Scholar

Crossref

5.

E.

Forest

and

R. D.

Ruth

,

Physica D

43

,

105

(

1990

).

https://doi.org/10.1016/0167-2789(90)90019-L

Google Scholar

Crossref

6.

P. J.

Channell

and

C.

Scovel

,

Nonlinearity

3

,

231

(

1990

).

https://doi.org/10.1088/0951-7715/3/2/001

Google Scholar

Crossref

7.

J.

Candy

and

W.

Rozmus

,

J. Comput. Phys.

92

,

230

(

1991

).

https://doi.org/10.1016/0021-9991(91)90299-Z

Google Scholar

Crossref

8.

J. E.

Marsden

and

M.

West

,

Acta Numer.

10

,

357

(

2001

).

https://doi.org/10.1017/S096249290100006X

Google Scholar

Crossref

9.

E.

Hairer

,

C.

Lubich

, and

G.

Wanner

,

Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations

(

Springer

,

New York

,

2002

).

Google Scholar

Crossref

10.

S. A.

Chin

,

Phys. Rev. E

77

,

066401

(

2008

).

https://doi.org/10.1103/PhysRevE.77.066401

Google Scholar

Crossref

11.

Y.

He

,

H.

Qin

,

Y.

Sun

,

J.

Xiao

,

R.

Zhang

, and

J.

Liu

, preprint arXiv:1505.06076 (

2015

).

12.

Y.

He

,

Y.

Sun

,

Z.

Zhou

,

J.

Liu

, and

H.

Qin

, preprint arXiv:1509.07794 (

2015

).

13.

H.

Qin

and

X.

Guan

,

Phys. Rev. Lett.

100

,

035006

(

2008

).

https://doi.org/10.1103/PhysRevLett.100.035006

Google Scholar

Crossref

PubMed

14.

H.

Qin

,

X.

Guan

, and

W. M.

Tang

,

Phys. Plasmas (1994-present)

16

,

042510

(

2009

).

https://doi.org/10.1063/1.3099055

Google Scholar

Crossref

15.

J.

Li

,

H.

Qin

,

Z.

Pu

,

L.

Xie

, and

S.

Fu

,

Phys. Plasmas (1994-present)

18

,

052902

(

2011

).

https://doi.org/10.1063/1.3589275

Google Scholar

Crossref

16.

J.

Squire

,

H.

Qin

, and

W. M.

Tang

,

Phys. Plasmas (1994-present)

19

,

052501

(

2012

).

https://doi.org/10.1063/1.4742985

Google Scholar

Crossref

17.

M.

Kraus

, preprint arXiv:1307.5665 (

2013

).

18.

R.

Zhang

,

J.

Liu

,

Y.

Tang

,

H.

Qin

,

J.

Xiao

, and

B.

Zhu

,

Phys. Plasmas (1994-present)

21

,

032504

(

2014

).

https://doi.org/10.1063/1.4867669

Google Scholar

Crossref

19.

C. L.

Ellison

,

J.

Finn

,

H.

Qin

, and

W. M.

Tang

,

Plasma Phys. Controlled Fusion

57

,

054007

(

2015

).

https://doi.org/10.1088/0741-3335/57/5/054007

Google Scholar

Crossref

20.

H.

Qin

,

S.

Zhang

,

J.

Xiao

,

J.

Liu

,

Y.

Sun

, and

W. M.

Tang

,

Phys. Plasmas (1994-present)

20

,

084503

(

2013

).

https://doi.org/10.1063/1.4818428

Google Scholar

Crossref

21.

S.

Zhang

,

Y.

Jia

, and

Q.

Sun

,

J. Comput. Phys.

282

,

43

(

2014

).

https://doi.org/10.1016/j.jcp.2014.10.062

Google Scholar

Crossref

22.

C.

Ellison

,

J.

Burby

, and

H.

Qin

,

J. Comput. Phys.

301

,

489

(

2015

).

https://doi.org/10.1016/j.jcp.2015.09.007

Google Scholar

Crossref

23.

Y.

He

,

Y.

Sun

,

J.

Liu

, and

H.

Qin

,

J. Comput. Phys.

281

,

135

(

2015

).

https://doi.org/10.1016/j.jcp.2014.10.032

Google Scholar

Crossref

24.

R.

Zhang

,

J.

Liu

,

H.

Qin

,

Y.

Wang

,

Y.

He

, and

Y.

Sun

,

Phys. Plasmas (1994-present)

22

,

044501

(

2015

).

https://doi.org/10.1063/1.4916570

Google Scholar

Crossref

25.

H.

Qin

,

R.

Cohen

,

W.

Nevins

, and

X.

Xu

,

Phys. Plasmas (1994-present)

14

,

056110

(

2007

).

https://doi.org/10.1063/1.2472596

Google Scholar

Crossref

26.

H.

Qin

,

J. W.

Burby

, and

R. C.

Davidson

,

Phys. Rev. E

90

,

043102

(

2014

).

https://doi.org/10.1103/PhysRevE.90.043102

Google Scholar

Crossref

27.

J.

Squire

,

H.

Qin

, and

W. M.

Tang

, “

Geometric integration of the Vlasov-Maxwell system with a variational particle-in-cell scheme

,”

Technical Report No. PPPL-4748

, Princeton Plasma Physics Laboratory,

2012

.

Google Scholar

28.

J.

Squire

,

H.

Qin

, and

W. M.

Tang

,

Phys. Plasmas (1994-present)

19

,

084501

(

2012

).

https://doi.org/10.1063/1.4714608

Google Scholar

Crossref

29.

J.

Xiao

,

J.

Liu

,

H.

Qin

, and

Z.

Yu

,

Phys. Plasmas

20

,

102517

(

2013

).

https://doi.org/10.1063/1.4826218

Google Scholar

Crossref

30.

J.

Xiao

,

J.

Liu

,

H.

Qin

,

Z.

Yu

, and

N.

Xiang

,

Phys. Plasmas (1994-present)

22

,

092305

(

2015

).

https://doi.org/10.1063/1.4930118

Google Scholar

Crossref

31.

E.

Evstatiev

and

B.

Shadwick

,

J. Comput. Phys.

245

,

376

(

2013

).

https://doi.org/10.1016/j.jcp.2013.03.006

Google Scholar

Crossref

32.

E.

Evstatiev

,

Comput. Phys. Commun.

185

,

2851

(

2014

).

https://doi.org/10.1016/j.cpc.2014.07.004

Google Scholar

Crossref

33.

B. A.

Shadwick

,

A. B.

Stamm

, and

E. G.

Evstatiev

,

Phys. Plasmas

21

,

055708

(

2014

).

https://doi.org/10.1063/1.4874338

Google Scholar

Crossref

34.

H.

Qin

,

J.

Liu

,

J.

Xiao

,

R.

Zhang

,

Y.

He

,

Y.

Wang

,

J. W.

Burby

,

L.

Ellison

, and

Y.

Zhou

, “

Canonical symplectic particle-in-cell method for long-term large-scale simulations of the Vlasov-Maxwell system

,”

Nucl. Fusion

(in press

2015

).

Google Scholar

35.

N.

Crouseilles

,

L.

Einkemmer

, and

E.

Faou

,

J. Comput. Phys.

283

,

224

(

2015

).

https://doi.org/10.1016/j.jcp.2014.11.029

Google Scholar

Crossref

36.

H.

Qin

,

Y.

He

,

R.

Zhang

,

J.

Liu

,

J.

Xiao

, and

Y.

Wang

,

J. Comput. Phys.

297

,

721

(

2015

).

https://doi.org/10.1016/j.jcp.2015.04.056

Google Scholar

Crossref

37.

P. J.

Morrison

,

Phys. Lett. A

80

,

383

(

1980

).

https://doi.org/10.1016/0375-9601(80)90776-8

Google Scholar

Crossref

38.

A.

Weinstein

and

P. J.

Morrison

,

Phys. Lett. A

86

,

235

(

1981

).

https://doi.org/10.1016/0375-9601(81)90496-5

Google Scholar

Crossref

39.

J. E.

Marsden

and

A.

Weinstein

,

Physica D

4

,

394

(

1982

).

https://doi.org/10.1016/0167-2789(82)90043-4

Google Scholar

Crossref

40.

J.

Burby

,

A.

Brizard

,

P.

Morrison

, and

H.

Qin

,

Phys. Lett. A

379

,

2073

(

2015

).

https://doi.org/10.1016/j.physleta.2015.06.051

Google Scholar

Crossref

41.

Y.

Zhou

,

H.

Qin

,

J.

Burby

, and

A.

Bhattacharjee

,

Phys. Plasmas (1994-present)

21

,

102109

(

2014

).

https://doi.org/10.1063/1.4897372

Google Scholar

Crossref

42.

Y.

Zhou

,

Y.-M.

Huang

,

H.

Qin

, and

A.

Bhattacharjee

, preprint arXiv:1509.08163 (

2015

).

43.

A. N.

Hirani

, “

Discrete exterior calculus

,” Ph.D. thesis (

California Institute of Technology

,

2003

).

Google Scholar

44.

A.

Stern

,

Y.

Tong

,

M.

Desbrun

, and

J. E.

Marsden

,

Geometry, Mechanics, and Dynamics

(

Springer

,

2015

), pp.

437

–

475

.

Google Scholar

Crossref

45.

H.

Whitney

,

Geometric Integration Theory

(

Princeton University Press

,

1957

).

Google Scholar

Crossref

46.

C. K.

Birdsall

and

A. B.

Langdon

,

Plasma Physics via Computer Simulation

(

IOP Publishing

,

1991

), p.

293

.

Google Scholar

Crossref

47.

R. W.

Hockney

and

J. W.

Eastwood

,

Computer Simulation Using Particles

(

CRC Press

,

1988

).

Google Scholar

Crossref

48.

C.

Nieter

and

J. R.

Cary

,

J. Comput. Phys.

196

,

448

(

2004

).

https://doi.org/10.1016/j.jcp.2003.11.004

Google Scholar

Crossref

49.

H.

Qin

,

Fields Inst. Commun.

46

,

171

(

2005

).

Google Scholar

50.

N.

Crouseilles

,

M.

Mehrenberger

, and

E.

Sonnendrucker

,

J. Comput. Phys.

229

,

1927

(

2010

).

https://doi.org/10.1016/j.jcp.2009.11.007

Google Scholar

Crossref

51.

J.

Dawson

,

Phys. Fluids (1958–1988)

4

,

869

(

1961

).

https://doi.org/10.1063/1.1706419

Google Scholar

Crossref

52.

T.

O'Neil

,

Phys. Fluids (1958–1988)

8

,

2255

(

1965

).

https://doi.org/10.1063/1.1761193

Google Scholar

Crossref

53.

C.

Mouhot

and

C.

Villani

,

Acta Math.

207

,

29

(

2011

).

https://doi.org/10.1007/s11511-011-0068-9

Google Scholar

Crossref

54.

G.

Manfredi

,

Phys. Rev. Lett.

79

,

2815

(

1997

).

https://doi.org/10.1103/PhysRevLett.79.2815

Google Scholar

Crossref

55.

T.

Zhou

,

Y.

Guo

, and

C.-W.

Shu

,

Physica D

157

,

322

(

2001

).

https://doi.org/10.1016/S0167-2789(01)00289-5

Google Scholar

Crossref

56.

T. H.

Stix

,

Waves in Plasmas

(

Springer

,

1992

), pp.

276

–

277

.

Google Scholar

57.

M.

Desbrun

,

E.

Kanso

, and

Y.

Tong

,

Discrete Differential Geometry

(

Springer

,

2008

), pp.

287

–

324

.

Google Scholar

Crossref

2015

AIP Publishing LLC

Explicit high-order non-canonical symplectic particle-in-cell algorithms for Vlasov-Maxwell systems

I. INTRODUCTION

II. NON-CANONICAL SYMPLECTIC PARTICLE-IN-CELL ALGORITHMS

III. NUMERICAL EXAMPLES

IV. SUMMARY AND DISCUSSION

ACKNOWLEDGMENTS

APPENDIX: HIGH ORDER INTERPOLATION FORMS FOR A CUBIC MESH

References

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Explicit high-order non-canonical symplectic particle-in-cell algorithms for Vlasov-Maxwell systems Free

I. INTRODUCTION

II. NON-CANONICAL SYMPLECTIC PARTICLE-IN-CELL ALGORITHMS

III. NUMERICAL EXAMPLES

IV. SUMMARY AND DISCUSSION

ACKNOWLEDGMENTS

APPENDIX: HIGH ORDER INTERPOLATION FORMS FOR A CUBIC MESH

References

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

Explicit high-order non-canonical symplectic particle-in-cell algorithms for Vlasov-Maxwell systems