# A Convolutional Method for Solving Non-quasi-neutral Hybrid-PIC Model of
Plasma

Jin-Long Jiao<EMAIL_ADDRESS>
###### Abstract

The hybrid model of particle-ion-fluid-electron is widely utilized in
astrophysical and laser-produced plasmas due to its kinetic characteristics
and remarkable computational efficiency. However, the current numerical
methods for the hybrid model have limitations when it comes to simulating non-
quasi-neutral plasmas. To address this issue, I proposed a novel method called
Convolutional Hybrid PIC (conv-HPIC). The conv-HPIC is based on a convolution
equation derived from the Poisson’s equation and the generalized Ohm’s law. I
performed numerical simulations on the collisionless electrostatic shock to
validate the effectiveness of this method. The results demonstrate that the
proposed method can accurately describe and predict the behavior of non-quasi-
neutral plasmas. This method can be applied to both stationary and non-
stationary electrostatic issues. It avoids the need for iterative
calculations, effectively suppresses numerical noise, exhibits good
computational efficiency, and can be easily parallelized.

###### keywords:

Plasma physics , Hybrid model , Convolutional method

organization=Zhejiang Institute of Modern Physics, Institute of Astronomy,
School of Physics, Zhejiang University,addressline=No. 866 Yuhangtang Road,
city=Hangzhou, postcode=310027, state=Zhejiang, country=China

##  1 Introduction

The hybrid particle-in-cell (HPIC) method is a simulation model used in plasma
physics, where some plasma species are treated as fluids while others are
treated as kinetic particles. The HPIC category discussed here specifically
refers to particle-ion-fluid-electron. In this model, ions are treated
kinetically, electrons are assumed to be a massless quasi-neutral fluid, and
the electromagnetic fields are solved using the low-frequency approximation
[1]. The HPIC offers a significant improvement in computational efficiency
compared to the kinetic PIC method, particularly for simulating plasmas within
ion response times [2]. Therefore, the HPIC is widely employed in numerical
studies involving various astrophysical and laser-produced plasmas on large
scales, e.g., magnetic reconnection, bow shock waves, and inertial confinement
fusion [3, 4, 5].

The quasi-neutral condition is commonly employed in the HPIC to simplify the
algorithm, but it eliminates the charge separation effect. Consequently,
simulating non-quasi-neutral plasmas becomes challenging using quasi-neutral
HPIC. To address this limitation, several methods have been proposed. One
method is replacing the quasi-neutral condition with the Poisson’s equation.
The Poisson’s equation involves two unknown physical quantities (electrostatic
potential and electron density), which need to be solved together with the
generalized Ohm’s law. The main approach for solving these equations is
iterative method [3, 6]. This method requires the manual selection of
iteration number (or tolerance) and initial value, and may lead to program
interruption or divergence. Another method is retaining the displacement
current term in the Ampere’s law. This method ensures that the Poisson’s
equation automatically holds in a charge conservation system [5]. It avoids
iterative calculations. However, this method is time-dependent and therefore
cannot be applied to stationary electrostatic issues.

Here, I introduced a novel method, called Convolutional Hybrid PIC (conv-
HPIC), for simulating non-quasi-neutral plasmas. This method is based on a
convolution equation derived from the Poisson’s equation and the generalized
Ohm’s law. The conv-HPIC method avoids iterative calculations and is suited
for stationary electrostatic plasmas. Furthermore, it can effectively suppress
the inherent numerical noise in particle methods.

##  2 The Convolutional Hybrid-PIC Method

Firstly, I review the quasi-neutral HPIC model [1]. In this model, ions are
regarded as pseudo-particles and their motion satisfy the Lorentz’s equation
(in c.g.s. unit),

$\displaystyle
m_{i}\frac{d\textbf{v}_{i}}{dt}=q_{i}\left(\textbf{E}+\frac{1}{c}\textbf{v}_{i}\times\textbf{B}\right),$
(1) $\displaystyle\frac{d\textbf{r}_{i}}{dt}=\textbf{v}_{i}.$ (2)

Where $\textbf{v}_{i}$ and $\textbf{r}_{i}$ are the velocity and position of
any pseudo-particle, respectively. The electric field E and magnetic field B
are determined by the Faraday’s law and the generalized Ohm’s law,

