2.1. Simplicial mesh

In this subsection, we build the geometric electrostatic PIC algorithm on an unstructured simplicial mesh. The model treats ions as fully kinetic six-dimensional particles. The response of the electrons is adiabatic (Horton Reference Horton1999; Weiland Reference Weiland2012; Sturdevant Reference Sturdevant2016; Hu et al. Reference Hu, Miecnikowski, Chen and Parker2018; Miecnikowski et al. Reference Miecnikowski, Sturdevant, Chen and Parker2018), which, with the quasi-neutrality condition, leads to

(2.1)\begin{equation} -\frac{q_{i}}{q_{e}}n_{i}=n_{e0}\exp \left(-\frac{q_{e}\phi}{T_{e}}\right),\end{equation}

where $n_{i}$ is the ion density, $n_{e0}$ is the background electron density, $\phi$ is the electric potential, $T_{e}$ is the electron temperature and $q_{e}$ and $q_{i}$ are electron and ion charges, respectively. The action integral of the system is (Xiao & Qin Reference Xiao and Qin2019)

(2.2)\begin{gather} S(\boldsymbol{x},\phi) =\int \, \textrm{d} tL(\boldsymbol{x},\phi), \end{gather}

(2.3)\begin{gather} L =\int \, \textrm{d} \boldsymbol{x}\, \textrm{d} \boldsymbol{v}\,f_{i}(\boldsymbol{x}, \boldsymbol{v})\left[\frac{1}{2}m_{i}\dot{\boldsymbol{x}}^{2}+q_{i} \dot{\boldsymbol{x}}\cdot\boldsymbol{A}_{0}(\boldsymbol{x},t)\right] \nonumber\\ \quad -q_{i}[\phi(\boldsymbol{x},t)+\phi_{0}(\boldsymbol{x},t)] +n_{e0}T_{e}\exp \left[-\frac{q_{e}\phi(\boldsymbol{x},t)}{T_{e}}\right], \end{gather}

where $\boldsymbol{x}=\boldsymbol{x}(\boldsymbol{x}_{0},\boldsymbol{v}_{0},t)$, $f_{i}$ is the ion distribution function, $m_{i}$ is ion mass and $\boldsymbol {A}_{0}$ and $\phi _{0}$ are the external vector and scalar potentials. The dynamics of the system is governed by the Euler–Lagrange equations:

(2.4)\begin{gather} \frac{\delta S}{\delta\phi} =0, \end{gather}

(2.5)\begin{gather}\frac{\delta S}{\delta\boldsymbol{x}} =0. \end{gather}

Equation (2.4) links the electric potential $\phi$ and the ion charge density $\rho _{i}$:

(2.6)\begin{equation} \phi={-}\frac{T_{e}}{q_{e}}\log \left(-\frac{\rho_{i}}{q_{e}n_{e0}}\right),\end{equation}

where $\rho _{i}=\int \textrm {d} \boldsymbol {v}q_{i}\,f_{i}(\boldsymbol {x},\boldsymbol {v})$. Equation (2.6) recovers the electron adiabatic response and charge-neutrality condition in (2.1). Equation (2.5) gives the equation of motion for particles:

(2.7)\begin{equation} \ddot{\boldsymbol{x}}=\frac{q_{i}}{m_{i}} [\boldsymbol{E}_{0}(\boldsymbol{x},t)+\boldsymbol{E} (\boldsymbol{x},t)+\dot{\boldsymbol{x}}\times\boldsymbol{B}_{0}(\boldsymbol{x},t)], \end{equation}

where $\boldsymbol {E}_{0}(\boldsymbol {x},t)$ and $\boldsymbol {B}_{0}(\boldsymbol {x},t)$ are the external electromagnetic field, and $\boldsymbol {E}(\boldsymbol {x},t)=-\boldsymbol {\nabla } \phi (\boldsymbol {x},t)$ is the perturbed electrostatic field.

We use the method introduced in Squire et al. (Reference Squire, Qin and Tang2012a) and Xiao & Qin (Reference Xiao and Qin2019) to discretize the action integral by particles and Whitney forms (Whitney Reference Whitney1957) on a 2-D unstructured triangular mesh. The discrete action integral $S_{d}$ can be written as

(2.8)\begin{gather} S_{d}(\boldsymbol{x}_{p},\phi_{I}) =\int L_{d}(\boldsymbol{x}_{p},\phi_{I})\, \textrm{d} t, \end{gather}

(2.9)\begin{gather} L_{d}(\boldsymbol{x}_{p},\phi_{I}) =\sum_{p}\left[\frac{1}{2}m_{i}\dot{\boldsymbol{x}}_{p}^{2}+q_{i} \dot{\boldsymbol{x}}_{p}\cdot\boldsymbol{A}_{0}(\boldsymbol{x}_{p})-q_{i} \sum_{I}W_{\sigma_{0},I}(\boldsymbol{x}_{p})\phi_{I}\right. \nonumber\\ \quad \left.\vphantom{\frac{1}{2}}-q_{i}\phi_{0}(\boldsymbol{x}_{p})\right]+ \sum_{I}n_{e0,I}T_{e,I}\exp \left(-\frac{q_{e}\phi_{I}}{T_{eI}}\right), \end{gather}

where $I$ is the triangular vertex index, $\boldsymbol {x}_{p}$ is the particle position of the $p$th particle, $\phi _{I}$ is the electric potential defined on the triangular vertex and $W_{\sigma _{0},I}$ is the Whitney 0-form interpolating the value of $\phi$ in continuous space using $\phi _{I}$. The discrete action integral is the same as in Xiao & Qin (Reference Xiao and Qin2019), except that here it is on a 2-D unstructured triangular mesh.

Variations of $S_{d}$ with respect to $\phi _{I}$ and $\boldsymbol {x}_{p}$ lead to the equations of motion of the electrostatic system:

(2.10)\begin{gather} \frac{\delta S_{d}}{\delta\phi_{I}}=0, \end{gather}

(2.11)\begin{gather}\frac{\delta S_{d}}{\delta\boldsymbol{x}_{p}}=0. \end{gather}

Equation (2.10) gives

(2.12)\begin{equation} \phi_{I}={-}\frac{T_{e,I}}{q_{e}}\log \left(-\frac{\rho_{I}}{q_{e}n_{e0,I}}\right),\end{equation}

where

(2.13)\begin{equation} \rho_{I}=q_{i}W_{\sigma_{0},I}(\boldsymbol{x}_{p})\end{equation}

is the charge density on the triangular vertex. Equation (2.11) is the governing equation for ion dynamics:

