Development and benchmark simulations of a hybrid fluid particle-in-cell code SHYPIC

1.   Introduction

Plasma-jet-driven magneto-inertial fusion (PJMIF) is an innovative approach to achieving controlled nuclear fusion, a long-sought-after goal for its promise of providing a virtually limitless and clean energy source [1]. In PJMIF, a high-velocity, dense plasma jet is employed to compress a magnetized target to fusion-relevant conditions. This approach exploits the adiabatic heating of the target plasma during compression, potentially leading to a more energy-efficient path to ignition compared to conventional inertial confinement fusion (ICF) [2]. This innovative fusion concept combines aspects of both magnetic confinement fusion (MCF) and ICF, offering the advantages of enhanced plasma confinement and reduced energy losses. The unique combination of features in PJMIF, such as the use of plasma jets for compression and the incorporation of magnetic fields in the target, opens up new possibilities for achieving fusion. However, it also presents a complex set of challenges that require sophisticated simulation tools for their understanding and optimization.

The hybrid fluid-kinetic approach is a numerical technique that combines the advantages of fluid and kinetic models. In this approach, electrons are treated as a fluid, thereby avoiding the need to resolve their fast timescales, while ions are treated kinetically, which allows the capture of non-Maxwellian effects that are crucial in many plasma phenomena [3].

The motivation for developing a hybrid fluid-kinetic code stems from the need for a computationally efficient yet physically accurate tool to simulate the complex dynamics in PJMIF. Kinetic particle-in-cell (PIC) simulations, while more accurate, are computationally intensive and not always feasible for large-scale simulations. On the other hand, fluid codes, though computationally less demanding, often lack the ability to accurately model kinetic effects. The hybrid fluid-kinetic approach strikes a balance between these extremes, providing a computationally manageable yet physically insightful tool [4].

This paper introduces the hybrid fluid-kinetic code called SHYPIC (ShanghaiTech HYbrid PIC code) developed as a baseline framework for investigating fundamental plasma dynamics relevant to PJMIF. While the full simulation of PJMIF processes (e.g., plasma gun operation, jet merging, and liner convergence) will require further development, this initial version focuses on capturing essential ion kinetic effects and validating benchmark plasma behaviors. In this version, electron behavior is modeled using a state equation without solving the electron energy equation, and collisional and non-ideal effects (e.g., resistivity, viscosity, and hyper-viscosity) are not yet incorporated. The following sections detail the principles of the hybrid PIC approach, the governing equations of the code, its implementation and operation, and the verification work carried out to validate its performance. This foundational work lays the groundwork for future enhancements that will progressively extend SHYPIC’s capabilities toward comprehensive PJMIF simulations.

2.   Fundamental principles and control equations

The hybrid fluid-kinetic approach subjects a portion of the plasma component to a relatively expensive kinetic treatment, while using a simplified fluid model for the rest. In our code, the electron component is treated as a fluid. This treatment is based on the assumption that electron motion is fast compared to the timescales of interest, allowing the electron distribution to be approximated as a Maxwellian distribution at each time step. This fluid treatment of the electrons avoids the need to resolve the fast electron timescale, thereby significantly reducing the computational cost. On the other hand, the ion component is treated kinetically, specifically, the PIC method is used to implement the ion evolution process. The PIC method involves tracking the motion of individual or a group of ions under the influence of electric and magnetic fields. This kinetic treatment of ions allows us to capture important non-Maxwellian effects that are crucial in many plasma phenomena, including the formation and compression of the plasma liner in PJMIF.

By combining the fluid treatment of electrons and the kinetic treatment of ions, the hybrid PIC method provides a computationally efficient yet physically accurate tool for simulating the complex dynamics in PJMIF. The governing equations of our code are described in detail below.

In the hybrid PIC model, ions are treated kinetically, following the standard PIC technique. The motion of ion particles is governed by the following equation:

