A High-Order WENO-based Staggered Godunov-type Scheme with Constrained Transport for Force-free Electrodynamics

The force-free (or low inertia) limit of magnetohydrodynamics (MHD) can be applied to many astrophysical objects, including black holes, neutron stars, and accretion disks, where the electromagnetic field is so strong that the inertia and pressure of the plasma can be ignored. This is difficult to achieve with the standard MHD numerical methods because they still have to deal with plasma inertial terms even when these terms are much smaller than the electromagnetic terms. Under the force free approximation, the plasma dynamics is entirely determined by the magnetic field. The plasma provides the currents and charge densities required by the dynamics of electromagnetic fields, but these currents carry no inertia. We present a high order Godunov scheme to study such force-free electrodynamics. We have implemented weighted essentially non-oscillatory (WENO) spatial interpolations in our scheme. An exact Riemann solver is implemented, which requires spectral decomposition into characteristic waves. We advance the magnetic field with the constrained transport (CT) scheme to preserve the divergence free condition to machine round-off error. We apply the third order total variation diminishing (TVD) Runge-Kutta scheme for the temporal integration. The mapping from face-centered variables to volume-centered variables is carefully considered. Extensive testing are performed to demonstrate the ability of our scheme to address force-free electrodynamics correctly. We finally apply the scheme to study relativistic magnetically dominated tearing instabilities and neutron star magnetospheres.

💡 Research Summary

The paper presents a comprehensive high‑order numerical framework for solving the force‑free electrodynamics (FFE) equations, which describe plasma dynamics in regimes where electromagnetic energy dominates over plasma inertia and pressure. After a concise introduction to the astrophysical relevance of FFE—such as pulsar magnetospheres, black‑hole jets, and magnetar coronae—the authors formulate the governing equations. The force‑free condition ρₑ E + J × B = 0 is combined with Maxwell’s equations, yielding an explicit expression for the current density J that contains a perpendicular term (E × B) and a parallel term proportional to B. To avoid the computational burden of the parallel term, the authors adopt a simplified current that retains only the perpendicular contribution while enforcing charge conservation through ρₑ = ∇·E.

The core of the algorithm is a Godunov‑type finite‑volume scheme. Primitive variables P = (Bₓ,Bᵧ,B_z,Eₓ,Eᵧ,E_z)ᵀ are evolved on a Cartesian grid. The Jacobian matrices for the x‑ and y‑directions are derived analytically, and their eigenvalues (±1, 0) and eigenvectors are listed explicitly. Using these spectral decompositions, an exact Roe Riemann solver is constructed, providing numerical fluxes at cell interfaces. High‑order spatial accuracy is achieved through a Weighted Essentially Non‑Oscillatory (WENO) reconstruction. The paper details the one‑dimensional WENO‑3 (third‑order) and WENO‑5 (fifth‑order) procedures, including the smoothness indicators, nonlinear weights, and the dimension‑by‑dimension extension to two dimensions.

A critical challenge in electromagnetic simulations is maintaining the divergence‑free condition ∇·B = 0. The authors adopt the Constrained Transport (CT) method on a staggered mesh: magnetic field components are stored at face centers, while electric fields and other variables reside at cell centers. The induction equation is integrated in its integral (area‑averaged) form, and the electric field circulation is computed at cell edges to update face‑centered B. To couple face‑centered and cell‑centered quantities without degrading the overall order, a third‑order mapping technique (following Li 2008) is employed.

Temporal integration uses a third‑order Total Variation Diminishing (TVD) Runge‑Kutta scheme, which combines the spatial operator H(Pⁿ) (including flux differences and source terms) in three stages to achieve strong stability and low numerical dissipation.

The authors validate the method with a suite of standard tests: propagation of Alfvén‑type waves, current‑sheet evolution, cylindrical wave propagation, and multidimensional shock‑tube problems. In all cases the scheme preserves ∇·B to machine precision, exhibits the expected convergence rates, and shows markedly reduced numerical diffusion compared with traditional FDTD implementations.

Finally, two astrophysical applications are demonstrated. First, a rotating neutron‑star magnetosphere is simulated, reproducing the formation of a Y‑point, equatorial current sheet, and realistic field line opening beyond the light cylinder. Second, the relativistic tearing (magnetic reconnection) instability in a highly magnetized plasma is investigated; the growth rates and nonlinear saturation agree with analytical predictions, and the high‑resolution WENO‑CT scheme captures fine current‑sheet structures and plasmoid formation.

In conclusion, the paper delivers a robust, high‑order Godunov‑WENO scheme combined with constrained transport that accurately solves force‑free electrodynamics while preserving the divergence‑free constraint. The methodology is poised to become a valuable tool for studying magnetically dominated astrophysical phenomena, and the authors suggest extensions to full general‑relativistic settings, coupling with particle dynamics, and inclusion of radiative processes.