Planar Collisionless Shock Simulations with Semi-Implicit Particle-in-Cell Model FLEKS

TL;DR

This work demonstrates that a refined semi-implicit PIC code (FLEKS) can robustly reproduce key kinetic features of collisionless planar shocks at sub-ion scales under heliospheric conditions, including shock structure, foreshock/magnetosonic waves, downstream reconnection, and jets. By coupling GL-ECSIM with diffusion and comoving-current techniques within SWMF, the authors validate 1D and 2D setups across quasi-perpendicular and quasi-parallel geometries, and show that multidimensionality and mass-ratio/resolution critically affect downstream wave physics. The results provide concrete guidance on choosing physical and numerical parameters for integrating kinetic shock processes into global MHD–AEPIC simulations, advancing the goal of bridging large-scale magnetospheric dynamics with small-scale kinetic phenomena. The study also outlines limitations and pathways for future work, including extended MHD anisotropic-pressure solvers and fully 3D global simulations to capture curved shocks and downstream jet behavior more realistically.

Abstract

This study investigates the applicability of the semi-implicit particle-in-cell code FLEKS to heliospheric shock simulations. We examine one- and two-dimensional local planar shock simulations, initialized using MHD states with upstream conditions representative of plasmas in the hypersonic, $β\sim 1$ regime, for both quasi-perpendicular and quasi-parallel configurations. The refined algorithm in FLEKS proves robust, enabling accurate shock simulations with a grid resolution on the order of the electron inertial length $d_e$. Our simulations successfully capture key shock features, including shock structures (foot, ramp, overshoot, and undershoot), upstream and downstream waves (fast magnetosonic, whistler, Alfvén ion-cyclotron, and mirror modes), and non-Maxwellian particle distributions. Crucially, we find that at least two spatial dimensions are critical for accurately reproducing downstream wave physics in quasi-perpendicular shocks and capturing the complex dynamics of quasi-parallel shocks, including surface rippling, shocklets, SLAMS, magnetic reconnection and jets. Furthermore, our parameter studies demonstrate the impact of mass ratio and grid resolution on shock physics. This work provides valuable guidance for selecting appropriate physical and numerical parameters for shock simulations using a semi-implicit PIC method, paving the way for incorporating kinetic shock processes into large-scale collisionless plasma simulations with the MHD-AEPIC model.

Figures (11)

• Figure 1: 1D simulation of a quasi-perpendicular shock with $\theta_{Bn} = 87^\circ$ and $m_i / m_e = 25$. The plasma moments and EM fields are normalized by the upstream quantities. Phase space distributions are normalized by the maximum phase space density in the displayed regions. (a-e) Close-up of plasma quantities and electromagnetic fields near the shock front at $t=20\,\omega_{ci,\text{up}}$. (f) $x-t$ stack plot of the total magnetic field, overlaid with ion density contours. (g,h) Ion and electron $x-v_{x}$ phase space plot, with the $B_z$ profile shown in blue for reference. Ion and electron velocity space distributions are displayed in the downstream region (i-l), downstream right behind the shock front (m-p), and upstream right next to the shock front (q-t).
• Figure 2: 1D simulation of a quasi-parallel shock with $\theta_{Bn} = 30^\circ$ and $m_i / m_e = 25$. The plasma moments and EM fields are normalized by the upstream quantities. Phase space distributions are normalized by the maximum phase space density in the displayed regions. (a-e) Close-up of plasma quantities and electromagnetic fields near the shock front at $t=19\,\omega_{ci,\text{up}}^{-1}$. (f) $x-t$ stack plot of the total magnetic field, overlaid with ion density contours. (g,h) Ion and electron $x-v_{x}$ phase space plot, with the $B_z$ profile shown in blue for reference. (i-l) Downstream ion and electron velocity space distributions. (m-p) Ion and electron velocity space distributions at the shock front. (q-t) Upstream ion and electron velocity space distributions.
• Figure 3: Influence of the ion-to-electron mass ratio $m_i/m_e$ on the normalized magnetic field profile across 1D collisionless shocks at $t=20\,\omega_{ci,\text{up}}^{-1}$. The grid resolution $\Delta x = 1\,d_e$. Panels (a-c) show quasi-perpendicular shocks ($\theta_{Bn}=87^\circ$), while panels (d-f) show quasi-parallel shocks ($\theta_{Bn}=30^\circ$). The x-axes are shifted to align the shock fronts for easier comparison. The spatial extent shown is $50\,d_{i,\text{up}}$ for quasi-perpendicular shocks and $200\,d_{i,\text{up}}$ for quasi-parallel shocks.
• Figure 4: Influence of the grid resolution $\Delta x$ on the normalized electromagnetic field profile across 1D collisionless shocks of $m_i / m_e =25$. Panels (a-h) show quasi-perpendicular shocks ($\theta_{Bn}=87^\circ$) at $t=20\,\omega_{ci,\text{up}}^{-1}$, while panels (i-p) show quasi-parallel shocks ($\theta_{Bn}=30^\circ$) at $t=20\,\omega_{ci,\text{up}}^{-1}$. The x-axes are shifted to align the shock fronts for easier comparison. The spatial extent shown is $50\,d_{i,\text{up}}$ for quasi-perpendicular shocks and $200\,d_{i,\text{up}}$ for quasi-parallel shocks.
• Figure 5: 2D $\theta_{Bn} = 87^\circ$ quasi-perpendicular shock run with $m_i / m_e = 25$ at $t=20\,\omega_{ci,\text{up}}^{-1}$. The plasma moments and EM fields are normalized by the upstream quantities. Phase space distributions are normalized by the maximum phase space density in the displayed regions. (a-c) Zoomed-in view of magnetic field components near the shock front. (d-h) Plasma quantities and electromagnetic fields along the y=0 dashed line indicated in panel a. (i,j) Ion and electron $x-v_{x}$ phase space plots, with the $B_y$ profile shown in blue for reference. (k-n) Ion and electron velocity space distributions in the downstream region. (o-r) Ion and electron velocity space distributions at the shock front. (s-v) Ion and electron velocity space distributions ahead of the shock ramp. The y-extent for taking the velocity distributions goes from $-2\,d_{i,\text{up}}$ to $2\,d_{i,\text{up}}$.