$\displaystyle\frac{\partial\textbf{B}}{\partial t}=-c\nabla\times\textbf{E},$
(3) $\displaystyle\textbf{E}=\frac{-1}{en_{e}}\nabla
p_{e}-\frac{1}{c}\textbf{u}_{\textbf{e}}\times\textbf{B}-\frac{1}{4\pi
en_{e}}\textbf{B}\times\left(\nabla\times\textbf{B}\right).$ (4)

The electron fluid velocity $\textbf{u}_{e}$ and electron density $n_{e}$ in
the generalized Ohm’s law (4) are obtained from the ion fluid velocity
$\textbf{u}_{i}$ and ion density $n_{i}$ by the quasi-neutral conditions of
current and charge (i.e. $\textbf{u}_{e}\cong\textbf{u}_{i}$ and $n_{e}\cong
Zn_{i}$). $\textbf{u}_{i}$ and $n_{i}$ can be obtained by mapping pseudo-
particles to grids.

To describe the non-quasi-neutral plasmas, the quasi-neutral condition of
charge $n_{e}\cong Zn_{i}$ should be replaced by the Poisson’s equation,

$\nabla\cdot\textbf{E}=4\pi e\left(n_{i}-n_{e}\right).$ (5)

The electric field E and electron density $n_{e}$ are two unknown variables,
and are obtained by solving the combination of the Poisson’s equation (5) and
the generalized Ohm’s law (4). The mainly solving method for the two equations
is the iterative calculation at present [6]. Below I introduce a new method to
solve the two equations.

The electron pressure $p_{e}$ is a function of electron density $n_{e}$ and
temperature $T_{e}$ in general. So, the gradient of the pressure can be
expanded as $\nabla p_{e}\left(n_{e},T_{e}\right)=\left.\frac{\partial
p_{e}}{\partial n_{e}}\right|_{T_{e}}\nabla n_{e}+\left.\frac{\partial
p_{e}}{\partial T_{e}}\right|_{n_{e}}\nabla T_{e}$. Here defines
$\alpha=\left.\frac{\partial p_{e}}{\partial n_{e}}\right|_{T_{e}}$,
$\beta=\left.\frac{\partial p_{e}}{\partial T_{e}}\right|_{n_{e}}$ and
$\varepsilon=\frac{1}{4\pi e^{2}n_{e}}$. The Eq. (4) can be written as,

$\textbf{E}=-4\pi e\varepsilon\left(\alpha\nabla n_{e}+\beta\nabla
T_{e}\right)-\textbf{F}.$ (6)

Where
$\textbf{F}=\frac{1}{c}\textbf{u}_{e}\times\textbf{B}+e\varepsilon\textbf{B}\times\left(\nabla\times\textbf{B}\right)$.

Figure 1: Division of simulation domain.

The simulation domain is divided into several subregions $\Omega_{N}$, as
shown in Fig. 1. For any subregion $\Omega_{i}$, its range is
$\textbf{r}\in\left[\textbf{r}_{i}-\delta\textbf{r}_{i}/2,\textbf{r}_{i}+\delta\textbf{r}_{i}/2\right]$.
The Eq. (6) on the subregion $\Omega_{i}$ can be written as,

$\displaystyle\Pi\left(\frac{\textbf{r}-\textbf{r}_{i}}{\delta\textbf{r}_{i}}\right)\textbf{E}=-\Pi\left(\frac{\textbf{r}-\textbf{r}_{i}}{\delta\textbf{r}_{i}}\right)\left[4\pi
e\varepsilon\left(\alpha\nabla n_{e}+\beta\nabla
T_{e}\right)+\textbf{F}\right].$ (7)

Where $\Pi\left(\frac{\textbf{r}-\textbf{r}_{i}}{\delta\textbf{r}_{i}}\right)$
is the rectangular function. Its value is 1 on the subregion $\Omega_{i}$,
otherwise 0. If the subregion $\Omega_{i}$ is small enough (i.e.
$\delta\textbf{r}_{i}\rightarrow 0$), one can assume that,

$\displaystyle\alpha\left(\textbf{r}\right)\approx\alpha\left(\textbf{r}_{i}\right),\
\beta\left(\textbf{r}\right)\approx\beta\left(\textbf{r}_{i}\right),\
\varepsilon\left(\textbf{r}\right)\approx\varepsilon\left(\textbf{r}_{i}\right).$
(8)