(2.14)\begin{equation} \ddot{\boldsymbol{x}}_{p}=\frac{q_{i}}{m_{i}} \left[\boldsymbol{E}_{0}(\boldsymbol{x}_{p},t)+ \dot{\boldsymbol{x}}_{p}\times\boldsymbol{B}_{0}(\boldsymbol{x}_{p},t)- \frac{\partial}{\partial\boldsymbol{x}_{p}}\sum_{I}W_{\sigma_{0},I}(\boldsymbol{x}_{p}) \phi_{I}\right].\end{equation}

The last term of (2.14) is the derivative of $W_{\sigma _{0},I}$ with respect to $\boldsymbol {x}_{p}$ in continuous space. According to the property of Whitney forms (Whitney Reference Whitney1957; Hirani Reference Hirani2003; Desbrun, Kanso & Tong Reference Desbrun, Kanso and Tong2008; Squire et al. Reference Squire, Qin and Tang2012a; Xiao et al. Reference Xiao, Qin, Liu, He, Zhang and Sun2015b),

(2.15)\begin{equation} \boldsymbol{\nabla} \sum_{I}W_{\sigma_{0},I}(\boldsymbol{x}_{p})\phi_{I}=\sum_{J}W_{\sigma_{1},J} (\boldsymbol{x}_{p})\sum_{I}\nabla_{dJ,I}\phi_{I},\end{equation}

where $\sum _{I}\nabla _{dJ,I}\phi _{I}$ is the discrete gradient of $\phi _{I}$, and $W_{\sigma _{1},J}$ is the Whitney 1-form that interpolates a continuous 1-form from the discrete 1-form defined on the triangular edge. The construction of $W_{\sigma _{0},I}, W_{\sigma _{1},J}$ and $\nabla _{dJ,I}$ is further discussed below. With property (2.15), (2.14) becomes

(2.16)\begin{equation} \ddot{\boldsymbol{x}}_{p}=\frac{q_{i}}{m_{i}}\left[\boldsymbol{E}_{0} (\boldsymbol{x}_{p},t)+\dot{\boldsymbol{x}}_{p}\times\boldsymbol{B}_{0} (\boldsymbol{x}_{p},t)+\sum_{J}W_{\sigma_{1},J}(\boldsymbol{x}_{p}) \boldsymbol{E}_{J}\right]\end{equation}

and

(2.17)\begin{equation} \boldsymbol{E}_{J}={-}\sum_{I}\nabla_{dJ,I}\phi_{I},\end{equation}

where $\boldsymbol {E}_{J}$ is the discrete electrical field defined on the triangular edge labelled by $J$. We want to emphasize again that, similar to the scenario in a cubic mesh (Xiao & Qin Reference Xiao and Qin2019), without Whitney forms and discrete exterior calculus, it is difficult to calculate on the electric field on an unstructured mesh to advance particle positions.

As is well known, the key parts of a PIC method include charge deposition, solving discrete field and field interpolation, which are encapsulated in a systematic way in (2.13), (2.17) and (2.16), respectively. Once $W_{\sigma _{0},I}, W_{\sigma _{1},J}$ and $\nabla _{dJ,I}$ are chosen, the PIC algorithm is defined.

Now we describe in detail the construction of $W_{\sigma _{0},I}, W_{\sigma _{1},J}$ and $\nabla _{dJ,I}$ on the triangular mesh. Figure 1(a) shows a particle in a triangle, where $(x,y)$ is the position of the particle $p$, and the three vertices of the $i$th triangle, $i_{1}, i_{2}, i_{3}$, have coordinates $(x_{i_{1}},y_{i_{1}}), (x_{i_{2}},y_{i_{2}}), (x_{i_{3}},y_{i_{3}})$, respectively. To deposit charge at each vertex according to (2.13), we need to specify the Whitney 0-forms, which are chosen to be linear barycentric functions. For $(x,y)$ inside the triangle, the Whitney 0-forms are

(2.18)\begin{align} W_{\sigma_{0},i_{1}}(x,y) & =\frac{(y_{i_{2}}-y_{i_{3}})(x-x_{i_{3}}) +(x_{i_{3}}-x_{i_{2}})(y-y_{i_{3}})}{(x_{i_{1}}-x_{i_{3}})(y_{i_{2}} -y_{i_{3}})+(x_{i_{3}}-x_{i_{2}})(y_{i_{1}}-y_{i_{3}})}, \end{align}

(2.19)\begin{align} W_{\sigma_{0},i_{2}}(x,y) & =\frac{(y_{i_{3}}-y_{i_{1}})(x-x_{i_{3}}) +(x_{i_{1}}-x_{i_{3}})(y-y_{i_{3}})}{(x_{i_{1}}-x_{i_{3}})(y_{i_{2}} -y_{i_{3}})+(x_{i_{3}}-x_{i_{2}})(y_{i_{1}}-y_{i_{3}})}, \end{align}

(2.20)\begin{align} W_{\sigma_{0},i_{3}}(x,y) & =1-W_{\sigma_{0},i_{1}}(x,y) -W_{\sigma_{0},i_{2}}(x,y). \end{align}

When $(x,y)$ is outside the triangle, all the Whitney 0-forms vanish. Note that $W_{\sigma _{0},i_{1}}(x,y)$, $W_{\sigma _{0},i_{2}}(x,y)$ and $W_{\sigma _{0},i_{3}}(x,y)$ are the areas of the triangle $\Delta pi_{2}i_{3}, \Delta pi_{3}i_{1}$ and $\Delta pi_{1}i_{2}$, respectively. Thus, the weighting of charge deposition with respect to a vertex is the area weighting of the triangle enclosed by the particle and the opposite edge of the vertex. The density plots of $W_{\sigma _{0},i_{1}}(x,y)$, $W_{\sigma _{0},i_{2}}(x,y)$ and $W_{\sigma _{0},i_{3}}(x,y)$ are shown in figure 2(a–c). Figure 2(a) shows $W_{\sigma _{0},i_{1}}(x,y)$ approaching 1 as the particle is close to the vertex $i_{1}$. The chosen Whitney 0-forms in (2.18)–(2.20) obviously satisfy the condition

(2.21)\begin{equation} \sum_{I=i_{1},i_{2},i_{3}} W_{\sigma_{0},I}(\boldsymbol{x}_{p}) = 1. \end{equation}

In the 2-D triangular mesh, 1-forms, such as $\boldsymbol {E}_{J}$, are defined on the triangular edges. And the index $J$ for the edges consists of an ordered pair of indices of the vertices. For example, $J=i_{1}i_{2}$ labels the oriented edge from $i_{1}$ to $i_{2}$, as shown in figure 1(b). The discrete gradient operator $\nabla _{dJ,I}$ consistent with (2.15) is

(2.22)\begin{equation} \nabla_{di_{1}i_{2},I} =\delta_{i_{2}I}-\delta_{i_{1}I}. \end{equation}

