NEOPIC: A Neural Operator Framework for Particle-based Kinetic Plasma Simulations

Sriramkrishnan Muralikrishnan; Chelsea Maria John; Jeffrey Kelling; Michael Bussmann; Robert Speck; Paul Gibbon

math.NAOct 29, 2026

NEOPIC: A Neural Operator Framework for Particle-based Kinetic Plasma Simulations

Authors: Sriramkrishnan Muralikrishnan, Chelsea Maria John, Jeffrey Kelling, Michael Bussmann, Robert Speck, Paul Gibbon

Organizations: Jülich Supercomputing Centre, Forschungszentrum Jülich GmbH, Jülich, Germany · Helmholtz-Zentrum Dresden-Rossendorf, Dresden, Germany · Chemnitz University of Technology, Chemnitz, Germany · Center for Advanced Systems Understanding (CASUS), Görlitz, Germany

Abstract

We introduce NEOPIC, a framework for learning the field operator in particle-based methods applied to kinetic plasma simulations. Its key novelty is a neural operator that acts directly on the empirical measure of the particles, rather than on a grid, and is evaluated at the particle positions. No mesh or tree is needed and, unlike grid-based neural operators, no deposition onto a mesh takes place, so the particle noise that dominates the inputs of particle simulations is not aliased. Since it replaces only the field solve, it serves as a plug-and-play field solver for particle-based methods: it can be trained on data from particle-in-cell, particle-in-Fourier or tree codes, and combined with any time integrator, including implicit schemes. We prove that the discrete evaluation coincides with the continuous operator at any number of particles (representation equivalence) and derive an error analysis that separates sampling, learning and time-stepping errors: up to the learning error, NEOPIC has the error structure of standard Monte Carlo particle methods. We show that the particle noise in the training data acts as an implicit regularizer, which gives guidance on the number of training particles in relation to the particle numbers used at inference. A complexity analysis shows that NEOPIC is an $O(N_p)$ method whose constants are set by the architecture alone, independently of any mesh or tree. Numerical experiments on electrostatic benchmarks in two and three dimensions, including Landau damping, two-stream and bump-on-tail instabilities, a cyclotron mini-app and a Penning trap, confirm the theoretical findings. A single model trained on several mini-apps generated by different reference methods captures all of them and remains stable up to twice the training time window, which points towards foundation models for fields in particle-based methods.

Figures & tables

Method	Kernel $K_{h}(\mathbf{x},\mathbf{y};X)$	Structure, cost
PIC (mesh)	$\sum_{i,j}\psi^{*}_{i}(\mathbf{x})\,(\mathsf{L}_{h}^{-1})_{ij}\,\psi_{j}(\mathbf{y})$	factorized ( 10 ); $\mathcal{O}(N_{p})+C_{h}$
PIF	$\frac{1}{\|\Omega\|}\sum_{\mathbf{k}\in\mathcal{K}\setminus\{0\}}\|\hat{S}_{\mathbf{k}}\|^{2}\,\frac{-\mathrm{i}\mathbf{k}}{\|\mathbf{k}\|^{2}}\,e^{\mathrm{i}\mathbf{k}\cdot(\mathbf{x}-\mathbf{y})}$	factorized ( 10 ), translation invariant; $\mathcal{O}\big((\|\log\varepsilon\|+1)^{d}N_{p}+\|\mathcal{K}\|\log\|\mathcal{K}\|\big)$
Direct sum	$-\nabla G_{\delta}(\mathbf{x}-\mathbf{y})$	translation invariant; $\mathcal{O}(N_{p}^{2})$
Barnes–Hut	$-\nabla G_{\delta}(\mathbf{x}-\mathbf{y})$ for $\mathbf{y}\in\mathcal{N}(\mathbf{x})$ ; $T_{P}[-\nabla G](\mathbf{x}-\mathbf{x}_{c};\mathbf{y}-\mathbf{x}_{c})$ for $\mathbf{y}\in c\in\mathcal{F}(\mathbf{x})$	depends on $X$ and on the target; $\mathcal{O}(N_{p}\log N_{p})$
FMM	near field as above; far field separable between well-separated cells (multipole-to-local translations)	hierarchical low rank, depends on $X$ ; $\mathcal{O}(N_{p})$