Thus, the Fourier transform of the first term on the right side of the Eq. (7)
is,

$\displaystyle\int_{\textbf{R}}{\Pi\left(\frac{\textbf{r}-\textbf{r}_{i}}{\delta\textbf{r}_{i}}\right)\varepsilon\left(\textbf{r}\right)\alpha\left(\textbf{r}\right)\nabla
n_{e}\left(\textbf{r}\right)e^{-\mathrm{i}\textbf{k}\cdot\textbf{r}}d\textbf{r}}\approx\varepsilon\left(\textbf{r}_{i}\right)\alpha\left(\textbf{r}_{i}\right)\int_{\textbf{R}}{\Pi\left(\frac{\textbf{r}-\textbf{r}_{i}}{\delta\textbf{r}_{i}}\right)\nabla
n_{e}\left(\textbf{r}\right)e^{-\mathrm{i}\textbf{k}\cdot\textbf{r}}d\textbf{r}}=\varepsilon\left(\textbf{r}_{i}\right)\alpha\left(\textbf{r}_{i}\right)\Theta\ast\left[\mathrm{i}\textbf{k}n_{e}\left(\textbf{k}\right)\right].$
(9)

Where $\Theta$ represents the Fourier transform of
$\Pi\left(\frac{\textbf{r}-\textbf{r}_{i}}{\delta\textbf{r}_{i}}\right)$. The
symbol [$\ast$] represents the convolution operator. The remaining terms in
Eq. (7) can be transformed in the similar procedure. Finally, the Fourier
transform of the Eq. (7) is,

$\displaystyle\textbf{E}\left(\textbf{k}\right)=-4\pi
e\mathrm{i}\textbf{k}\left[\varepsilon\left(\textbf{r}_{i}\right)\alpha\left(\textbf{r}_{i}\right)n_{e}\left(\textbf{k}\right)+\varepsilon\left(\textbf{r}_{i}\right)\beta\left(\textbf{r}_{i}\right)T_{e}\left(\textbf{k}\right)\right]-\textbf{F}\left(\textbf{k}\right).$
(10)

Similarly, the Fourier transform of the Poisson’s equation (5) is,

$\mathrm{i}\textbf{k}\cdot\textbf{E}\left(\textbf{k}\right)=4\pi
e\left[n_{i}\left(\textbf{k}\right)-n_{e}\left(\textbf{k}\right)\right].$ (11)

The electric field $\textbf{E}\left(\textbf{k}\right)$ can be eliminated by
substituting the Eq. (10) into (11). Then one can obtain,

$\displaystyle
n_{e}\left(\textbf{k}\right)=\left[n_{i}\left(\textbf{k}\right)-\varepsilon\left(\textbf{r}_{i}\right)\beta\left(\textbf{r}_{i}\right)k^{2}T_{e}\left(\textbf{k}\right)+\frac{1}{4\pi
e}\rm{i}\textbf{k}\cdot\textbf{F}\left(\textbf{k}\right)\right]f\left(k\right).$
(12)

Where
$f\left(k\right)=\frac{1}{1+{\varepsilon\left(\textbf{r}_{i}\right)\alpha\left(\textbf{r}_{i}\right)k}^{2}}$,
$k=\left|\textbf{k}\right|=\sqrt{k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}$. Applying the
inverse Fourier transform for the Eq. (12), the electron density in r-space
can be written as a convolution form,

$\displaystyle
n_{e}\left(\textbf{r},\textbf{r}_{i}\right)=\left[n_{i}\left(\textbf{r}\right)+\varepsilon\left(\textbf{r}_{i}\right)\beta\left(\textbf{r}_{i}\right)\nabla^{2}T_{e}\left(\textbf{r}\right)+\frac{1}{4\pi
e}\nabla\cdot\textbf{F}\left(\textbf{r}\right)\right]\ast f\left(r\right)$
$\displaystyle=\int_{\textbf{R}}\left[n_{i}\left(\textbf{s}\right)+\varepsilon\left(\textbf{r}_{i}\right)\beta\left(\textbf{r}_{i}\right)\nabla^{2}T_{e}\left(\textbf{s}\right)+\frac{1}{4\pi
e}\nabla\cdot\textbf{F}\left(\textbf{s}\right)\right]f\left(\left|\textbf{r}-\textbf{s}\right|\right)d\textbf{s}.$
(13)