where ${\boldsymbol{E}}_{\mathrm{p}}$​ and ${\boldsymbol{B}}_{\mathrm{p}}$​ are the electric and magnetic fields gathered at the particle, respectively. The ion density $n\mathrm{_i}$ and ion flow velocity ${\boldsymbol{U}}_{\mathrm{i}}$​ on the grid can be obtained by interpolating the particle velocity ${\boldsymbol{v}}_{\mathrm{p}}$​ and position ${\boldsymbol{x}}_{\mathrm{p}}$​. For multiple ion species, the density and flow velocity of each ion species need to be collected separately on the grid.

To eliminate the kinetic effects of electrons, they are treated as an inertialess fluid [5]. The electron momentum equation then becomes:

where ${\boldsymbol{U}}_{\mathrm{e}}$​ represents the electron fluid velocity, ${\boldsymbol{P}}_{\mathrm{e}}$​ represents the electron pressure tensor, $\boldsymbol{J}$ represents the total current, and $\boldsymbol{\eta }$ represents the resistivity. In general, the pressure tensor ${\boldsymbol{P}}_{\mathrm{e}}$ can be simplified to a scalar ${p}_{\mathrm{e}}$​. For simplicity, we ignore the friction between electrons and ions (i.e., the resistive term) in this preliminary version, which is adopted under collisionless or weakly collisional conditions. Future versions will incorporate the resistive term to account for collisional effects expected in PJMIF plasmas.

The quasi-neutrality condition is adopted:

where $e$ is the unit charge and ${n}_{\mathrm{e}}$​ is the electron number density.

The electromagnetic fields are described using the low-frequency approximation of Ampere’s law

and Faraday’s law

In the hybrid PIC model, Faraday’s law (6) is used for the time advancement of the magnetic field [6]. Similar to the two-fluid model, Ampere’s law (5) is substituted into the electron momentum equation (3) and the electron fluid velocity ${\boldsymbol{U}}_{\mathrm{e}}$ is eliminated, yielding

In the above formula, the terms on the right side of the equal sign are the convective term, the electron pressure term, and the Hall term, respectively. The electric field $\boldsymbol{E}$ is calculated using formula (7). It can be seen that the electric field $\boldsymbol{E}$ can be directly calculated from the magnetic field $\boldsymbol{B}$, ion number density ${n}_{\mathrm{i}}$​, and ion flow velocity ${\boldsymbol{U}}_{\mathrm{i}}$​, without the need for time advancement. Finally, we use the adiabatic relation of electrons

to close the above set of equations. The adiabatic relation in equation (8) is used as a closure for the electron fluid model in this baseline framework. For strongly collisional plasmas, however, an electron energy equation that accounts for heat conduction and dissipative processes would be more appropriate, and such extensions will be pursued in future code developments.

It is worth noting that no linearization was used in the derivation of the governing equations. Therefore, our simulation code covers all nonlinear effects of the plasma.

3.   Program workflow and implementation

The following is a brief description of the time advancement scheme used in our SHYPIC code. In the PIC part of our program, the time advancement of ions is performed using the Boris algorithm. The Boris algorithm preserves the phase space volume, thus providing excellent long-term accuracy [7]. Using the Boris algorithm implies that the position of the particle is defined at time $N$, while the velocity of the particle is defined at time $N+1/2$:

The ion number density ${n}_{\mathrm{i}}$ and flow velocity ${\boldsymbol{U}}_{\mathrm{i}}$ on the grid are obtained by linear interpolation of the particle velocity ${\boldsymbol{v}}_{\mathrm{p}}$ and position ${\boldsymbol{x}}_{\mathrm{p}}$. The time advancement of the magnetic field is performed using a Runge–Kutta sub-stepping scheme [8]. Substituting the electric field expression into Faraday’s law, we obtain:

