3.1 Discretization of $f$ and $\phi$

Here, we focus on the one-dimensional case with a periodic boundary condition, but higher dimensional cases can be treated similarly. The distribution function $f$ is approximated as

(3.1)\begin{equation} f_h({x}, {v}, t) = \sum_{k=1}^{N_p} w_k S ({x} - {x}_k) \delta ({v} - {v}_k), \end{equation}

where $N_p$ is the total particle number, and $w_k$, ${x}_k$ and ${v}_k$ denote the weight, position and velocity of the $k$th particle. Additionally, $S$ is the shape function of the particle, typically chosen as a B-spline. We use the vector ${\boldsymbol {X}}$ to denote $(x_1, \ldots, x_{N_p})^\top$ and vector ${\boldsymbol {V}}$ to denote $(v_1, \ldots, v_{N_p})^\top$.

The electrostatic potential $\phi$ is discretized using the finite difference method, i.e.

(3.2)\begin{equation} \phi_j \approx\phi({x}_j), \quad j = 1, \ldots, N, \end{equation}

with a set of uniform grids $\{ {x}_j \}$, $N$ is the number of grids, $(\phi _1, \ldots, \phi _N)^\top$ is denoted as ${\boldsymbol {\phi }}$ and $(\phi _0(x_1), \ldots, \phi _0(x_N))^\top$ is denoted as ${\boldsymbol {\phi }}_0$.

3.2 Phase-space discretization

3.2.1 Discretization of action integral

We approximate the variational action integral (2.6) as

(3.3)\begin{align} \mathcal{A}_h({\boldsymbol{X}}, {\boldsymbol{\phi}}) & = \sum_{k=1}^{N_p} w_k \left( \frac{1}{2} \dot{x}_k^2 - Z \sum_{j=1}^N \Delta x S({x}_j - {x}_k) \phi_j \right) + \frac{1}{2} {\boldsymbol{\phi}}^\top \boldsymbol{\mathsf{A}} {\boldsymbol{\phi}} \Delta x\nonumber\\ & \quad + \sum_{j=1}^N \Delta x n_0({x}_j) T_e({x}_j) \exp\left(\frac{\phi_j - \phi_0(x_j)}{T_e(x_j)}\right), \end{align}

where the matrix $\boldsymbol{\mathsf{A}}$ of size $N \times N$ is

(3.4)\begin{equation} \boldsymbol{\mathsf{A}} = \frac{1}{\Delta x^2} \left(\begin{matrix} -2 & 1 & 0 & \cdots & 0 & 0 & 1 \\ 1 & -2 & 1 & 0 & \cdots & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & 1 & -2 & 1 & 0\\ 0 & \cdots & 0 & 0 & 1 & -2 & 1 \\ 1 & 0 & 0 & \cdots & 0 & 1 & -2 \end{matrix} \right). \end{equation}

By calculating the variations about ${x}_k$ and ${\boldsymbol {\phi }}$, we have

(3.5)\begin{align} \ddot{x}_k & = Z \sum_j \Delta x \partial_x S({x}_j - {x}_k) \phi_j, \quad k = 1, \ldots, N_p,\nonumber\\ & \quad -Z \sum_{k=1}^{N_p} w_k S({x}_j - {x}_k) + (\boldsymbol{\mathsf{A}} {\boldsymbol{\phi}})_j + n_0({x}_j) \exp\left(\frac{\phi_j - \phi_0(x_j)}{T_e({x}_j)}\right) = 0, \quad j = 1, \ldots, N. \end{align}

Remark 3.1 The fixed point iteration method of Cohen et al. (Reference Cohen, Lasinski, Langdon and Williams1997) is used for solving the above discretized Poisson–Boltzmann equation, i.e. the second equation of (3.5),

(3.6)\begin{equation} \boldsymbol{\mathsf{A}} {\boldsymbol{\phi}}^{m+1} - c{\boldsymbol{\phi}}^{m+1} ={-}{\boldsymbol{n}}_i - c{\boldsymbol{\phi}}^{m} + {\boldsymbol{n}}_{e}^m, \end{equation}

where $m$ is the iteration index, the linear system in iteration $m \rightarrow m+1$ is solved by the conjugate gradient method Kelley (Reference Kelley1995), ${\boldsymbol {n}}_i$ is the ion density at grids and the electron density at the $j$-grid is

(3.7a,b)\begin{align} n_{e,j}^m = n_0({x}_j) \exp\left({\frac{\phi_j^m}{T_e({x}_j)}}\right)\quad \text{and}\quad c = \max \left\{\frac{1}{T_e(x_j)} \exp\left({{{\phi_j^m/T_e({x}_j), \, j = 1, \ldots, N}}}\right)\right\}. \end{align}

