Explicit structure-preserving geometric particle-in-cell algorithm in curvilinear orthogonal coordinate systems and its applications to whole-device 6D kinetic simulations of tokamak physics

The procedure of designing a structure-preserving geometric PIC algorithm starts from a field theory, or variational principle, for the system of charged particles and electromagnetic field. Instead of discretizing the corresponding Vlasov–Maxwell equations, the variational principle is discretized.

For spatial discretization of the electromagnetic field, discrete exterior calculus [71, 72] is adopted. As indicated by its name, a PIC algorithm contains two components not found in other simulation methods: charge and current deposition from particles to grid points, and field interpolation from grid points to particles. Note that these two components are independent from the field solver, e.g. the Yee-FDTD method, and the particle pusher, e.g. the Boris algorithm. In conventional PIC algorithms [1–3], the function of charge and current deposition and the function of field interpolation are implemented using intuitive techniques without fundamental guiding principles other than the consistency requirement. It was found [15, 24] that a systematic method to realize these two functions is to apply the Whitney interpolation (differential) forms [73] and their high-order generalizations [17] on the discrete Lagrangian.

The application of Whitney forms is a key technology for achieving the goal of preserving geometric structures and conservation laws. It stipulates from first principles how charge and current deposition and field interpolation should be performed [15–25]. It also enabled the derivation of the discrete non-canonical symplectic structure and Poisson bracket for the charged particle-electromagnetic field system [17]. At the continuous limit when the size of the space-time grid approaches zero, the discrete non-canonical Poisson bracket reduces to the Morrison–Marsden–Weinstein (MMW) bracket for the Vlasov–Maxwell equations [74–77, 26]. (The MMW bracket was also independently discovered by Iwinski and Turski [78].) To derive the discrete Poisson bracket, Whitney forms and their high-order generalizations are adopted for the purpose of geometric spatial discretization of the field Lagrangian density as an 1-form in the phase space of the charged particle-electromagnetic field system, whose exterior derivative automatically generates a discrete non-canonical symplectic structure and consequently a Poisson bracket. As a different approach, He et al [21] and Kraus et al [22] used the method of finite element exterior calculus to discretize the MMW bracket directly, and verified explicitly that the discretized bracket satisfies the Jacobi identity through lengthy calculations, in order for the discretized bracket to be a legitimate Poisson bracket. If one chooses to discretize the MMW bracket directly [32, 31], such verification is necessary because there are other discretizations of the MMW bracket which satisfy the numerical consistency requirement but not the Jacobi identity. On the other hand, the discrete Poisson bracket in [17] was not a discretization of the MMW bracket. It was derived from the spatially discretized 1-form using Whitney forms and is naturally endowed with a symplectic structure and Poisson bracket, and reduces to the MMW bracket in the continuous limit. The advantage of this discrete Poisson bracket in this respect demonstrates the power of Whitney forms in the design of structure-preserving geometric PIC algorithms.

To numerically integrate a Hamiltonian or Poisson system for the purpose of studying multiscale physics, symplectic integrators are necessities. Without an effective symplectic integrator, a symplectic or Poisson structure is not beneficial in terms of providing a better algorithm with long-term accuracy and fidelity. However, essentially all known symplectic integrators are designed for canonical Hamiltonian systems [43–47, 79–98] except for recent investigations on non-canonical symplectic integrators [99–112] for charged particle dynamics in a given electromagnetic field. Generic symplectic integrators for non-canonical symplectic systems are not known to exist. Fortunately, a high-order explicit symplectic integrator for the MMW bracket was discovered [19] using a splitting technique in the Cartesian coordinate system. The validity of this explicit non-canonical symplectic splitting method is built upon the serendipitous one-way decoupling of the dynamic components of the particle-field system in the Cartesian coordinate system. It is immediately realized [17] that this explicit non-canonical symplectic splitting applies to the discrete non-canonical Poisson bracket without modification for the charged particle-electromagnetic field system. This method was also adopted in [21, 22] to integrate the discretized MMW bracket using finite element exterior calculus. See also [113, 114] for early effort for developing symplectic splitting method for the Vlasov–Maxwell system. It was recently proven [24] that the discrete non-canonical Poisson bracket [17] and the explicit non-canonical symplectic splitting algorithm [19] can be equivalently formulated as a discrete field theory [48] using a Lagrangian discretized in both space and time.

The geometric discretization of the field theory for the system of charged particles and electromagnetic field using Whitney forms [73, 15, 17], discrete exterior calculus [71, 72], and explicit non-canonical symplectic integration [19, 17] results in a PIC algorithm preserving many important physical symmetries and conservation laws. In addition to the truncated infinitely dimensional symplectic structure, the algorithms preserves exactly the local conservation laws for charge [15, 24] and energy [23]. It was shown that the discrete local charge conservation law is a natural consequence of the discrete gauge symmetry admitted by the system [15, 24, 115]. The correspondence between discrete space-time translation symmetry [116, 117] and discrete local energy-momentum conservation law has also been established for the Maxwell system [118].

As mentioned in section 1, the goal of this section is to extend the explicit structure-preserving geometric PIC algorithm, especially the discrete non-canonical Poisson bracket [17] and the explicit non-canonical symplectic splitting algorithm [19], from the Cartesian coordinate systems to curvilinear orthogonal coordinate systems that are suitable for the tokamak geometry.

We start from the action integral of the system of charged particles and electromagnetic field. In a curvilinear orthogonal coordinate system $\left({x}_{1},{x}_{2},{x}_{3}\right)$ with line element

$$\begin{eqnarray}\begin{array}{rcl}{\rm{d}}{s}^{2} & = & {h}_{1}{\left({x}_{1},{x}_{2},{x}_{3}\right)}^{2}{\rm{d}}{x}_{1}^{2}+{h}_{2}{\left({x}_{1},{x}_{2},{x}_{3}\right)}^{2}{\rm{d}}{x}_{2}^{2}\\ & & +{h}_{3}{\left({x}_{1},{x}_{2},{x}_{3}\right)}^{2}{\rm{d}}{x}_{3}^{2},\end{array}\end{eqnarray} \tag{ 1 }$$

the action integral of the system is

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal A }[{{\bf{x}}}_{{sp}},{\dot{{\bf{x}}}}_{{sp}},{\bf{A}},\phi ]\\ \,=\displaystyle \int {\rm{d}}t\displaystyle \sum _{s,p}\left({L}_{{sp}}\left({m}_{s},{{\bf{v}}}_{{sp}}\right)+{q}_{s}\left({{\bf{v}}}_{{sp}}\cdot {\bf{A}}\left({{\bf{x}}}_{{sp}},t\right)-\phi \left({{\bf{x}}}_{{sp}},t\right)\right)\right)\\ \,+\,\displaystyle \int {\rm{d}}V{\rm{d}}t\displaystyle \frac{1}{2}\left({\left(-\dot{{\bf{A}}}\left({\bf{x}},t\right)-{\rm{\nabla }}\phi \left({\bf{x}},t\right)\right)}^{2}-{\left({\rm{\nabla }}\times {\bf{A}}\left({\bf{x}},t\right)\right)}^{2}\right),\end{array}\end{eqnarray} \tag{ 2 }$$

where x = [x₁, x₂, x₃] and ${\rm{d}}V=| {h}_{1}\left({\bf{x}}\right){h}_{2}\left({\bf{x}}\right){h}_{3}\left({\bf{x}}\right)| {\rm{d}}{x}_{1}{\rm{d}}{x}_{2}{\rm{d}}{x}_{3}$. In this coordinate system, the position and velocity of the pth particle of species s are ${{\bf{x}}}_{{sp}}=[{x}_{{{sp}}_{1}},{x}_{{{sp}}_{2}},{x}_{{{sp}}_{3}}]$ and ${{\bf{v}}}_{{sp}}=[{\dot{x}}_{{{sp}}_{1}}{h}_{1}\left({{\bf{x}}}_{{sp}}\right),{\dot{x}}_{{{sp}}_{2}}{h}_{2}\left({{\bf{x}}}_{{sp}}\right),{\dot{x}}_{{{sp}}_{3}}{h}_{3}\left({{\bf{x}}}_{{sp}}\right)]$. L_sp is the free Lagrangian for the pth particle of species s. For the non-relativistic case L_sp = m_s ∣v_sp ∣²/2, and for the relativistic case ${L}_{{sp}}=-{m}_{s}\sqrt{1-| {{\bf{v}}}_{{sp}}{| }^{2}}$. Here, we set both permittivity ₀ and permeability μ₀ in the vacuum to 1 to simplify the notation. The evolution of this system is governed by the variational principle,

$$\begin{eqnarray}&&\displaystyle \frac{\delta { \mathcal A }}{\delta {\bf{A}}}=0,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\delta { \mathcal A }}{\delta \phi }=0,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\delta { \mathcal A }}{\delta {{\bf{x}}}_{{sp}}}=0.\end{eqnarray} \tag{ 5 }$$

It can be verified that equations (3) and (4) are the Maxwell equations, and equation (5) is Newton's equation for the pth particle of species s.

According to the general principle of structure-preserving geometric algorithm, the first step is to discretize the field theory [15, 17, 24, 25]. For the electromagnetic field, we use the technique of Whitney forms and discrete manifold [73, 71]. For example, in a tetrahedron mesh the magnetic field B as a 2-form field can be discretized on a 2-simplex σ₂ (a side of a 3D tetrahedron mesh) as ${\int }_{{\sigma }_{2}}{\bf{B}}\cdot {\rm{d}}{\bf{S}}$, where S is the unit normal vector of σ₂. Note that σ₂ is a common side of two tetrahedron cells. In a mesh constructed using a curvilinear orthogonal coordinate system, which will be referred to as a Curvilinear Orthogonal Mesh (COM), such discretization can be also performed. Let