Using a 4th order Runge–Kutta scheme, the magnetic field $\boldsymbol{B}$ can be advanced from time $N$ to time $N+1$ (the latter two parameters are obtained from PIC calculation, and remain constant during the magnetic field sub-steps):

Here, the sub-step time is $\mathrm{\Delta }{t}'=\mathrm{\Delta }t/L$ [9, 10]. Finally, the electric field $\boldsymbol{E}$ at time $N+1$ is calculated from ${\boldsymbol{B}}^{N+1}$, ${\boldsymbol{U}}_{\mathrm{i}}^{N+1}$, and ${n}_{\mathrm{i}}^{N+1}$:

In the above equation, ${\boldsymbol{U}}_{\mathrm{i}}^{N+1}$ can be extrapolated from ${\boldsymbol{U}}_{\mathrm{i}}^{N+{1}/{2}}$ and ${\boldsymbol{U}}_{\mathrm{i}}^{N-{1}/{2}}$.

It is worth noting that the scattering calculation of the ion flow velocity ${\boldsymbol{U}}_{\mathrm{i}}$ involves not only the particle velocity ${\boldsymbol{v}}_{\mathrm{p}}$, but also the particle position ${\boldsymbol{x}}_{\mathrm{p}}$. Due to the Boris algorithm, this means that the position and velocity of the particle are staggered in time. Traditionally, ${\boldsymbol{v}}_{\mathrm{p}}^{N+{1}/{2}}$ and ${\boldsymbol{x}}_{\mathrm{p}}^{N}$ are used to calculate the ion flow velocity at time $N+1/2$, i.e., ${\boldsymbol{v}}_{\mathrm{p}}^{N+{1}/{2}},{\boldsymbol{x}}_{\mathrm{p}}^{N}\to {\boldsymbol{U}}_{\mathrm{i}}^{N+{1}/{2}}$, or ${\boldsymbol{v}}_{\mathrm{p}}^{N+{1}/{2}}$ and ${\boldsymbol{x}}_{\mathrm{p}}^{N+1}$ are used, i.e., ${\boldsymbol{v}}_{\mathrm{p}}^{N+{1}/{2}}, {\boldsymbol{x}}_{\mathrm{p}}^{N+1}\to {\boldsymbol{U}}_{\mathrm{i}}^{N+{1}/{2}}$. The advantages and disadvantages of these two methods have been discussed in reference [11]. For the PJMIF jet convergence and compression experiments, the plasma jet velocity is very fast, reaching about 100 km/s. At this time, particle position errors will have a significant impact on the simulation. In view of this feature, the SHYPIC code has created a high-precision ion current collection algorithm. After calculating the particle position ${\boldsymbol{x}}_{\mathrm{p}}^{N+1}$, we use interpolation to calculate the particle position at the half time step ${\boldsymbol{x}}_{\mathrm{p}}^{N+1/2}$. Therefore, we can use ${\boldsymbol{v}}_{\mathrm{p}}^{N+{1}/{2}}$ and ${\boldsymbol{x}}_{\mathrm{p}}^{N+{1}/{2}}$ to calculate a more accurate ion flow velocity ${\boldsymbol{v}}_{\mathrm{p}}^{N+{1}/{2}}, {\boldsymbol{x}}_{\mathrm{p}}^{N+{1}/{2}}\to {\boldsymbol{U}}_{\mathrm{i}}^{N+{1}/{2}}$. Moreover, with ${\boldsymbol{x}}_{\mathrm{p}}^{N+1/2}$, we can use scatter to calculate the ion density ${n}_{\mathrm{i}}^{N+{1}/{2}}$ at time $N+1/2$ more accurately, rather than using the average of the macroscopic quantities ${n}_{\mathrm{i}}^{N}$ and ${n}_{\mathrm{i}}^{N+1}$ to calculate ${n}_{\mathrm{i}}^{N+{1}/{2}}$. This special treatment of the ion flow velocity calculation in the SHYPIC code allows for a more accurate description of the plasma dynamics and effectively avoids artificial heating of particles, especially for the simulation of high-speed plasma jets in PJMIF [12].

