Kinetic Plasma Simulations Research Articles

Electrostatic Bursts Generated by the Ion–Ion Acoustic Instability with Solar Wind Plasma Parameters

This study is motivated by recent observations from the Parker Solar Probe (PSP) mission, which have been identified as ion-acoustic waves from 15 to 25 solar radii. These observations reveal characteristic sequences of narrowband, high-frequency bursts exceeding 100 Hz embedded into a slower evolution around 1 Hz, persisting for several hours. To explore the potential role of the ion-acoustic instability (IAI) in these phenomena, we begin by reviewing classical findings on the IAI within the framework of linear kinetic theory. Focusing on proton distributions comprising both a core and a beam component, we analyze the IAI instability range and growth rates within the parameter regime relevant to PSP observations. Our findings indicate that the IAI can indeed occur in this regime, albeit requiring electron-to-core and beam-to-core temperature ratios slightly different from reported values during electrostatic burst detection. Furthermore, employing one-dimensional kinetic plasma simulations, we validate the growth rates predicted by linear theory and observe the saturation behavior of the instability. The resultant nonlinear structures exhibit trapped proton beam populations and oscillatory signatures comparable to those observed, both in terms of timescales and amplitude.

Reconstruction and interpretation of ionization asymmetry in magnetic confinement via synthetic diagnostics

Strong poloidal refueling asymmetry in the DIII-D tokamak is inferred from line radiation measurements. Synthetic diagnostics in neutral transport modeling coupled to gyrokinetic simulations illuminate implications for the plasma flow profile in the scrape-off layer of single-null beam-driven discharges. Recycling occurs primarily either on the inner or outer divertor legs, depending on the toroidal magnetic field direction. By reversing the toroidal magnetic field, the observed line radiation asymmetry is nearly eliminated or reversed. It is determined that, while relatively simple physics can describe the observed ionization asymmetry, predicting the overall brightness of the hydrogenic Lyman-α signal requires detailed simulation of the plasma and resulting turbulence. To this end, kinetic plasma simulations fully coupled to comprehensive neutral transport calculations—a novel capability—provide first-principles reproduction of Lyman-α observations on DIII-D.

Accelerating particle-in-cell kinetic plasma simulations via reduced-order modeling of space-charge dynamics using dynamic mode decomposition.

We present a data-driven reduced-order modeling of the space-charge dynamics for electromagnetic particle-in-cell (EMPIC) plasma simulations based on dynamic mode decomposition (DMD). The dynamics of the charged particles in kinetic plasma simulations such as EMPIC is manifested through the plasma current density defined along the edges of the spatial mesh. We showcase the efficacy of DMD in modeling the time evolution of current density through a low-dimensional feature space. Not only do such DMD-based predictive reduced-order models help accelerate EMPIC simulations, they also have the potential to facilitate investigative analysis and control applications. We demonstrate the proposed DMD-EMPIC scheme for reduced-order modeling of current density and speedup in EMPIC simulations involving electron beam under the influence of magnetic field, virtual cathode oscillations, and backward wave oscillator.

Interplay between accelerated protons, x rays and neutrinos in the corona of NGC 1068: Constraints from kinetic plasma simulations

Interplay between accelerated protons, x rays and neutrinos in the corona of NGC 1068: Constraints from kinetic plasma simulations

An unstructured body-of-revolution electromagnetic particle-in-cell algorithm with radial perfectly matched layers and dual polarizations

A novel electromagnetic particle-in-cell algorithm has been developed for fully kinetic plasma simulations on unstructured (irregular) meshes in complex body-of-revolution geometries. The algorithm, implemented in the BORPIC++ code, utilizes a set of field scalings and a coordinate mapping, reducing the Maxwell field problem in a cylindrical system to a Cartesian finite element Maxwell solver in the meridian plane. The latter obviates the cylindrical coordinate singularity in the symmetry axis. The choice of an unstructured finite element discretization enhances the geometrical flexibility of the BORPIC++ solver compared to the more traditional finite difference solvers. Symmetries in Maxwell's equations are explored to decompose the problem into two dual polarization states with isomorphic representations that enable code reuse. The particle-in-cell scatter and gather steps preserve charge-conservation at the discrete level. Our previous algorithm (BORPIC+) discretized the E and B field components of TEϕ and TMϕ polarizations on the finite element (primal) mesh [1,2]. Here, we employ a new field-update scheme. Using the same finite element (primal) mesh, this scheme advances two sets of field components independently: (1) E and B of TEϕ polarized fields, (Ez,Eρ,Bϕ) and (2) D and H of TMϕ polarized fields, (Dϕ,Hz,Hρ). Since these field updates are not explicitly coupled, the new field solver obviates the coordinate singularity, which otherwise arises at the cylindrical symmetric axis, ρ=0 when defining the discrete Hodge matrices (generalized finite element mass matrices). A cylindrical perfectly matched layer is implemented as a boundary condition in the radial direction to simulate open space problems, with periodic boundary conditions in the axial direction. We investigate effects of charged particles moving next to the cylindrical perfectly matched layer. We model azimuthal currents arising from rotational motion of charged rings, which produce TMϕ polarized fields. Several numerical examples are provided to illustrate the first application of the algorithm.

