Relativistic Semi-Implicit Methods (RelSIM)

• RelSIM is a family of semi-implicit schemes for relativistic plasma and radiation dynamics that handle stiff couplings implicitly while preserving explicit advances in less critical parts.
• One formulation combines leap-frog Maxwell updates with Crank-Nicolson particle pushes and Picard iterations to achieve exact discrete energy and local charge conservation.
• Alternative formulations use mass-matrix based linear solves to enable unconditional stability in underresolved astrophysical regimes, reducing numerical dispersion and computational cost.

Searching arXiv for the cited RelSIM papers to ground the article in current records. Relativistic Semi-Implicit Method (RelSIM) denotes a family of semi-implicit discretizations for relativistic plasma and radiation dynamics in which stiff couplings are handled implicitly while retaining explicit structure in the rest of the update. In the particle-in-cell (PIC) literature, the name is associated most directly with relativistic Vlasov-Maxwell solvers designed to relax explicit stability constraints while controlling dispersion and conservation error. One formulation selectively combines a leap-frog advance for Maxwell’s equations with a Crank-Nicolson advance for particle motion and a small number of Picard iterations, yielding exact discrete energy and local charge conservation (Chen et al., 2019). A second formulation expresses the half-step current through particle-dependent mass matrices and reduces each step to one global linear field solve, with no nonlinear iteration and improved behavior in underresolved astrophysical regimes (Bacchini, 2023). The same label also appears, in separate numerical contexts, in semi-implicit schemes for resistive relativistic magnetohydrodynamics and general-relativistic radiation hydrodynamics [(Cordero-Carrión et al., 19 Jan 2025); (Sadowski et al., 2012)].

1. Scope of the term and major formulations

In PIC usage, RelSIM refers to algorithmic hybrids that depart from fully explicit electromagnetic PIC by introducing implicit field-particle coupling only where needed. The 2019 formulation is a second-order, semi-implicit PIC algorithm for the relativistic Vlasov-Maxwell system that simultaneously conserves the discrete total energy exactly, enforces exact local charge conservation, preserves the low numerical dispersion of the standard leap-frog solver for light waves, and supports large time steps up to the CFL limit $c\Delta t/\Delta x \le 1$ with only mild iteration (Chen et al., 2019). It was implemented in a code named iVPIC, based on the Los Alamos National Laboratory VPIC code.

The 2023 formulation presented by Bacchini is also called RelSIM, but its mechanics are different. It is described as relatively simple, free of nonlinear iterations, and only requiring a global linear solve of the field equations; it targets astrophysical plasma simulations and is presented as a code framework ready for production runs (Bacchini, 2023). In that line of work, RelSIM is explicitly contrasted with both standard explicit PIC and the Relativistic Implicit-Moment Method (RelIMM).

A compact comparison clarifies the distinction between the two PIC formulations.

Formulation	Core update	Stated properties
2019 exact-conservation PIC variant	Leap-frog Maxwell update + Crank-Nicolson particle push + Picard iteration	Exact global energy conservation, exact local charge conservation, leap-frog light-wave dispersion (Chen et al., 2019)
2023 mass-matrix PIC variant	$\theta$-scheme Maxwell update + staggered particle update + one global linear field solve	Unconditionally stable, no nonlinear iteration, smaller numerical errors in total energy than explicit PIC in underresolved regimes (Bacchini, 2023)

This multiplicity of usage is important: RelSIM is not a single canonical algorithm but a label applied to several relativistic semi-implicit schemes sharing the goal of reducing stiffness without reverting to fully explicit time stepping.

2. Discrete structure in relativistic PIC RelSIM

The 2019 PIC formulation places the fields on a Yee mesh and advances Maxwell’s equations with a leap-frog discretization. Faraday’s law is split into two half-steps,