The initial value of the fixed point iteration is set as zero in § 6.

Equation (3.5) can be formulated as a Hamiltonian system by the Legendre transformation (Marsden & Ratiu Reference Marsden and Ratiu2013). In the following, we present another way to derive a finite dimensional Hamiltonian system.

3.2.2 Discretization of Poisson bracket

Here, we discretize the Poisson bracket (2.8) according to Qin et al. (Reference Qin, Liu, Xiao, Zhang, He, Wang, Sun, Burby, Ellison and Zhou2015b), Kraus et al. (Reference Kraus, Kormann, Morrison and Sonnendrücker2017) and obtain

(3.8)\begin{equation} \{F, G\}_h = \sum_{k=1}^{N_p} \frac{1}{w_k} \left(\partial_{x_k}F \partial_{v_k}G - \partial_{x_k}G \partial_{v_k}F \right). \end{equation}

The discrete Hamiltonian is

(3.9)\begin{align} H & ={-} \frac{1}{2} {\boldsymbol{\phi}}^\top \boldsymbol{\mathsf{A}} {\boldsymbol{\phi}} \Delta x + Z \sum_{j=1}^N \Delta x \sum_{k=1}^{N_p} w_k S({x}_j - {x}_k) \phi_j + \frac{1}{2} \sum_{k=1}^{N_p} w_k {v}_k^2\nonumber\\ & \quad - \sum_{j=1}^N \Delta x T_e({x}_j) n_0({x}_j) \exp\left({\frac{\phi_j - \phi_0(x_j)}{T_e({x}_j)}}\right), \end{align}

where the ${\boldsymbol {\phi }}$ is determined by the particles via the second equation in (3.5). Then, we can obtain the following finite dimensional Hamiltonian system after phase-space discretization:

(3.10)\begin{equation} \dot{x}_k = \frac{1}{w_k} \partial_{v_k}H, \quad \dot{v}_k ={-} \frac{1}{w_k} \partial_{x_k}H, \quad k = 1, \ldots, N_p. \end{equation}

The following theorem shows that the discretizations of the action principle and the Poisson bracket as above are equivalent.

Theorem 3.2 The Hamiltonian system (3.10) is equivalent to (3.5).

Proof. As we know that $\partial _{v_k}H = w_k {v}_k$, we have $\dot {x}_k = v_k$. To obtain (3.5), the only thing we need to prove is that $\partial _{x_k}H = - Z w_k \sum _{j=1}^N \Delta x \partial _x S({x}_j - {x}_k) \phi _j$, which is obtained by calculating $\partial _{x_k}H$ with the discrete Poisson–Boltzmann equation (the second equation in (3.5)),

(3.11)\begin{align} \partial_{x_k}H & ={-} Z w_k \sum_{j=1}^N \Delta x \partial_x S({x}_j - {x}_k) \phi_j - \boldsymbol{\mathsf{A}} {\boldsymbol{\phi}} \boldsymbol{\cdot} \frac{\partial \boldsymbol{\phi}}{\partial {x_k}} \Delta x\nonumber\\ & \quad + Z \sum_{j=1}^N \Delta x \sum_{k'=1}^{N_p} w_k S({x}_j - {x}_k') \frac{\partial \phi_j}{\partial {x_k}} - \sum_{j=1}^N \Delta x_j n_0({x}_j) \exp\left({\frac{\phi_j - \phi_0(x_j)}{T_e({x}_j)}}\right) \boldsymbol{\cdot} \frac{\partial \phi_j}{\partial {x_k}}, \end{align}

where the sum of the last three terms is zero because of the discrete Poisson–Boltzmann equation (the second equation in (3.5)).

Theorem 3.3 The discrete neutrality is conserved by the discretizations (3.5) and (3.10).

Proof. Taking the sum over $j$ in the discrete Poisson–Boltzmann equation (the second equation in (3.5)) gives

(3.12)\begin{equation} -Z \sum_{j=1}^N \sum_{k=1}^{N_p} w_k S({x}_j - {x}_k) \Delta x + \sum_{j=1}^N n_0({x}_j) \exp\left(\frac{\phi_j - \phi_0(x_j)}{T_e({x}_j)}\right) \Delta x = 0, \end{equation}

where we use that $\sum _j (\boldsymbol{\mathsf{A}} {\boldsymbol {\phi }})_j = 0$. This proves the discrete neutrality.

3.3 Time discretization

