M3D-C1: High-Order Extended-MHD Code

• M3D-C1 is a high-order extended-MHD code that models magnetically confined plasmas using implicit time-stepping and C1-continuous finite elements.
• It supports both axisymmetric and fully 3D simulations, enabling studies of vertical displacement events, sawtooth oscillations, and energetic-particle interactions.
• The platform integrates advanced neoclassical bootstrap-current closures and hybrid kinetic modules to enhance predictive capabilities in fusion research.

M3D-C$^1$ is a high-order finite-element extended-magnetohydrodynamics code for global, nonlinear simulation of magnetically confined plasmas in realistic toroidal geometry. Across the cited literature it appears as a fully 3D-capable code used in axisymmetric and non-axisymmetric settings, with fully implicit or split-implicit time advance, high-order $C^1$-continuous finite elements, and applications spanning vertical displacement events, energetic-particle–MHD interaction, bootstrap-current closure, and internal-kink and sawtooth dynamics in SPARC-like plasmas (Krebs et al., 2019, Liu et al., 2021, Saxena et al., 7 Jul 2025, Wang et al., 2 Apr 2026).

1. Definition, scope, and code family

M3D-C$^1$ is described as a high-fidelity extended-MHD code designed for global, nonlinear simulations of magnetically confined plasmas on transport as well as MHD timescales. In the 2025 bootstrap-current study it is characterized technically as solving single-fluid, and more general extended-MHD, equations on a high-order finite-element grid with $C^1$ continuity, with fully implicit or split-implicit time advance so that the timestep is not limited by the Alfvén Courant condition and can reach transport timescales (Saxena et al., 7 Jul 2025). In the 2019 vertical-displacement-event benchmark it is likewise described as a high-order finite element code that solves the nonlinear time-dependent extended MHD equations and uses a split-implicit time advance to enable simulations over transport time scales (Krebs et al., 2019).

The code family is explicitly broader than a single equilibrium class. One source states that M3D-C$^1$ was originally developed for axisymmetric tokamaks and has been generalized to fully 3D domains so it can evolve strongly non-axisymmetric stellarators (Saxena et al., 7 Jul 2025). Another shows the code run in strictly axisymmetric mode for an $n=0$ benchmark even though it is fully 3D capable (Krebs et al., 2019). A further extension, M3D-C1-K, is developed on top of the M3D-C1 finite-element code to add a kinetic energetic-particle module while retaining the MHD solver (Liu et al., 2021).

A concise view of representative uses is given below.

Study	Problem class	Role of M3D-C$^1$
(Krebs et al., 2019)	Axisymmetric VDE benchmark	One of three large-scale nonlinear MHD codes
(Liu et al., 2021)	EP–MHD hybrid simulation	Base code for M3D-C1-K
(Saxena et al., 7 Jul 2025)	Bootstrap-current closure	Self-consistent Sauter/Redl implementation
(Wang et al., 2 Apr 2026)	SPARC sawtooth modeling	High-fidelity 3D extended-MHD solver

This suggests that M3D-C$^1$ is best understood not as a single-purpose disruption code, but as a general extended-MHD platform onto which additional closures and kinetic couplings can be grafted.

2. Numerical architecture and discretization

The numerical core is consistently described as finite-element based. In the VDE benchmark, M3D-C$^1$ uses triangular wedge finite elements; in axisymmetry these reduce to high-order triangular finite elements in the poloidal $R$–$C^1$0 plane (Krebs et al., 2019). In the bootstrap-current and hybrid-kinetic descriptions, the code is associated with high-order $C^1$1-continuous finite elements and realistic toroidal geometry, including unstructured meshes and implicit time integration (Saxena et al., 7 Jul 2025, Liu et al., 2021).

The time-advance strategy is central to its operating regime. The split-implicit formulation is used to treat stiff MHD components implicitly, permitting long-time evolution across resistive and transport timescales rather than only Alfvénic timescales (Krebs et al., 2019). In the hybrid code, the MHD system remains implicit or semi-implicit, while energetic particles are advanced explicitly in subcycles between MHD updates (Liu et al., 2021). In SPARC sawtooth simulations, implicit time stepping is highlighted as enabling long multi-timescale runs from Alfvénic to resistive dynamics (Wang et al., 2 Apr 2026).

The mesh and domain structure depend on application. For axisymmetric VDEs, a single unified mesh spans a central plasma region, a finite-thickness resistive wall region, and an outer vacuum region, so the wall is modeled as a true volumetric conductor rather than a thin shell (Krebs et al., 2019). For SPARC baseline-like scenarios, two meshes are used: a whole-device mesh including plasma, first wall, conducting vacuum vessel, and coils, and a cheaper plasma-only mesh with the boundary close to the last closed flux surface (Wang et al., 2 Apr 2026). In the QA stellarator nonlinear runs reported in the bootstrap-current paper, the 3D grid uses 36 toroidal planes and about $C^1$2 3D elements (Saxena et al., 7 Jul 2025).

A plausible implication is that M3D-C$C^1$3 derives much of its versatility from the combination of global geometry, high-order continuity, and implicit advancement rather than from any single plasma model.