$$\frac{B^{n+1/2}-B^n}{\Delta t/2}=-\nabla_h\times E^n,\qquad \frac{B^{n+1}-B^{n+1/2}}{\Delta t/2}=-\nabla_h\times E^{n+1},$$

while Ampère’s law is centered in time,

$$\epsilon_0\frac{E^{n+1}-E^n}{\Delta t}=\frac{1}{\mu_0}\nabla_h\times B^{n+1/2}-j^{n+1/2}.$$

Particles are advanced over the full step with a Crank-Nicolson discretization of the relativistic Lorentz force, using $u\equiv \gamma v$ and the centered velocities

$$v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{\gamma_p^{n+1}+\gamma_p^n},\qquad \bar v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{2\sqrt{1+\left|(u_p^{n+1}+u_p^n)/(2c)\right|^2}}.$$

Because $E^{n+1}$ enters the particle push through $E_p^{n+1/2}$ and the particles determine $j^{n+1/2}$, the method uses a fixed-point Picard iteration, typically with $K=4$–$\theta$0 iterations to reach round-off precision (Chen et al., 2019).

The Bacchini formulation starts from a $\theta$1-scheme for Maxwell’s equations in CGS units,

$\theta$2

combined with staggered particle positions,

$\theta$3

Its key step is to linearize the half-step current as

$\theta$4

with the particle-dependent mass matrix

$\theta$5

Substitution into the field equations yields the global linear system

$\theta$6

Under the usual boundary conditions on $\theta$7, this becomes a single sparse linear system; $\theta$8 follows from Faraday’s law, and particles are then pushed (Bacchini, 2023).

A later reconnection study by Schoeffler et al. presents RelSIM as built on the energy-conserving semi-implicit moment method (ECSIM) of Lapenta et al., extended to relativistic pair and electron-ion plasmas. There the semi-implicit field solve is again expressed as a matrix problem for $\theta$9, followed by a relativistic Boris-type particle update (Schoeffler et al., 14 Sep 2025).

3. Conservation laws, charge continuity, and Gauss’s law

The most distinctive feature of the 2019 formulation is its exact conservation structure. Its particle pusher is constructed so that the discretized Lorentz force does no work through the magnetic rotation term,

$$\frac{B^{n+1/2}-B^n}{\Delta t/2}=-\nabla_h\times E^n,\qquad \frac{B^{n+1}-B^{n+1/2}}{\Delta t/2}=-\nabla_h\times E^{n+1},$$0

which is identified as crucial for exact energy balance. Discrete field and particle energies are defined as

$$\frac{B^{n+1/2}-B^n}{\Delta t/2}=-\nabla_h\times E^n,\qquad \frac{B^{n+1}-B^{n+1/2}}{\Delta t/2}=-\nabla_h\times E^{n+1},$$1

and satisfy

$$\frac{B^{n+1/2}-B^n}{\Delta t/2}=-\nabla_h\times E^n,\qquad \frac{B^{n+1}-B^{n+1/2}}{\Delta t/2}=-\nabla_h\times E^{n+1},$$2

so total energy is conserved exactly to round-off (Chen et al., 2019).

The same algorithm also enforces exact local charge conservation. For each particle,

$$\frac{B^{n+1/2}-B^n}{\Delta t/2}=-\nabla_h\times E^n,\qquad \frac{B^{n+1}-B^{n+1/2}}{\Delta t/2}=-\nabla_h\times E^{n+1},$$3