Figure 1. The PIC algorithm on a triangle. (a) Depositing a particle's charge to the triangular vertices using Whitney 0-forms $W_{\sigma _{0},I}$. (b) Computing $\boldsymbol {E}_{J}$ on the edge with the discrete gradient operator $\nabla _{dJ,I}$. (c) Interpolating $\boldsymbol {E}_{J}$ from edges to the particle's location through Whitney 1-forms $W_{\sigma _{1},J}$.

Figure 2. The values of the Whitney 0-forms and 1-forms. The density plot of the value of (a) $W_{\sigma _{0},i_{1}}(x,y)$, (b) $W_{\sigma _{0},i_{2}}(x,y)$ and (c) $W_{\sigma _{0},i_{3}}(x,y)$. The amplitude density plot and the quiver plot of (d) $W_{\sigma _{1},i_{2}i_{3}}(x,y)$, (e) $W_{\sigma _{1},i_{3}i_{1}}(x,y)$ and (f) $W_{\sigma _{1},i_{1}i_{2}}(x,y)$.

According to the definition of Whitney forms (Whitney Reference Whitney1957), the Whitney 1-form on an unstructured triangular mesh is

(2.23)\begin{equation} W_{\sigma_{1},j^{\prime}j}=W_{\sigma_{0},j}\boldsymbol{\nabla} W_{\sigma_{0},j^{\prime}}-W_{\sigma_{0},j^{\prime}}\boldsymbol{\nabla} W_{\sigma_{0},j}. \end{equation}

For the triangular mesh, the expressions of $W_{\sigma _{1},j^{\prime }j}(x,y)$ for $j^{\prime },j=i_{1},i_{2},i_{3}$ and $j^{\prime }\neq j$ are

(2.24)\begin{align} W_{\sigma_{1},i_{1}i_{2}}(x,y) & =\left[\frac{(y_{i_{3}}-y_{i_{1}}) W_{\sigma_{0},i_{1}}(x,y)-(y_{i_{2}}-y_{i_{3}})W_{\sigma_{0},i_{2}} (x,y)}{(x_{i_{1}}-x_{i_{3}})(y_{i_{2}}-y_{i_{3}})+(x_{i_{3}} -x_{i_{2}})(y_{i_{1}}-y_{i_{3}})},\right. \nonumber\\ & \quad \left.\frac{(x_{i_{1}}-x_{i_{3}})W_{\sigma_{0},i_{1}}(x,y) -(x_{i_{3}}-x_{i_{2}})W_{\sigma_{0},i_{2}}(x,y)}{(x_{i_{1}} -x_{i_{3}})(y_{i_{2}}-y_{i_{3}})+(x_{i_{3}}-x_{i_{2}}) (y_{i_{1}}-y_{i_{3}})}\right], \end{align}

(2.25)\begin{align} W_{\sigma_{1},i_{2}i_{3}}(x,y) & =\left[\frac{(y_{i_{1}}-y_{i_{2}})W_{\sigma_{0}, i_{2}}(x,y)-(y_{i_{3}}-y_{i_{1}})W_{\sigma_{0},i_{3}}(x,y)}{(x_{i_{1}} -x_{i_{3}})(y_{i_{2}}-y_{i_{3}})+(x_{i_{3}}-x_{i_{2}})(y_{i_{1}}-y_{i_{3}})},\right. \nonumber\\ & \quad \left.\frac{(x_{i_{2}}-x_{i_{1}})W_{\sigma_{0},i_{2}}(x,y)-(x_{i_{1}} -x_{i_{3}})W_{\sigma_{0},i_{3}}(x,y)}{(x_{i_{1}}-x_{i_{3}})(y_{i_{2}} -y_{i_{3}})+(x_{i_{3}}-x_{i_{2}})(y_{i_{1}}-y_{i_{3}})}\right], \end{align}

(2.26)\begin{align} W_{\sigma_{1},i_{3}i_{1}}(x,y) & =\left[\frac{(y_{i_{2}}-y_{i_{3}}) W_{\sigma_{0},i_{3}}(x,y)-(y_{i_{1}}-y_{i_{2}})W_{\sigma_{0},i_{1}} (x,y)}{(x_{i_{1}}-x_{i_{3}})(y_{i_{2}}-y_{i_{3}})+(x_{i_{3}}-x_{i_{2}}) (y_{i_{1}}-y_{i_{3}})},\right. \nonumber\\ & \quad \left.\frac{(x_{i_{3}}-x_{i_{2}})W_{\sigma_{0},i_{3}}(x,y) -(x_{i_{2}}-x_{i_{1}})W_{\sigma_{0},i_{1}}(x,y)}{(x_{i_{1}} -x_{i_{3}})(y_{i_{2}}-y_{i_{3}})+(x_{i_{3}}-x_{i_{2}}) (y_{i_{1}}-y_{i_{3}})}\right], \end{align}

for $(x,y)$ inside the triangle. All $W_{\sigma _{1},j^{\prime }j}(x,y)$ vanish when $(x,y)$ is outside the triangle. The amplitude density plot and the quiver plot of $W_{\sigma _{1},i_{2}i_{3}}(x,y)$ are shown in figure 2(d), where the amplitude approaches zero when a particle is close to the opposite vertex of the edge and maximizes when the particle approaches the edge. Figures 2(e) and 2(f) plot the value and direction of $W_{\sigma _{1},i_{3}i_{1}}(x,y)$ and $W_{\sigma _{1},i_{1}i_{2}}(x,y)$.

After $W_{\sigma _{0},I}$, $W_{\sigma _{1},J}$ and $\nabla _{dJ,I}$ are chosen, the discrete electric potential $\phi _{I}$ at each vertex can be calculated by (2.13) and (2.12), and the electric field on the edges according to (2.17) is

(2.27)\begin{equation} E_{i_{1}i_{2}}=\phi_{i_{2}}-\phi_{i_{1}}, \end{equation}

where $E_{i_{1}i_{2}}$ is a discrete 1-form, denoted by a boldface symbol in figure 1 following the convention of physicists. Particles’ positions and velocities are advanced according to (2.16), which interpolates the electrical field at $\boldsymbol {x}_{p}$ as ${\boldsymbol {E}}(\boldsymbol {x}_{p})=\sum _{J}W_{\sigma _{1},J} (\boldsymbol {x}_{p})\boldsymbol {E}_{J}$ using Whitney 1-forms. This process is illustrated in figure 1(c). In the current implementation, (2.16) is integrated by the Boris algorithm (Boris Reference Boris1970), which preserves the phase space volume (Qin et al. Reference Qin, Zhang, Xiao, Liu, Sun and Tang2013; Ellison, Burby & Qin Reference Ellison, Burby and Qin2015; He et al. Reference He, Sun, Liu and Qin2015b, Reference He, Sun, Liu and Qin2016a; Zhang et al. Reference Zhang, Liu, Qin, Wang, He and Sun2015; He et al. Reference He, Sun, Zhang, Wang, Liu and Qin2016c) although is not symplectic. The algorithm on a three-dimensional tetrahedral mesh can also be constructed in the same way.