The SHYPIC code is designed as an object-oriented application implemented in C++. It adopts a modular structure to enhance flexibility, extensibility, and computational efficiency. The program consists of several key modules, each responsible for a specific part of the simulation process. The Particle Solver Module is responsible for the time advancement of the PIC macroparticles. The Boris algorithm is used to advance the ion velocities and positions, thereby preserving the phase space volume and providing excellent long-term accuracy. This module is similar to the corresponding part in the PPSC PIC code [13]. The Electromagnetic Field Solver Module calculates the magnetic and electric fields at the grid points. The magnetic field is time-advanced using the ion flow velocity and ion density on the grid, by using the Runge–Kutta sub-stepping technique. In addition, the electric field on the grid is calculated using the updated magnetic field and the ion flow velocity and density information. The Interpolation Module converts the discrete particle velocity, position, and weight information into macro quantities on the grid; it also interpolates the magnetic and electric field information on the grid to the position of the particles. In addition, this module is used to calculate the ion velocity distribution function, etc. These modules work together to simulate the complex dynamics in plasmas. Although the present implementation focuses on a collisionless framework for establishing a robust baseline, future versions will incorporate collisional effects to further enhance the physical fidelity, particularly for PJMIF applications.

4.   Verification work

A series of tests have been conducted to verify the implementation of the hybrid PIC algorithm in our SHYPIC code. The simulation results are compared with those of other researchers in the literature and with known theoretical solutions. In this section, we describe the evolution of the total energy in the system under the simulation of a quiet plasma. In addition, we will introduce the simulation of the dispersion relation of plasma waves, as well as the verification simulation of its excitation and propagation. In the following simulation studies, the simulations are carried out in 1D ($x$-axis) or 2D ($x{\text{-}}y$ plane) slab geometry. All boundaries in the simulations are set to periodic boundary conditions. The ions in the simulations are set to be spatially uniform hydrogen ions with a Maxwellian distribution. We will introduce the simulation results of the three cases in detail in the following sections. All simulation results mentioned below have passed the convergence test; that is, when decreasing the time step and the grid width, or increasing the number of macroparticles, the same simulation results will be obtained.

4.1   Quiet plasma test

The first benchmark case is the quiet plasma test, which involves simulating the time evolution of a uniform, undisturbed plasma. This test is an important benchmark to verify the accuracy and reliability of our SHYPIC code, and is a necessary validation step before simulating the PJMIF experiments. The results of this test should show only small statistical fluctuations and should preserve energy conservation over long simulation times. In the simulations, we use the same simulation parameters as in references [14, 15], i.e., 16 cells with 16 particles per cell for 1D simulations, and 64×64 cells with 32 particles per cell for 2D simulations. The grid width $\mathrm{\Delta }x=0.5{d}_{\mathrm{i}}$, where ${d}_{\mathrm{i}}=\dfrac{{v}_{\mathrm{A}}}{{\mathrm{\Omega }}_{\mathrm{i}}}$ is the ion inertial length. The time step $\mathrm{\Delta }t=0.1{\mathrm{\Omega }}_{\mathrm{i}}^{-1}$, where ${\mathrm{\Omega }}_{\mathrm{i}}=\dfrac{{q}_{\mathrm{i}}B}{{m}_{\mathrm{i}}}$ is the ion cyclotron frequency. The ion and electron temperatures are given by ${\beta }_{\mathrm{i}}=1$ and ${\beta }_{\mathrm{e}}=0$, where $\beta =\dfrac{n{k}_{\mathrm{B}}T}{{B}^{2}/2{\mu }_{0}}$. In our simulation, ${\beta }_{\mathrm{e}}=0$ means that the initial electron temperature ${T}_{\mathrm{e}0}$ is set to 0, which follows the parameter settings in references [14, 15] to facilitate direct comparison.