holds at every cell when first-order (trilinear) splines $$\frac{B^{n+1/2}-B^n}{\Delta t/2}=-\nabla_h\times E^n,\qquad \frac{B^{n+1}-B^{n+1/2}}{\Delta t/2}=-\nabla_h\times E^{n+1},$$4 are used for $$\frac{B^{n+1/2}-B^n}{\Delta t/2}=-\nabla_h\times E^n,\qquad \frac{B^{n+1}-B^{n+1/2}}{\Delta t/2}=-\nabla_h\times E^{n+1},$$5 and a mixed $$\frac{B^{n+1/2}-B^n}{\Delta t/2}=-\nabla_h\times E^n,\qquad \frac{B^{n+1}-B^{n+1/2}}{\Delta t/2}=-\nabla_h\times E^{n+1},$$6–$$\frac{B^{n+1/2}-B^n}{\Delta t/2}=-\nabla_h\times E^n,\qquad \frac{B^{n+1}-B^{n+1/2}}{\Delta t/2}=-\nabla_h\times E^{n+1},$$7 scheme is used for $$\frac{B^{n+1/2}-B^n}{\Delta t/2}=-\nabla_h\times E^n,\qquad \frac{B^{n+1}-B^{n+1/2}}{\Delta t/2}=-\nabla_h\times E^{n+1},$$8, evaluated at the trajectory midpoint. The proof is exact so long as each sub-trajectory segment remains within one cell; multi-cell crossings are handled by splitting the trajectory. Summing over particles yields the discrete continuity equation, and Gauss’s law $$\frac{B^{n+1/2}-B^n}{\Delta t/2}=-\nabla_h\times E^n,\qquad \frac{B^{n+1}-B^{n+1/2}}{\Delta t/2}=-\nabla_h\times E^{n+1},$$9 is therefore preserved exactly if it is true initially (Chen et al., 2019).

The Bacchini formulation makes a different trade-off. It avoids nonlinear iteration by approximating the particle response through the rotation matrix $$\epsilon_0\frac{E^{n+1}-E^n}{\Delta t}=\frac{1}{\mu_0}\nabla_h\times B^{n+1/2}-j^{n+1/2}.$$0 and an approximation $$\epsilon_0\frac{E^{n+1}-E^n}{\Delta t}=\frac{1}{\mu_0}\nabla_h\times B^{n+1/2}-j^{n+1/2}.$$1 extracted from the chosen pusher. Exact energy conservation is therefore sacrificed, but the residual non-conservation is stated to scale as

$$\epsilon_0\frac{E^{n+1}-E^n}{\Delta t}=\frac{1}{\mu_0}\nabla_h\times B^{n+1/2}-j^{n+1/2}.$$2

and in the reported tests the relative energy error remains below $$\epsilon_0\frac{E^{n+1}-E^n}{\Delta t}=\frac{1}{\mu_0}\nabla_h\times B^{n+1/2}-j^{n+1/2}.$$3 over $$\epsilon_0\frac{E^{n+1}-E^n}{\Delta t}=\frac{1}{\mu_0}\nabla_h\times B^{n+1/2}-j^{n+1/2}.$$4 in the worst case (Bacchini, 2023). Charge conservation is not implicit in that form; the implementation guidance is to add a divergence-cleaning step, such as a Hodge projection or Marder corrector, or to implement the charge-conserving variant of Esirkepov (2001) adapted to the semi-implicit moment (Bacchini, 2023).

4. Dispersion, stability, and accuracy regimes

RelSIM was developed partly to address a central drawback of conventional explicit PIC: explicit electromagnetic schemes do not conserve discrete energy exactly and, when heavily underresolved, can accumulate numerical heating or dispersion error. The 2019 formulation retains the leap-frog light-wave dispersion relation,

$$\epsilon_0\frac{E^{n+1}-E^n}{\Delta t}=\frac{1}{\mu_0}\nabla_h\times B^{n+1/2}-j^{n+1/2}.$$5

which is identical to that of classical explicit PIC. As $$\epsilon_0\frac{E^{n+1}-E^n}{\Delta t}=\frac{1}{\mu_0}\nabla_h\times B^{n+1/2}-j^{n+1/2}.$$6, the numerical phase speed tends to unity and no spurious slow-light “Cerenkov” emission arises. Its stability remains tied to the familiar CFL condition $$\epsilon_0\frac{E^{n+1}-E^n}{\Delta t}=\frac{1}{\mu_0}\nabla_h\times B^{n+1/2}-j^{n+1/2}.$$7, but because only Ampère’s law and the particle push are iterated, the field-particle coupling remains only mildly stiff and typically requires about five iterations for convergence (Chen et al., 2019).