In this subsection, we introduce the time discretizations for (3.5) and (3.10). The first method is the Hamiltonian splitting method (Hairer et al. Reference Hairer, Lubich and Wanner2006), which is explicit and symplectic for (3.5) and (3.10). This method was applied by Crouseilles et al. (Reference Crouseilles, Einkemmer and Faou2015), Qin et al. (Reference Qin, He, Zhang, Liu, Xiao and Wang2015a) and He et al. (Reference He, Qin, Sun, Xiao, Zhang and Liu2015) for the Vlasov–Maxwell equations, and was used by Xiao et al. (Reference Xiao, Qin, Liu, He, Zhang and Sun2015), He et al. (Reference He, Sun, Qin and Liu2016) and Kraus et al. (Reference Kraus, Kormann, Morrison and Sonnendrücker2017) for the construction of the fully discrete structure-preserving methods. The other time discretization is the implicit discrete gradient method (McLachlan et al. Reference McLachlan, Quispel and Robidoux1999), which preserves the energy exactly.

Hamiltonian splitting method. We split the Hamiltonian (3.9) as $H = H_1 + H_2$, where

(3.13)\begin{equation} \left. \begin{aligned} H_1 & = \frac{1}{2} \sum_{k=1}^{N_p} w_k {v}_k^2,\\ H_2 & ={-} \frac{1}{2} {\boldsymbol{\phi}}^\top \boldsymbol{\mathsf{A}} {\boldsymbol{\phi}} \Delta x + Z \sum_{j=1}^N \Delta x \sum_{k=1}^{N_p} w_k S({x}_j - {x}_k) \phi_j\\ & \quad - \sum_{j=1}^N \Delta x T_e({x}_j) n_0({x}_j) \exp\left({\frac{\phi_j - \phi_0(x_j)}{T_e({x}_j)}}\right), \end{aligned} \right\} \end{equation}

which give the following two corresponding subsystems,

(3.14)\begin{equation} \left. \begin{aligned} & \text{sub-step I}: \dot{x}_k = {v}_k, \quad \dot{v}_k = 0,\\ & \text{sub-step II}: \dot{x}_k = 0, \quad \dot{v}_k = Z \sum_{j=1}^{N} \Delta x \partial_x S({x}_j - {x}_k) \phi_j, \end{aligned} \right\} \end{equation}

where $\phi _j (j = 1, \ldots, N)$ are given by the discrete Poisson–Boltzmann equation (the second equation in (3.5)). Both sub-steps can be solved exactly. Here, we present the first- and second-order methods by the composition method (Hairer et al. Reference Hairer, Lubich and Wanner2006),

(3.15)\begin{equation} \left. \begin{aligned} & \text{First-order Lie splitting}: \varPhi^{1}_{\Delta t} \circ \varPhi^{2}_{\Delta t},\\ & \text{Second-order Strang splitting}: \varPhi^{2}_{\Delta t/2} \circ \varPhi^{1}_{\Delta t} \circ \varPhi^{2}_{\Delta t/2}, \end{aligned} \right\} \end{equation}

where $\varPhi ^{1}_{\Delta t}$ and $\varPhi^2_{\Delta t}$ are solution maps of sub-steps I and II, respectively. Higher order structure-preserving schemes can be constructed by composition methods (Hairer et al. Reference Hairer, Lubich and Wanner2006).

Discrete gradient method. We use the second-order discrete gradient method proposed by Gonzalez (Reference Gonzalez1996) to conserve energy exactly,

(3.16a,b)\begin{equation} \frac{{\boldsymbol{X}}^{n+1} - {\boldsymbol{X}}^{n}}{\Delta t} = {\boldsymbol{\mathsf{W}}}^{{-}1} \bar{\boldsymbol{\nabla}}_{\boldsymbol{V}} H,\quad \frac{{\boldsymbol{V}}^{n+1} - {\boldsymbol{V}}^{n}}{\Delta t} ={-} \boldsymbol{\mathsf{W}}^{{-}1} \bar{\boldsymbol{\nabla}}_{\boldsymbol{X}} H, \end{equation}

where