Figures 1 and 2 show the variation of particle energy, electromagnetic field energy, and total energy with time in the 1D and 2D simulation cases, respectively. In the figure, the red dashed line and green dashed line represent the particle energy and field energy, respectively, and the blue solid line represents total energy. It can be seen that the total energy of the system changes very little during the entire simulation process, which verifies the accuracy of the program from the perspective of energy conservation. We calculated the error of the total energy of the plasma and electromagnetic field relative to the initial moment of the simulation when the simulation time is $t=100{\mathrm{\Omega }}_{\mathrm{i}}^{-1}$ and $t=300{\mathrm{\Omega }}_{\mathrm{i}}^{-1}$, and compared the results with those of references [14, 15], see table 1. It can be seen that our results are better than those in reference [14] and comparable to those in reference [15].

Figure  1.  1D quiet plasma test. In the figure, the red and the green dashed lines represent particle energy and field energy, respectively, and the blue solid line represents the total energy.

DownLoad: Full-Size Img PowerPoint

Figure  2.  2D quiet plasma test. Again, the red and the green dashed lines represent the particle energy and the field energy, respectively, and the blue solid line represents the total energy.

DownLoad: Full-Size Img PowerPoint

This quiet plasma test is a basic verification of our SHYPIC code, showing that it can accurately simulate the behavior of a homogeneous plasma evolving over time. It should be noted that the simulation parameters here (including the time step, the grid width, and the number of macroparticles, etc.) are relatively stringent. The parameters used in actual simulations will produce better results. More benchmark cases will be introduced in the following parts.

Table  1.  Total energy errors for quiet plasma simulations.

	Errors in total energy [14]	Errors in total energy [15]	Errors in total energy here
1D, t $=100{\mathrm{\Omega }}_{\mathrm{i}}^{-1}$	9%	0.9%	0.1%
1D, t $=300{\mathrm{\Omega }}_{\mathrm{i}}^{-1}$	47%	3%	0.2%
2D, t $=100{\mathrm{\Omega }}_{\mathrm{i}}^{-1}$	2.6%	0.9%	−0.5%
2D, t $=300{\mathrm{\Omega }}_{\mathrm{i}}^{-1}$	14%	3%	−1.5%

 | Show Table

DownLoad: CSV

4.2   Dispersion relation test

The second benchmark case we proposed is a study of the dispersion relation of plasma waves. In this case, we obtained the dispersion relations of magnetosonic waves and whistler/ion cyclotron waves through simulation, and analyzed the consistency of the simulation results with theoretical expectations. The simulations were carried out in a homogeneous plasma. There was no external perturbation in the simulation, and the eigenmodes in the plasma were excited only by the noise of the PIC model. The parameters used in the simulation are as follows: the number of cells is 512, with 200 particles per cell; the grid width ${\mathrm{\Delta }x=0.1d}_{\mathrm{i}}$, where ${d}_{\mathrm{i}}$ is the ion inertial length; the time step $\mathrm{\Delta }t=0.005{\mathrm{\Omega }}_{\mathrm{i}}^{-1}$, where ${\mathrm{\Omega }}_{\mathrm{i}}$ is the ion cyclotron frequency; the ion temperature is the same as the electron temperature, $T_{\mathrm{e}}=T_{\mathrm{i}}=100\mathrm{\ e}\mathrm{V}$; the plasma density is $1.0\times10^{19\ }\mathrm{m}^{-3}$; and the external magnetic field is $B_0=1.0\mathrm{\ T}$.