$$\begin{eqnarray}&&{\phi }_{i,j,k,l}=\phi \left({x}_{J,l}^{4}\right),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{\bf{A}}}_{i,j,k,l}={\left[\begin{array}{c}{A}_{{x}_{1}}\left({x}_{J,l}^{4}\right){h}_{1}\left({x}_{J,l}^{4}\right),\\ {A}_{{x}_{2}}\left({x}_{J,l}^{4}\right){h}_{2}\left({x}_{J,l}^{4}\right),\\ {A}_{{x}_{3}}\left({x}_{J,l}^{4}\right){h}_{3}\left({x}_{J,l}^{4}\right)\end{array}\right]}^{{\rm{T}}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{{\bf{B}}}_{i,j,k,l}={\left[\begin{array}{c}{B}_{{x}_{1}}\left({x}_{J,l}^{4}\right){h}_{2}\left({x}_{J,l}^{4}\right){h}_{3}\left({x}_{J,l}^{4}\right),\\ {B}_{{x}_{2}}\left({x}_{J,l}^{4}\right){h}_{3}\left({x}_{J,l}^{4}\right){h}_{1}\left({x}_{J,l}^{4}\right),\\ {B}_{{x}_{3}}\left({x}_{J,l}^{4}\right){h}_{1}\left({x}_{J,l}^{4}\right){h}_{2}\left({x}_{J,l}^{4}\right)\end{array}\right]}^{{\rm{T}}},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\rho }_{i,j,k,l}=\rho \left({x}_{J,l}^{4}\right){h}_{1}\left({x}_{J,l}^{4}\right){h}_{2}\left({x}_{J,l}^{4}\right){h}_{3}\left({x}_{J,l}^{4}\right),\end{eqnarray} \tag{ 9 }$$

where $\left({x}_{J,l}^{4}\right)$ denotes $\left(i{\rm{\Delta }}{x}_{1},j{\rm{\Delta }}{x}_{2},k{\rm{\Delta }}{x}_{3},l{\rm{\Delta }}t\right)$, ϕ_i,j,k,l , A_i,j,k,l , B_i,j,k,l and ρ_i,j,k,l are discrete 0-, 1-, 2- and 3-forms, respectively. In this discretization, the discrete gradient, curl and divergence operators can simply be finite difference operators. In the present study, the following difference operators are adopted,

$$\begin{eqnarray}&&{\left({{\rm{\nabla }}}_{{\rm{d}}}\phi \right)}_{i,j,k}=[{\phi }_{i+1,j,k}-{\phi }_{i,j,k},{\phi }_{i,j+1,k}-{\phi }_{i,j,k},{\phi }_{i,j,k+1}-{\phi }_{i,j,k}],\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\left({{\rm{\nabla }}}_{{\rm{d}}}\times {\bf{A}}\right)}_{i,j,k}={\left[\begin{array}{c}\left({{A}_{{x}_{3}}}_{i,j+1,k}-{{A}_{{x}_{3}}}_{i,j,k}\right)-\left({{A}_{{x}_{2}}}_{i,j,k+1}-{{A}_{{x}_{2}}}_{i,j,k}\right)\\ \left({{A}_{{x}_{1}}}_{i,j,k+1}-{{A}_{{x}_{1}}}_{i,j,k}\right)-\left({{A}_{{x}_{3}}}_{i+1,j,k}-{{A}_{{x}_{3}}}_{i,j,k}\right)\\ \left({{A}_{{x}_{2}}}_{i+1,j,k}-{{A}_{{x}_{2}}}_{i,j,k}\right)-\left({{A}_{{x}_{1}}}_{i,j+1,k}-{{A}_{{x}_{1}}}_{i,j,k}\right)\end{array}\right]}^{{\rm{T}}},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\left({{\rm{\nabla }}}_{{\rm{d}}}{\bf{B}}\right)}_{i,j,k} & = & \left({{B}_{{x}_{1}}}_{i+1,j,k}-{{B}_{{x}_{1}}}_{i,j,k}\right)+\left({{B}_{{x}_{2}}}_{i,j+1,k}-{{B}_{{x}_{2}}}_{i,j,k}\right)\\ & & +\left({{B}_{z}}_{i,j,k+1}-{{B}_{z}}_{i,j,k}\right).\end{array}\end{eqnarray} \tag{ 12 }$$

To discretize the particle-field interaction while preserving geometric structures and symmetries, Whitney interpolating forms [73] and their high-order generalizations [17] are used. Akin to previous results on Whitney interpolating forms in a cubic mesh [17], interpolating forms ${W}_{{\rm{O}}{\sigma }_{0}J}\left({\bf{x}}\right),{W}_{{\rm{O}}{\sigma }_{1}J}\left({\bf{x}}\right),{W}_{{\rm{O}}{\sigma }_{2}J}\left({\bf{x}}\right)$ and ${W}_{{\rm{O}}{\sigma }_{3}J}\left({\bf{x}}\right)$ for 0-, 1-, 2- and 3-forms can be constructed on a COM as follows,

$$\begin{eqnarray}&&{W}_{{\rm{O}}{\sigma }_{0}J}\left({\bf{x}}\right)={W}_{{\sigma }_{0}J}\left({\bf{x}}\right),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{W}_{{\rm{O}}{\sigma }_{1}J}\left({\bf{x}}\right)={W}_{{\sigma }_{1}J}\left({\bf{x}}\right)/[{h}_{1}\left({\bf{x}}\right),{h}_{2}\left({\bf{x}}\right),{h}_{3}\left({\bf{x}}\right)],\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{W}_{{\rm{O}}{\sigma }_{2}J}\left({\bf{x}}\right)={W}_{{\sigma }_{2}J}\left({\bf{x}}\right)/[{h}_{2}\left({\bf{x}}\right){h}_{3}\left({\bf{x}}\right),{h}_{3}\left({\bf{x}}\right){h}_{1}\left({\bf{x}}\right),{h}_{1}\left({\bf{x}}\right){h}_{2}\left({\bf{x}}\right)],\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{W}_{{\rm{O}}{\sigma }_{3}J}\left({\bf{x}}\right)={W}_{{\sigma }_{3}J}\left({\bf{x}}\right)/{h}_{1}\left({\bf{x}}\right){h}_{2}\left({\bf{x}}\right){h}_{3}\left({\bf{x}}\right),\end{eqnarray} \tag{ 16 }$$

where ${W}_{\sigma {iJ}}\left({\bf{x}}\right),0\leqslant i\leqslant 3$ are the interpolating forms in a cubic mesh defined in [17] and the quotient (product) of vectors means component-wise quotient (product), i.e.

$$\begin{eqnarray}&&[{A}_{{x}_{1}},{A}_{{x}_{2}},{A}_{{x}_{3}}]/[{B}_{{x}_{1}},{B}_{{x}_{2}},{B}_{{x}_{3}}]=\left[\displaystyle \frac{{A}_{{x}_{1}}}{{B}_{{x}_{1}}},\displaystyle \frac{{A}_{{x}_{2}}}{{B}_{{x}_{2}}},\displaystyle \frac{{A}_{{x}_{3}}}{{B}_{{x}_{3}}}\right].\end{eqnarray} \tag{ 17 }$$

With these operators and interpolating forms, we discretize the action integral as

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal A }}_{{\rm{d}}} & = & \displaystyle \sum _{s,p,l}\left({L}_{{sp}}\left({m}_{s},{{\bf{v}}}_{{sp},l}\right)+{q}_{s}\left(\overline{{{\bf{v}}}_{{sp}}\cdot {\bf{A}}}\left({{\bf{x}}}_{{sp},l},{{\bf{x}}}_{{sp},l+1},l\right)\right.\right.\\ & & \left.\left.-\phi \left({{\bf{x}}}_{{sp},l+1},l+1\right)\right)\right)\\ & & +\displaystyle \frac{{h}_{\sigma 3}\left({{\bf{x}}}_{J}\right)}{2}\displaystyle \sum _{J}\left({\left(-\displaystyle \frac{{{\bf{A}}}_{J,l+1}-{{\bf{A}}}_{J,l}}{{\rm{\Delta }}t{{\bf{h}}}_{\sigma 1}\left({{\bf{x}}}_{J}\right)}-\displaystyle \frac{{\left({{\rm{\nabla }}}_{{\rm{d}}}\phi \right)}_{J,l+1}}{{{\bf{h}}}_{\sigma 1}\left({{\bf{x}}}_{J}\right)}\right)}^{2}\right.\\ & & -\left.{\left(\displaystyle \frac{{\left({{\rm{\nabla }}}_{{\rm{d}}}\times {\bf{A}}\right)}_{K,l}}{{{\bf{h}}}_{\sigma 2}\left({{\bf{x}}}_{K}\right)}\right)}^{2}\right),\end{array}\end{eqnarray} \tag{ 18 }$$

where l is the index for the temporal grid, J and K are the index vectors for the spatial grid, i.e., J = i, j, k and K = i', j', k', and