Where $r=\left|\textbf{r}\right|=\sqrt{x^{2}+y^{2}+z^{2}}$. The Fourier
transform of $f\left(k\right)$ can be converted to the $n$-dimensional Hankel
transform due to $k^{2}$ is symmetrical [7],

$f\left(r\right)=\frac{\left(2\pi\right)^{-\frac{n}{2}}}{r^{\frac{n}{2}-1}}\int_{0}^{\infty}{f\left(k\right)J_{\frac{n}{2}-1}\left(kr\right)k^{n/2}dk}.$
(14)

Where $J$ is the Bessel function of the first kind. Here defines
$a^{2}=\frac{1}{\varepsilon\left(\textbf{r}_{i}\right)\alpha\left(\textbf{r}_{i}\right)}$.

(1) For 1D case ($n$=1, $k=|kx|$, $r=|x|$),

$f\left(r\right)=\frac{a}{2}e^{-ar}.$ (15)

One can prove that $\int_{-\infty}^{\infty}f\left(r\right)dr=1$, so
$f\left(r\right)$ is normalized.

(2) For 2D case ($n$=2, $k=\sqrt{k_{x}^{2}+k_{y}^{2}}$,
$r=\sqrt{x^{2}+y^{2}}$),

$f\left(r\right)=\frac{a^{2}}{2\pi}K_{0}\left(ar\right).$ (16)

Where $K_{0}$ is the modified Bessel functions of second kind of orders zero.
One can prove that $\lim\limits_{\delta r\rightarrow 0}{\int_{\delta
r}^{\infty}f\left(r\right)2\pi rdr}=1$, so $f\left(r\right)$ is normalized.

(3) For 3D case ($n$=3, $k=\sqrt{k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}$,
$r=\sqrt{x^{2}+y^{2}+z^{2}}$),

$f\left(r\right)=\frac{a^{2}}{4\pi r}e^{-ar}.$ (17)

One can prove that $\int_{0}^{\infty}{f\left(r\right)4\pi r^{2}dr}=1$, so
$f\left(r\right)$ is normalized. Note that the kernel function
$f\left(r\right)$ has a same form with the Debye screening potential in
plasma. That means the Eq. (13) contains the charge separation effect self-
consistently.

The Eq. (13) holds only in the subregion $\Omega_{i}$ because the assumption
(8) is used. The electron density in the whole region can be obtained by
integrating all the subregions,

$\displaystyle
n_{e}\left(\textbf{r}\right)=\sum_{i}\left[\Pi\left(\frac{\textbf{r}-\textbf{r}_{i}}{\delta\textbf{r}_{i}}\right)n_{e}\left(\textbf{r},\textbf{r}_{i}\right)\right]\approx\int_{\textbf{R}}{\delta\left(\textbf{r}-\textbf{r}_{i}\right)n_{e}\left(\textbf{r},\textbf{r}_{i}\right)d\textbf{r}_{i}}$
$\displaystyle=\int_{\textbf{R}}\left[n_{i}\left(\textbf{s}\right)+\varepsilon\left(\textbf{r}\right)\beta\left(\textbf{r}\right)\nabla^{2}T_{e}\left(\textbf{s}\right)+\frac{1}{4\pi
e}\nabla\cdot\textbf{F}\left(\textbf{s}\right)\right]f\left(\left|\textbf{r}-\textbf{s}\right|\right)d\textbf{s}.$
(18)

Where the coefficient $a^{2}$ in the kernel function $f\left(r\right)$ is
$a^{2}=\frac{1}{\varepsilon\left(\textbf{r}\right)\alpha\left(\textbf{r}\right)}$.
Here we use a relationship $\lim\limits_{\delta\textbf{r}_{i}\rightarrow
0}{\frac{1}{\delta\textbf{r}_{i}}\Pi\left(\frac{\textbf{r}-\textbf{r}_{i}}{\delta\textbf{r}_{i}}\right)}=\delta\left(\textbf{r}-\textbf{r}_{i}\right)$.

