Particle-in-Cell (PIC) methods have achieved widespread recognition as simple and flexible approaches to model collisionless plasma physics in fully kinetic simulations of astrophysical environments. However, in many situations the standard PIC algorithm must be extended to include macroscopic effects in microscale simulations. For plasmas subjected to shearing or expansion, shearing-box and expanding-box methods can be incorporated into PIC to account for these global effects. For plasmas subjected to local acceleration in confined regions of space, a leaky-box method can allow closed-box PIC simulations to account for particle escape from the accelerator region. In this work, we review and improve methods to include shearing, expansion, and escape in PIC simulations. We provide the numerical details of how Maxwell's equations and the particle equations of motion are solved in each case, and introduce generalized Boris-like particle pushers to solve the momentum equation in the presence of extra forces. This work is intended to serve as a comprehensive reference for the implementation of shearing-box, expanding-box, and leaky-box algorithms in PIC.
Fully kinetic models of astrophysical plasmas have gained widespread recognition for their capability of reproducing collisionless plasma dynamics. Particle-in-Cell (PIC) methods [Dawson, 1983, Birdsall andLangdon, 1991], in particular, are the most commonly employed numerical tool to conduct fully kinetic simulations. The simplicity and reliability of PIC methods allow for studies of a wide variety of (astrophysical) plasma phenomena, both in the nonrelativistic and relativistic regimes [Lapenta, 2012, Nishikawa et al., 2021].
Standard PIC methods include a Maxwell solver to advance the electric and magnetic fields E, B in time as well as a particle pusher to propagate particles in position-velocity phase space (x, u) according to the charged-particle equations of motion. Moments of the velocity distribution function are collected from the particles onto a numerical grid, to construct the sources of Maxwell’s equations; the electromagnetic fields, in turn, are interpolated from the grid to the particles to determine the particle motion through the Lorentz force.
For plasmas embedded within realistic astrophysical environments, however, the relatively simple algorithm described above may not suffice. Particularly in local boxes, there is a need to couple local plasma evolution to the global-scale configuration of the system and its surrounding medium. This includes, for instance, the effects of background shear in local accretion-disk simulations; the modeling of local plasma expansion in magnetized winds launched from a central source; or the inclusion of diffusive cosmic-ray escape in local boxes. The PIC method needs modifications to be able to model kinetic plasmas in such scenarios, which are the subject of ongoing investigations.
For plasmas that experience shearing forces, the shearing-box method has been used in kinetic simulations primarily to model plasma turbulence in accretion disks. The general approach is to simulate an initially thermal plasma threaded by a weak magnetic field and subjected to background shear mimicking the differential rotation of the disk. This allows the development of the magnetorotational instability (MRI) [Balbus and Hawley, 1991], which nonlinearly develops into sustained turbulence. The shearing, and potentially Coriolis and gravitational forces, are incorporated into Maxwell’s equations and into the particle equations of motion to augment standard PIC algorithms. Initial works using the shearing box focused on 2D systems with or without “shearing coordinates” [Riquelme et al., 2012, Hoshino, 2013, Inchingolo et al., 2018]. Pioneering 3D simulations were also conducted with relatively small system sizes [Hoshino, 2015]. However, only hybrid simulations (discarding electron physics) employed large enough 3D system sizes that could retrieve the expected sustained-turbulence nonlinear state [Kunz et al., 2016]. In fully kinetic runs, this expected behavior was first reported by us [Bacchini et al., 2022, Bacchini et al., 2024] for the pair-plasma case and later extended to incorporate radiation and ion-electron physics [Gorbunov et al., 2025a]. The stratified-density case was also recently simulated for the first time with fully kinetic methods in 2D and 3D [Sandoval et al., 2024].
For plasmas that are subjected to expansion or compression, the expanding-box method has been designed to enable local kinetic simulations of e.g. expanding solar-wind plasma or compressing plasma in accretion disks or the intracluster medium. In PIC runs, expansion/contraction can be included by modifying Maxwell’s equations to take length dilation/contraction into account in the curl terms, and by adding extra forces on the particle equations of motion. Most previous works employing the expanding box have focused on hybrid simulations (treating electrons as a massless fluid) [Hellinger et al., 2019, Bott et al., 2021]; fully kinetic runs, including electron physics, have rarely been explored, with notable exceptions for the expanding solar wind [Elena Innocenti et al., 2019, Innocenti et al., 2019, Micera et al., 2021] and for compressing plasma in accretion disks and galaxy clusters [Sironi andNarayan, 2015, Tran et al., 2023].
In plasmas where particles are subjected to efficient acceleration, e.g. in the surroundings of compact astrophysical objects, a typical approach is to consider a PIC model of a local turbulent box. In such models, turbulence is usually driven by some external force, and the plasma continuously heats up due to the turbulent cascade. Accelerated particles are expected to leave the turbulent region by diffusing out of the source, carrying away energy. To model these dynamics and attain a steady state without secular, unbound energy pile-up, a leaky-box approach can be used [Gorbunov et al., 2025b, Groselj et al., 2026]. In this case, particles can leave the simulation box according to some prescription, potentially reaching a balance between energy injection and e
This content is AI-processed based on open access ArXiv data.