$$\begin{eqnarray}&&{{\bf{h}}}_{\sigma 1}\left({\bf{x}}\right)=[{h}_{1}\left({\bf{x}}\right),{h}_{2}\left({\bf{x}}\right),{h}_{3}\left({\bf{x}}\right)],\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{{\bf{h}}}_{\sigma 2}\left({\bf{x}}\right)=[{h}_{2}\left({\bf{x}}\right){h}_{3}\left({\bf{x}}\right),{h}_{3}\left({\bf{x}}\right){h}_{1}\left({\bf{x}}\right),{h}_{1}\left({\bf{x}}\right){h}_{2}\left({\bf{x}}\right)],\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{h}_{\sigma 3}\left({\bf{x}}\right)={h}_{1}\left({\bf{x}}\right){h}_{2}\left({\bf{x}}\right){h}_{3}\left({\bf{x}}\right),\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{{\bf{v}}}_{{sp},l}=\displaystyle \frac{{{\bf{x}}}_{{sp},l+1}-{{\bf{x}}}_{{sp},l}}{{\rm{\Delta }}t}{{\bf{h}}}_{\sigma 1}\left(\displaystyle \frac{{{\bf{x}}}_{{sp},l+1}+{{\bf{x}}}_{{sp},l}}{2}\right),\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&\phi \left({\bf{x}},l\right)=\displaystyle \sum _{J}{\phi }_{J,l}{W}_{{\sigma }_{0}J}\left({\bf{x}}\right),\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}\begin{array}{rcl}\overline{{{\bf{v}}}_{{sp}}\cdot {\bf{A}}}\left({{\bf{x}}}_{{sp},l},{{\bf{x}}}_{{sp},l+1},l\right) & = & {\displaystyle \int }_{0}^{1}{\rm{d}}\tau \left[\begin{array}{c}\displaystyle \frac{{x}_{1,{sp},l+1}-{x}_{1,{sp},l}}{{\rm{\Delta }}t}{h}_{1}\left({x}_{1,{sp},l}+\tau \left({x}_{1,{sp},l+1}-{x}_{1,{sp},l},{x}_{2,{sp},l},{x}_{3,{sp},l}\right)\right),\\ \displaystyle \frac{{x}_{2,{sp},l+1}-{x}_{2,{sp},l}}{{\rm{\Delta }}t}{h}_{2}\left({x}_{1,{sp},l+1},{x}_{2,{sp},l}+\tau \left({x}_{2,{sp},l+1}-{x}_{2,{sp},l},{x}_{3,{sp},l}\right)\right),\\ \displaystyle \frac{{x}_{3,{sp},l+1}-{x}_{3,{sp},l}}{{\rm{\Delta }}t}{h}_{3}\left({x}_{1,{sp},l+1},{x}_{2,{sp},l+1},{x}_{3,{sp},l}+\tau \left({x}_{3,{sp},l+1}-{x}_{3,{sp},l}\right)\right)\end{array}\right]\\ & & \cdot \left[\begin{array}{l}{\tilde{A}}_{{x}_{1},l}\left({x}_{1,{sp},l}+\tau \left({x}_{1,{sp},l+1}-{x}_{1,{sp},l}\right),{x}_{2,{sp},l},{x}_{3,{sp},l}\right),\\ {\tilde{A}}_{{x}_{2},l}\left({x}_{1,{sp},l+1},{x}_{2,{sp},l}+\tau \left({x}_{2,{sp},l+1}-{x}_{2,{sp},l}\right),{x}_{3,{sp},l}\right),\\ {\tilde{A}}_{{x}_{3},l}\left({x}_{1,{sp},l+1},{x}_{2,{sp},l+1},{x}_{3,{sp},l}+\tau \left({x}_{3,{sp},l+1}-{x}_{3,{sp},l}\right)\right)\end{array}\right]\end{array}\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}\begin{array}{rcl} & = & {\displaystyle \int }_{0}^{1}{\rm{d}}\tau \left[\displaystyle \frac{{x}_{1,{sp},l+1}-{x}_{1,{sp},l}}{{\rm{\Delta }}t},\displaystyle \frac{{x}_{2,{sp},l+1}-{x}_{2,{sp},l}}{{\rm{\Delta }}t},\displaystyle \frac{{x}_{3,{sp},l+1}-{x}_{3,{sp},l}}{{\rm{\Delta }}t}\right]\\ & & \cdot \left[\begin{array}{l}{A}_{{x}_{1},l}\left({x}_{1,{sp},l}+\tau \left({x}_{1,{sp},l+1}-{x}_{1,{sp},l}\right),{x}_{2,{sp},l},{x}_{3,{sp},l}\right),\\ {A}_{{x}_{2},l}\left({x}_{1,{sp},l+1},{x}_{2,{sp},l}+\tau \left({x}_{2,{sp},l+1}-{x}_{2,{sp},l}\right),{x}_{3,{sp},l}\right),\\ {A}_{{x}_{3},l}\left({x}_{1,{sp},l+1},{x}_{2,{sp},l+1},{x}_{3,{sp},l}+\tau \left({x}_{3,{sp},l+1}-{x}_{3,{sp},l}\right)\right)\end{array}\right],\end{array}\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&\left[\begin{array}{c}{\tilde{A}}_{{x}_{1},l}\left({\bf{x}}\right),\\ {\tilde{A}}_{{x}_{2},l}\left({\bf{x}}\right),\\ {\tilde{A}}_{{x}_{3},l}\left({\bf{x}}\right)\end{array}\right]\equiv \displaystyle \sum _{J}{{\bf{A}}}_{J,l}{W}_{{\rm{O}}{\sigma }_{1}J}\left({\bf{x}}\right),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\left[\begin{array}{c}{A}_{{x}_{1},l}\left({\bf{x}}\right),\\ {A}_{{x}_{2},l}\left({\bf{x}}\right),\\ {A}_{{x}_{3},l}\left({\bf{x}}\right)\end{array}\right]\equiv \displaystyle \sum _{J}{{\bf{A}}}_{J,l}{W}_{{\sigma }_{1}J}\left({\bf{x}}\right).\end{eqnarray} \tag{ 27 }$$

Finally, the time advance rule is given by the variation of the action integral with respect to the discrete fields,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{ \mathcal A }}_{{\rm{d}}}}{\partial {{\bf{x}}}_{{sp},l}}=0,\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{ \mathcal A }}_{{\rm{d}}}}{\partial {{\bf{A}}}_{J,l}}=0,\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{ \mathcal A }}_{{\rm{d}}}}{\partial {\phi }_{J,l}}=0.\end{eqnarray} \tag{ 30 }$$

From equation (28), the variation with respect to x_sp,l leads to

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial {{\bf{x}}}_{l}}\left({L}_{{sp}}\left({m}_{s},{{\bf{v}}}_{{sp},l-1}\right)+{L}_{{sp}}\left({m}_{s},{{\bf{v}}}_{{sp},l}\right)\right)\\ \quad ={q}_{s}\left({{\bf{E}}}_{l}\left({{\bf{x}}}_{{sp},l}\right)+\overline{{{\bf{v}}}_{{sp}}\times {\bf{B}}}\left(l-1,l,l+1\right)\right),\end{array}\end{eqnarray} \tag{ 31 }$$

where

$$\begin{eqnarray}&&{{\bf{E}}}_{l}\left({\bf{x}}\right)=\displaystyle \sum _{J}{{\bf{E}}}_{J,l}{W}_{{\sigma }_{1}J}\left({\bf{x}}\right),\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{{\bf{E}}}_{J,l}=-\displaystyle \frac{{{\bf{A}}}_{J,l}-{{\bf{A}}}_{J,l-1}}{{\rm{\Delta }}t},\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&\begin{array}{l}\overline{{{\bf{v}}}_{{sp}}\times {\bf{B}}}\left(l-1,l,l+1\right)\\ \,=\,\displaystyle \frac{1}{{\rm{\Delta }}t}\left[\begin{array}{l}{\displaystyle \int }_{{x}_{2,{sp},l-1}}^{{x}_{2,{sp},l}}{\rm{d}}{x}_{2}{B}_{{x}_{3},l-1}\left({x}_{1,{sp},l},{x}_{2},{x}_{3,{sp},l-1}\right),\\ {\displaystyle \int }_{{x}_{3,{sp},l-1}}^{{x}_{3,{sp},l}}{\rm{d}}{x}_{3}{B}_{{x}_{1},l-1}\left({x}_{1,{sp},l},{x}_{2,{sp},l},{x}_{3}\right),\\ {\displaystyle \int }_{{x}_{1,{sp},l}}^{{x}_{1,{sp},l+1}}{\rm{d}}{x}_{1}{B}_{{x}_{2},l}\left({x}_{1},{x}_{2,{sp},l},{x}_{3,{sp},l}\right)\end{array}\right]\\ \,-\,\displaystyle \frac{1}{{\rm{\Delta }}t}\left[\begin{array}{l}{\displaystyle \int }_{{x}_{3,{sp},l-1}}^{{x}_{3,{sp},l}}{\rm{d}}{x}_{3}{B}_{{x}_{2},l-1}\left({x}_{1,{sp},l},{x}_{2,{sp},l},{x}_{3}\right),\\ {\displaystyle \int }_{{x}_{1,{sp},l}}^{{x}_{1,{sp},l+1}}{\rm{d}}{x}_{1}{B}_{{x}_{3},l}\left({x}_{1},{x}_{2,{sp},l},{x}_{3,{sp},l}\right),\\ {\displaystyle \int }_{{x}_{2,{sp},l}}^{{x}_{2,{sp},l+1}}{\rm{d}}{x}_{2}{B}_{{x}_{1},l}\left({x}_{1,{sp},l+1},{x}_{2},{x}_{3,{sp},l}\right)\end{array}\right],\end{array}\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&\left[\begin{array}{c}{B}_{{x}_{1},l}\left({\bf{x}}\right),\\ {B}_{{x}_{2},l}\left({\bf{x}}\right),\\ {B}_{{x}_{3},l}\left({\bf{x}}\right)\end{array}\right]=\displaystyle \sum _{J}{{\bf{B}}}_{J,l}{W}_{{\sigma }_{2}J}\left({\bf{x}}\right),\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{{\bf{B}}}_{K,l}=\displaystyle \sum _{J}{{\rm{\nabla }}}_{{\rm{d}}}{\times }_{K,J}{{\bf{A}}}_{J,l}.\end{eqnarray} \tag{ 36 }$$

To reduce simulation noise, we use 2nd-order Whitney forms for field interpolation. The concept of 2nd-order Whitney forms and their constructions were systematically developed in [17]. In general, only piece-wise polynomials are used to construct high-order Whitney forms, and it is straightforward to calculate integrals along particles' trajectories. These integrals can be calculated explicitly even with more complex Whitney interpolating forms, because these interpolating forms contain derivatives that are easy to integrate.

From equation (29), the discrete Ampere's law is

$$\begin{eqnarray*}&&\begin{array}{l}\displaystyle \frac{1}{{{\bf{h}}}_{\sigma 1}\left({{\bf{x}}}_{J}\right){{\bf{h}}}_{\sigma 1}\left({{\bf{x}}}_{J}\right)}\displaystyle \frac{{{\bf{E}}}_{J,l+1}-{{\bf{E}}}_{J,l}}{{\rm{\Delta }}t}{h}_{\sigma 3}\left({{\bf{x}}}_{J}\right)\\ \quad =\displaystyle \sum _{K}{{\rm{\nabla }}}_{{\rm{d}}}^{{\rm{T}}}{\times }_{J,K}\displaystyle \frac{{h}_{\sigma 3}\left({{\bf{x}}}_{J}\right)}{{{\bf{h}}}_{\sigma 2}\left({{\bf{x}}}_{K}\right){{\bf{h}}}_{\sigma 2}\left({{\bf{x}}}_{K}\right)}{{\bf{B}}}_{K,l}-{{\bf{j}}}_{J,l},\end{array}\end{eqnarray*}$$

where

$$\begin{eqnarray}&&{{\bf{j}}}_{J,l}=\displaystyle \frac{1}{{\rm{\Delta }}t}\displaystyle \sum _{s,p}{q}_{s}{\int }_{{C}_{{sp},l,l+1}}{W}_{{\sigma }_{1}J}\left({\bf{x}}\right){\rm{d}}{\bf{x}},\end{eqnarray} \tag{ 37 }$$