Hybrid particle-spectral method for kinetic plasma simulations

A hybrid model for numerical solutions of the Vlasov–Poisson equations is presented, which blends spectral and particle approaches. The model splits the distribution function for plasma species into both spectral and particle representations in the velocity space to combine the advantages of each approach. The spectral representation leverages asymmetrically weighted Hermite basis, whereas the particle representation leverages the particle-in-cell method. Configuration phase space is decomposed with the Fourier method, which is well suited for periodic problems. We derive conservation equations for mass, momentum, and energy for the proposed combined method. It is shown that the coupling error between the two methods is absent in the semi-discrete setting (not taking into account time discretization). Finally, numerical test cases are presented simulating a weak electron beam interaction with plasma, leading to beam–plasma instability. The initially localized electron beam evolved into a highly non-equilibrium distribution function in the velocity space. A small growth rate and the resonance nature of instability make it difficult to obtain accurate solutions for purely particle methods due to noise, which falls as ∼1/Np with a number of particles. At the same time, purely spectral methods may require a large number of modes to capture the highly non-equilibrium state of the evolved beam. We show that the hybrid method is well suited for such problems: it reproduces the linear stage as well as nonlinear dynamics with sufficient accuracy using a highly non-equilibrium distribution function.

Electron Injection via Modified Diffusive Shock Acceleration in High-Mach-number Collisionless Shocks

The ability of collisionless shocks to efficiently accelerate nonthermal electrons via diffusive shock acceleration (DSA) is thought to require an injection mechanism capable of preaccelerating electrons to high enough energy where they can start crossing the shock front potential. We propose, and show via fully kinetic plasma simulations, that in high-Mach-number shocks electrons can be effectively injected by scattering in kinetic-scale magnetic turbulence produced near the shock transition by the ion Weibel, or current filamentation, instability. We describe this process as a modified DSA mechanism where initially thermal electrons experience the flow velocity gradient in the shock transition and are accelerated via a first-order Fermi process as they scatter back and forth. The electron energization rate, diffusion coefficient, and acceleration time obtained in the model are consistent with particle-in-cell simulations and with the results of recent laboratory experiments where nonthermal electron acceleration was observed. This injection model represents a natural extension of DSA and could account for electron injection in high-Mach-number astrophysical shocks, such as those associated with young supernova remnants and accretion shocks in galaxy clusters.

Radio Emission and Electric Gaps in Pulsar Magnetospheres

The origin of pulsar radio emission is one of the old puzzles in theoretical astrophysics. In this Letter, we present a global kinetic plasma simulation that shows from first principles how and where radio emission can be produced in pulsar magnetospheres. We observe the self-consistent formation of electric gaps that periodically ignite electron-positron discharge. The gaps form above the polar cap and in the bulk return current. Discharge of the gaps excites electromagnetic modes, which share several features with the radio emission of real pulsars. We also observe the excitation of plasma waves and charge bunches by beam instabilities in the outer magnetosphere. Our numerical experiment demonstrates that global kinetic models can provide deep insight into the emission physics of pulsars and may help interpret their multiwavelength observations.

Phase space dynamics of unmagnetized plasmas: Collisionless and collisional regimes

Eulerian electrostatic kinetic simulations of unmagnetized plasmas (kinetic electrons and motionless protons) with high-frequency equilibrium perturbations have been employed to investigate the phase space free energy transfer across spatial and velocity scales, associated with the resonant interaction of electrons with the self-induced electric field. Numerical runs cover a wide range of collisionless and weakly collisional plasma regimes. An analysis technique based on the Fourier–Hermite transform of the particle distribution function allows to point out how kinetic processes trigger the free energy cascade, which is instead inhibited at finer scales when collisions are turned on. Numerical results are presented and discussed for the cases of linear wave Landau damping, nonlinear electron trapping, and bump-on-tail and two-stream instabilities. A more realistic situation of turbulent Langmuir fluctuations is also discussed in detail. Fourier–Hermite transform shows a free energy spread, highly conditioned by collisions, which involves velocity scales more quickly than the spatial scales, even when nonlinear effects are dominant. This results in anisotropic spectra whose slopes are compatible with theoretical expectations. Finally, an exact conservation law has been derived, which describes the time evolution of the free energy of the system, taking into account the collisional dissipation.