(3.17)\begin{equation} \left. \begin{aligned} \boldsymbol{\mathsf{W}} & = \text{diag}(w_1, \ldots, w_{N_p}),\\ \bar{\boldsymbol{\nabla}}_{\boldsymbol{X}} H & = \boldsymbol{\nabla}_{\boldsymbol{X}} H \left(\frac{{\boldsymbol{X}}^n + {\boldsymbol{X}}^{n+1} }{2}\right) + d_{c} \left( {\boldsymbol{X}}^{n+1} - {\boldsymbol{X}}^n \right),\\ \bar{\boldsymbol{\nabla}}_{\boldsymbol{V}} H & = \boldsymbol{\nabla}_{\boldsymbol{V}} H \left(\frac{{\boldsymbol{V}}^n + {\boldsymbol{V}}^{n+1} }{2}\right) + d_{c} \left( {\boldsymbol{V}}^{n+1} - {\boldsymbol{V}}^n \right),\\ d_{c} & = \frac{H_d- \boldsymbol{\nabla} H ({\boldsymbol{X}}^{n+{1}/{2}}, {\boldsymbol{V}}^{n+{1}/{2}}) \boldsymbol{\cdot} ( ({\boldsymbol{X}}^{n+1} - {\boldsymbol{X}}^n)^\top, ({\boldsymbol{V}}^{n+1} - {\boldsymbol{V}}^n)^\top)^\top }{|{\boldsymbol{X}}^{n+1} - {\boldsymbol{X}}^n|^2 + |{\boldsymbol{V}}^{n+1} - {\boldsymbol{V}}^n|^2},\\ H_d & = H({\boldsymbol{H}}^{n+1}, {\boldsymbol{V}}^{n+1}) - H({\boldsymbol{X}}^n, {\boldsymbol{V}}^{n}). \end{aligned} \right\} \end{equation}

This discrete gradient method is implicit, for which the fixed-point iteration method is used. The degree of shape function $S$ is chosen to be at least two for the convergence of iterations.

4.1 Quasi-neutral limit

We do the normalization as

(4.1a-h)\begin{align} & \tilde{\boldsymbol{x}} = \frac{\boldsymbol{x}}{x^*}, \quad \tilde{\boldsymbol{v}} = \frac{\boldsymbol{v}}{C_0}, \quad \tilde{t} = t \omega_i, \quad \tilde{f} = \frac{C_0^3}{\bar{n}}f, \nonumber\\ & \tilde{n}_0 = \frac{n_0}{\bar{n}}, \quad \tilde{\phi} = \frac{e\phi}{k_B \bar{T}_e}, \quad \tilde{\phi}_0 = \frac{e\phi_0}{k_B \bar{T}_e}, \quad \tilde{T}_e = \frac{T_e}{\bar{T}_e}, \end{align}

where $x^*$ is the space scale of interest, $Z\omega _i = \sqrt {{\bar {n} Z^2 e^2}/{\epsilon _0 m_i}}$ is the ion plasma frequency, $\lambda _D = \sqrt {{\epsilon _0 k_B \bar {T}_e}/{\bar {n} e^2}}$ is the electron Debye length, $C_0 = \lambda _D \omega _i = \sqrt {{k_B \bar {T}_e}/{m_i}}$, $\bar {n}$ is the characteristic ion density and $\bar {T}_e$ is the characteristic electron temperature. Then, we get the normalized HBS model (tilde symbol is omitted for convenience)

(4.2)\begin{equation} \left. \begin{gathered} \frac{\partial f}{\partial t} + {\boldsymbol{v}} \boldsymbol{\cdot} \boldsymbol{\nabla} f + Z{\boldsymbol{E}} \boldsymbol{\cdot} \boldsymbol{\nabla}_v f = 0,\\ {\boldsymbol{E}} ={-} \boldsymbol{\nabla} \phi,\\ - \lambda^2 \Delta \phi = Z \int f\, \mathrm{d}{\boldsymbol{v}} - n_0 \exp\left(\frac{\phi-\phi_0}{T_e}\right), \quad \text{Poisson}{-}\text{Boltzmann}, \end{gathered} \right\} \end{equation}

where $\lambda = \lambda _D / x^*$.

When we take the quasi-neutral limit $\lambda \rightarrow 0$ for the normalized HBS model (4.2), we get the following hybrid model with quasi-neutrality and Boltzmann electrons (Rambo Reference Rambo1995; Xiao & Qin Reference Xiao and Qin2019a):

(4.3)\begin{equation} \left. \begin{gathered} \frac{\partial f}{\partial t} + {\boldsymbol{v}} \boldsymbol{\cdot} \boldsymbol{\nabla} f + Z{\boldsymbol{E}} \boldsymbol{\cdot} \boldsymbol{\nabla}_v f = 0,\\ {\boldsymbol{E}} ={-} \boldsymbol{\nabla} \phi,\\ 0 = Z \int f \,\mathrm{d}{\boldsymbol{v}} - n_0 \exp\left(\frac{\phi - \phi_0}{T_e}\right). \end{gathered} \right\} \end{equation}