We simulated the external magnetic field along the $z$ direction ($\boldsymbol{k}\perp \boldsymbol{B}$) and along the $x$ direction ($\boldsymbol{k}\parallel \boldsymbol{B}$), respectively. By performing Fourier transform on the simulated electromagnetic field data, the distribution of physical quantities was transformed from $t{\text{-}}x$ space to $\omega {\text{-}}k$ space, thereby obtaining the dispersion relation of the simulated plasma waves. Figure 3 shows the dispersion relation curve obtained by simulation when $\boldsymbol{k}\perp \boldsymbol{B}$, and the red dashed line in the figure represents the theoretical solution of the dispersion relations of magnetosonic waves. Figure 4 shows the simulated dispersion relation curve when $\boldsymbol{k}\parallel \boldsymbol{B}$, and the red dashed line and the black dashed line represent the theoretical solutions of the dispersion relations of whistler waves and ion cyclotron waves, respectively. It can be seen that the dispersion relation curve obtained by SHYPIC simulation is very consistent with the theoretical solution. The theoretical dispersion relations in the figure can be obtained from references [16, 17].

Figure  3.  The dispersion relation curve obtained by simulation when $\boldsymbol{k}\perp \boldsymbol{B}$. In the figure, the red dashed line represents the theoretical solution of the magnetosonic wave dispersion relation.

DownLoad: Full-Size Img PowerPoint

Figure  4.  The simulated dispersion relation curve when $\boldsymbol{k}\parallel \boldsymbol{B}$. In the figure, the red dashed line and the black dashed line represent the theoretical solutions of the whistler wave and ion cyclotron wave dispersion relations, respectively.

DownLoad: Full-Size Img PowerPoint

4.3   External perturbation test

The third benchmark case involves the excitation and propagation of plasma waves by introducing a periodic external perturbation. This simulation is also performed in a homogeneous plasma with the same plasma parameters as in the previous case. The main difference here is that we introduce a small periodic perturbation of the form $\delta =$ $A\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{0}t\right)\mathrm{exp}\left[-\dfrac{{\left(x-{p}_{x}\right)}^{2}}{{2w}_{x}^{2}}\right]\mathrm{exp}\left[-\dfrac{{\left(y-{p}_{y}\right)}^{2}}{{2w}_{y}^{2}}\right]\left[1-\mathrm{exp}\left(-\dfrac{{t}^{2}}{{t}_{\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{x}}^{2}}\right)\right]$, where A is the perturbation amplitude, ${\omega }_{0}$ is the frequency of the perturbation source, ${p}_{x}$ and ${p}_{y}$ are the positions of the perturbation source in the $x$ and $y$ directions, ${w}_{x}$ and ${w}_{y}$ are the widths of the source in the $x$ and $y$ directions, respectively, and ${t}_{\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{x}}$ is the ramp-up factor [18]. In the following simulation, the frequency of the driving source is set to ${\omega }_{0}=3{\mathrm{\Omega }}_{\mathrm{i}}$, where ${\mathrm{\Omega }}_{\mathrm{i}}$ is the ion cyclotron frequency.

It is important to note that our dispersion relation tests are performed using two distinct methods. In the dispersion relation test, the intrinsic numerical noise of the PIC method excites a broad spectrum of frequencies and wave numbers, leading to a continuous dispersion relation curve that agrees with theoretical expectations. In contrast, in the external perturbation test, a periodic external perturbation with a fixed frequency (${\omega }_{0}=3{\mathrm{\Omega }}_{\mathrm{i}}$) is applied deliberately, so that only the plasma mode corresponding to ${\omega }_{0}$ is predominantly excited. This results in a single, distinct peak in the $\omega {\text{-}}k$ spectrum. Similar observations have been reported in references [13, 19].