An extended Vlasov–Fokker–Planck approach for kinetic simulations of laser plasmas

Vlasov–Fokker–Planck simulation codes occupy an important niche in modeling laser-produced plasmas, since they are well suited to studying the effect of collisions on electron kinetic phenomena, especially energy transport. One of the most important elements of energy transport is the absorption of laser light by the plasma; however, simulating this in detail requires resolving oscillations of the laser light, whose characteristic timescale is orders of magnitude shorter than the simulation time needed to study transport physics. For this reason, most Vlasov–Fokker–Planck codes used to study electron transport in laser plasmas rely on simplified models of the laser–plasma coupling. Their underlying assumptions nominally preclude their use for modeling laser light having short-scale structure in space or time, such as broadband lasers. In this work, we derive a more general computational framework suitable for arbitrarily structured laser fields. Our approach is based on an extended set of Vlasov–Fokker–Planck equations that separately solve for the low- and high-frequency plasma response. We implement these extended Vlasov–Fokker–Planck equations in the spherical harmonic code K2 and demonstrate the performance of the method on several laser absorption test problems, with particular attention to the judicious selection of time steps, time integrators, and spherical harmonic truncation, according to the intensity and spectrum of the laser light under consideration. Comparison with the widely used Langdon absorption operator shows the Langdon operator performs remarkably well for predicting laser heating in the simple cases considered here, even in situations that would seem to violate its underlying assumptions.

Plasma‐enhanced evaporation and its impact on plasma properties and gaseous chemistry in a pin‐to‐water pulsed discharge

AbstractThe plasma in contact with liquids has led to various novel applications such as plasma biomedicine, material synthesis, and so on. However, the phenomenon of evaporation under plasma treatment and its impact on plasma–liquid interactions has a limited understanding. In this study, the spatially and temporally resolved behavior of water vapor production and its induced influences on plasma properties and gaseous chemistry were studied in detail in an atmospheric pressure pin‐to‐water pulsed He discharge. Diagnostic methods such as laser‐induced fluorescence (LIF) and high‐resolution optical emission spectroscopy (OES) were applied to determine the water vapor and OH radical densities, as well as key plasma parameters such as the gas temperature and electron density. It shows that the physicochemical properties of plasma vary among different discharge regions due to evaporation behavior stimulated during the pulsed discharge‐on phase. In addition, using simulation based on the experimental data, the mechanisms of how water vapor affects the observed spatiotemporal behaviors of OH radicals in different discharge regions are understood. Compared to the pin‐anode and liquid‐cathode sheath regions, proper electron parameters such as density and temperature, as well as water vapor density in the plasma‐positive column, significantly enhance the production of reactive OH radical through the dominant path of electron‐stimulated H2O dissociation. However, higher levels of electron parameters in the intense discharge region near the positive‐pin boundary enhance OH dissociation and finally result in the hollow distribution of OH density. From the global kinetic plasma simulation, the production of reactive hydroxide species playing key roles in plasma medicine treatments, such as O, H, HO2, H2O2, and hydrated ions including H+(H2O)4 and H+(H2O)5, are promoted noticeably as a result of the enhanced water evaporation process.

Efficient parallelization for 3d-3v sparse grid Particle-In-Cell: Shared memory architectures

Particle-In-Cell (PIC) schemes are ones of the most broadly used numerical methods in kinetic simulation of plasmas. The contribution of the present paper is dedicated to the introduction of novel parallelization strategies specific to shared memory architectures tailored for Particle-In-Cell methods implementing sparse grid reconstructions. These strategies operate the different parallelisms genuine to Sparse-PIC methods to obtain speed-up exceeding 100 on 128 cores using a strong scaling study of a 3d-3v simulation. On top of that, substantial gains (two orders of magnitude) are introduced for sequential as well as parallel computation of the electric field thanks to the hierarchization procedure. It consists in decomposing the information carried by sparse grids onto hierarchical basis functions, entailing a significantly reduced number of operations.