The Bacchini formulation instead emphasizes unconditional stability. Explicit PIC must resolve both the Courant limit $$\epsilon_0\frac{E^{n+1}-E^n}{\Delta t}=\frac{1}{\mu_0}\nabla_h\times B^{n+1/2}-j^{n+1/2}.$$8 and the plasma-frequency limit $$\epsilon_0\frac{E^{n+1}-E^n}{\Delta t}=\frac{1}{\mu_0}\nabla_h\times B^{n+1/2}-j^{n+1/2}.$$9, whereas the mass-matrix RelSIM is described as unconditionally stable, with no discrete unstable modes even when $u\equiv \gamma v$0 or $u\equiv \gamma v$1 (Bacchini, 2023). That does not remove accuracy considerations: the same source states that in practice $u\equiv \gamma v$2 is chosen to balance the $u\equiv \gamma v$3 truncation error against linear-solver cost.

The reconnection study by Schoeffler et al. sharpens this point. There, a von Neumann analysis of the semi-implicit update indicates that no spurious high-frequency Langmuir or light-wave modes appear, and that the numerical dispersion relation approaches the exact continuum dispersion for $u\equiv \gamma v$4, even for $u\equiv \gamma v$5. Reported runs remain stable up to $u\equiv \gamma v$6, while measured tearing-mode growth rates remain within approximately $u\equiv \gamma v$7 of theory for $u\equiv \gamma v$8; nearly correct growth rates are emphasized up to $u\equiv \gamma v$9 (Schoeffler et al., 14 Sep 2025). A plausible implication is that the practical operating window is set not by linear instability but by the onset of unacceptable phase and growth-rate error.

5. Implementations, benchmarks, and representative applications

The 2019 exact-conservation RelSIM was implemented in iVPIC. Its efficiency analysis states that each Picard iteration carries one particle push and one sparse Maxwell solve of cost comparable to explicit PIC, that no global nonlinear solve or Jacobian construction is required, and that the iteration is matrix-free and trivially parallel. With $$v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{\gamma_p^{n+1}+\gamma_p^n},\qquad \bar v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{2\sqrt{1+\left|(u_p^{n+1}+u_p^n)/(2c)\right|^2}}.$$0–$$v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{\gamma_p^{n+1}+\gamma_p^n},\qquad \bar v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{2\sqrt{1+\left|(u_p^{n+1}+u_p^n)/(2c)\right|^2}}.$$1, energy errors of approximately $$v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{\gamma_p^{n+1}+\gamma_p^n},\qquad \bar v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{2\sqrt{1+\left|(u_p^{n+1}+u_p^n)/(2c)\right|^2}}.$$2 in single precision are reported, for an overall cost per $$v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{\gamma_p^{n+1}+\gamma_p^n},\qquad \bar v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{2\sqrt{1+\left|(u_p^{n+1}+u_p^n)/(2c)\right|^2}}.$$3 of roughly five times explicit PIC. Test problems include the relativistic two-stream instability, Weibel instability, and laser-plasma instabilities. In the Weibel case, over $$v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{\gamma_p^{n+1}+\gamma_p^n},\qquad \bar v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{2\sqrt{1+\left|(u_p^{n+1}+u_p^n)/(2c)\right|^2}}.$$4 time steps with five Picard iterations per step, the total energy error remains bounded at approximately $$v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{\gamma_p^{n+1}+\gamma_p^n},\qquad \bar v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{2\sqrt{1+\left|(u_p^{n+1}+u_p^n)/(2c)\right|^2}}.$$5 and the discrete Gauss-law error at approximately $$v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{\gamma_p^{n+1}+\gamma_p^n},\qquad \bar v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{2\sqrt{1+\left|(u_p^{n+1}+u_p^n)/(2c)\right|^2}}.$$6; in quasi-1D laser-plasma simulations, instantaneous and time-averaged reflectivities agree closely with explicit VPIC, including the transition from stimulated Raman scattering to Brillouin scattering (Chen et al., 2019).