Numerically solving the Eq. (18) requires two approximations. (1) Electron
density approximation: The functions $\varepsilon\left(\textbf{r}\right)$,
$\beta\left(\textbf{r}\right)$, $\textbf{F}\left(\textbf{r}\right)$ and
$f\left(r\right)$ on the right side of the Eq. (18) are associated with the
electron density $n_{e}$. The four functions need to be calculated before
solving the Eq. (18) to obtain the electron density at time $t+\Delta t$. I
propose two schemes for computing these four functions. (i) Calculate the four
functions using the electron density at time $t$. (ii) First, compute the
electron density at time $t+\Delta t$ using the scheme (i), and then
recalculate the four functions using the new electron density. Scheme (ii)
offers higher accuracy but requires more computational resources. (2) Kernel
function approximation: The kernel function’s domain spans entire r-space,
implying that the convolution should be performed across the entire r-space.
However, one can observe that the kernel function’s value exponentially
decays, allowing for acceptance of a finite domain. I recommend a domain size
$r=4/a$, where the value of the kernel function has decreased by two orders of
magnitude ($\sim e^{-4}$).

The Eqs. (1-4) and (18) serve as the controlled equations of the conv-HPIC
method, enabling simulation of non-quasi-neutral plasmas. The calculation
procedure for a given time step is as follows. (a) Pushing ions: Solve the
motion equations of ions, the Eqs. (1,2), to obtain their new positions and
velocities under the influence of electric and magnetic fields. (b) Computing
ion density and fluid velocity: Map the ion pseudo-particles onto grids to
obtain the ion density and fluid velocity. (c) Calculating electron density:
Solve the convolution equation (18) to obtain the electron density. (d)
Determining electric and magnetic fields: Solve the combination of the
Faraday’s law (3) and the generalized Ohm’s law (4) to compute the electric
and magnetic fields.

##  3 Test Examples

To verify the effectiveness of the conv-HPIC method, I have selected a
classical 1D example, the collisionless electrostatic shock (CES). I compare
the simulation results obtained using the quasi-neutral HPIC, conv-HPIC, and
kinetic PIC methods.

The CES arises from the electrostatic ion-ion acoustic instability in a
supersonic counter-stream plasma system [8, 9, 10, 11], and has been observed
in various plasma experiments [12, 13]. The width of the shock region is
comparable to the Debye length, leading to the breakdown of the quasi-neutral
condition within this region. Hence, the CES provides an ideal test
environment for the non-quasi-neutral plasmas.

I utilize colliding plasmas to generate the CES. The two colliding plasmas
consist of fully ionized hydrogen and have identical parameters. The electron
density is set to $10^{20}$ cm-3, while the electron and ion temperatures are
40 keV and 1 keV, respectively. The colliding speed is 1$\%$ of the speed of
light. The Mach number is
$Ma=\frac{u}{c_{s}}=u/\sqrt{\frac{ZT_{e}}{Am_{p}}}=1.53>1$, confirming the
shock condition. For the simulation, the domain length is set to 100 $\mu$m,
and the colliding plasmas initially make contact at $x$=50 $\mu$m. The total
simulation time is 5 ps. I apply an injector boundary for the particles and an
open boundary for electric and magnetic fields. Note that no additional
heating is applied during the simulation. Therefore, I adopt an equation of
state (EOS) with adiabatic approximation,

$p_{e}n_{e}^{-\gamma}=const.$ (19)

Where the specific heat ratio is $\gamma=\frac{f+2}{f}$, and $f$ is the degree
of freedom. For 1D case, $f$=1, corresponds to $\gamma$=3.

I provide a brief introduction to the quasi-neutral HPIC, conv-HPIC, and
kinetic PIC methods in the context of electrostatic simulations.

(a) Electrostatic quasi-neutral HPIC

In electrostatic case, the magnetic field is ignored, so the generalized Ohm’s
law (4) is simplified as,

$\textbf{E}=\frac{-1}{en_{e}}\nabla p_{e}.$ (20)

Using the EOS (19) and the Clapeyron equation of the ideal gas
$p_{e}=n_{e}k_{B}T_{e}$, the gradient of the electron pressure is,

$\nabla p_{e}=\gamma k_{B}T_{e}\nabla n_{e}.$ (21)

