New FDTD scheme tames numerical errors in plasma simulations

The method, presented by Harune Sekido and Takayuki Umeda in the journal Earth, Planets and Space, corrects the current density so that the discrete form of the charge conservation law is satisfied exactly on the computational grid, while also relaxing the restrictive Courant stability condition that normally limits the size of the time step in explicit electromagnetic field solvers.

The FDTD approach, first introduced by Yee in 1966, remains one of the most widely used techniques for electromagnetic analysis. It discretizes Maxwell’s equations with central finite differences in both time and space, and it places the electric and magnetic field components on a staggered grid, known as the Yee grid, so that both fields satisfy Gauss’s law. This structural advantage has kept the method in service for more than half a century, and it has been adopted for plasma kinetic simulations, where the motion of charged particles and the electromagnetic fields they generate must be advanced self-consistently.

The standard second-order scheme, denoted FDTD(2,2), suffers from numerical dispersion: the numerical phase velocity of waves on the grid deviates from the physical speed of light, producing spurious oscillations. Higher-order spatial discretizations such as FDTD(2,4), which uses fourth-order accuracy in space, reduce these phase velocity errors. However, the Courant condition, which bounds the time step for stability, generally becomes more severe as the spatial order and the number of dimensions increase. Implicit schemes such as Crank-Nicolson and ADI-FDTD can sidestep this limit, but they demand larger computational costs because of iterative convergence or matrix inversions.

Building on their earlier work, Sekido and Umeda had previously developed higher-order FDTD schemes that relax the Courant condition up to a Courant number of one, allowing longer time steps and reducing total computational time. In plasma kinetic simulations, where particle computations dominate the cost, a longer time step directly cuts the number of loops for both particles and fields. Yet these previous methods carried a hidden flaw: numerical errors appeared in the charge conservation law, which prevented them from being straightforwardly implemented in plasma kinetic simulations.

The charge conservation problem is not academic. In Particle-In-Cell (PIC) codes, the current density deposited by moving particles must satisfy the continuity equation relating charge and current; if it does not, non-physical electrostatic fields are generated that can contaminate the entire simulation. The charge conservation methods conventionally used in PIC simulations assume second-order spatial accuracy, so when a fourth-order FDTD solver is paired with them, the mismatch shows up immediately. In test simulations with a loop current injected at the center of a two- or three-dimensional domain, the researchers demonstrated that FDTD(2,4) fails to satisfy the charge conservation law, producing visible errors in the electric field in the regions where current is injected, in both two and three dimensions.

The team’s remedy for FDTD(2,4) is a corrected current density. By replacing the raw current density components with corrected versions built from second-degree difference operators, the divergence of the current vanishes identically in the discrete sense, so the charge conservation law is satisfied numerically. Test simulations in two and three dimensions with the corrected current densities showed that the numerical errors in the charge conservation law were suppressed, restoring clean electric field profiles.

Correcting the current density alone was not enough for the more ambitious scheme from the researchers’ previous study, referred to as Sekido24, which adds third-degree difference terms derived from a Taylor expansion of Maxwell’s equations to relax the Courant condition while keeping fourth-order spatial accuracy. Because of those additional third-degree terms, a current density correction for Sekido24 is difficult. The new scheme therefore restructures the third-degree difference terms: instead of combining first-degree differences with second-order accuracy and second-degree differences, the new formulation combines the first-degree difference with fourth-order accuracy and the second-degree difference with second-order accuracy. This preserves the same spatial accuracy as Sekido24 but makes the current density correction tractable, and the corrected current densities are obtained by substituting the FDTD(2,4) correction into the new time-development equations.

The coefficients of the additional difference terms, alpha one and alpha two, were determined by brute-force searches over the dispersion relation for each Courant number, excluding any coefficient set for which the dispersion relation becomes unstable. Two variants emerged. Scheme 1 minimizes the phase velocity error, while Scheme 2 adds the requirement that the numerical frequency increase monotonically across the entire wavenumber space, which guarantees positive group velocities and prevents non-physical wave modes. This monotonicity condition could only be satisfied for Courant numbers up to 0.99 in two dimensions and 0.81 in three dimensions. Comparisons showed that the phase velocity errors of the new schemes follow trends similar to Sekido24, and that at a Courant number of one the numerical errors are smaller than those of the standard FDTD(2,2) at its usual Courant numbers, except at large wavenumbers.

The decisive test came when the new scheme was implemented in a two-and-a-half-dimensional PIC code with fourth-degree particle shape, a relativistic particle integrator, and charge-conserving current deposition. In non-relativistic tests, results with a Courant number of 1.0 were almost consistent with those at 0.5, and the elapsed time for the entire computation dropped from about 221 minutes to about 119 minutes. In relativistic tests with strong electron temperature anisotropy, which excites electromagnetic whistler mode waves, the team found that large Courant numbers triggered strong numerical Cherenkov radiation at wavenumbers where the dispersion relation has a negative group velocity, causing electric field and total energy densities to grow rapidly. With Scheme 2, whose numerical frequency increases monotonically, the growth rate of this spurious radiation was lower, and the increase in electric field energy density was markedly reduced compared with runs using schemes lacking this property or lacking the current density correction altogether.

The study demonstrates that high-order FDTD methods can now be made compatible with the strict charge conservation demands of plasma kinetic simulations, combining relaxed time-step limits, reduced phase velocity errors, and controlled numerical Cherenkov radiation. The authors note that the optimal coefficients were derived under the assumption of isotropic grid sizes, and that anisotropic grids require a new coefficient search using the dispersion relation. For the plasma modeling community, the result opens a path to faster kinetic simulations without sacrificing the physical fidelity that charge conservation protects.

Subject of Research: A numerical method that suppresses charge conservation errors in higher-order FDTD solvers for plasma kinetic simulations

Article Title: Suppression of numerical errors in the charge conservation law for finite-difference time-domain method with higher-order accuracy

Article References: Sekido, H., & Umeda, T. (2026). Suppression of numerical errors in the charge conservation law for finite-difference time-domain method with higher-order accuracy. Earth, Planets and Space, 78(1), Article 228. https://doi.org/10.1186/s40623-026-02540-6

Image Credits: AI Generated

DOI: 10.1186/s40623-026-02540-6

Keywords: FDTD method, charge conservation law, plasma kinetic simulations, Particle-In-Cell, numerical dispersion, Courant condition, current density correction, numerical Cherenkov radiation, Maxwell's equations, fourth-order accuracy, group velocity, computational electromagnetics