Axisymmetric Pulsar Magnetosphere Revisited

We present a global kinetic plasma simulation of an axisymmetric pulsar magnetosphere with self-consistent e ± pair production. We use the particle-in-cell method and log-spherical coordinates with a grid size 4096 × 4096. This allows us to achieve a high voltage induced by the pulsar rotation and investigate pair creation in a young pulsar far from the death line. We find the following: (1) The energy release and e ± creation are strongly concentrated in the thin, Y-shaped current sheet, with a peak localized in a small volume at the Y-point. (2) The Y-point is shifted inward from the light cylinder by ∼15% and “breathes” with a small amplitude. (3) The dense e ± cloud at the Y-point is in ultrarelativistic rotation, which we call superrotation, because it exceeds corotation with the star. The cloud receives angular momentum flowing from the star along the poloidal magnetic field lines. (4) Gamma-ray emission peaks at the Y-point and is collimated in the azimuthal direction, tangent to the Y-point circle. (5) The separatrix current sheet between the closed magnetosphere and the open magnetic field lines is sustained by the electron backflow from the Y-point cloud. Its thickness is self-regulated to marginal charge starvation. (6) Only a small fraction of dissipation occurs in the separatrix inward of the Y-point. A much higher power is released in the equatorial plane, including the Y-point where the created dense e ± plasma is spun up and intermittently ejected through the nozzle between the two open magnetic fluxes.

Quadrature-based moment methods for kinetic plasma simulations

Quadrature-based moment methods (QBMM) are applied to rarefied gas flow and plasma physics problems represented by the Boltzmann equation or the Vlasov-Poisson system of equations. QBMM are a computationally advantageous alternative to both direct Eulerian and Lagrangian particle solvers, able to provide a noise-free solution to the Boltzmann equation and capture non-equilibrium velocity distribution functions (VDF) without excessive computational cost. To provide a closure of the moment equations, the VDF is assumed to be represented by the sum of weighted kernel functions that are placed at given velocity abscissas. Three QBMM approaches, QMOM (Dirac delta kernels), HyQMOM (hyperbolic QMOM), and EQMOM (Gaussian kernels), are applied to canonical one-dimensional rarefied gas flow and plasma physics problems. A new algorithm for EQMOM that is able to treat pathological VDFs and provide control over large abscissas for more than two nodes is developed. QMOM and EQMOM are suitable for problems where the flow does not have shock-like features (for QMOM) or a discontinuity in the VDF (for EQMOM) and where the VDF is far from non-equilibrium and features strong phase-mixing. HyQMOM is generally applicable to all test problems and yields good representations of both the moments and the VDF. HyQMOM also provides significant wall-time improvement (factor of 10–100) compared to other direct Eulerian solvers.

A charge-momentum-energy-conserving 1D3V hybrid Lagrangian–Eulerian method for Vlasov–Maxwell system

Here we propose hybrid Lagrangian–Eulerian (HLE) method for kinetic simulations of plasmas. The HLE method solves advection equations of shape functions unlike particle-in-cell (PIC) method. Although the PIC method cannot preserve the conservation laws of momentum and energy simultaneously, the HLE method can resolve this issue while the particle velocity is formulated by the Lagrangian description. The HLE method discretizes the advection equations and Maxwell's equations with the same scheme, so the momentum and energy exchange between charged particles and electromagnetic field exactly balances in discrete level. Numerical experiments show that linear growth rates of two-stream and Weibel instabilities are reproduced well, and global conservation is exactly preserved except for round-off errors although the velocity space is discretized by the Lagrangian description like the PIC method.

Slowdown of interpenetration of two counterpropagating plasma slabs due to collective effects.

The nonlinear evolution of electromagnetic instabilities driven by the interpenetration of two e^{-},e^{+} plasma clouds is explored using ab initio kinetic plasma simulations. We show that the plasma clouds slow down due to both oblique and Weibel generated electromagnetic fields, which deflect the particle trajectories, transferring bulk forward momentum into transverse momentum and thermal velocity spread. This process causes the flow velocity v_{inst} to decrease approximately by a factor of sqrt[1/3] in a time interval Δt_{αB}ω_{p}∼c/(v_{fl}sqrt[α_{B}]), where α_{B} is the magnetic equipartition parameter determined by the nonlinear saturation of the instabilities, v_{fl} is the initial flow speed, and ω_{p} is the plasma frequency. For the α_{B} measured in our simulations, Δt_{αB} is close to 10 times the instability growth time. We show that as long as the plasma slab length L>v_{fl}Δt_{αB}, the plasma flow is expected to slow down by a factor close to sqrt[1/3].