Substituting the Eq. (21) and the quasi-neutral condition $n_{e}\cong n_{i}$
into the Eq. (20), the electrostatic field is $\textbf{E}=-4\pi
e\gamma\lambda_{D}^{2}\nabla n_{i}$, where
$\lambda_{D}=\left(\frac{k_{B}T_{e}}{4\pi n_{e}e^{2}}\right)^{1/2}$ is the
Debye length. Therefore, the controlled equations of the quasi-neutral HPIC
are,

$\displaystyle m_{i}\frac{d\textbf{v}_{i}}{dt}=q_{i}\textbf{E},$ (22)
$\displaystyle\textbf{E}=-4\pi e\gamma\lambda_{D}^{2}\nabla n_{i}.$ (23)

The source code sees the attached MATLAB script file HPIC1D.m.

(b) Electrostatic conv-HPIC

Using the EOS (19), the Eq. (18) is simplified as $n_{e}=n_{i}\ast
f\left(r\right)$, where the coefficient $a^{2}$ in the kernel function
$f\left(r\right)$ is $a^{2}=\frac{1}{\gamma\lambda_{D}^{2}}$. Therefore, the
controlled equations of the conv-HPIC are,

$\displaystyle m_{i}\frac{d\textbf{v}_{i}}{dt}=q_{i}\textbf{E},$ (24)
$\displaystyle\textbf{E}=-4\pi e\gamma\lambda_{D}^{2}\nabla n_{e},$ (25)
$\displaystyle n_{e}=n_{i}\ast f\left(r\right).$ (26)

The source code sees the attached MATLAB script file HPIC1D.m.

(c) Electrostatic kinetic PIC

The kinetic PIC method is a first-principal method, providing highly accurate
simulation results, but requiring expensive computations. Here it is used as a
benchmark to validate other simulation methods. In the kinetic PIC, both
electrons and ions are treated as the pseudo-particles. The charge densities
are obtained by mapping these particles to grids. The electrostatic potential
is then calculated by solving the Poisson’s equation [14]. The controlled
equations of the kinetic HPIC are,

$\displaystyle m_{e}\frac{d\textbf{v}_{e}}{dt}=-e\textbf{E},$ (27)
$\displaystyle m_{i}\frac{d\textbf{v}_{i}}{dt}=q_{i}\textbf{E},$ (28)
$\displaystyle\nabla^{2}\phi=4\pi e\left({n_{e}-n}_{i}\right),$ (29)
$\displaystyle\textbf{E}=-\nabla\phi.$ (30)

The source code sees the attached MATLAB script file ESPIC1D.m.

Figure 2: Electric fields of collisionless electrostatic shock in quasi-
neutral HPIC and kinetic PIC simulations.

Figure 3: Electric fields of collisionless electrostatic shock in conv-HPIC
and kinetic PIC simulations.

Firstly, I compared the quasi-neutral HPIC with the kinetic PIC. The two
results revealed notable differences, as shown in Fig. 2. In the quasi-neutral
HPIC, we observed that the amplitude and frequency of the CES electric field
increased as the grid size decreased, indicating a grid-dependent result. This
is incorrect because the physical result must be grid-independent. This result
indicates that the quasi-neutral HPIC is unable to simulate the non-quasi-
neutral plasmas.

Next, I compared the conv-HPIC with the kinetic PIC. The CES electric field
obtained from the conv-HPIC is agree with that of the kinetic PIC and displays
grid-independent, as shown in Fig. 3. Additionally, I examined the ion
densities in the $x-v_{x}$ space for both methods, and they exhibited
agreement, as illustrated in Fig. 4. These results indicate that the conv-HPIC
method is capable of simulating the non-quasi-neutral plasmas.

Figure 4: Ion densities (arbitrary unit) in $x-v_{x}$ space of kinetic PIC (a)
and conv-HPIC (b) simulations.

The conv-HPIC method offers several advantages over other methods when
simulating non-quasi-neutral plasmas. (1) Avoiding iteration: The convolution
equation (18) can be explicit solved, avoids iterative calculation. (2)
Working on stationary electrostatic issues: The convolution equation (18) is
time-independent, so it can be applied to the stationary issues. Note that the
Eq. (18) can also be applied to the non-stationary issues by considering the
time-dependence of ion density, velocity, electric and magnetic fields. (3)
Suppressing numerical noise: Particle methods often introduce strong numerical
noise during mapping particles to grids. When numerical noise grows rapidly,
it may result in non-physical simulation results. The convolution operation in
the Eq. (18) enables the smoothing of grid quantities, significantly reducing
numerical noise.