Table 1: Discrete kernels

K_{h}

in ( 9 ) (up to the factor

Q_{\mathrm{tot}}/\varepsilon_{0}

\psi_{j},\psi^{*}_{i}

: deposition and interpolation functions;

\mathsf{L}_{h}

: discrete field operator;

\hat{S}_{\mathbf{k}}

: Fourier transform of the shape function;

\mathcal{K}

: retained Fourier modes;

G_{\delta}

: softened Green’s function;

T_{P}

: multipole expansion of order

P

about the cell centre

\mathbf{x}_{c}

;

\mathcal{F}(\mathbf{x})

\mathcal{N}(\mathbf{x})

: far-field cells and near-field particles of target

\mathbf{x}

for opening angle

\vartheta

;

\varepsilon

: NUFFT tolerance. Structure: whether the kernel factorizes through a space that does not depend on the particles, as in ( 10 ), whether it is translation invariant, and whether it depends on the particle configuration

X

. Costs per field evaluation, with

C_{h}

the cost of the mesh solve.

Figure 1: (a) The particle–field cycle of a particle method for the Vlasov–Poisson system. The particles are advanced by the time integrator (green); the field solver (grey box) evaluates the field at the particles with a discrete kernel

K_{h}

( 9 ), whose form depends on the method (Table 1 ): a mesh or basis (

\mathscr{R}^{*}\mathscr{L}_{h}^{-1}\mathscr{R}

), a tree, or direct summation. (b) NEOPIC keeps the particles and the time integrator and replaces the field solver by a neural operator

\mathbf{E}_{\theta}

, which maps the charge measure of the particles directly to the field at the particles. It is trained offline on data from particle methods; in this work PIC, PIF and Barnes–Hut, and the neural operator is the NUFNO of Section 4 . With an implicit time integrator ( E ) the cycle is the same, but the field is evaluated several times per step within the nonlinear solve.

Figure 2: The nonuniform Fourier neural operator (NUFNO) of NEOPIC. (a) The training data come from particle methods with different kernels

K_{h}

(Table 1 ): PIC uses a mesh, PIF a truncated Fourier basis, and Barnes–Hut multipole expansions on the cells of a tree. The examples illustrate their complementary strengths: PIF for smooth densities on periodic domains, Barnes–Hut for clustered particles and open boundaries (particles sampled from the density of the cyclotron mini-app at

t=12

), and PIC for complex geometries (a tokamak cross-section, for illustration only; the PIC teachers of this work use uniform periodic meshes, Section 8 ). All of them provide the particles

\{\mathbf{x}_{p},w_{p}\}

and the field at the particles. (b) The NUFNO maps the charge measure

\mu_{X}

, normalized from its bounding box to the torus by

\mathcal{N}_{X}

, Eq. ( 13 ), through a pointwise lifting

\mathcal{P}

L

Fourier layers and a pointwise projection

\mathcal{Q}

to the field at the particles, followed by the scaling and mean subtraction of ( 19 ). (c) A Fourier layer: a type-1 nonuniform transform (NUDFT or NUFFT) from the particle features

\mathbf{U}^{\ell}_{p}

to the coefficients on the modes

\mathcal{K}_{m}

( 22 ), multiplication by the spectral weights, a type-2 nonuniform transform back to the particles ( 23 ), and the pointwise update ( 24 ). No mesh is used at any stage. The point clouds and the spectra are illustrations; in the PIF panel the colour shows the smooth density from which the particles are sampled.

Operation	NUDFT (DSE)	NUFFT
Lifting $\mathcal{P}$	$N_{p}C_{\mathcal{P}}$	$N_{p}C_{\mathcal{P}}$
Forward transform (per layer)	$N_{p}Md_{v}$	$(N_{p}w^{d}+\sigma^{d}M\log M)d_{v}$
Spectral mixing (per layer)	$Md_{v}^{2}$	$Md_{v}^{2}$
Inverse transform (per layer)	$N_{p}Md_{v}$	$(N_{p}w^{d}+\sigma^{d}M\log M)d_{v}$
Skip $\mathcal{W}_{\ell}$ , activation (per layer)	$N_{p}(C_{\mathcal{W}}+d_{v})$	$N_{p}(C_{\mathcal{W}}+d_{v})$
Projection $\mathcal{Q}$	$N_{p}C_{\mathcal{Q}}$	$N_{p}C_{\mathcal{Q}}$

Table 2: Leading-order cost of one evaluation of the NEOPIC field operator.

L

: number of Fourier layers;

d_{v}

: lifting (channel) width;

M=(2m)^{d}

: number of retained modes;

w

: NUFFT kernel width;

\sigma

: upsampling factor.

Figure 3: 2D-2V weak Landau damping, teacher-forced field error $\langle e_{E}\rangle_{t}$ of the five models (training particle number $N_{\mathrm{tr}}$ ) evaluated with $N_{p}$ particles. (a) All pairs, in units of $10^{-3}$ ; diagonal entries ( $N_{p}=N_{\mathrm{tr}}$ ) outlined. (b) Error against $N_{p}$ , with the diagonal (stars) and a power-law fit to it.

Figure 4: 2D-2V weak Landau damping, in-loop errors ( 66 ) averaged over

t\leq 20

, against the inference particle number

N_{p}

for the five trained models. (a)–(c) Reference: the CIC teacher with the same particles; (d)–(f) reference: high-fidelity PIF (fourth-order shape). Stars mark

N_{p}=N_{\mathrm{tr}}

; dashed lines are power-law fits to the diagonal (a) and to the model with

N_{\mathrm{tr}}=8\times 10^{5}

(d)–(f).

Figure 5: 2D-2V weak Landau damping, time histories of the in-loop errors against the CIC teacher at

N_{p}=1.6\times 10^{6}

for the five trained models.

Figure 6: Generalization of the weak-Landau model (

k=0.5

N_{\mathrm{tr}}=8\times 10^{5}

) at

N_{p}=1.6\times 10^{6}

, errors against the CIC teacher. (a) Other wave numbers; the dotted line marks the training wave number. (b) Other mini-apps: time-averaged errors (bars) and time medians (black marks).

Figure 7: 2D-2V multi-app model at

N_{p}=10^{6}

, errors ( 66 ) against each mini-app’s teacher (colour: teacher). (a)–(c) Time histories with time normalized by the training horizon

T_{\mathrm{train}}

; (d)–(f) time averages within (plain) and beyond (hatched) the training horizon.

Figure 8: 2D-2V multi-app model, field energy per component for the four periodic mini-apps: teacher (wide lines, colour as in Fig. 7 ) and NEOPIC (red).

W_{y}

(the growing component of the instabilities) is solid and

W_{x}

dashed. Dash-dotted: linear theory; dotted: end of training.

Figure 9: 2D-2V multi-app model, relative change of the total energy (solid) and of the momentum (dashed; for the cyclotron the canonical angular momentum

L_{\mathrm{can}}

( 68 )) for the teachers (wide lines) and NEOPIC (red).

Figure 10: Cyclotron mini-app (2D magnetized electron cloud, open boundaries): particle density $n/N_{p}$ per unit area for Barnes–Hut (top) and NEOPIC (bottom) at $N_{p}=10^{6}$ , binned on a fixed grid about the centre of the cloud. Training ends at $t=6$ .

Figure 11: Cyclotron mini-app: local difference between the field of each method and a free-space PIC field evaluated at the same particle positions, averaged per bin and normalized by the maximum PIC field.

Figure 12: Cyclotron mini-app. (a) Relative

\ell^{2}

difference between the field of each method and a free-space PIC evaluation at the same positions. (b) Winding of the two spiral arms between

r=0.09

and

r=0.29

, from the phase

\varphi(r)

of the

m=2

azimuthal mode of the density. (c) Arm shape

|\varphi(r)-\varphi(0.09)|

t=11.98

Figure 13: 3D-3V multi-app models at

N_{p}=163\,840

, field energy per component: reference (wide lines; grey: Barnes–Hut, orange: PIF), model trained on four mini-apps (red) and model trained on five mini-apps including the Penning trap (purple).

W_{z}

(the growing component of the instabilities) is solid,

W_{x}

dashed and

W_{y}

dotted. For the bump-on-tail instability the linear phase is below the noise level and no theory line is shown.

Figure 14: Field fluctuation level of the 3D-3V models: runs from the same initial particles (Barnes–Hut) compared with both Barnes–Hut and PIF. (a) Weak Landau damping, mean of the three components; (b) bump-on-tail instability, mean of the two non-growing components; (c) bump-on-tail instability, growing component.

	reference		NEOPIC, 4 apps		NEOPIC, 5 apps
mini-app (teacher)	$t\leq T_{\mathrm{tr}}$	$t>T_{\mathrm{tr}}$	$t\leq T_{\mathrm{tr}}$	$t>T_{\mathrm{tr}}$	$t\leq T_{\mathrm{tr}}$	$t>T_{\mathrm{tr}}$
weak Landau (BH)	$3.5{\cdot}10^{-4}$	$3.6{\cdot}10^{-4}$	$3.8{\cdot}10^{-4}$	$4.1{\cdot}10^{-4}$	$3.1{\cdot}10^{-4}$	$3.2{\cdot}10^{-4}$
strong Landau (BH)	$8.3{\cdot}10^{-2}$	$7.9{\cdot}10^{-2}$	$7.2{\cdot}10^{-2}$	$6.8{\cdot}10^{-2}$	$6.9{\cdot}10^{-2}$	$6.6{\cdot}10^{-2}$
two-stream (PIF)	$1.9{\cdot}10^{-3}$	$6.7{\cdot}10^{-5}$	$2.1{\cdot}10^{-2}$	$3.5{\cdot}10^{-3}$	$8.9{\cdot}10^{-3}$	$5.1{\cdot}10^{-3}$
bump-on-tail (PIF)	$3.9{\cdot}10^{-4}$	$3.3{\cdot}10^{-4}$	$3.1{\cdot}10^{-3}$	$3.6{\cdot}10^{-3}$	$9.5{\cdot}10^{-4}$	$8.4{\cdot}10^{-4}$
Penning trap (PIC)	$3.0{\cdot}10^{-3}$	$1.6{\cdot}10^{-3}$	–	–	$6.6{\cdot}10^{-2}$	$5.0{\cdot}10^{-2}$

Table 3: 3D-3V multi-app models: maximum relative change of the total energy,

\max_{t}|\mathcal{E}(t)-\mathcal{E}(0)|/|\mathcal{E}(0)|

, within (

t\leq T_{\mathrm{tr}}

) and beyond (

t>T_{\mathrm{tr}}

) the training horizon, for the reference and the models trained on four and five mini-apps. The smaller of the two model errors is shown in bold.

Figure 15: Penning trap (3D, five-app model,

N_{p}=163\,840

). (a) Kinetic, self-field and external potential energy; (b) relative change of the total energy.

Figure 16: Time per field evaluation (one time step) against the number of particles for (a) 2D and (b) 3D Landau damping, on one GPU. Colours identify the method, line styles the implementation: Python/PyTorch (solid, filled markers) and the C++ codes IPPL and Ryoanji (dashed, open markers). Dotted lines: fits

T=\alpha N_{p}+\beta

to the NEOPIC timings.

Appendix figures & tables3 assets

Supplementary material from the paper’s appendix.

Appendix

Figure 17: Cyclotron mini-app, 2D multi-app model,

N_{p}=10^{6}

T=3

. (a) Mean distance of the guiding centres ( 75 ) at

t=3

between runs with time step

\Delta t

and the finest run (

\Delta t_{\min}=3.125\times 10^{-4}

), for PIF and NEOPIC, and between NEOPIC and PIF at the same and at the finest time step; the dash-dotted line has slope

2

. (b) Maximum relative change of the total energy (solid) and of the canonical angular momentum

L_{\mathrm{can}}

( 68 ) (dashed) for PIF and NEOPIC. The grey line marks the time step of the training data.

Figure 18: 1D-1V weak Landau damping with leapfrog PIF, leapfrog NEOPIC and implicit-midpoint NEOPIC (

N_{p}=10^{5}

T=60

; the training data end at

t=30

, dotted). (a)–(c) Field energy

W_{x}

for

\Delta t=0.05

1

and

2

; in (a) the dash-dotted line is the linear damping rate. (d) Maximum relative change of the total energy over the run against

\Delta t

; the leapfrog curves of PIF and NEOPIC lie on top of each other. The shaded region marks

\omega\Delta t>2

, where the leapfrog scheme is unstable for the Langmuir mode.

Figure 19: Cost of NEOPIC (NUFFT, 3D) for sub-networks of the trained model with

L

layers,

d_{v}

channels and

2m

modes per direction. (a) Cost per particle

\alpha

, from timings at

N_{p}=10^{6}

and

4\times 10^{6}

, against

Ld_{v}

; the line is the least-squares fit. (b) Time per step at

N_{p}=4\times 10^{6}

against the number of parameters. Colours:

d_{v}

; markers:

L

; filled:

m=16

, open:

m=8

Physics-Based Simulation

Physics-based simulation models physical systems with equations and numerical solvers, increasingly augmented by learned models for speed and control.

See all 231 papers

A Neural Score-Based Particle Method for the Vlasov-Maxwell-Landau System
Mar 26, 2026Vasily Ilin, Jingwei HuNeural PDE Solvers
Department of Mathematics, University of Washington, Seattle, WA, USA. · Department of Applied Mathematics, University of Washington, Seattle, WA, USA.
Learn the Dynamics, Enforce the Physics: Generative Latent Operators for High-Dimensional Beam Evolution
Oct 10, 2026Haohan Zou, Alexander Scheinker, Yuanyuan ShiLatent Dynamics Modeling Physics-Informed Generative Modeling
University of California San Diego · Los Alamos National Laboratory
DeepSPoC: A Deep Learning Based Sequential Propagation of Chaos
Aug 29, 2024Kai Du, Yongle Xie, Tao Zhou +1Stochastic Differential Equations Neural PDE Solvers
Shanghai Center for Mathematical Sciences, Fudan University, Shanghai, 200438, China. · Shanghai Center for Mathematical Sciences, Fudan University, Shanghai, 200433, China. · LSEC, Institute of Computational Mathematics and Scientific/Engineering Computing, AMSS, Chinese Academy of Sciences, Beijing, China.

Mar 26, 2026J/K move · Enter open · S save

A Neural Score-Based Particle Method for the Vlasov-Maxwell-Landau System

Vasily Ilin, Jingwei Hu

Department of Mathematics, University of Washington, Seattle, WA, USA. · Department of Applied Mathematics, University of Washington, Seattle, WA, USA.

Neural PDE Solvers

Open paper AskPDF arXiv

Plasma modeling is central to the design of nuclear fusion reactors, yet simulating collisional plasma kinetics from first principles remains a formidable computational challenge: the Vlasov-Maxwell-Landau (VML) system describes six-dimensional phase-space transport under self-consistent electromagnetic fields together with the nonlinear, nonlocal Landau collision operator. A recent deterministic particle method for the full VML system estimates the velocity score function via the blob method, a kernel-based approximation with

O(n^2)

cost. In this work, we replace the blob score estimator with score-based transport modeling (SBTM), in which a neural network is trained on-the-fly via implicit score matching at

O(n)

cost. We prove that the approximated collision operator preserves momentum and kinetic energy, and dissipates an estimated entropy. We also characterize the unique global steady state of the VML system and its electrostatic reduction, providing the ground truth for numerical validation. On three canonical benchmarks -- Landau damping, two-stream instability, and Weibel instability -- SBTM is more accurate than the blob method, achieves correct long-time relaxation to Maxwellian equilibrium where the blob method fails, and delivers

50\%

faster runtime with

4\times

lower peak memory.

Abstract

Physics-Based Simulation

Explore similar work

A Neural Score-Based Particle Method for the Vlasov-Maxwell-Landau System

Learn the Dynamics, Enforce the Physics: Generative Latent Operators for High-Dimensional Beam Evolution

DeepSPoC: A Deep Learning Based Sequential Propagation of Chaos