and the integral path C_sp,l,l+1 is defined as a zigzag path from x_sp,l to x_sp,l+1, i.e.

$$\begin{eqnarray}\begin{array}{rcl}{C}_{{sp},l,l+1} & = & \left\{\left({x}_{1,{sp},l}+\tau \left({x}_{1,{sp},l+1}-{x}_{1,{sp},l}\right),{x}_{2,{sp},l},{x}_{3,{sp},l}\right)| \tau \in \left[0,1\right)\right\}\bigcup \\ & & \left\{\left({x}_{1,{sp},l+1},{x}_{2,{sp},l}+\tau \left({x}_{2,{sp},l+1}-{x}_{2,{sp},l}\right),{x}_{3,{sp},l}\right)| \tau \in \left[0,1\right)\right\}\bigcup \\ & & \left\{\left({x}_{1,{sp},l+1},{x}_{2,{sp},l+1},{x}_{3,{sp},l}+\tau \left({x}_{3,{sp},l+1}-{x}_{3,{sp},l}\right)\right)| \tau \in \left[0,1\right)\right\}.\end{array}\end{eqnarray} \tag{ 38 }$$

Using the definition of B_K,l , i.e. equation (36), we can obtain its discrete time evolution,

$$\begin{eqnarray}&&\displaystyle \frac{{{\bf{B}}}_{K,l}-{{\bf{B}}}_{K,l-1}}{{\rm{\Delta }}t}=-\displaystyle \sum _{J}{{\rm{\nabla }}}_{{\rm{d}}}{\times }_{K,J}{{\bf{E}}}_{J,l},\end{eqnarray} \tag{ 39 }$$

which is the discrete version of Faraday's law.

Generally the above scheme is implicit. However, if particles are non-relativistic and the line element vector h of a curvilinear orthogonal coordinate system satisfies the following condition

$$\begin{eqnarray}&&\displaystyle \frac{\partial {h}_{1}}{\partial {x}_{1}}=\displaystyle \frac{\partial {h}_{2}}{\partial {x}_{2}}=\displaystyle \frac{\partial {h}_{3}}{\partial {x}_{3}}=0,\end{eqnarray} \tag{ 40 }$$

then high-order explicit schemes exist in the COM constructed using this coordinate system. The cylindrical coordinate is such a case. In section 2.5 we derive the 2nd-order explicit scheme for the non-relativistic Vlasov–Maxwell system in the cylindrical mesh.

For the structure-preserving geometric PIC algorithm presented in section 2.3, there exists a corresponding discrete Poisson bracket. When particles are non-relativistic and condition (40) is satisfied, an associate splitting algorithm, which is explicit and symplectic, can be constructed. The algorithm is similar to and generalizes the explicit non-canonical symplectic splitting algorithm in the Cartesian coordinate system designed in [19]. Since the algorithm formulated in section 2.3 is independent from these constructions, we only list the results here without detailed derivations.

In a COM built on a curvilinear orthogonal coordinate system, we can choose to discretize the space only using the same method described above to obtain

$$\begin{eqnarray}\begin{array}{rcl}{L}_{{sd}} & = & \displaystyle \sum _{s,p}\left({L}_{{sp}}\left({m}_{s},{\dot{{\bf{x}}}}_{{sp}}\right)+{q}_{s}{\dot{{\bf{x}}}}_{{sp}}\cdot \displaystyle \sum _{J}{{\bf{A}}}_{J}{W}_{{\sigma }_{1}J}\left({{\bf{x}}}_{{sp}}\right)\right)\\ & & +\displaystyle \frac{1}{2}\displaystyle \sum _{J}\left({\left(\displaystyle \frac{-{\dot{{\bf{A}}}}_{J}}{{{\bf{h}}}_{\sigma 1}\left({{\bf{x}}}_{J}\right)}\right)}^{2}-{\left(\displaystyle \frac{{\left({{\rm{\nabla }}}_{{\rm{d}}}\times {\bf{A}}\right)}_{J,l}}{{{\bf{h}}}_{\sigma 2}\left({{\bf{x}}}_{J}\right)}\right)}^{2}\right){h}_{\sigma 3}\left({{\bf{x}}}_{J}\right),\end{array}\end{eqnarray} \tag{ 41 }$$

where the temporal gauge (ϕ = 0) has been adopted. Following the procedure in [17], a non-canonical symplectic structure can be constructed from this Lagrangian, and the associated discrete Poisson bracket is

$$\begin{eqnarray}\begin{array}{rcl}\left\{F,G\right\} & = & \displaystyle \sum _{J}\left(\displaystyle \frac{\partial F}{\partial {{\bf{E}}}_{J}}\cdot \displaystyle \frac{\mathrm{diag}\left({{\bf{h}}}_{\sigma 1}{\left({{\bf{x}}}_{J}\right)}^{2}\right)}{{h}_{\sigma 3}\left({{\bf{x}}}_{J}\right)}\cdot \displaystyle \sum _{K}\displaystyle \frac{\partial G}{\partial {{\bf{B}}}_{K}}{{\rm{\nabla }}}_{{\rm{d}}}{\times }_{{KJ}}\right.\\ & & -\,\left.\displaystyle \sum _{K}\displaystyle \frac{\partial F}{\partial {{\bf{B}}}_{K}}{{\rm{\nabla }}}_{{\rm{d}}}{\times }_{{KJ}}\cdot \displaystyle \frac{\mathrm{diag}\left({{\bf{h}}}_{\sigma 1}{\left({{\bf{x}}}_{J}\right)}^{2}\right)}{{h}_{\sigma 3}\left({{\bf{x}}}_{J}\right)}\cdot \displaystyle \frac{\partial G}{\partial {{\bf{E}}}_{J}}\right)\\ & & +\,\displaystyle \sum _{s,p}\displaystyle \frac{1}{{m}_{s}}\displaystyle \sum _{i=1}^{3}\displaystyle \frac{1}{{h}_{i}{\left({{\bf{x}}}_{{sp}}\right)}^{2}}\left(\displaystyle \frac{\partial F}{\partial {x}_{i,{sp}}}\cdot \displaystyle \frac{\partial G}{\partial {\dot{x}}_{i,{sp}}}-\displaystyle \frac{\partial F}{\partial {\dot{x}}_{i,{sp}}}\cdot \displaystyle \frac{\partial G}{\partial {x}_{i,{sp}}}\right)\\ & & +\displaystyle \sum _{s,p}\displaystyle \sum _{J}\displaystyle \sum _{i=1}^{3}\displaystyle \frac{{h}_{\sigma 1,{x}_{i}}{\left({{\bf{x}}}_{J}\right)}^{2}{q}_{s}}{{m}_{s}{h}_{i}{\left({{\bf{x}}}_{{sp}}\right)}^{2}{h}_{\sigma 3,{x}_{i}}\left({{\bf{x}}}_{J}\right)}\\ & & \times \left(\displaystyle \frac{\partial F}{\partial {\dot{x}}_{i,{sp}}}{W}_{{\sigma }_{1J}}\left({x}_{i,{sp}}\right)\displaystyle \frac{\partial G}{\partial {E}_{{x}_{i},J}}-\displaystyle \frac{\partial G}{\partial {\dot{x}}_{i,{sp}}}{W}_{{\sigma }_{1J}}\left({x}_{i,{sp}}\right)\displaystyle \frac{\partial F}{\partial {E}_{{x}_{i},J}}\right)\\ & & +\displaystyle \sum _{s}\displaystyle \sum _{K}\displaystyle \frac{1}{{m}_{s}^{2}}\left(\mathrm{diag}\left(\displaystyle \frac{1}{{{\bf{h}}}_{\sigma 1}{\left({{\bf{x}}}_{{sp}}\right)}^{2}}\right)\cdot \displaystyle \frac{\partial F}{\partial {\dot{{\bf{x}}}}_{{sp}}}\right)\\ & & \times \left(\left({{\rm{\nabla }}}_{{{\bf{x}}}_{{sp}}}\times \left({m}_{s}\mathrm{diag}\left({{\bf{h}}}_{\sigma 1}{\left({{\bf{x}}}_{{sp}}\right)}^{2}\right)\cdot {\dot{{\bf{x}}}}_{{sp}}\right)+{q}_{s}{W}_{{\sigma }_{2K}}\left({{\bf{x}}}_{{sp}}\right){{\bf{B}}}_{K}\right)\right.\\ & & \left.\times \left(\displaystyle \frac{\partial G}{\partial {\dot{{\bf{x}}}}_{{sp}}}\cdot \mathrm{diag}\left(\displaystyle \frac{1}{{{\bf{h}}}_{\sigma 1}{\left({{\bf{x}}}_{{sp}}\right)}^{2}}\right)\right)\right),\end{array}\end{eqnarray} \tag{ 42 }$$

where $\mathrm{diag}\left({\bf{v}}\right)$ means the diagonal matrix

$$\begin{eqnarray*}\mathrm{diag}\left({\bf{v}}\right)=\left[\begin{array}{ccc}{v}_{x} & 0 & 0\\ 0 & {v}_{y} & 0\\ 0 & 0 & {v}_{z}\end{array}\right],\end{eqnarray*}$$

and the corresponding Hamiltonian is

$$\begin{eqnarray}\begin{array}{rcl}{H}_{\mathrm{sd}} & = & \displaystyle \frac{1}{2}\displaystyle \sum _{J}\left({\left|\displaystyle \frac{{{\bf{E}}}_{J}}{{{\bf{h}}}_{\sigma 1}\left({{\bf{x}}}_{J}\right)}\right|}^{2}+{\left|\displaystyle \frac{{{\bf{B}}}_{J,l}}{{{\bf{h}}}_{\sigma 2}\left({{\bf{x}}}_{J}\right)}\right|}^{2}\right){h}_{\sigma 3}\left({{\bf{x}}}_{J}\right)\\ & & +\displaystyle \sum _{{sp}}\displaystyle \sum _{i}\displaystyle \frac{{m}_{s}}{2}{\dot{x}}_{i,{sp}}^{2}{h}_{i}{\left({{\bf{x}}}_{{sp}}\right)}^{2}.\end{array}\end{eqnarray} \tag{ 43 }$$

Here particles are assumed to be non-relativistic and E_J , B_K are spatially discretized electromagnetic fields defined as

