FPIC: a new Particle-In-Cell code for stationary and axisymmetric black-hole spacetimes

TL;DR

FPIC delivers a comprehensive general-relativistic PIC framework for stationary and axisymmetric black-hole spacetimes, combining a Maxwell solver on a 2D spherical Yee grid with multiple particle pushers, including a novel hybrid integrator that switches between RK4 and a Hamiltonian scheme based on energy-violation criteria. The code, implemented in Fortran with MPI, leverages Kerr-Schild coordinates to regularize the horizon and achieves reproducible, high-precision simulations of both neutral and charged particle dynamics, as well as electrovacuum and plasma-filled magnetospheres. Validation spans neutral geodesics, vacuum Wald fields, and self-consistent plasma dynamics in Wald and split-monopole topologies, yielding physically consistent Meissner effects, Penrose-energy extraction, and Blandford-Znajek power in agreement with analytic expectations. The work demonstrates FPIC’s potential as a robust, extensible tool for kinetic studies near black holes, while outlining future upgrades to 3D, radiative processes, and GPU-accelerated architectures. Overall, FPIC advances reproducible, high-fidelity GRPIC simulations and provides a solid platform for exploring microphysical plasma processes in strong gravity.

Abstract

In this paper we present a newly developed GRPIC code framework called FPIC, providing a detailed description of the Maxwell-equations solver, of the particle ``pushers'', and of the other algorithms that are needed in this approach. We describe in detail the code, which is written in Fortran and exploits parallel architectures using MPI directives both for the fields and particles. FPIC adopts spherical Kerr-Schild coordinates, which encode the overall spherical topology of the problem while remaining regular at the event horizon. The Maxwell equations are evolved using a finite-difference time-domain solver with a leapfrog scheme, while multiple particle ``pushers'' are implemented for the evolution of the particles. In addition to well-known algorithms, we introduce a novel hybrid method that dynamically switches between the most appropriate scheme based on the violation of the Hamiltonian energy. We first present results for neutral particles orbiting around black holes, both in the Schwarzschild and Kerr metrics, monitoring the evolution of the Hamiltonian error across different integration schemes. We apply our hybrid approach, showing that it is capable of achieving improved energy conservation at reduced computational cost. We apply FPIC to investigate the Wald solution, first in electrovacuum and subsequently in plasma-filled configurations. In the latter case, particles with negative energy at infinity are present inside the ergosphere, indicating that the Penrose process is active. Finally, we present the split-monopole solution in a plasma-filled environment and successfully reproduce the Blandford-Znajek luminosity, finding very good agreement with analytical predictions.

Figures (10)

• Figure 1: Spherical Yee grid adopted in FPIC. The red arrows denote the components of the electric fields $\boldsymbol{D}, \boldsymbol{E}$, and the current density $\boldsymbol{J}$, the blue ones denote the magnetic fields $\boldsymbol{B}, \boldsymbol{H}$. The positions on the mesh are labelled by the indices $ir$ in the radial direction and $i\theta$ in the radial direction. Note that at $(ir, i\theta)$ and ($ir+1/2, i\theta+1/2$) the vectors are pointing outwards.
• Figure 2: Volume weighting procedure utilised in FPIC. Reported in the diagram is the geometry of a single cell in a 2D axisymmetric spherical mesh, with a particle located in the cell at position $P(r, \theta)$ (blue point). The volumes $\mathcal{V}$ involved in the interpolation scheme are also reported, accordingly.
• Figure 3: MPI-domain decomposition scheme employed in FPIC. The diagram shows a portion of the computational domain with 9 CPUs, where solid black lines delimit the physical sub-domains assigned to each CPU. Each computational sub-domain extends beyond the physical one by $\Delta r$ (radial direction) and $\Delta \theta$ (angular direction) to allow MPI communications, and CPUs exchange information with the neighbours across the light-blue shaded regions.
• Figure 4: Top row: orbits in the $(r,\varphi)$ plane for neutral particles around black holes. The black disk at the center of each panel represents the horizon; Bottom row: deviation of energy from its original value, $|H-H_0|/H_0$, vs time for each case. From left to right: a few precessing orbits for the "four-leaf" orbit with $E=0.987649$ and $L=3.9$ around a Schwarzschild black hole; orbits around Kerr black hole with $a_*=0.995$, and $\{E, L\}=\{0 .916235, 2\}, \{0 .841423, 1.82\}$, respectively (see the Tab. \ref{['tab:params']} for details). We report the energy errors for the RK4 (blue), the IMR (teal), and the Hamiltonian (magenta) schemes, up to $t=10^4 \, M.$
• Figure 5: Same trajectories for neutral particles as reported in Fig. \ref{['part0']}, computed using our hybrid approaches. Top row: particle trajectories colour-coded according to the integrator adopted for the "Hyb2" scheme (see text for details and Tab. \ref{['tab:params_hyb']} for the numerical parameters). The RK4 integrator (blue) is employed when the particle is far from the black hole, while the Hamiltonian scheme (magenta) is primarily used as the particle enters stronger gravitational-field regions. Bottom row: energy violations obtained using fixed integrator schemes (solid lines) and our hybrid approaches (dashed lines). The total runtime, in seconds, is also reported for each case.