By the similar calculations for the Vlasov–Poisson equation as Sonnendrücker (Reference Sonnendrücker2017), we get the dispersion relation of the linear Landau damping of (4.2) with $Z = 1, n_0=1, {\phi _0 = 0}$,

(4.4)\begin{equation} 1 + \lambda^2{k^2 T_e} = \frac{T_e}{T_i} \mathcal{Z}'\left(\frac{\omega}{kv_T} \right), \end{equation}

where $\mathcal {Z}$ is the plasma dispersion function and $T_i = v^2_T/2$. By $\lambda \rightarrow 0$, we get the dispersion relation of (4.3) (Kunz, Stone & Bai Reference Kunz, Stone and Bai2014),

(4.5)\begin{equation} 1 = \frac{T_e}{T_i} \mathcal{Z}'\left(\frac{\omega}{kv_T} \right). \end{equation}

Under this normalization (4.1a–h), the discrete Poisson–Boltzmann equation in our scheme (3.14) becomes

(4.6)\begin{equation} -Z \sum_{k=1}^{N_p} w_k S({x}_j - {x}_k) + \lambda^2(\boldsymbol{\mathsf{A}}{\boldsymbol{\phi}})_j + n_0({x}_j) \exp\left(\frac{\phi_j - \phi_0(x_j)}{T_e({x}_j)}\right) = 0, \quad j = 1, \ldots, N. \end{equation}

When taking the quasi-neutral limit for the scheme (3.14) with (4.6), we get the following structure-preserving scheme similar to the scheme proposed by Xiao & Qin (Reference Xiao and Qin2019a) derived using discrete exterior calculus and Whitney form,

(4.7)\begin{equation} \left. \begin{aligned} & \text{sub-step I}: \dot{x}_k= {v}_k, \quad \dot{v}_k= 0,\\ & \text{sub-step II}: \dot{x}_k = 0, \quad \dot{v}_k = Z \sum_{j=1}^N \Delta x \partial_x S({x}_j - {x}_k) \phi_j, \end{aligned} \right\} \end{equation}

where $\phi _j, j = 1, \ldots, N$ are determined by the following discrete equation about ${\boldsymbol {\phi }}$,

(4.8)\begin{equation} -Z \sum_{k=1}^{N_p} w_k S({x}_j - {x}_k) + n_0({x}_j) \exp\left(\frac{\phi_j - \phi_0(x_j)}{T_e}\right) = 0, \quad j = 1, \ldots, N. \end{equation}

Then the limiting scheme (4.7) is the Hamiltonian splitting method for the quasi-neutral limit model (4.3), i.e. scheme (3.14), and is asymptotic preserving (Jin Reference Jin1999) and structure-preserving at the same time. Similarly, the discrete gradient method (3.16a,b) for the HBS model becomes a discrete gradient method for the hybrid model with quasi-neutrality and Boltzmann electrons. Note that quasi-neutral limit is not a singular asymptotic limit as explained by Degond et al. (Reference Degond, Liu, Savelief and Vignal2012) for the case of the Euler–Poisson–Boltzmann model.

4.2 Large $T_e$ limit

Here, we adopt the normalization (2.3a–h). By taking the $T_e \rightarrow +\infty$ in (2.4), we get the equations

(4.9)\begin{equation} \left. \begin{gathered} \frac{\partial f}{\partial t} + {\boldsymbol{v}} \boldsymbol{\cdot} \boldsymbol{\nabla} f + Z{\boldsymbol{E}} \boldsymbol{\cdot} \boldsymbol{\nabla}_v f = 0,\\ {\boldsymbol{E}} ={-} \boldsymbol{\nabla} \phi,\\ -\Delta \phi = Z \int f\, \mathrm{d}{\boldsymbol{v}} - n_0. \end{gathered} \right\} \end{equation}

Dividing by $k^2T_e$ on both sides of the dispersion relation (4.4) of (2.4), we get the the following dispersion relation with the current normalization:

(4.10)\begin{equation} 1 + \frac{1}{k^2 T_e} = \frac{1}{k^2 T_i} \mathcal{Z}'\left(\frac{\omega}{kv_T} \right). \end{equation}

By $T_e \rightarrow +\infty$, we get the dispersion relation of model (4.9) (Sonnendrücker Reference Sonnendrücker2017),

(4.11)\begin{equation} 1 = \frac{1}{k^2 T_i} \mathcal{Z}'\left(\frac{\omega}{kv_T} \right). \end{equation}

Similar to the quasi-neutral limit, when $T_e \rightarrow +\infty$, the limiting schemes of the Hamiltonian splitting method (3.14) and the discrete gradient method (3.16a,b) become the Hamiltonian splitting method and the discrete gradient method for model (4.9).