$$\begin{eqnarray}&&{{\bf{E}}}_{J}=-{\dot{{\bf{A}}}}_{J},\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}&&{{\bf{B}}}_{K}=\displaystyle \sum _{J}{{\rm{\nabla }}}_{{\rm{d}}}{\times }_{{KJ}}{{\bf{A}}}_{J}.\end{eqnarray} \tag{ 45 }$$

The Poisson bracket given by equation (42) generalizes the previous Cartesian version [17] to arbitrary COMs. It automatically satisfies the Jacobi identity because it is derived from a Lagrangian 1-form. See [17] for detailed geometric constructions.

The dynamics equation, i.e. the Hamiltonian equation, is

$$\begin{eqnarray}&&\dot{F}=\{F,{H}_{\mathrm{sd}}\},\end{eqnarray} \tag{ 46 }$$

where

$$\begin{eqnarray}&&F=[{{\bf{E}}}_{J},{{\bf{B}}}_{J},{{\bf{x}}}_{{sp}},{\dot{{\bf{x}}}}_{{sp}}].\end{eqnarray} \tag{ 47 }$$

To build an explicit symplectic algorithm, we adopt the splitting method [19]. The Hamiltonian H_sd is divided into 5 parts,

$$\begin{eqnarray}&&{H}_{\mathrm{sd}}={H}_{E}+{H}_{B}+{H}_{1}+{H}_{2}+{H}_{3},\end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray}&&{H}_{E}=\displaystyle \frac{1}{2}\displaystyle \sum _{J}{\left|\displaystyle \frac{{{\bf{E}}}_{J}}{{{\bf{h}}}_{\sigma 1}\left({{\bf{x}}}_{J}\right)}\right|}^{2}{h}_{\sigma 3}\left({{\bf{x}}}_{J}\right),\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&{H}_{B}=\displaystyle \frac{1}{2}\displaystyle \sum _{J}{\left|\displaystyle \frac{{{\bf{B}}}_{J,l}}{{{\bf{h}}}_{\sigma 2}\left({{\bf{x}}}_{J}\right)}\right|}^{2}{h}_{\sigma 3}\left({{\bf{x}}}_{J}\right),\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&{H}_{i}=\displaystyle \sum _{{sp}}\displaystyle \frac{{m}_{s}}{2}{\dot{x}}_{i,{sp}}^{2}{h}_{i}{\left({{\bf{x}}}_{{sp}}\right)}^{2}\quad \mathrm{for}\ i\ \mathrm{in}\,\left\{1,2,3\right\}.\end{eqnarray} \tag{ 51 }$$

Each part defines a sub-system with the same Poisson bracket (42). It turns out that when condition (40) is satisfied the exact solution of each subsystem can be written down in a closed form, and explicit high-order symplectic algorithms for the entire system can be constructed by compositions using the exact solutions of the sub-systems. For H_E and H_B , the corresponding Hamiltonian equations are

$$\begin{eqnarray}&&\dot{F}=\left\{F,{H}_{E}\right\},\end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray}&&\dot{F}=\left\{F,{H}_{B}\right\},\end{eqnarray} \tag{ 53 }$$

i.e.