For computational efficiency, the conv-HPIC is comparable with the quasi-
neutral HPIC. The total computation of the HPIC methods can be divided into
particle computation and grid computation. Among them, particle computation
dominates. Therefore, although the convolution operation in the conv-HPIC
increases the grid computation, it has little impact on the total computation.
As a result, the computational efficiency of the conv-HPIC far exceeds that of
the kinetic PIC.

For parallelization, both the conv-HPIC and the quasi-neutral HPIC employ the
same scheme. Both methods can be easily implemented using the MPI (Message
Passing Interface) standard. The conv-HPIC requires additional guardian grids
to store grid information from adjacent parallel regions. The size of these
guardian grids depends on the length of the kernel function $f\left(r\right)$.
Although the conv-HPIC involves more grid communication, the main
communication is often focused on the particles. As a result, the total
communications required by both methods are comparable.

##  4 Conclusion

I introduced a method called the conv-HPIC, which is capable of simulating
non-quasi-neutral plasmas. This method utilizes a kernel function in the form
of the Debye screening potential to compute the electron density through
convolution. Subsequently, the electrostatic field is calculated by
substituting the electron density into the generalized Ohm’s law. The method
offers the benefits of avoiding iteration, suppressing numerical noise, and is
suitable for both stationary and non-stationary electrostatic issues. It
demonstrates strong computational efficiency and is easy to parallelize.

##  5 Acknowledgements

This work is supported by the National Science Foundation of China (No.
12147103) and the Fundamental Research Funds for the Central Universities (No.
226-2022-00216).

## References

  * [1] D. Winske, L. Yin, N. Omidi, Hybrid simulation codes: Past, present and future—a tutorial, Space plasma simulation (2003) 136–165.
  * [2] J. L. Jiao, Ion current screening modeling of the ion–weibel instability, The Astrophysical Journal 924(2) (2022) 89.
  * [3] A. S. Lipatov, The hybrid multiscale simulation technology: an introduction with application to astrophysical and laboratory plasmas, Springer Science and Business Media, 2013.
  * [4] D. R. Welch, D. V. Rose, M. E. Cuneo, Integrated simulation of the generation and transport of proton beams from laser-target interaction, Physics of Plasmas 13(6) (2006).
  * [5] H. Cai, X. Yan, P. Yao, Hybrid fluid–particle modeling of shock-driven hydrodynamic instabilities in a plasma, Matter and Radiation at Extremes 6(3) (2021).
  * [6] H. X. Vu, An adiabatic fluid electron particle-in-cell code for simulating ion-driven parametric instabilities, Journal of Computational Physics 124(2) (1996) 417–430.
  * [7] R. Bracewell, P. B. Kahn, The fourier transform and its applications, American Journal of Physics 34(8) (1966) 712.
  * [8] S. S. Moiseev, R. Z. Sagdeev, Collisionless shock waves in a plasma in a weak magnetic field, Journal of Nuclear Energy 5(1) (1963) 43.
  * [9] D. W. Forslund, C. R. Shonk, Formation and structure of electrostatic collisionless shocks, Physical Review Letters 25(25) (1970) 1699.
  * [10] T. N. Kato, H. Takabe, Electrostatic and electromagnetic instabilities associated with electrostatic shocks: Two-dimensional particle-in-cell simulation, Physics of Plasmas 17(3) (2010).
  * [11] W. Zhang, H. Cai, S. Zhu, The formation and dissipation of electrostatic shock waves: the role of ion–ion acoustic instabilities, Plasma Physics and Controlled Fusion 60(5) (2018) 055001.
  * [12] R. J. Taylor, D. R. Baker, H. Ikezi, Observation of collisionless electrostatic shocks, Physical Review Letters 24(5) (1970) 206.
  * [13] J. L. Jiao, S. K. He, H. B. Zhuo, Experimental observation of ion–ion acoustic instability associated with collisionless shocks in laser-produced plasmas, The Astrophysical Journal Letters 883(2) (2019) L37.
  * [14] C. K. Birdsall, A. B. Langdon, Plasma physics via computer simulation, CRC press, 2018.