2.2. Three-dimensional non-simplicial mesh

In this subsection, we develop the algorithm in a special non-simplicial mesh made of prisms. The mesh is constructed from a 2-D unstructured simplicial mesh and a structured grid in the third perpendicular dimension. This type of unstructured mesh has been adopted by the XGC code (Chang et al. Reference Chang, Ku and Weitzner2004, Reference Chang, Ku, Diamond, Lin, Parker, Hahm and Samatova2009, Reference Chang, Ku, Tynan, Hager, Churchill, Cziegler, Greenwald, Hubbard and Hughes2017; Ku et al. Reference Ku, Chang, Adams, Cummings, Hinton, Keyes, Klasky, Lee, Lin and Parker2006, Reference Ku, Chang, Hager, Churchill, Tynan, Cziegler, Greenwald, Hughes, Parker, Adams, D'Azevedo and Worley2018), where the structured direction is the toroidal direction, and the poloidal plane is covered by a 2-D unstructured simplicial mesh. The discrete variational principle, (2.4), gives the electric potential on a prism:

(2.28)\begin{equation} \phi_{I}={-}\frac{T_{e,I}}{q_{e}}\log \left(-\frac{\rho_{I}}{q_{e}n_{e0,I}}\right),\end{equation}

where

(2.29)\begin{equation} \rho_{I}=q_{i}f_{I}(\boldsymbol{x}_{p}).\end{equation}

Equation (2.29) describes charge deposition on a prism with the function $f_{I}$ defined as

(2.30)\begin{equation} \sum_{I=i_{1},\ldots,i_{6}}f_{I}(\boldsymbol{x}_{p})=f_{i_{1}} (\boldsymbol{x}_{p})+f_{i_{2}}(\boldsymbol{x}_{p})+f_{i_{3}} (\boldsymbol{x}_{p})+f_{i_{4}}(\boldsymbol{x}_{p})+f_{i_{5}} (\boldsymbol{x}_{p})+f_{i_{6}}(\boldsymbol{x}_{p}),\end{equation}

where $i_{1},\ldots ,i_{6}$ label vertices 1 to 6 of the $i$th prism. Define $f_{i_{1}}(\boldsymbol {x}_{p})=f_{\varDelta _{1},i_{1}}(x,y)\,f_{\varDelta _{1}}(z)$, where $\varDelta _{1}$ is the back triangle of the prism, $(x,y)$ is the position of the particle in the plane containing the triangle, $f_{\varDelta _{1},i_{1}}(x,y)$ gives the weighting of $(x,y)$ with respect to $i_{1}$ on the triangle $\varDelta _{1}$ and $z$ is the particle's position along the direction perpendicular to the triangle plane. Similarly, we define $f_{i_{4}}(\boldsymbol {x}_{p})=f_{\varDelta _{2},i_{4}}(x,y)\,f_{\varDelta _{2}}(z)$, where $\varDelta _{2}$ is the front triangle of the prism. Functions $f_{i_{2}}(\boldsymbol {x}_{p}), f_{i_{3}}(\boldsymbol {x}_{p}), f_{i_{5}}(\boldsymbol {x}_{p})$ and $f_{i_{6}}(\boldsymbol {x}_{p})$ are defined in a similar way (see figure 3a). Since $\varDelta _{1}$ and $\varDelta _{2}$ are identical, the weightings of $(x,y)$ with respect to the corresponding vertices $i_{1}$ and $i_{4}$ are identical, i.e.

(2.31)\begin{equation} f_{\varDelta_{1},i_{1}}(x,y)=f_{\varDelta_{2},i_{4}}(x,y),\end{equation}

and

(2.32)\begin{gather} f_{\varDelta_{1},i_{2}}(x,y)=f_{\varDelta_{2},i_{5}}(x,y), \end{gather}

(2.33)\begin{gather}f_{\varDelta_{1},i_{3}}(x,y)=f_{\varDelta_{2},i_{6}}(x,y). \end{gather}

Hence, (2.30) can be written as

(2.34)\begin{align} &f_{i_{1}}(\boldsymbol{x}_{p})+f_{i_{2}}(\boldsymbol{x}_{p}) +f_{i_{3}}(\boldsymbol{x}_{p})+f_{i_{4}}(\boldsymbol{x}_{p}) +f_{i_{5}}(\boldsymbol{x}_{p})+f_{i_{6}}(\boldsymbol{x}_{p})\nonumber\\ &\quad =[\,f_{\varDelta_{1},i_{1}}(x,y)+f_{\varDelta_{1},i_{2}}(x,y) +f_{\varDelta_{1},i_{3}}(x,y)]f_{\varDelta_{1}}(z)\nonumber\\ &\qquad +[\,f_{\varDelta_{2},i_{4}}(x,y)+f_{\varDelta_{2},i_{5}}(x,y) +f_{\varDelta_{2},i_{6}}(x,y)]f_{\varDelta_{2}}(z)\nonumber\\ &\quad =[\,f_{\varDelta_{1},i_{1}}(x,y)+f_{\varDelta_{1},i_{2}}(x,y) +f_{\varDelta_{1},i_{3}}(x,y)](\,f_{\varDelta_{1}}(z)+f_{\varDelta_{2}}(z)). \end{align}

By choosing $f_{\varDelta _{1},i_{1}}(x,y)+f_{\varDelta _{1},i_{2}}(x,y)+f_{\varDelta _{1},i_{3}}(x,y)=1$ and $f_{\varDelta _{1}}(z)+f_{\varDelta _{2}}(z)=1$, we have ${\sum _{I-i_{1},\ldots ,i_{6}}}f(\boldsymbol {x}_{p})=1$. Equation (2.11) gives the governing equation for ion dynamics:

(2.35)\begin{equation} \ddot{\boldsymbol{x}}_{p}=\frac{q_{i}}{m_{i}}\left[\boldsymbol{E}_{0} (\boldsymbol{x}_{p},t)+\dot{\boldsymbol{x}}_{p}\times\boldsymbol{B}_{0} (\boldsymbol{x}_{p},t)-\frac{\partial}{\partial\boldsymbol{x}_{p}} \sum_{I}f_{I}(\boldsymbol{x}_{p})\phi_{I}\right].\end{equation}

The last term of (2.35) is

(2.36)\begin{align} -\frac{\partial}{\partial\boldsymbol{x}_{p}}\sum_{I=i_{1},\ldots,i_{6}}f_{I} (\boldsymbol{x}_{p})\phi_{I}={-}\phi_{i_{1}}\boldsymbol{\nabla}\,f_{i_{1}}-\phi_{i_{2}}\boldsymbol{\nabla}\,f_{i_{2}}-\phi_{i_{3}}\boldsymbol{\nabla}\,f_{i_{3}}-\phi_{i_{4}}\boldsymbol{\nabla}\,f_{i_{4}}-\phi_{i_{5}}\boldsymbol{\nabla}\,f_{i_{5}}-\phi_{i_{6}}\boldsymbol{\nabla}\,f_{i_{6}},\end{align}

which can be further expressed as

(2.37)\begin{align} &-\phi_{i_{1}}\boldsymbol{\nabla}\,f_{i_{1}}-\phi_{i_{2}}\boldsymbol{\nabla}\,f_{i_{2}}-\phi_{i_{3}}\boldsymbol{\nabla}\,f_{i_{3}}-\phi_{i_{4}}\boldsymbol{\nabla}\,f_{i_{4}}-\phi_{i_{5}}\boldsymbol{\nabla}\,f_{i_{5}}-\phi_{i_{6}}\boldsymbol{\nabla}\,f_{i_{6}}\nonumber\\ &\quad =f_{\varDelta_{1}}(z)[-\phi_{i_{1}}\nabla_{{\perp}}\,f_{\varDelta_{1},i_{1}} (x,y)-\phi_{i_{2}}\nabla_{{\perp}}\,f_{\varDelta_{1},i_{2}}(x,y)-\phi_{i_{3}} \nabla_{{\perp}}\,f_{\varDelta_{1},i_{3}}(x,y)]\nonumber\\ &\qquad +\,f_{\varDelta_{2}}(z)[-\phi_{i_{4}}\nabla_{{\perp}}\,f_{\varDelta_{2}, i_{4}}(x,y)-\phi_{i_{5}}\nabla_{{\perp}}\,f_{\varDelta_{2},i_{5}}(x,y)- \phi_{i_{6}}\nabla_{{\perp}}\,f_{\varDelta_{2},i_{6}}(x,y)]\nonumber\\ &\qquad -\phi_{i_{1}}f_{\varDelta_{1},i_{1}}(x,y)\frac{\partial}{\partial z}\,f_{\varDelta_{1}}(z)-\phi_{i_{4}}f_{\varDelta_{2},i_{4}}(x,y)\frac{\partial}{\partial z}\,f_{\varDelta_{2}}(z)\nonumber\\ &\qquad -\phi_{i_{2}}f_{\varDelta_{1},i_{2}}(x,y)\frac{\partial}{\partial z}\,f_{\varDelta_{1}}(z)-\phi_{i_{5}}f_{\varDelta_{2},i_{5}}(x,y)\frac{\partial}{\partial z}\,f_{\varDelta_{2}}(z)\nonumber\\ &\qquad -\phi_{i_{3}}f_{\varDelta_{1},i_{3}}(x,y)\frac{\partial}{\partial z}\,f_{\varDelta_{1}}(z)-\phi_{i_{6}}f_{\varDelta_{2},i_{6}}(x,y)\frac{\partial}{\partial z}\,f_{\varDelta_{2}}(z), \end{align}

where $\nabla _{\perp }=({{\partial }/{\partial x}+{\partial }/{\partial y}})$. The first term of the right-hand side of (2.37) can be expressed as

(2.38)\begin{align} &f_{\varDelta_{1}} (z)[-\phi_{i_{1}}\nabla_{{\perp}}\,f_{\varDelta_{1},i_{1}} (x,y)-\phi_{i_{2}}\nabla_{{\perp}}\,f_{\varDelta_{1},i_{2}}(x,y)- \phi_{i_{3}}\nabla_{{\perp}}\,f_{\varDelta_{1},i_{3}}(x,y)]\nonumber\\ &\quad =f_{\varDelta_{1}}(z)[(\phi_{i_{2}}-\phi_{i_{1}})(\,f_{\varDelta_{1}, i_{2}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{1},i_{1}}(x,y)-f_{\varDelta_{1}, i_{1}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{1},i_{2}}(x,y))\nonumber\\ &\qquad +(\phi_{i_{3}}-\phi_{i_{2}})(\,f_{\varDelta_{1},i_{3}}(x,y) \nabla_{{\perp}}\,f_{\varDelta_{1},i_{2}}(x,y)-f_{\varDelta_{1},i_{2}} (x,y)\nabla_{{\perp}}\,f_{\varDelta_{1},i_{3}}(x,y))\nonumber\\ &\qquad +(\phi_{i_{1}}-\phi_{i_{3}})(\,f_{\varDelta_{1},i_{1}}(x,y) \nabla_{{\perp}}\,f_{\varDelta_{1},i_{3}}(x,y)-f_{\varDelta_{1},i_{3}} (x,y)\nabla_{{\perp}}\,f_{\varDelta_{1},i_{1}}(x,y))]. \end{align}

Similarly, the second term of the right-hand side of (2.37) is

(2.39)\begin{align} &f_{\varDelta_{2}}(z)[-\phi_{i_{4}}\nabla_{{\perp}}\,f_{\varDelta_{2}, i_{4}}(x,y)-\phi_{i_{5}}\nabla_{{\perp}}\,f_{\varDelta_{2},i_{5}}(x,y) -\phi_{i_{6}}\nabla_{{\perp}}\,f_{\varDelta_{2},i_{6}}(x,y)]\nonumber\\ &\quad =f_{\varDelta_{2}}(z)[(\phi_{i_{5}}-\phi_{i_{4}})(\,f_{\varDelta_{2}, i_{5}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{2},i_{4}}(x,y)-f_{\varDelta_{2}, i_{4}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{2},i_{5}}(x,y))\nonumber\\ &\qquad +(\phi_{i_{6}}-\phi_{i_{5}})(\,f_{\varDelta_{2},i_{6}}(x,y) \nabla_{{\perp}}\,f_{\varDelta_{2},i_{5}}(x,y)-f_{\varDelta_{2},i_{5}}(x,y) \nabla_{{\perp}}\,f_{\varDelta_{2},i_{6}}(x,y))\nonumber\\ &\qquad +(\phi_{i_{4}}-\phi_{i_{6}})(\,f_{\varDelta_{2},i_{4}}(x,y) \nabla_{{\perp}}\,f_{\varDelta_{2},i_{6}}(x,y)-f_{\varDelta_{2},i_{6}}(x,y) \nabla_{{\perp}}\,f_{\varDelta_{2},i_{4}}(x,y))]. \end{align}

Based on (2.31), the third term of (2.37) can be written as

(2.40)\begin{align} &-\phi_{i_{1}}f_{\varDelta_{1},i_{1}}(x,y)\frac{\partial}{\partial z}\,f_{\varDelta_{1}}(z)-\phi_{i_{4}}f_{\varDelta_{2},i_{4}}(x,y) \frac{\partial}{\partial z}\,f_{\varDelta_{2}}(z)\nonumber\\ &\quad =f_{\varDelta_{1},i_{1}}(x,y)\left[-\phi_{i_{1}} \frac{\partial}{\partial z}\,f_{\varDelta_{1}}(z)-\phi_{i_{4}} \frac{\partial}{\partial z}\,f_{\varDelta_{2}}(z)\right]\nonumber\\ &\quad =(\phi_{i_{4}}-\phi_{i_{1}})\,f_{\varDelta_{1},i_{1}}(x,y) \left[\,f_{\varDelta_{2}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{1}} (z)-f_{\varDelta_{1}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{2}}(z)\right]. \end{align}

Similarly, the last two terms of (2.37) can be written as

(2.41)\begin{align} &-\phi_{i_{2}}f_{\varDelta_{1},i_{2}}(x,y)\frac{\partial}{\partial z}\,f_{\varDelta_{1}}(z) -\phi_{i_{5}}f_{\varDelta_{2},i_{5}}(x,y)\frac{\partial}{\partial z}\,f_{\varDelta_{2}}(z)\nonumber\\ &\quad =(\phi_{i_{5}}-\phi_{i_{2}})\,f_{\varDelta_{1},i_{2}}(x,y) \left[\,f_{\varDelta_{2}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{1}}(z) -f_{\varDelta_{1}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{2}}(z)\right] \end{align}

and

(2.42)\begin{align} &-\phi_{i_{3}}f_{\varDelta_{1},i_{3}}(x,y)\frac{\partial}{\partial z}\,f_{\varDelta_{1}}(z)-\phi_{i_{6}}f_{\varDelta_{2},i_{6}}(x,y) \frac{\partial}{\partial z}\,f_{\varDelta_{2}}(z)\nonumber\\ &\quad =(\phi_{i_{6}}-\phi_{i_{3}})\,f_{\varDelta_{1},i_{3}}(x,y) \left[\,f_{\varDelta_{2}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{1}}(z) -f_{\varDelta_{1}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{2}}(z)\right]. \end{align}

Finally, (2.36) is written as

(2.43)\begin{align} &-\phi_{i_{1}}\boldsymbol{\nabla}\,f_{i_{1}}-\phi_{i_{2}}\boldsymbol{\nabla}\,f_{i_{2}}-\phi_{i_{3}} \boldsymbol{\nabla}\,f_{i_{3}}-\phi_{i_{4}}\boldsymbol{\nabla}\,f_{i_{4}}-\phi_{i_{5}}\boldsymbol{\nabla}\,f_{i_{5}}-\phi_{i_{6}}\boldsymbol{\nabla}\,f_{i_{6}}\nonumber\\ &\quad =f_{\varDelta_{1}}(z)[(\phi_{i_{2}}-\phi_{i_{1}})(\,f_{\varDelta_{1}, i_{2}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{1},i_{1}}(x,y)-f_{\varDelta_{1}, i_{1}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{1},i_{2}}(x,y))\nonumber\\ &\qquad +(\phi_{i_{3}}-\phi_{i_{2}})(\,f_{\varDelta_{1},i_{3}}(x,y) \nabla_{{\perp}}\,f_{\varDelta_{1},i_{2}}(x,y)-f_{\varDelta_{1},i_{2}}(x,y) \nabla_{{\perp}}\,f_{\varDelta_{1},i_{3}}(x,y))\nonumber\\ &\qquad +(\phi_{i_{1}}-\phi_{i_{3}})(\,f_{\varDelta_{1},i_{1}}(x,y) \nabla_{{\perp}}\,f_{\varDelta_{1},i_{3}}(x,y)-f_{\varDelta_{1},i_{3}}(x,y) \nabla_{{\perp}}\,f_{\varDelta_{1},i_{1}}(x,y))]\nonumber\\ &\qquad +f_{\varDelta_{2}}(z)[(\phi_{i_{5}}-\phi_{i_{4}})(\,f_{\varDelta_{2}, i_{5}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{2},i_{4}}(x,y)-f_{\varDelta_{2}, i_{4}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{2},i_{5}}(x,y))\nonumber\\ &\qquad +(\phi_{i_{6}}-\phi_{i_{5}})(\,f_{\varDelta_{2},i_{6}}(x,y) \nabla_{{\perp}}\,f_{\varDelta_{2},i_{5}}(x,y)-f_{\varDelta_{2},i_{5}}(x,y) \nabla_{{\perp}}\,f_{\varDelta_{2},i_{6}}(x,y))\nonumber\\ &\qquad +(\phi_{i_{4}}-\phi_{i_{6}})(\,f_{\varDelta_{2},i_{4}}(x,y) \nabla_{{\perp}}\,f_{\varDelta_{2},i_{6}}(x,y)-f_{\varDelta_{2},i_{6}}(x,y) \nabla_{{\perp}}\,f_{\varDelta_{2},i_{4}}(x,y))]\nonumber\\ &\qquad +(\phi_{i_{4}}-\phi_{i_{1}})\,f_{\varDelta_{1},i_{1}}(x,y) \left[\,f_{\varDelta_{2}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{1}}(z)- f_{\varDelta_{1}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{2}}(z)\right]\nonumber\\ &\qquad +(\phi_{i_{5}}-\phi_{i_{2}})\,f_{\varDelta_{1},i_{2}}(x,y) \left[\,f_{\varDelta_{2}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{1}}(z) -f_{\varDelta_{1}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{2}}(z)\right]\nonumber\\ &\qquad +(\phi_{i_{6}}-\phi_{i_{3}})\,f_{\varDelta_{1},i_{3}}(x,y) \left[\,f_{\varDelta_{2}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{1}}(z) -f_{\varDelta_{1}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{2}}(z)\right]. \end{align}

Inserting (2.43) into (2.35), the equation of motion for particles on the prism mesh is

(2.44)\begin{align} \ddot{\boldsymbol{x}}_{p}&=\frac{q_{i}}{m_{i}}[\boldsymbol{E}_{0}(\boldsymbol{x}_{p},t) +\dot{\boldsymbol{x}}_{p}\times\boldsymbol{B}_{0}(\boldsymbol{x}_{p},t)\nonumber\\ &\quad +g_{i_{1}i_{2}}(\boldsymbol{x}_{p})\boldsymbol{E}_{i_{1}i_{2}} +g_{i_{2}i_{3}}(\boldsymbol{x}_{p})\boldsymbol{E}_{i_{2}i_{3}} +g_{i_{3}i_{1}}(\boldsymbol{x}_{p})\boldsymbol{E}_{i_{3}i_{1}}\nonumber\\ &\quad +g_{i_{4}i_{5}}(\boldsymbol{x}_{p})\boldsymbol{E}_{i_{4}i_{5}}+g_{i_{5}i_{6}}(\boldsymbol{x}_{p})\boldsymbol{E}_{i_{5}i_{6}}+g_{i_{6}i_{4}}(\boldsymbol{x}_{p})\boldsymbol{E}_{i_{6}i_{4}}\nonumber\\ &\quad +g_{i_{1}i_{4}}(\boldsymbol{x}_{p})\boldsymbol{E}_{i_{1}i_{4}}+g_{i_{2}i_{5}}(\boldsymbol{x}_{p})\boldsymbol{E}_{i_{2}i_{5}}+g_{i_{3}i_{6}}(\boldsymbol{x}_{p})\boldsymbol{E}_{i_{3}i_{6}}], \end{align}

where the discrete electric field defined on edge $i_{1}i_{2}$ is

(2.45)\begin{equation} \boldsymbol{E}_{i_{1}i_{2}}=\phi_{i_{2}}-\phi_{i_{1}},\end{equation}

and $\boldsymbol {E}_{i_{2}i_{3}}, \boldsymbol {E}_{i_{3}i_{1}}, \boldsymbol {E}_{i_{4}i_{5}}, \boldsymbol {E}_{i_{5}i_{6}}, \boldsymbol {E}_{i_{6}i_{4}}, \boldsymbol {E}_{i_{1}i_{4}}, \boldsymbol {E}_{i_{2}i_{5}}$ and $\boldsymbol {E}_{i_{3}i_{6}}$ are defined in a similar way. The interpolation from discrete electric field to particle is defined as

(2.46)\begin{gather} g_{i_{1}i_{2}}=f_{\varDelta_{1}}(z)[\,f_{\varDelta_{1},i_{2}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{1},i_{1}}(x,y)-f_{\varDelta_{1},i_{1}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{1},i_{2}}(x,y)], \end{gather}

(2.47)\begin{gather}g_{i_{2}i_{3}}=f_{\varDelta_{1}}(z)[\,f_{\varDelta_{1},i_{3}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{1},i_{2}}(x,y)-f_{\varDelta_{1},i_{2}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{1},i_{3}}(x,y)], \end{gather}

(2.48)\begin{gather}g_{i_{3}i_{1}}=f_{\varDelta_{1}}(z)[\,f_{\varDelta_{1},i_{1}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{1},i_{3}}(x,y)-f_{\varDelta_{1},i_{3}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{1},i_{1}}(x,y)], \end{gather}

(2.49)\begin{gather}g_{i_{4}i_{5}}=f_{\varDelta_{2}}(z)[\,f_{\varDelta_{2},i_{5}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{2},i_{4}}(x,y)-f_{\varDelta_{2},i_{4}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{2},i_{5}}(x,y)], \end{gather}

(2.50)\begin{gather}g_{i_{5}i_{6}}=f_{\varDelta_{2}}(z)[\,f_{\varDelta_{2},i_{6}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{2},i_{5}}(x,y)-f_{\varDelta_{2},i_{5}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{2},i_{6}}(x,y)], \end{gather}

(2.51)\begin{gather}g_{i_{6}i_{4}}=f_{\varDelta_{2}}(z)[\,f_{\varDelta_{2},i_{4}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{2},i_{6}}(x,y)-f_{\varDelta_{2},i_{6}}(x,y)\nabla_{{\perp}}\,f_{\varDelta_{2},i_{4}}(x,y)], \end{gather}

and

(2.52)\begin{gather} g_{i_{1}i_{4}}=f_{\varDelta_{1},i_{1}}(x,y)\left[\,f_{\varDelta_{2}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{1}}(z)-f_{\varDelta_{1}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{2}}(z)\right], \end{gather}

(2.53)\begin{gather}g_{i_{2}i_{5}}=f_{\varDelta_{1},i_{2}}(x,y)\left[\,f_{\varDelta_{2}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{1}}(z)-f_{\varDelta_{1}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{2}}(z)\right], \end{gather}

(2.54)\begin{gather}g_{i_{3}i_{6}}=f_{\varDelta_{1},i_{3}}(x,y)\left[\,f_{\varDelta_{2}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{1}}(z)-f_{\varDelta_{1}}(z)\frac{\partial}{\partial z}\,f_{\varDelta_{2}}(z)\right]. \end{gather}

Figure 3. The PIC algorithm on prism mesh. (a) Charge-deposition algorithm. Here $\rho _{i_{1}},\ldots ,\rho _{i_{6}}$ are the charge density on the vertices deposited from a particle at ($x, y, z$). (b) Potentials $\phi _{i_{1}},\ldots ,\phi _{i_{6}}$ are electric potentials on the vertices, and $\boldsymbol {E}_{i_{2}i_{3}}, \boldsymbol {E}_{i_{3}i_{1}}, \boldsymbol {E}_{i_{4}i_{5}}, \boldsymbol {E}_{i_{5}i_{6}}, \boldsymbol {E}_{i_{6}i_{4}}, \boldsymbol {E}_{i_{1}i_{4}}, \boldsymbol {E}_{i_{2}i_{5}}, \boldsymbol {E}_{i_{3}i_{6}}$ are discrete electric field on the edges. (c) Field-interpolation algorithm. Note that $\boldsymbol {E}_{i_{3}i_{1}}$, $\boldsymbol {E}_{i_{5}i_{6}}$ and $\boldsymbol {E}_{i_{1}i_{4}}$ are not labelled for clarity.

The diagram of the algorithm is shown in figure 3. The charge-deposition algorithm is shown in figure 3(a). Figure 3(b) shows the electric field on each edge calculated as (2.45). The electric field is interpolated to a particle's position by (2.46)–(2.54) as shown in figure 3(c). Note that (2.46)–(2.54) are consistent with the formalism of Whitney 1-forms on a prism mesh (Nedelec Reference Nedelec1980; Lohi & Kettunen Reference Lohi and Kettunen2021). In the PIC code, we take the linear barycentric function for $f_{\varDelta _{1},i_{1}}(x,y)$, $f_{\varDelta _{1},i_{2}}(x,y)$ and $f_{\varDelta _{1},i_{3}}(x,y)$, and we take the tent function for $f_{\varDelta _{1}}(z)$ and $f_{\varDelta _{2}}(z)$.

3.1. The IBW with periodic boundary conditions in a rectangular region

First, the IBW in a rectangular region of a uniform plasma is examined on an unstructured mesh with periodic boundary conditions. The simulation domain is shown in figure 4(a), and figure 4(b) shows a zoom-in of the red region. The simulation domain has 20 201 vertices and 40 000 triangles. To implement the periodic boundary conditions, each vertex at the left-hand boundary is identified with the corresponding vertex at the right-hand boundary. Similar identification is imposed for the top and bottom boundaries. The ion mass is $m_{i}=1.67 \times 10^{-27}$ kg and charge is $q_{i}=1.6 \times 10^{-19}$ C, the initial ion density is $n_{i0}=10^{20}~\textrm {m}^{-3}$, the electron and ion temperatures are 1000 eV and the out-of-plane background magnetic field is $B_{0}=2$ T. The simulation time step is $\Delta t=0.01/\varOmega _{i}$ and the total number of simulation particles is $5.12 \times 10^{7}$.

Figure 4. (a) The rectangular simulation domain. (b) A zoom-in of the red rectangular region in (a).

During the simulation, the electric potential $\phi$ on each vertex is recorded. The potentials $\phi$ on the vertices are interpolated to a rectangular mesh; after interpolation, data are analysed using the fast Fourier transform. The dispersion relations of the IBW can be inferred from $\tilde {\phi }(k_{x}, k_{y}, \omega )$, where $k_{x}$ and $k_{y}$ are the wave numbers in the $x$ and $y$ directions and $\omega$ is the angular frequency. Figure 5 plots the contours of the spectral power of $\tilde {\phi }$, and the contour peaks show the dispersion relation $k_{x}$–$\omega$ at a fixed $k_{y}$, where $k_{x}$ and $k_{y}$ are normalized to the ion gyroradius $\rho _{i}$ and $\omega$ to the ion gyrofrequency $\varOmega _{i}$. The contour plots of $\tilde {\phi }$ at $k_{y}\rho _{i}=$ 0, 1.0, 2.0 and 3.0 are shown in figure 5. The contour plots are compared with the theoretical dispersion relation with kinetic ions and adiabatic electrons (Sturdevant Reference Sturdevant2016):

(3.1)\begin{equation} 1+\theta\sum_{n={-}\infty}^{\infty}\frac{n\varOmega_{i} \varGamma_{n}(\textit{b})}{\omega+n\varOmega_{i}}=0.\end{equation}

Here, $\theta =q_{i}T_{e}/q_{e}T_{i}$, $b\equiv (k_{\perp }\rho _{i})^{2}$, $\varGamma _{n}\equiv I_{n}(\textit {b})e^{-b}$, $k_{\perp }$ is the perpendicular wave number and $I_{n}$ is the $n$th modified Bessel function of the first kind. For the present 2-D simulation, $k_{\perp }=\sqrt {k_{x}^{2}+k_{y}^{2}}$. The dispersion relation in terms of $(k_{x},k_{y},\omega )$ can be directly compared with the dispersion relation from the PIC simulation. The red dashed lines in figure 5 are the dispersion relation curves at $k_{y}\rho _{i}=$ 0, 1.0, 2.0, 3.0 from (3.1). The dispersion relations from the PIC simulations agree well with the theory. We also carry out the simulations on prism and tetrahedral meshes and obtain the dispersion relations of IBW for $k_{||}=0$; the dispersion relations are consistent with the theory results.

Figure 5. Dispersion relation of the IBW in a 2-D rectangular domain. The contour plot of $\tilde {\phi }(k_{x},\omega )$ at (a) $k_{y}\rho _{i}=0$, (b) $k_{y}\rho _{i}=1.0$, (c) $k_{y}\rho _{i}=2.0$ and (d) $k_{y}\rho _{i}=3.0$. The red dashed lines represent the theoretical dispersion relation.

3.2. The IBW in a 2-D circular region with fixed boundary conditions

In this subsection, we simulate the IBW in a 2-D circular domain using an unstructured mesh. The simulation domain is shown in figure 6(a), and figure 6(b) shows a zoom-in of the red rectangular region. The simulation domain has 7477 vertices and 14 646 triangles. The physical parameters are the same as in § 3.1. The total number of simulation particles is $4.9 \times 10^{10}$ in order to reduce noise and obtain eigenmode structures by the nonlinear PIC simulation, and the time step is $\Delta t=0.02/\varOmega _{i}$. The boundary condition is that $\phi =0$ at the outermost vertices, and particles are reflected when entering the outermost triangles.

Figure 6. (a) The 2-D circular simulation domain. (b) A zoom-in of the red rectangular region in (a).

For the fast Fourier transform analysis, $\phi$ is interpolated to a diagnostic circular mesh which has 101 co-concentric circles with equal intervals, and 61 grid points are distributed on each circle with the same angular interval. On the diagnostic mesh, $\phi$ is a function of time $t$ and the radial coordinates $(r,\theta )$. Performing the fast Fourier transform analysis of $\phi (r,\theta ,t)$ along the $\theta$ and $t$ directions, we discover that for each azimuthal mode number $m$, the spectrum of the system has a rich structure that can be labelled by two integer indices, $n$ and $l$. In the neighbourhood of each integer harmonics of $\varOmega _{i}$ labelled by $n$, there exists a family of eigenmodes labelled by $l$. The value of $l$ indicates the number of oscillations of the eigenmode in the radial direction. Thus, the eigenmode expansion for $\phi$ is

(3.2)\begin{equation} \phi=\sum_{n,m,l}\tilde{\phi}_{nml}(r)\exp (im\theta-i\omega_{nml}t), \end{equation}

where $\tilde {\phi }_{nml}(r)$ is the eigenfunction of the mode at $\omega =\omega _{nml}$. Figure 7(a) plots the spectrum of $\tilde {\phi }$ for $m=0,1,2,3$, where the frequency $\omega$ is normalized to ion gyrofrequency $\varOmega _{i}$. Figure 7(b) shows a zoom-in of the spectrum in the range $0\leq \omega /\varOmega _{i}\leq 3$. The eigenmode structures $\tilde {\phi }_{nml}(r)$ for $m=0,1,2,3$, $l=1,2,3$ in the neighbourhood of the first harmonic $(n=1)$ are shown in figure 8. We are not aware of any previous study of these eigenmodes of the IBW in a circular domain.

Figure 7. (a) Spectrum of the IBW for $m=0,1,2,3$. (b) Zoom-in of (a) for $0\leq \omega /\varOmega _{i}\leq 3$.

Figure 8. Eigenmode structures $\tilde {\phi }_{nml}(r)$ for $m=0,1,2,3$ and $l=1,2,3$ in the neighbourhood of the first harmonic ($n=1$). The black solid line is the real part of $\tilde {\phi }_{nml}(r)$ and the red dashed line is the imaginary part. The eigenfrequency $\omega _{nml}$ for each eigenmode is listed.