$$\begin{eqnarray}\left\{\begin{array}{rcl}{\dot{{\bf{E}}}}_{J} & = & 0,\\ {\dot{{\bf{B}}}}_{K} & = & -\displaystyle \sum _{J}{{\rm{\nabla }}}_{{\rm{d}}}{\times }_{{KJ}}{{\bf{E}}}_{J},\\ {\dot{{\bf{x}}}}_{{sp}} & = & 0,\\ {\ddot{{\bf{x}}}}_{{sp}} & = & \displaystyle \frac{{q}_{s}}{{m}_{s}}\mathrm{diag}\left(\displaystyle \frac{1}{{{\bf{h}}}_{\sigma 1}{\left({{\bf{x}}}_{{sp}}\right)}^{2}}\right)\cdot \displaystyle \sum _{J}{W}_{{\sigma }_{1J}}\left({{\bf{x}}}_{{sp}}\right){{\bf{E}}}_{J},\end{array}\right.\end{eqnarray} \tag{ 54 }$$

and

$$\begin{eqnarray}\left\{\begin{array}{rcl}{\dot{{\bf{E}}}}_{J} & = & \displaystyle \frac{\mathrm{diag}\left({{\bf{h}}}_{\sigma 1}{\left({{\bf{x}}}_{J}\right)}^{2}\right)}{{h}_{\sigma 3}\left({{\bf{x}}}_{J}\right)}\cdot \displaystyle \sum _{K}{h}_{\sigma 3}\left({{\bf{x}}}_{K}\right)\mathrm{diag}\left(\displaystyle \frac{1}{{{\bf{h}}}_{\sigma 2}{\left({{\bf{x}}}_{K}\right)}^{2}}\right)\cdot {{\rm{\nabla }}}_{{\rm{d}}}{\times }_{{KJ}}{{\bf{B}}}_{K},\\ {\dot{{\bf{B}}}}_{K} & = & 0,\\ {\dot{{\bf{x}}}}_{{sp}} & = & 0,\\ {\ddot{{\bf{x}}}}_{{sp}} & = & 0.\end{array}\right.\end{eqnarray} \tag{ 55 }$$

Their analytical solutions are

$$\begin{eqnarray}{{\rm{\Theta }}}_{E}:\left\{\begin{array}{rcl}{{\bf{E}}}_{J}\left(t+{\rm{\Delta }}t\right) & = & {{\bf{E}}}_{J}\left(t\right),\\ {{\bf{B}}}_{K}\left(t+{\rm{\Delta }}t\right) & = & {{\bf{B}}}_{K}\left(t\right)-{\rm{\Delta }}t\displaystyle \sum _{J}{{\rm{\nabla }}}_{{\rm{d}}}{\times }_{{KJ}}{{\bf{E}}}_{J}(t),\\ {{\bf{x}}}_{{sp}}\left(t+{\rm{\Delta }}t\right) & = & {{\bf{x}}}_{s}\left(t\right),\\ {\dot{{\bf{x}}}}_{{sp}}\left(t+{\rm{\Delta }}t\right) & = & {\dot{{\bf{x}}}}_{s}\left(t\right)+{\rm{\Delta }}t\displaystyle \frac{{q}_{s}}{{m}_{s}}\mathrm{diag}\left(\displaystyle \frac{1}{{{\bf{h}}}_{\sigma 1}{\left({{\bf{x}}}_{{sp}}\right)}^{2}}\right)\cdot \displaystyle \sum _{J}{W}_{{\sigma }_{1J}}\left({{\bf{x}}}_{{sp}}(t)\right){{\bf{E}}}_{J}(t),\end{array}\right.\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}{{\rm{\Theta }}}_{B}:\left\{\begin{array}{rcl}{{\bf{E}}}_{J}\left(t+{\rm{\Delta }}t\right) & = & {{\bf{E}}}_{J}\left(t\right)+{\rm{\Delta }}t\displaystyle \frac{\mathrm{diag}\left({{\bf{h}}}_{\sigma 1}{\left({{\bf{x}}}_{J}\right)}^{2}\right)}{{h}_{\sigma 3}\left({{\bf{x}}}_{J}\right)}\cdot \displaystyle \sum _{K}{{\rm{\nabla }}}_{{\rm{d}}}{\times }_{{KJ}}{h}_{\sigma 3}\left({{\bf{x}}}_{K}\right)\mathrm{diag}\left(\displaystyle \frac{1}{{{\bf{h}}}_{\sigma 2}{\left({{\bf{x}}}_{K}\right)}^{2}}\right)\cdot {{\bf{B}}}_{K}(t),\\ {{\bf{B}}}_{K}\left(t+{\rm{\Delta }}t\right) & = & {{\bf{B}}}_{K}\left(t\right),\\ {{\bf{x}}}_{{sp}}\left(t+{\rm{\Delta }}t\right) & = & {{\bf{x}}}_{{sp}}\left(t\right),\\ {\dot{{\bf{x}}}}_{{sp}}\left(t+{\rm{\Delta }}t\right) & = & {\dot{{\bf{x}}}}_{{sp}}\left(t\right).\end{array}\right.\end{eqnarray} \tag{ 57 }$$

For H₁, the dynamic equation is $\dot{F}=\{F,{H}_{1}\}$, or more specifically,

$$\begin{eqnarray}\left\{\begin{array}{rcl}{\dot{{\bf{E}}}}_{J} & = & -\displaystyle \sum _{{sp}}\displaystyle \frac{{q}_{s}}{{h}_{\sigma 3}\left({{\bf{x}}}_{J}\right)}\mathrm{diag}\left({{\bf{h}}}_{\sigma 1}{\left({{\bf{x}}}_{J}\right)}^{2}\right)\cdot {\dot{x}}_{1,{sp}}{{\bf{e}}}_{1}{W}_{{\sigma }_{1J}}\left({{\bf{x}}}_{{sp}}\right),\\ {\dot{{\bf{B}}}}_{K} & = & 0,\\ {\dot{{\bf{x}}}}_{{sp}} & = & {\dot{x}}_{1,{sp}}{{\bf{e}}}_{1},\\ {\ddot{{\bf{x}}}}_{{sp}} & = & -\mathrm{diag}\left(\displaystyle \frac{1}{2{{\bf{h}}}_{\sigma 1}{\left({{\bf{x}}}_{{sp}}\right)}^{2}}\right)\cdot {{\rm{\nabla }}}_{{{\bf{x}}}_{{sp}}}{h}_{1}{\left({{\bf{x}}}_{{sp}}\right)}^{2}{\dot{x}}_{1,{sp}}^{2}\\ & & +\mathrm{diag}\left(\displaystyle \frac{1}{{{\bf{h}}}_{\sigma 1}{\left({{\bf{x}}}_{{sp}}\right)}^{2}}\right){\dot{x}}_{1,{sp}}{{\bf{e}}}_{1}\times \left({{\rm{\nabla }}}_{{{\bf{x}}}_{{sp}}}\times \left(\mathrm{diag}\left({{\bf{h}}}_{\sigma 1}{\left({{\bf{x}}}_{{sp}}\right)}^{2}\cdot {\dot{{\bf{x}}}}_{{sp}}\right)\right)\right)\\ & & +\mathrm{diag}\left(\displaystyle \frac{{q}_{s}}{{{\bf{h}}}_{\sigma 1}{\left({{\bf{x}}}_{{sp}}\right)}^{2}{m}_{s}}\right){\dot{x}}_{1,{sp}}{{\bf{e}}}_{1}\times \displaystyle \sum _{K}{W}_{{\sigma }_{2K}}\left({{\bf{x}}}_{s}\right){{\bf{B}}}_{K}.\end{array}\right.\end{eqnarray} \tag{ 58 }$$

Because the equation for x_1,sp contains both ${\dot{x}}_{1,{sp}}$ and ${\ddot{x}}_{1,{sp}}$ explicitly, equation (58) is difficult to solve in general. However, when $\partial {h}_{1}\left({\bf{x}}\right)/\partial {x}_{1}=0$, i.e.

$$\begin{eqnarray}&&{h}_{1}\left({\bf{x}}\right)={h}_{1}\left({x}_{2},{x}_{3}\right),\end{eqnarray} \tag{ 59 }$$

the dynamics equation for particles in equation (58) reduces to

$$\begin{eqnarray}\left\{\begin{array}{rcl}{\dot{x}}_{1,{sp}} & = & {\dot{x}}_{1,{sp}},\\ {\dot{x}}_{2,{sp}} & = & 0,\\ {\dot{x}}_{3,{sp}} & = & 0,\\ {\ddot{x}}_{1,{sp}}{h}_{1}{\left({{\bf{x}}}_{{sp}}\right)}^{2} & = & 0,\\ \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\left({\dot{x}}_{2,{sp}}{h}_{2}{\left({{\bf{x}}}_{{sp}}\right)}^{2}\right) & = & 2{\dot{x}}_{1,{sp}}^{2}{h}_{1}\left({{\bf{x}}}_{{sp}}\right)\displaystyle \frac{\partial {h}_{1}\left({{\bf{x}}}_{{sp}}\right)}{\partial {x}_{2,{sp}}}-\displaystyle \frac{{q}_{s}}{{m}_{s}}{\dot{x}}_{1,{sp}}\displaystyle \sum _{K}{W}_{{\sigma }_{2K},{x}_{3}}\left({{\bf{x}}}_{{sp}}\right){B}_{{x}_{3},K}\\ & = & {\tilde{F}}_{{x}_{2}}\left({{\bf{x}}}_{{sp}},{\dot{x}}_{1,{sp}}\right),\\ \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\left({\dot{x}}_{3,{sp}}{h}_{3}{\left({{\bf{x}}}_{{sp}}\right)}^{2}\right) & = & 2{\dot{x}}_{1,{sp}}^{2}{h}_{1}\left({{\bf{x}}}_{{sp}}\right)\displaystyle \frac{\partial {h}_{1}\left({{\bf{x}}}_{{sp}}\right)}{\partial {x}_{3,{sp}}}+\displaystyle \frac{{q}_{s}}{{m}_{s}}{\dot{x}}_{1,{sp}}\displaystyle \sum _{K}{W}_{{\sigma }_{2K},{x}_{2}}\left({{\bf{x}}}_{{sp}}\right){B}_{{x}_{2},K}\\ & = & {\tilde{F}}_{{x}_{3}}\left({{\bf{x}}}_{{sp}},{\dot{x}}_{1,{sp}}\right).\end{array}\right.\end{eqnarray} \tag{ 60 }$$

In this case, equation (58) admits an analytical solution,

$$\begin{eqnarray}{{\rm{\Theta }}}_{1}\left({\rm{\Delta }}t\right):\left\{\begin{array}{rcl}{{\bf{E}}}_{J}\left(t+{\rm{\Delta }}t\right) & = & -\displaystyle \sum _{{sp}}\displaystyle \frac{{q}_{s}}{{h}_{\sigma 3}\left({{\bf{x}}}_{J}\right)}\mathrm{diag}\left({{\bf{h}}}_{\sigma 1}{\left({{\bf{x}}}_{J}\right)}^{2}\right)\cdot {\displaystyle \int }_{0}^{{\rm{\Delta }}t}{\rm{d}}t^{\prime} {\dot{x}}_{1,{sp}}{{\bf{e}}}_{1}{W}_{{\sigma }_{1J}}\left({{\bf{x}}}_{{sp}}+{{\bf{e}}}_{1}{\dot{x}}_{1,{sp}}t^{\prime} \right),\\ {{\bf{B}}}_{K}\left(t+{\rm{\Delta }}t\right) & = & {{\bf{B}}}_{K}\left(t\right),\\ {x}_{1,{sp}}\left(t+{\rm{\Delta }}t\right) & = & {x}_{1,{sp}}\left(t\right)+{\dot{x}}_{1,{sp}}\left(t\right){\rm{\Delta }}t,\\ {x}_{2,{sp}}\left(t+{\rm{\Delta }}t\right) & = & {x}_{2,{sp}}\left(t\right),\\ {x}_{3,{sp}}\left(t+{\rm{\Delta }}t\right) & = & {x}_{3,{sp}}\left(t\right),\\ {\dot{x}}_{1,{sp}}\left(t+{\rm{\Delta }}t\right) & = & {\dot{x}}_{1,{sp}}\left(t\right),\\ {\dot{x}}_{2,{sp}}\left(t+{\rm{\Delta }}t\right) & = & \displaystyle \frac{{h}_{2}{\left({{\bf{x}}}_{{sp}}\left(t\right)\right)}^{2}}{{h}_{2}{\left({{\bf{x}}}_{{sp}}\left(t\right)+{{\bf{e}}}_{1}{\dot{x}}_{1,{sp}}{\rm{\Delta }}t\right)}^{2}}{\dot{x}}_{2,{sp}}\left(t\right)+\\ & & \displaystyle \frac{1}{{h}_{2}{\left({{\bf{x}}}_{{sp}}\left(t\right)+{{\bf{e}}}_{1}{\dot{x}}_{1,{sp}}{\rm{\Delta }}t\right)}^{2}}{\displaystyle \int }_{0}^{{\rm{\Delta }}t}{\rm{d}}t^{\prime} {\tilde{F}}_{{x}_{2}}\left({{\bf{x}}}_{{sp}}+{{\bf{e}}}_{1}{\dot{x}}_{1,{sp}}t^{\prime} ,{\dot{x}}_{1,{sp}}\right),\\ {\dot{x}}_{3,{sp}}\left(t+{\rm{\Delta }}t\right) & = & \displaystyle \frac{{h}_{3}{\left({{\bf{x}}}_{{sp}}\left(t\right)\right)}^{2}}{{h}_{3}{\left({{\bf{x}}}_{{sp}}\left(t\right)+{{\bf{e}}}_{1}{\dot{x}}_{1,{sp}}{\rm{\Delta }}t\right)}^{2}}{\dot{x}}_{3,{sp}}\left(t\right)+\\ & & \displaystyle \frac{1}{{h}_{3}{\left({{\bf{x}}}_{{sp}}+{{\bf{e}}}_{1}{\dot{x}}_{1,{sp}}{\rm{\Delta }}t\right)}^{2}}{\displaystyle \int }_{0}^{{\rm{\Delta }}t}{\rm{d}}t^{\prime} {\tilde{F}}_{{x}_{3}}\left({{\bf{x}}}_{{sp}}+{{\bf{e}}}_{1}{\dot{x}}_{1,{sp}}t^{\prime} ,{\dot{x}}_{1,{sp}}\right).\end{array}\right.\end{eqnarray} \tag{ 61 }$$

Similarly, analytical solutions Θ₂ (Θ₃) for H₂ (H₃) can be also derived when $\partial {h}_{2}\left({\bf{x}}\right)/\partial {x}_{2}=\partial {h}_{3}\left({\bf{x}}\right)/\partial {x}_{3}=0$. Finally, we can compose these analytical solutions to obtain explicit symplectic integration algorithms for the entire system. For example, a first order scheme can be constructed as

$$\begin{eqnarray}&&{{\rm{\Theta }}}_{1}\left({\rm{\Delta }}t\right)={{\rm{\Theta }}}_{E}\left({\rm{\Delta }}t\right){{\rm{\Theta }}}_{B}\left({\rm{\Delta }}t\right){{\rm{\Theta }}}_{x}\left({\rm{\Delta }}t\right){{\rm{\Theta }}}_{y}\left({\rm{\Delta }}t\right){{\rm{\Theta }}}_{z}\left({\rm{\Delta }}t\right),\end{eqnarray} \tag{ 62 }$$

and a 2nd-order symmetric scheme is

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Theta }}}_{2}\left({\rm{\Delta }}t\right) & = & {{\rm{\Theta }}}_{x}\left({\rm{\Delta }}t/2\right){{\rm{\Theta }}}_{y}\left({\rm{\Delta }}t/2\right){{\rm{\Theta }}}_{z}\left({\rm{\Delta }}t/2\right){{\rm{\Theta }}}_{B}\left({\rm{\Delta }}t/2\right){{\rm{\Theta }}}_{E}\left({\rm{\Delta }}t\right)\\ & & \times {{\rm{\Theta }}}_{B}\left({\rm{\Delta }}t/2\right){{\rm{\Theta }}}_{z}\left({\rm{\Delta }}t/2\right){{\rm{\Theta }}}_{y}\left({\rm{\Delta }}t/2\right){{\rm{\Theta }}}_{x}\left({\rm{\Delta }}t/2\right).\end{array}\end{eqnarray} \tag{ 63 }$$

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

$$\begin{eqnarray}&&{{\rm{\Theta }}}_{2(l+1)}({\rm{\Delta }}t)={{\rm{\Theta }}}_{2l}({\alpha }_{l}{\rm{\Delta }}t){{\rm{\Theta }}}_{2l}({\beta }_{l}{\rm{\Delta }}t){{\rm{\Theta }}}_{2l}({\alpha }_{l}{\rm{\Delta }}t),\end{eqnarray} \tag{ 64 }$$

$$\begin{eqnarray}&&{\alpha }_{l}=1/(2-{2}^{1/(2l+1)}),\end{eqnarray} \tag{ 65 }$$

$$\begin{eqnarray}&&{\beta }_{l}=1-2{\alpha }_{l}.\end{eqnarray} \tag{ 66 }$$

Magnetic fusion plasmas are often confined in the toroidal geometry, for which the cylindrical coordinate system is convenient. We now present the high-order explicit structure-preserving geometric PIC algorithm in a cylindrical mesh. In this coordinate system, the line element is

$$\begin{eqnarray}&&{\rm{d}}{s}^{2}={\left({\rm{d}}r\right)}^{2}+{\left(\displaystyle \frac{r+{R}_{0}}{{R}_{0}}{\rm{d}}({R}_{0}\xi )\right)}^{2}+{\left({\rm{d}}z\right)}^{2},\end{eqnarray} \tag{ 67 }$$

where R₀ is a fixed radial length, r + R₀ is the radius in the standard cylindrical coordinate system, and R₀ ξ is the polar angle coordinate normalized by 1/R₀. For typical applications in tokamak physics, R₀ is the major radius and ξ is called toroidal angle. To simplify the notation, we also refer to (r, R₀ ξ, z) as (x, y, z).

For non-relativistic particles in the cylindrical mesh, if the discrete velocity in equation (22) is changed to

$$\begin{eqnarray}&&{{\bf{v}}}_{{sp},l}=\displaystyle \frac{{{\bf{x}}}_{{sp},l+1}-{{\bf{x}}}_{{sp},l}}{{\rm{\Delta }}t}{{\bf{h}}}_{\sigma 1}\left({{\bf{x}}}_{{sp},l+1}\right),\end{eqnarray} \tag{ 68 }$$

then the 1st-order scheme given by equation (28) will be explicit. To construct an explicit 2nd-order scheme, the 2nd-order action integral can be chosen as

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal A }}_{{\rm{d}}2} & = & \displaystyle \sum _{s,p,l}\displaystyle \frac{1}{2}\left({L}_{{sp}}\left({m}_{s},{{\bf{v}}}_{{sp},2l}\right)+{L}_{{sp}}\left({m}_{s},{{\bf{v}}}_{{sp},2l+1}^{* }\right)\right.\\ & & +{q}_{s}\left(2\overline{{{\bf{v}}}_{{sp}}\cdot {\bf{A}}}\left({{\bf{x}}}_{{sp},2l},{{\bf{x}}}_{{sp},2l+1},l\right)+2\overline{{{\bf{v}}}_{{sp}}^{* }\cdot {\bf{A}}}\left({{\bf{x}}}_{{sp},2l+1},{{\bf{x}}}_{{sp},2l+2},l\right)\right.\,\left.\left.-2\phi \left({{\bf{x}}}_{{sp},2l+2},l+1\right)\right)\right)\\ & & +\displaystyle \frac{1}{2}\displaystyle \sum _{J}\left({\left(\left(-\displaystyle \frac{{{\bf{A}}}_{J,l+1}-{{\bf{A}}}_{J,l}}{{\rm{\Delta }}t}-{\left({{\rm{\nabla }}}_{{\rm{d}}}\phi \right)}_{J,l+1}\right)\displaystyle \frac{1}{{{\bf{h}}}_{\sigma 1}\left({{\bf{x}}}_{J}\right)}\right)}^{2}\right.-\left.{\left(\displaystyle \frac{{\left({{\rm{\nabla }}}_{{\rm{d}}}\times {\bf{A}}\right)}_{K,l}}{{{\bf{h}}}_{\sigma 2}\left({{\bf{x}}}_{K}\right)}\right)}^{2}\right){h}_{\sigma 3}\left({{\bf{x}}}_{J}\right),\end{array}\end{eqnarray} \tag{ 69 }$$

where

$$\begin{eqnarray*}\begin{array}{rcl} & & {{\bf{v}}}_{{sp},2l}=2\displaystyle \frac{{{\bf{x}}}_{2l+1}-{{\bf{x}}}_{2l}}{{\rm{\Delta }}t}{{\bf{h}}}_{\sigma 1}\left({{\bf{x}}}_{{sp},2l+1}\right),\\ & & {{\bf{v}}}_{{sp},2l+1}^{* }=2\displaystyle \frac{{{\bf{x}}}_{2l+2}-{{\bf{x}}}_{2l+1}}{{\rm{\Delta }}t}{{\bf{h}}}_{\sigma 1}\left({{\bf{x}}}_{{sp},2l+1}\right),\\ & & \overline{{{\bf{v}}}_{{sp}}^{* }\cdot {\bf{A}}}\left({{\bf{x}}}_{{sp},2l+1},{{\bf{x}}}_{{sp},2l+2},l\right)={\displaystyle \int }_{0}^{1}{\rm{d}}\tau \\ & & \left[\displaystyle \frac{{x}_{1,{sp},2l+2}-{x}_{1,{sp},2l+1}}{{\rm{\Delta }}t},\displaystyle \frac{{x}_{2,{sp},2l+2}-{x}_{2,{sp},2l+1}}{{\rm{\Delta }}t},\displaystyle \frac{{x}_{3,{sp},2l+2}-{x}_{3,{sp},2l+1}}{{\rm{\Delta }}t}\right]\\ & & \cdot \left[\begin{array}{l}{A}_{{x}_{1},l}\left({x}_{1,{sp},2l+1}+\tau \left({x}_{1,{sp},2l+2}-{x}_{1,{sp},2l+1}\right),{x}_{2,{sp},2l+2},{x}_{3,{sp},2l+2}\right),\\ {A}_{{x}_{2},l}\left({x}_{1,{sp},2l+1},{x}_{2,{sp},2l+1}+\tau \left({x}_{2,{sp},2l+2}-{x}_{2,{sp},2l+1}\right),{x}_{3,{sp},2l+2}\right),\\ {A}_{{x}_{3},l}\left({x}_{1,{sp},2l+1},{x}_{2,{sp},2l+1},{x}_{3,{sp},2l+1}+\tau \left({x}_{3,{sp},2l+2}-{x}_{3,{sp},2l+1}\right)\right)\end{array}\right].\end{array}\end{eqnarray*}$$

Taking the following discrete variation yields discrete time-advance rules,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{ \mathcal A }}_{{\rm{d}}2}}{\partial {{\bf{x}}}_{2l}}=0,\end{eqnarray} \tag{ 70 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{ \mathcal A }}_{{\rm{d}}2}}{\partial {{\bf{x}}}_{2l+1}}=0,\end{eqnarray} \tag{ 71 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{ \mathcal A }}_{{\rm{d}}2}}{\partial {{\bf{A}}}_{J,l}}=0,\end{eqnarray} \tag{ 72 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{ \mathcal A }}_{{\rm{d}}2}}{\partial {\phi }_{J,l}}=0.\end{eqnarray} \tag{ 73 }$$