The simulations are performed with external magnetic fields applied along the $z$ direction ($\boldsymbol{k}\perp \boldsymbol{B}$) and along the $x$ direction ($\boldsymbol{k}\parallel \boldsymbol{B}$). Figures 5 and 6 show the wave field distribution in the $x{\text{-}}y$ space for $\boldsymbol{k}\perp \boldsymbol{B}$ and $\boldsymbol{k}\parallel \boldsymbol{B}$, respectively. The wave field distribution data along the $x$ direction at different times are plotted together, to obtain figures 7 ($\boldsymbol{k}\perp \boldsymbol{B}$) and 8 ($\boldsymbol{k}\parallel \boldsymbol{B}$). By performing Fourier transform, the wave field distribution in the $\omega {\text{-}}k$ space under two different magnetic field configurations is obtained, as shown in figures 9 and 10. As in the previous case, the dashed lines in the figure represent the theoretical dispersion relations.

Figure  5.  Wave field distribution in $x{\text{-}}y$ real space when $\boldsymbol{k}\perp \boldsymbol{B}$.

DownLoad: Full-Size Img PowerPoint

Figure  6.  Wave field distribution in $x{\text{-}}y$ real space when $\boldsymbol{k}\parallel \boldsymbol{B}$.

DownLoad: Full-Size Img PowerPoint

As can be seen from figures 5–8, the plasma wave is excited from the position of the perturbation source and gradually propagates outward. In the $\omega {\text{-}}k$ space shown in figures 9 and 10, the distribution of the wave field presents a peak: its frequency $\omega$ corresponds to the frequency of the driving source ${\omega }_{0}$, and its wave number $k$ satisfies the position of the theoretical dispersion relation curve corresponding to ${\omega }_{0}$. These simulation results confirm the correctness of our SHYPIC code in modeling the excitation and propagation of plasma waves in response to external perturbations.

Figure  7.  Wave field distribution in $x{\text{-}}t$ space when $\boldsymbol{k}\perp \boldsymbol{B}$.

DownLoad: Full-Size Img PowerPoint

Figure  8.  Wave field distribution in $x{\text{-}}t$ space when $\boldsymbol{k}\parallel \boldsymbol{B}$.

DownLoad: Full-Size Img PowerPoint

Figure  9.  Wave field distribution in $\omega {\text{-}}k$ space when $\boldsymbol{k}\perp \boldsymbol{B}$. In the figure, the red dashed line represents the theoretical dispersion relation of the magnetosonic wave.

DownLoad: Full-Size Img PowerPoint

Figure  10.  Wave field distribution in $\omega {\text{-}}k$ space when $\boldsymbol{k}\parallel \boldsymbol{B}$. In the figure, the red dashed line and the black dashed line represent the theoretical solutions of the whistler wave and ion cyclotron wave dispersion relations, respectively.

DownLoad: Full-Size Img PowerPoint

5.   Conclusion and future work

In this paper, a hybrid fluid-kinetic code SHYPIC is presented as a baseline framework for simulating plasma dynamics relevant to PJMIF. Our code combines the kinetic PIC technique with a fluid model for electrons, offering a significant reduction in computational cost compared to pure PIC codes, while capturing essential non-Maxwellian effects that are crucial in many plasma phenomena, in contrast to MHD codes. A key innovation of our SHYPIC code is the adoption of a high-precision ion flow velocity collection algorithm, which improves the accuracy of high-speed plasma jet simulations, a critical aspect of PJMIF simulations. A series of validation tests—including quiet plasma tests, dispersion relation simulations, and simulations of plasma wave excitation and propagation via periodic external perturbations—verify the accuracy and reliability of our code by comparison with theoretical solutions and previous simulation works.

It should be noted that the current version of SHYPIC is a preliminary framework that does not yet incorporate collisional effects, resistive or viscous terms in the electron fluid, or other non-ideal physics necessary for a fully comprehensive simulation of PJMIF. Future developments will focus on extending SHYPIC to include these additional physical effects, thereby enabling simulations of key phenomena in PJMIF such as plasma gun operation, plasma liner formation through jet convergence, and subsequent compression of the magnetized target. These future enhancements will provide more accurate support for experimental studies and deeper insights into the complex physical processes in PJMIF, ultimately contributing to the development of fusion as a viable energy source.