The Bacchini formulation is explicitly benchmarked against standard explicit PIC and RelIMM. For the test parameters summarized in its Table 1, explicit PIC has maximum stable $$v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{\gamma_p^{n+1}+\gamma_p^n},\qquad \bar v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{2\sqrt{1+\left|(u_p^{n+1}+u_p^n)/(2c)\right|^2}}.$$7, RelIMM $$v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{\gamma_p^{n+1}+\gamma_p^n},\qquad \bar v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{2\sqrt{1+\left|(u_p^{n+1}+u_p^n)/(2c)\right|^2}}.$$8, and RelSIM $$v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{\gamma_p^{n+1}+\gamma_p^n},\qquad \bar v_p^{n+1/2}=\frac{u_p^{n+1}+u_p^n}{2\sqrt{1+\left|(u_p^{n+1}+u_p^n)/(2c)\right|^2}}.$$9; the corresponding per-cycle costs are $E^{n+1}$0 field $E^{n+1}$1 push for explicit, $E^{n+1}$2 linear solve $E^{n+1}$3–$E^{n+1}$4 push iterates for RelIMM, and $E^{n+1}$5 linear solve $E^{n+1}$6 push for RelSIM; and the reported total-energy errors at $E^{n+1}$7 are $E^{n+1}$8–$E^{n+1}$9 for explicit PIC in underresolved conditions, $E_p^{n+1/2}$0–$E_p^{n+1/2}$1 for RelIMM, and $E_p^{n+1/2}$2–$E_p^{n+1/2}$3 for RelSIM (Bacchini, 2023). The same work reports one-dimensional beam-instability tests, an ion-electron shock with underresolved electron scales by factors $E_p^{n+1/2}$4–$E_p^{n+1/2}$5, and two-dimensional relativistic reconnection on a double Harris sheet. In the coarse-resolution reconnection cases, explicit PIC fails for $E_p^{n+1/2}$6, RelIMM exhibits unphysical magnetic-energy depletion, and RelSIM converges in reconnection rate and energy partition by $E_p^{n+1/2}$7 (Bacchini, 2023).

Schoeffler et al. extend the method to large-scale reconnection and particle acceleration. In their linear tearing benchmarks, OSIRIS uses $E_p^{n+1/2}$8 and $E_p^{n+1/2}$9, while RelSIM uses $j^{n+1/2}$0 at the same time step. Under those conditions, OSIRIS requires the finer mesh to reproduce the theoretical growth rate $j^{n+1/2}$1, with energy error around $j^{n+1/2}$2, whereas RelSIM on the eight-times coarser mesh reproduces $j^{n+1/2}$3–$j^{n+1/2}$4 with energy error around $j^{n+1/2}$5 (Schoeffler et al., 14 Sep 2025). Although the implicit field solve is stated to be about $j^{n+1/2}$6 slower per step than an explicit push, coarsening the mesh by a factor of eight in each transverse direction and increasing $j^{n+1/2}$7 by a factor of four yields a nominal $j^{n+1/2}$8 net speed-up, and with parameter tuning the authors report an overall computational advantage of approximately $j^{n+1/2}$9 over an equivalent high-resolution OSIRIS run. In nonlinear two- and three-dimensional reconnection, the simulations produce hard power laws with cutoffs in the energetic particle spectra, with representative fits

$K=4$0

and spectral indices and cutoffs consistent with previous explicit-PIC studies (Schoeffler et al., 14 Sep 2025).

6. Related semi-implicit relativistic solvers beyond PIC