The explicit expressions for these time-advance are listed in appendix B. Higher-order discrete action integral for building explicit schemes can be also derived using a similar technique, or using the splitting method described in section 2.4.

In addition, we have also implemented the 1st-order relativistic charge-conserving geometric PIC in the cylindrical mesh. However the particle pusher in the relativistic scheme is implicit, and it is about 4 times slower than the explicit scheme. When studying relativistic effects, the implicit relativistic algorithm can be applied without significantly increasing the computational cost.

Kinetic equilibrium is the starting point for analytical and numerical studies of kinetic instabilities and associated transport phenomena. Because no self-consistent kinetic equilibrium is known in the tokamak geometry, most previous studies adopted non-self-consistent distributions as assumed kinetic equilibria, especially for simulations based on the δ f-method. Here, we numerically obtain an axisymmetric self-consistent kinetic steady state for a small tokamak using parameters similar to those of the Alcator C-Mod tokamak [49, 50]. The machine parameters are tabulated in table 1. The device is numerically constructed with poloidal field (PF) coils displayed in figure 1(a).

Zoom In Zoom Out Reset image size

At the beginning of the simulation, non-equilibrium distributions for deuterium ions and electrons are loaded into the device. The initial density profile in the poloidal plane is shown in figure 1(b). The 2D profiles of the external field and initial velocity and temperature profiles for both species are plotted in figure 2. Simulation parameters are chosen as

$$\begin{eqnarray*}\begin{array}{rcl}{\rm{\Delta }}x & = & 8\times {10}^{-4}\,{\rm{m}},\,{\rm{\Delta }}t=0.1{\rm{\Delta }}x/{\rm{c}},\,{R}_{0}=500{\rm{\Delta }}x,\\ {T}_{{\rm{e}},0} & = & {T}_{{\rm{i}},0}={T}_{0}=800\,\mathrm{eV},\,{n}_{{\rm{e}},0}={n}_{{\rm{i}},0}=5.0\times {10}^{19}\,{{\rm{m}}}^{-3},\\ {m}_{{\rm{i}}} & = & 3672{m}_{{\rm{e}}}=3.34\times {10}^{-27}\mathrm{kg},\\ {J}_{\mathrm{PF},i} & = & 67.9\,\mathrm{kA},\quad \mathrm{for}\,1\leqslant i\leqslant 6,\\ {J}_{\mathrm{PF},7} & = & {J}_{\mathrm{PF},8}=211.1\,\mathrm{kA},\,{J}_{\mathrm{PF},9}={J}_{\mathrm{PF},10}=238.1\,\mathrm{kA},\\ {J}_{\mathrm{PF},11} & = & -189.9\,\mathrm{kA},\,{J}_{\mathrm{PF},12}=-211.1\,\mathrm{kA}.\end{array}\end{eqnarray*}$$

Detailed calculation of the external magnetic field, initial particle distributions, and boundary setup are outlined in appendix A. Simulations show that the system reaches a steady state after 1.6 × 10⁷ time-steps. At this time, the amplitude of magnetic perturbation at the middle plane of the edge (x = 0.48 m and z = 0.41 m) is smaller than B₀/200, and the oscillation of position of plasma core is very small. Therefore, we can treat this state as a steady state. One such calculation requires about 1.3 × 10⁵ core-hours on the Tianhe 3 prototype cluster. Profiles of flow velocities, densities, temperatures of electrons and ions as well as electromagnetic field at the steady state are shown in figures 3–6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From figure 3, the flow velocity of ions at the steady state is in the range of 10 km s⁻¹, which is consistent with experimental observations [51, 52] and theoretical calculations [53, 54]. This ion flow is much slower than the thermal velocity of ions and thus is negligible in the force balance for the steady state, which can be written as

$$\begin{eqnarray}&&{\bf{j}}\times {\bf{B}}-{\rm{\nabla }}\cdot {\bf{p}}=0.\end{eqnarray} \tag{ 74 }$$