3. Governing equations and physical closures

The code is used with several related MHD models. In the VDE study, the plasma region employs single-fluid, resistive MHD with mass continuity, a momentum equation including viscosity, an induction equation via Ohm’s law with Spitzer resistivity, an electron temperature equation with anisotropic heat conduction and Ohmic heating, and particle diffusion (Krebs et al., 2019). In the SPARC study, the governing model is presented as a single-fluid extended MHD with two temperatures, $C^1$4 and $C^1$5, including ion continuity, momentum balance, induction, and separate electron and ion energy equations (Wang et al., 2 Apr 2026).

The SPARC formulation writes

$C^1$6

with Ohm’s law

$C^1$7

and total pressure $C^1$8 under quasineutrality $C^1$9 (Wang et al., 2 Apr 2026). The same study states explicitly that Hall, electron-inertia, and FLR terms are not used there, so the operative model is effectively single-fluid resistive or ideal MHD (Wang et al., 2 Apr 2026).

Geometry-dependent field representations are also specified. In the SPARC work the magnetic field and velocity are written in Clebsch-like form,

$^1$0

$^1$1

with primary unknowns $^1$2 (Wang et al., 2 Apr 2026). In the VDE benchmark, the wall region evolves according to

$^1$3

while the vacuum region is treated with $^1$4 (Krebs et al., 2019).

Boundary conditions are problem-specific and physically consequential. In the VDE benchmark, no-slip velocity is imposed at the plasma–wall interface, temperature and density are fixed there, and no boundary condition is imposed on the magnetic field at the resistive wall within the mesh, allowing halo currents to flow into and out of the wall (Krebs et al., 2019). In SPARC-like plasma-only runs, perturbations vanish at a boundary placed near the last closed flux surface (Wang et al., 2 Apr 2026).

4. Bootstrap-current modeling and quasisymmetric stellarators

A major 2025 extension adds self-consistent bootstrap-current models to M3D-C$^1$5. Two analytical closures are implemented: a generalized Sauter model and a revised Sauter-like model due to Redl et al. (Saxena et al., 7 Jul 2025). In both cases the quantity supplied by neoclassical theory is the flux-surface-averaged parallel current density $^1$6.

For the original Sauter model, the implementation uses

$^1$7

while the Redl model is implemented as

$^1$8

(Saxena et al., 7 Jul 2025). The same source notes that Sauter’s model works well at low collisionality but becomes inaccurate for $^1$9, which motivates inclusion of Redl’s revision (Saxena et al., 7 Jul 2025).

The bootstrap-current source enters the MHD system through a modified Ohm’s law,

$C^1$0

so the induction equation becomes

$C^1$1

The model current is taken as purely parallel and divergence-free, with

$C^1$2

(Saxena et al., 7 Jul 2025).

The same study generalizes the closure to quasisymmetric stellarators using the Landreman isomorphism, which maps neoclassical transport in a QS stellarator to an effective axisymmetric problem (Saxena et al., 7 Jul 2025). In this setting M3D-C$C^1$3 uses effective tokamak-like quantities such as $C^1$4, $C^1$5, and $C^1$6, derived from stellarator geometry, and applies a Landreman-modified Sauter-Redl formula to compute bootstrap current in QA configurations (Saxena et al., 7 Jul 2025). The implementation is benchmarked against NEO, XGCa, and SFINCS, with percent-level agreement reported for the benchmark cases, and nonlinear QA simulations show that enabling the bootstrap model largely maintains the toroidal current density whereas disabling it allows strong decay (Saxena et al., 7 Jul 2025).

This suggests that M3D-C$C^1$7 has evolved from an extended-MHD solver with prescribed transport coefficients into a platform that can incorporate neoclassical current closures in both axisymmetric and quasisymmetric geometries.

5. Hybrid energetic-particle extension: M3D-C1-K

M3D-C1-K is a hybrid PIC–MHD extension built on M3D-C$C^1$8 for simulations of energetic particles interacting with global plasma modes (Liu et al., 2021). The bulk plasma remains fluid, while energetic particles are advanced kinetically. The kinetic module uses a slow-manifold “classical Pauli particle” Boris pusher, a $C^1$9 formulation for low-noise moment calculation, and two coupling schemes, pressure coupling and current coupling (Liu et al., 2021).

The slow-manifold algorithm is designed to preserve long-term orbit fidelity. Instead of explicit full-orbit integration or conventional guiding-center RK4 alone, it advances a fictitious particle whose orbit lies on the slow manifold of the Lorentz dynamics by using the Boris update with an effective field $^1$0 (Liu et al., 2021). The hybrid paper reports that errors in toroidal canonical momentum and energy remain bounded under the slow-manifold Boris method, whereas RK4 guiding-center integration exhibits secular drift, especially for passing particles with large $^1$1 (Liu et al., 2021).