6.1 Finite grid instability

Finite grid instability in the context of hybrid simulations was first reported by Rambo (Reference Rambo1995) for the hybrid model with quasi-neutrality and Boltzmann electrons. This numerical phenomenon typically arises in standard particle-in-cell methods when the temperature ratio $T_e/T_i \gg 1$, and ions are heated until the ion thermal speed becomes comparable to the ion acoustic speed (and therefore, $T_e/T_i \approx 1$). Rambo (Reference Rambo1997) noted that finite grid instability also occurs when using traditional particle-in-cell methods for the HBS model, although it is weaker than the hybrid model with quasi-neutrality and Boltzmann electrons. Stanier, Chacón & Chen (Reference Stanier, Chacón and Chen2019) and Li et al. (Reference Li, Campos-Pinto, Holderied, Possanner and Sonnendrücker2024a,Reference Li, Holderied, Possanner and Sonnendrückerb) numerically reduced the finite grid instability by using the conservative or bracket-based particle-in-cell methods for the hybrid model with quasi-neutrality and massless electrons. The finite grid instability of particle-in-cell methods has also been studied by Huang et al. (Reference Huang, Zeng, Wang, Meyers, Yi and Albright2016) and Xiao & Qin (Reference Xiao and Qin2019b), which reveals that the aliased spatial modes are the major cause of the finite grid instability in the particle-in-cell methods, and geometric particle-in-cell methods are able to suppress the finite grid instability.

In this test, we investigate the finite grid instability using the following initial condition (an equilibrium of (2.4)) by the numerical simulations conducted with schemes (3.14) and (3.16a,b),

(6.1)\begin{equation} f = \frac{n_i}{{\rm \pi}^{{1}/{2}} v_T^{{1}/{2}}} \exp\left({\frac{|v - 0.1|^2}{v_T^2} }\right), \quad T_e = 0.08, \quad v_T = 0.1, \quad n_i = 1. \end{equation}

Computational parameters include: domain $[0, 5{\rm \pi} ]$, time step size $\Delta t = 0.05$ and particle number per cell $100$. We run the simulations with the numerical methods (3.14) and (3.16a,b) with different cell sizes, i.e. $\Delta x = 5{\rm \pi} /12, 5{\rm \pi} /25, 5{\rm \pi} /50, 5{\rm \pi} /100$, and the results are presented in figures 1 and 2. We can see $k(t)/k(0)$ ($k(t) = \tfrac {1}{2} \sum _{k=1}^{N_p} w_k v_k^2$ the ion kinetic energy) oscillates with time without exhibiting rapid linear growth, as observed in figure 3(a) of Rambo (Reference Rambo1997). This indicates that the finite grid instability is reduced numerically. As the cell size decreases and the electron Debye length is resolved with higher resolution, the oscillating level of $k(t)/k(0)$ becomes closer to 1 and the momentum error also becomes smaller. The momentum error also depends on the particle number. When there are 100 cells and 2000 particles in each, the Hamiltonian splitting method with quadratic weighting gives the momentum error at the level of $10^{-4}$.

Figure 1.

Finite grid instability of the HBS model by Hamiltonian splitting method. Time evolution of $k(t)/k(0)$ with $k$ denoting the ion kinetic energy, relative energy error and momentum error.

Figure 2.

Finite grid instability of the HBS model by discrete gradient method with quadratic weighting. Time evolution of $k(t)/k(0)$ with $k$ denoting the ion kinetic energy, relative energy error and momentum error.

As the derivatives of B-splines appear in the schemes (3.14) and (3.16a,b), second-order at least B-spline shape functions should be used. As the Hamiltonian splitting method (3.14) (Strang splitting) is symplectic, it has superior long-time numerical behaviour, although energy is not conserved exactly (with a relative error of approximately $10^{-5}$) with quadratic weighting. In figure 1, the relative energy error is also very small ($10^{-4}$) even when the first order B-spine is used (the derivative of the B-spine is taken as its right derivative). Relative energy error of the discrete gradient method with quadratic weighting is approximately $10^{-13}$.

Here, we discuss the time step size of the Hamiltonian splitting methods. We consider the case with 100 cells in space, time interval $[0,500]$ and quadratic weighting. To achieve satisfying results, when $n_i =1$ and the number of particles per cell is 100, the maximum time step size is 0.5 approximately with energy error at the level of $10^{-4}$. For higher initial densities, such as $n_i = 4, 16$ with $400, 1600$ particles per cell, the maximum time step sizes giving satisfying results are approximately $0.3,0.3$ with energy errors at the level of $10^{-4}, 10^{-4}$, respectively.