Here, p is the pressure tensor. Equation (74) is obtained by the familiar procedure of taking the second moment of the Vlasov equation, subtracting the flow velocity and summing over species. From the simulation data, the p_ij component of the pressure tensor is calculated as

$$\begin{eqnarray}&&{p}_{{ij},J}=\displaystyle \sum _{{{\bf{x}}}_{{sp}}\in \mathrm{grid}J}{m}_{s}v{{\prime} }_{i,{sp}}v{{\prime} }_{j,{sp}},\quad {\rm{for}}\,i,j\,\mathrm{in}\left\{x,y,z\right\},\end{eqnarray} \tag{ 75 }$$

$$\begin{eqnarray}&&{\bf{v}}{{\prime} }_{{sp}}\equiv {{\bf{v}}}_{{sp}}-\left(\displaystyle \sum _{{{\bf{x}}}_{{sp}}\in \mathrm{grid}J}{{\bf{v}}}_{{sp}}\right)/\left(\displaystyle \sum _{{{\bf{x}}}_{{sp}}\in \mathrm{grid}J}1\right),\end{eqnarray} \tag{ 76 }$$

where ${\sum }_{{{\bf{x}}}_{{sp}}\in \mathrm{grid}J}$ means that the summary is over all particles satisfying iΔx < x_sp ≤ (i + 1)Δx, jΔy < y_sp ≤ (j + 1)Δy and kΔz < z_sp ≤ (k + 1)Δz. The profile of pressure tensor at the steady state at z = 0.41 m is shown in figure 7. Clearly, the pressure tensor is predominately diagonal and anisotropic with p_yy > p_xx ≈ p_zz . Note that the pressure is almost isotropic in the poloidal plane, which indicates that it is valid to adopt the ideal magnetohydrodynamics (MHD) model with a scalar pressure for force balance in the 2D tokamak equilibrium. However, for 3D physics, the effect of pressure anisotropy needs to be considered. Since the steady state is 2D in space for the present case, the force balance equation reduces to

$$\begin{eqnarray}&&\displaystyle \frac{\partial {p}_{{xx}}}{\partial x}={\left({\bf{j}}\times {\bf{B}}\right)}_{x},\end{eqnarray} \tag{ 77 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {p}_{{zz}}}{\partial z}={\left({\bf{j}}\times {\bf{B}}\right)}_{z}.\end{eqnarray} \tag{ 78 }$$

To verify the force balance of the steady state, the four terms in equations (77) and (78) are plotted in figure 8. For comparison, the components of ∂p_yy /∂x and ∂p_yy /∂z are also plotted. Figure 8 shows that the force balance is approximately satisfied for the numerically obtained steady state, which can be viewed as a self-consistent kinetic steady state.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To verify the energy conservation in the simulation, we have recorded the time-history of the total energy in figure 9. The total energy drops a little, because some particles outside the last closed magnetic surface hit the simulation boundary, and these particles are removed from the simulation.

Zoom In Zoom Out Reset image size

Before reaching the kinetic steady state, the plasma oscillates in the poloidal plane. It is expected that this oscillation can be described as an MHD process whose characteristic velocity is the Alfvn velocity ${v}_{{\rm{A}}}={B}_{0}/\sqrt{{\mu }_{0}{n}_{{\rm{i}}}{m}_{{\rm{i}}}}$. To observe this oscillation, we plot in figure 10 the evolution of the magnetic field at x = 600Δx and z = 512Δx. From the parameters of the simulation, the characteristic frequency of the oscillation is

$$\begin{eqnarray}&&{\omega }_{{\rm{A}}}=\displaystyle \frac{{v}_{{\rm{A}}}}{{{qR}}_{0}}\sim 2.26\times {10}^{-2}{\omega }_{{\rm{c}},{\rm{i}}},\end{eqnarray} \tag{ 79 }$$

which agrees with the frequency of the B_x oscillation in figure 10.

Zoom In Zoom Out Reset image size

KBM [119], characterized by both electromagnetic perturbations of the MHD type and nontrivial kinetic effects, plays an important role in tokamak edge physics. Traditionally, it has been simulated using electromagnetic gyrokinetic codes such as the Kin-2DEM [55, 56], LIGKA [120], GTC [37, 57] and GEM [42]. However, for edge plasmas, the gyrokinetic ordering [56, 121–124] may become invalid under certain parameter regimes for modern tokamaks. For instance, the characteristic length in the edge of the H-mode plasma simulated by Wan et al [133, 134] can be as short as about 5 times of the gyroradius of thermal ions, and in this situation, the gyrokinetic density calculation may be inaccurate. We have applied SymPIC code developed to carry out the first whole-device 6D kinetic simulations of the KBM in a tokamak geometry.

The machine parameters are the same as in section 3.1. To trigger the KBM instability, we increase the plasma density to n₀ = 1 × 10²⁰ m⁻³, and the rest of parameters are

$$\begin{eqnarray*}\begin{array}{rcl}{I}_{{\rm{p}}} & = & 0.858\,\mathrm{MA},\,{T}_{{\rm{e}},0}={T}_{{\rm{i}},0}={T}_{0}=800\,\mathrm{eV},\\ {\rm{\Delta }}x & = & 3.2\times {10}^{-3}\,{\rm{m}},\,{\rm{\Delta }}t=0.5{\rm{\Delta }}x/{c}_{{\rm{n}}},\,{R}_{0}=125{\rm{\Delta }}x,\\ {c}_{{\rm{n}}} & = & {r}_{c}{c},\,{m}_{{\rm{i}}}={r}_{{\rm{m}}}{m}_{{\rm{e}}}=3.34\times {10}^{-27}\,\mathrm{kg},\\ {J}_{\mathrm{PF},i} & = & 144.6\,\mathrm{kA},\quad \mathrm{for}\,1\leqslant i\leqslant 6,\\ {J}_{\mathrm{PF},7} & = & {J}_{\mathrm{PF},8}=448.3\,\mathrm{kA},\,{J}_{\mathrm{PF},9}={J}_{\mathrm{PF},10}=506.2\,\mathrm{kA},\\ {J}_{\mathrm{PF},11} & = & -367.6\,\mathrm{kA},\,{J}_{\mathrm{PF},12}=-448.3\,\mathrm{kA}.\end{array}\end{eqnarray*}$$

Here r_m is the mass ratio between the deuterium and electron, c_n is the speed of light in the simulation and r_c is the ratio between c_n and the real speed of light in the vacuum c. For real plasmas, r_m ≈ 3672 and r_c = 1. Limited by available computation power, we reduce r_m and r_c in some of the simulations. Such an approximation is valid because the low frequency ion motion is relatively independent from the mass of electron and the speed of light, as long as r_m ≫ 1 and ${\left(c{r}_{{\rm{c}}}/{v}_{{\rm{A}}}\right)}^{2}\gg 1.$ In the present work, we take r_m = 100 and r_c = 0.16 to obtain long-term simulation results. Short-term results for r_m = 300, r_c = 0.5 and r_m = 3672, r_c = 1 (in this case Δt = 0.1Δx/c) are also obtained for comparison. The simulation domain is a 192 × 64 × 256 mesh, where perfect electric conductor is assumed at the boundaries in the x- and z-directions, and the periodic boundary is selected in y direction.

Because of the steep pressure gradient in the edge of the plasma, the threshold β_crit for ballooning mode instability is low. An estimated scaling for β_crit is [135]

$$\begin{eqnarray*}&&{\beta }_{\mathrm{crit}}=0.6\hat{s}/\left(\displaystyle \frac{2{q}^{2}{R}_{0}}{{L}_{{\rm{p}}}}\right),\end{eqnarray*}$$

where $\hat{s}$ is the magnetic shear and L_p is the pressure scaling length. For our simulated plasma, β_crit ∼ 1 × 10⁻³ and β ≈ 3 × 10⁻³. We except to observe unstable KBM.

To obtain a self-consistent kinetic steady state, we first perform a simulation as described in section 3.1 and obtain the 2D kinetic steady state. The profiles of temperature, safety factor, number density, pressure, toroidal current and bulk velocity of this steady state at z = 0.41 m are shown in figure 11.

Zoom In Zoom Out Reset image size

A random perturbation is then added in as the initial condition of the 6D simulation. The total simulation time is ${t}_{{\rm{a}}}=1.5\times {10}^{6}{\rm{\Delta }}t\approx 10000{\omega }_{{\rm{c}},{\rm{i}}}^{-1}\approx 5\times {10}^{-5}$ s. For one such simulation it takes about 3 × 10⁵ core-hours on the Tianhe 3 prototype cluster. The resulting mode structures of ion density for toroidal mode number n = 1, 2, 3, 6, 10, 14, 18 are shown in figures 12, 13 and 14 for t ω_c,i = 637, 1241 and 9693. It is clear that the unstable modes are triggered at the edge of the plasma, and the ballooning structure can be observed for modes with large n. The growth rate as a function of n is plotted in figure 15. For comparison, the grow rate obtained using r_c = 0.5, r_m = 300 and r_c = 1, r_m = 3672 are also plotted. It is clear that the growth rate has little correlation with the reduction of r_c and r_m. Figure 15 shows that the growth rate increases with n, consistent with the early gyrokinetic simulation results obtained using the Kin-2DEM eigenvalue code [55, 56]. Because the number of grids in the toroidal direction is 64 and the width of interpolating function is 4 times the grid size, the results for modes with n > 16 may not be accurate. The results displayed here are thus preliminary. In the next step, we plan to perform a larger scale simulation with more realistic 2D equilibria to obtain improved accuracy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The time-history of the mode amplitude is shown in figure 16. The unstable mode saturates approximately at t ≈ 5a/v_t,i, and the saturation level is in the range of 2%. Recent nonlinear gyrokinetic simulation [57] suggested that the instability is saturated by the E × B zonal flow generated by the instability. To verify this mechanism in our 6D fully kinetic simulation, the toroidally averaged E, E × B velocity and the measured phase velocity of the n = 12 mode at z = 128Δx and t = 8a/v_t,i are compared in figure 17. The E × B velocity at the edge correlates strongly with the phase velocity of the perturbation in terms of amplitude and profile. As a result, the E × B flow for the background plasma generated by instability interacts coherently with the mode structure, significantly modifies the space-time structure of the perturbation relative to the background plasma and reduces the drive of the instability. For this case simulated, the nonlinear saturation mechanism agrees qualitatively with the nonlinear gyrokinetic simulation [57].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size