Energetic-particle moments are coupled back into the MHD solver either through $^1$2 in pressure coupling or through $^1$3 in current coupling (Liu et al., 2021). The code also supports FLR effects through 4-point gyro-averaging of fields and deposition (Liu et al., 2021). GPU acceleration is implemented for particle pushing and weight evolution with OpenACC. On the cited Summit benchmark, 4 million particles pushed for 50 time steps gave 134.22 s for RK4 guiding center on CPUs, 12.65 s for RK4 guiding center on GPUs, and 10.32 s for slow-manifold Boris on GPUs (Liu et al., 2021).

Validation is performed on linear fishbone, toroidal Alfvén eigenmode, and reversed-shear Alfvén eigenmode problems, with good agreement reported against M3D-K, NIMROD, and other eigenvalue, kinetic, and hybrid codes (Liu et al., 2021). The paper states that differences between pressure and current coupling are negligible in the benchmarked linear problems (Liu et al., 2021).

6. Benchmarks, applications, and scientific role

A defining use of M3D-C$^1$4 is benchmark-oriented cross-code comparison. In the 2019 VDE study it is one of three 3D nonlinear MHD codes, together with JOREK and NIMROD, used for axisymmetric simulations of vertical displacement events in a vertically unstable NSTX-like equilibrium (Krebs et al., 2019). That benchmark reports good agreement in linear growth rates between linear NIMROD and M3D-C$^1$5, and excellent agreement in full nonlinear VDE evolution regarding plasma location, plasma currents, and eddy and halo currents in the wall (Krebs et al., 2019). The same work identifies several physical regimes: for sufficiently small wall resistivity and large edge resistivity, the linear VDE growth rate scales as $^1$6, whereas at higher edge temperature or larger wall resistivity the growth rate becomes limited by the open-field-line plasma response rather than the wall alone (Krebs et al., 2019).

The SPARC study uses M3D-C$^1$7 as a high-fidelity 3D extended-MHD code to investigate low-$^1$8 instabilities and sawtooth crashes in relaxed baseline cases based on the SPARC Primary Reference Discharge (Wang et al., 2 Apr 2026). The linear simulations identify a dominant $^1$9 internal kink at the $n=0$0 surface, with growth rates that are strongly sensitive to the temperature profile, plasma $n=0$1, and $n=0$2 near unity (Wang et al., 2 Apr 2026). Nonlinear runs show two qualitatively distinct regimes: low-$n=0$3 current-driven sawtooth-like oscillations with moderate pressure flattening, and a baseline high-$n=0$4 case in which combined current and pressure drive yields a strong sawtooth with magnetic reconnection and a hollowed pressure profile (Wang et al., 2 Apr 2026). The paper interprets the baseline event as a mixture of Kadomtsev and Wesson mechanisms rather than a purely one-model crash (Wang et al., 2 Apr 2026).

Taken together, these studies place M3D-C$n=0$5 in a particular niche. It is repeatedly used where global geometry, current redistribution, wall coupling, or profile-driven mode structure matter enough that reduced or local models are inadequate. This suggests that its principal scientific role is as a predictive extended-MHD framework for problems in which equilibrium evolution, stability, and nonlinear topology change must be handled self-consistently.

7. Limits, assumptions, and nomenclature

The cited literature is explicit about modeling limits. In the VDE benchmark all calculations are axisymmetric, so inherently non-axisymmetric stages of VDEs are not simulated there, and the thermal quench is imposed artificially through transport coefficients (Krebs et al., 2019). In the bootstrap-current work, Sauter and Redl are local closures with finite validity ranges, and the Landreman extension is expected to fail once strong non-QS structure or magnetic stochasticity develops; in such cases the paper advocates post-processing M3D-C$n=0$6 equilibria with SFINCS or NEO 3D to quantify error (Saxena et al., 7 Jul 2025). In the SPARC sawtooth study, external heating, energetic particles, Hall physics, electron inertia, and FLR effects are omitted, so the reported crash thresholds and recovery times are intentionally an MHD-only baseline rather than a full burning-plasma prediction (Wang et al., 2 Apr 2026). In the hybrid energetic-particle extension, the published benchmarks focus on linear physics, with nonlinear EP–MHD interaction identified as a target for future work (Liu et al., 2021).

A recurring source of ambiguity is the string “M3D” itself. Several unrelated arXiv works use “M3D” for topics such as 3D medical image analysis, dataset condensation, methylation-profile testing, and monolithic 3D DRAM, and those papers explicitly do not define “M3D-C1” as a concept within their own nomenclature (Bai et al., 2024, Zhang et al., 2023, Mayo et al., 2014, Huang et al., 2020). In the plasma-physics literature, by contrast, M3D-C$n=0$7 consistently denotes the finite-element MHD code and its extensions (Krebs et al., 2019, Liu et al., 2021, Saxena et al., 7 Jul 2025, Wang et al., 2 Apr 2026).

In this sense, M3D-C$n=0$8 is not a generic “M3D” variant but a specific computational framework in fusion-plasma simulation: a high-order, implicit, geometry-resolving extended-MHD code that has been extended toward neoclassical bootstrap-current closure and hybrid energetic-particle kinetics while retaining its role as a benchmark-grade nonlinear MHD solver.