6.2 Landau damping

First, we simulate the linear ion Landau damping by one-dimensional simulations. The initial distribution function is

(6.2)\begin{equation} f = \frac{n_i}{{\rm \pi}^{{1}/{2}} v_T^{{1}/{2}}} (1 + 0.02\cos(0.25x)) \exp\left({-\frac{v^2}{v_T^2} }\right). \end{equation}

The computational parameters are as follows: grid number 64, domain size $[0, 8{\rm \pi} ]$, time step size $\Delta t = 0.05$, final computation time $20$, $v_T = 1.4142$, $n_i = 1$ and total particle number $10^7$. In this test, the quadratic weighting is used. See the numerical results with $T_e=5$ in figure 3 by the Hamiltonian splitting method (3.14) and discrete gradient method (3.16a,b). Solving the dispersion relation mentioned in § 4, $1 + {k^2 T_e} = ({T_e}/{T_i}) \mathcal {Z}'({\omega }/{kv_T} ),$ we find $\omega = 0.6986 - 0.0810\textrm {i}$ when $k = 0.25$. Methods (3.14)–(3.16a,b) give an accurate damping rate of the electric energy $\tfrac {1}{2}\int |\boldsymbol {\nabla } \phi |^2 \,\mathrm {d} x$. The dispersion relation $1 = ({T_e}/{T_i}) \mathcal {Z}'({\omega }/{kv_T} )$ of the model with quasi-neutrality and Boltzmann electrons (Xiao & Qin Reference Xiao and Qin2019a) in § 4 yields $\omega = 0.7528 - 0.05806\textrm {i}$, i.e. a slower damping rate. The energy errors of the schemes (3.14)–(3.16a,b) are approximately $10^{-4}$ and $10^{-13}$, respectively.

Figure 3.

Linear Landau damping of the HBS model with $T_e=5$ by Hamiltonian splitting and discrete gradient methods. Time evolution of the electric energy $\tfrac {1}{2}\int |\boldsymbol {\nabla } \phi |^2 \,\mathrm {d} x$ and total energy error.

Then, we simulate nonlinear ion Landau damping. The initial distribution function is

(6.3)\begin{equation} f = \frac{n_i}{{\rm \pi}^{{1}/{2}} v_T^{{1}/{2}}} (1 + 0.5\cos(0.5x)) \exp\left({-\frac{v^2}{v_T^2} }\right). \end{equation}

The computational parameters are as follows: grid number 65, domain size $[0, 4{\rm \pi} ]$, time step size $\Delta t = 0.05$, final computation time $40$, $v_T = 1.4142, n_i = 1$ and total particle number $10^5$. In this test, the quadratic weighting is used. See the numerical results with large $T_e=100$ in figure 4 by the Hamiltonian splitting method (3.14) and discrete gradient method (3.16a,b). Due to the large $T_e$, the term $\exp (\phi /T_e)$ approximates 1, making the solution of the HBS model (2.4) approximate the solution of the Vlasov–Poisson system (with static electron density as 1). In figure 4, we observe the nonlinear Landau damping. The time evolution of energy component $\tfrac {1}{2}\int |\boldsymbol {\nabla } \phi |^2 \,\mathrm {d}{x}$ decays exponentially initially, and the decay rate is very close to the decay rate 0.2854 of the Vlasov–Poisson system (Kraus et al. Reference Kraus, Kormann, Morrison and Sonnendrücker2017) before time $T=10$. Here, $\tfrac {1}{2}\int |\boldsymbol {\nabla } \phi |^2 \,\mathrm {d}{x}$ oscillates when $t \in [10, 30]$, then grows exponentially with time when $t \in [30, 40]$ with a rate close to 0.086671 (Kraus et al. Reference Kraus, Kormann, Morrison and Sonnendrücker2017) of the Vlasov–Poisson system. For the Hamiltonian splitting method, the energy error is approximately $10^{-2}$. The discrete gradient method gives a smaller total energy error of approximately $10^{-12}$, and a similar behaviour of the electric energy is presented. The errors of neutrality given by both numerical methods are at the level of iteration tolerance. For the time step size of the Hamiltonian splitting methods, we consider the case with 65 cells in space, time interval $[0,40]$, $10^5$ particles and quadratic weighting. As $T_e=100$ is large, $\exp (\phi /T_e)$ is close to 1, the numerical stability property is close to the result of Kormann & Sonnendrücker (Reference Kormann and Sonnendrücker2021), i.e. the stability condition is approximately $\Delta t \omega _i < 2$. To get the acceptable accuracy, the time step size is usually chosen smaller than 2. When $n_i = 1$, the maximum time step size yielding good numerical behaviour is $\Delta t = 0.4$, resulting in an energy error of approximately $0.45$ after saturation; for $n_i = 4$, $\Delta t = 0.2$ gives an energy error of approximately $0.4$ after saturation, and for $n_i = 16$, $\Delta t = 0.12$ results in an energy error of approximately $0.3$ after saturation. We also consider the case with a small electron temperature $T_e = 1$. In this case, when $n_i = 1$, $\Delta t = 0.5$ gives an energy error of approximately 0.05 after saturation; when $n_i = 4$, $\Delta t = 0.2$ gives an energy error of approximately 0.002 after saturation; when $n_i = 16$, $\Delta t = 0.1$ gives an energy error of approximately 0.004.