Kinetic simulations of collision-less plasmas in open magnetic geometries *

Laboratory plasmas in open magnetic geometries can be found in many different applications such as (a) scrape-of-layer (SOL) and divertor regions in toroidal confinement fusion devices, (b) linear divertor simulators, (c) plasma-based thrusters and (d) magnetic mirrors etc. A common feature of these plasma systems is the need to resolve, in addition to velocity space, at least one physical dimension (e.g. along flux lines) to capture the relevant physics. In general, this requires a kinetic treatment. Fully kinetic particle-in-cell (PIC) simulations can be applied but at the expense of large computational effort. A common way to resolve this is to use a hybrid approach: kinetic ions and fluid electrons. In the present work, the development of a hybrid PIC computational tool suitable for open magnetic geometries is described which includes (a) the effect of non-uniform magnetic fields, (b) finite fully-absorbing boundaries for the particles and (c) volumetric particle sources. Analytical expressions for the momentum transport in the paraxial limit are presented with their underlying assumptions and are used to validate the results from the PIC simulations. A general method is described to construct discrete particle distribution functions in a state of mirror-equilibrium. This method is used to obtain the initial state for the PIC simulation. Collisionless simulations in a mirror geometry are performed. The results show that the effect of magnetic compression is correctly described and momentum is conserved. The self-consistent electric field is calculated and is shown to modify the ion velocity distribution function in a manner consistent with analytic theory. Based on this analysis, the ion distribution function is understood in terms of a loss-cone distribution and an isotropic Maxwell-Boltzmann distribution driven by a volumetric plasma source. Finally, the inclusion of a Monte Carlo based Fokker-Planck collision operator is discussed in the context of future work.

Simulation investigations on stepwise penetration of virtual anode in pseudospark discharge

Pseudospark discharge is a special low pressure discharge and has been widely used in the gas switch and electron beam sources. From experiments, when the work pressure is relatively low, the electron beam current or loop current generated by a pseudospark discharge usually has two or more peaks, which has not been fully explained. In this paper, a single-gap pseudospark discharge model is established using 2D kinetic plasma simulation code VSim to study this phenomenon. According to the simulation results, when the anode voltage is 20 kV and the helium pressure is 100 Pa, the current has two peaks, which is similar to the experimental results, accompanied by the stepwise penetration of the virtual anode. This is mainly related to the formation and disappearance of the potential barrier in the cathode hole region. The formation of the potential barrier is caused by the consumption of ions at the cathode hole, and the disappearance is caused by the increase in electrons in the cathode cavity. By classifying the electrons, it is found that the increased electrons are generated by secondary emission caused by ion bombardment on the wall of the cathode. The simulation results also show that the stepwise penetration of the virtual anode can be suppressed or eliminated by increasing the working gas pressure, the secondary electron yield of the cathode material, or the trigger intensity.

Detection and prediction of equilibrium states in kinetic plasma simulations via mode tracking using reduced-order dynamic mode decomposition

A dynamic mode decomposition (DMD) based reduced-order model (ROM) is developed for tracking, detection, and prediction of kinetic plasma behavior. DMD is applied to the high-fidelity kinetic plasma model based on the electromagnetic particle-in-cell (EMPIC) algorithm to extract the underlying dynamics and key features of the model. In particular, the ability of DMD to reconstruct the spatial pattern of the self electric field from high-fidelity data and the effect of DMD extrapolated self-fields on charged particle dynamics are investigated. An in-line sliding-window DMD method is presented for identifying the transition from transient to equilibrium state based on the loci of DMD eigenvalues in the complex plane. The in-line detection of equilibrium state combined with time extrapolation ability of DMD has the potential to effectively expedite the simulation. Case studies involving electron beams and plasma ball are presented to assess the strengths and limitations of the proposed method.

Koopman Autoencoder Architecture for Current Density Modeling in Kinetic Plasma Simulations

We present a reduced-order model (ROM) based on Koopman autoencoder architecture for prediction of current density in kinetic plasma simulations. We apply the proposed method to forecast the time evolution of the current density derived from high-fidelity electromagnetic particle-in-cell (EMPIC) simulation based on the Maxwell- Vlasov system of equations. We test the Koopman auto encoder architecture for a two-dimensional oscillating electron beam inside a square cavity. Such analysis can be crucial in modeling of complex nonlinear plasma system as well as expediting time-consuming particle-in-cell simulations.