The name RelSIM also appears in other relativistic numerical settings, though these are methodologically distinct from the Vlasov-Maxwell PIC algorithms. In resistive relativistic magnetohydrodynamics, the 2025 minimally implicit Runge-Kutta work presents a semi-implicit treatment of the stiff conductivity term in the electric-field equation. The conserved variables are split as $K=4$1, and only the $K=4$2-dependent part of the $K=4$3 update is treated implicitly and linearly in the conserved electric field. The resulting PIRK(1) and PIRK(2) schemes require one primitive recovery per stage and only a diagonal inversion for the implicit part; no nonlinear solve and no extra recovery loops are needed. The stated computational complexity is therefore $K=4$4, comparable to explicit Runge-Kutta and substantially cheaper than IMEX schemes requiring repeated primitive recoveries inside Newton or fixed-point iterations (Cordero-Carrión et al., 19 Jan 2025).

In general-relativistic radiation hydrodynamics, the KORAL implementation combines M1 closure with an IMEX splitting of radiation-matter coupling. The explicit step advances advection and geometric terms, while the stiff source term $K=4$5 is handled by a local backward-Euler implicit solve for the radiation stress-energy components. Because the source depends nonlinearly on the comoving radiation moments and temperature, each cell requires a $K=4$6 nonlinear system solved by Newton-Raphson, with a fallback analytical approximate implicit update if the solve fails or produces unphysical $K=4$7 (Sadowski et al., 2012). This usage shares the semi-implicit philosophy of handling optically thin and optically thick regimes within one update, but it belongs to radiation hydrodynamics rather than kinetic PIC.

A related antecedent in relativistic PIC is the adjustable-damping direct implicit method implemented in the ELIXIRS code. That method uses a time-centered, semi-implicit discretization with a damping parameter $K=4$8 controlling high-frequency numerical dissipation, Newton linearization of the coupled field-particle update, and high-order particle weighting to mitigate finite-grid heating when $K=4$9 (0906.2912). It is not named RelSIM in the paper title, but it occupies the same methodological lineage of relativistic semi-implicit plasma solvers.

7. Limitations, misconceptions, and ongoing directions

A persistent misconception is that semi-implicit relativistic methods form a single algorithmic template. The record is more heterogeneous. The 2019 exact-conservation PIC method preserves leap-frog light-wave dispersion and exact conservation laws, but it still retains the CFL stability limit and requires Picard iteration (Chen et al., 2019). The Bacchini formulation removes nonlinear iteration and is unconditionally stable, but exact charge conservation is not intrinsic in its current form and exact total-energy conservation is replaced by controlled residual error (Bacchini, 2023). The reconnection study further indicates that unconditional stability does not imply arbitrary accuracy: stable operation is reported up to $\theta$00, but nearly correct tearing growth rates are emphasized only up to about $\theta$01 (Schoeffler et al., 14 Sep 2025).

The limitations and forward directions stated for the Bacchini formulation include recovering exact total-energy conservation via a few Newton-Krylov field-particle nonlinear iterations, subcycling particles in strong-field regions, combining RelSIM with adaptive mesh refinement, and GPU acceleration of mass-matrix assembly and the linear solver (Bacchini, 2023). The same source recommends $\theta$02 up to $\theta$03 to retain $\theta$04 accuracy while benefiting from relaxed $\theta$05 constraints, and states that in fully underresolved electron-scale scenarios $\theta$06–$\theta$07 is acceptable provided it is understood that plasma physics is averaged at sub-cell scales. This suggests that RelSIM is best viewed not as a universal replacement for explicit PIC, but as a controlled underresolution strategy for regimes where macroscopic extent and long duration are otherwise prohibitive.

Across these variants, the common numerical objective is consistent: to bridge relativistic multi-scale dynamics by treating the stiffest couplings implicitly while preserving as much of the explicit algorithm’s structure, locality, and physical fidelity as possible.