Figure 4.

Nonlinear Landau damping of the HBS model with $T_e=100$ by Hamiltonian splitting and discrete gradient methods. Time evolution of total energy error and electric energy $\tfrac {1}{2}\int |\boldsymbol {\nabla } \phi |^2 \,\mathrm {d} x$.

6.3 Simulations with the ponderomotive driving term

Cohen et al. (Reference Cohen, Lasinski, Langdon and Williams1997) numerically solved the HBS model with a ponderomotive driving term to study nonlinear ion acoustic waves. Here, following Cohen et al. (Reference Cohen, Lasinski, Langdon and Williams1997), we conduct a simulation with a non-zero given time-dependent function $\phi _0$ in (2.4). Specifically, the initial condition and $\phi _0$ are given by

(6.4)\begin{equation} f = \frac{n_i}{{\rm \pi}^{{1}/{2}} v_T^{{1}/{2}}} \exp\left({-\frac{v^2}{v_T^2} }\right), \quad \phi_0 = \tilde{\phi}_0 \cos(\varOmega t - k x), \end{equation}

where $n_i = 1, v_T = \tfrac {\sqrt {2}}{10}$, $\tilde {\phi }_0 = 0.05T_e$, $\varOmega = 0.4472$, $k = 1.49$. Other computational parameters are: grid number 64, domain size $[0, {10{\rm \pi} }/{k}]$, time step size $\Delta t = 0.1$, final computation time $600$, $T_e = 0.1125$ and total particle number $10^6$. In this test, the quadratic weighting is used. Since $\phi _0$ is time dependent, the Hamiltonian system (3.10) is a non-autonomous Hamiltonian system, for which we use the the technique of extending the dimension (Zhou et al. Reference Zhou, He, Sun, Liu and Qin2017). From figures 5 and 6, we can see that the peak value of the response function $R(t) = \max _{x}({\phi }/{\phi _0})$ is approximately 5, which is consistent with the result of Cohen et al. (Reference Cohen, Lasinski, Langdon and Williams1997). There are a rapid oscillation at $2.1\varOmega$ and a slow modulation at $0.04\varOmega$ in the fourth figure obtained by the fast Fourier transformation of $\max _x({\phi }/{\phi _0})$ in time, which are close to the results of Cohen et al. (Reference Cohen, Lasinski, Langdon and Williams1997) with $1.85 \varOmega$ (fast) and $0.15\varOmega$ (slow). The Hamiltonian splitting method and discrete gradient method give the energy errors of approximately $10^{-5}$ and $10^{-12}$, respectively. Additionally, when $t=400$, both methods give similar five vortices in phase-space contour plots, due to the driving force being the fifth mode, i.e. $k = 5 ({2{\rm \pi} }/{L})$, where $L$ is the domain size.

Figure 5.

Simulations with the ponderomotive driving term by Hamiltonian splitting method. Time evolutions of $R(t) = \max ({\phi }/{\tilde {\phi }_0})$ and energy error, the contour plot of the distribution function at time $t=400$, and the fast Fourier transformation of $R(t)$.

Figure 6.

Simulations with the ponderomotive driving term by discrete gradient method. Time evolutions of $R(t) = \max ({\phi }/{\tilde {\phi }_0})$ and energy error, the contour plot of the distribution function at time $t=400$, and the fast Fourier transformation of $R(t)$.

Regarding the time step size for the Hamiltonian splitting methods, we use the same parameters as above but vary the initial density and time step size. When the initial density $n_i = 1$, the maximum step size for giving good numerical behaviour is $\Delta t = 1.5$, which gives an energy error of approximately 0.003; When the initial density $n_i = 4$, the maximum step size $\Delta t = 1$ gives an energy error of approximately 0.02; when the initial density $n_i = 16$, the maximum step size $\Delta t = 0.5$ gives an energy error of approximately 0.008.