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

1. Introduction

• Standard PIC’s grid can alias unresolved modes, while its noise and Debye-length resolution requirements increase simulation cost; more advanced schemes can be harder to adopt.PIF can be expensive when many Fourier modes are needed and cannot handle some boundary conditions; tree methods avoid a grid but have larger cost constants and harder parallelization.Source: §1 · pp. 4–5
• NEOPIC learns the field operator shared by particle methods from data that can combine methods in their respective strengths, then returns fields at particle positions for plug-and-play use.The framework focuses on collisionless electrostatic simulations and can be combined with explicit or implicit time integrators.Source: §1 · p. 5
• NEOPIC’s NUFNO acts on empirical particle measures and evaluates fields at those particles, with representation equivalence preventing particle noise from aliasing onto retained Fourier modes.Without this property, nonlinear activations could worsen aliasing and contribute to instabilities like finite-grid instability and grid heating.Source: §1 · p. 6
• The analysis separates Monte Carlo sampling, time-stepping, and model errors; training noise acts as an implicit regularizer, and the method has O(Np) complexity with an architecture-determined constant.The training-noise analysis gives guidance on training particle counts relative to expected inference counts.Source: §1 · p. 6
• A single model trained on several mini-apps from different methods remains stable up to twice the training horizon, supporting exploration of foundation models for particle-method fields.The paper leaves fine-tuning such a model with application-specific samples to future work.Source: §1 · p. 6

1.1. Related work

• Earlier learned plasma field solvers either operated on deposited grid sources or learned a different part of the simulation; NEOPIC instead learns the field directly from particles.Prior grid-based surrogates require deposition and may require a dedicated compression stage, while graph approaches can require many layers or large neighborhoods for long-range interactions.Source: §1.1 · p. 7
• NEOPIC shares direct nonuniform spectral evaluation with Lingsch et al. [49] but treats the particle measure as the state, so its discrete evaluation is the continuous operator applied to that measure.Unlike the single linear convolutional layer of Peng et al. [67], NEOPIC composes Fourier layers and nonlinearities to learn nonlinear and non-translation-invariant field operators.Source: §1.1 · p. 8

2.1. The Vlasov–Poisson system

• The Vlasov–Poisson setting evolves the electron distribution under self-consistent and possibly external fields, with periodic or open boundary conditions and conservation properties depending on the setting.The space-charge field is a linear, nonlocal functional of charge density; periodic domains include a neutralizing ion background, while open domains use a decaying potential without that background.Source: §2.1 · pp. 8–9

2.2. Particle discretization

• Particle discretization samples the distribution with weighted macro-particles that follow characteristics, conserving total charge exactly while approximating the field of a smoothed charge distribution.The point-measure field is singular at particle locations, so particle methods regularize it through shape functions, Fourier truncation, or softened kernels.Source: §2.2 · pp. 9–10
• Explicit time integrators evaluate fields once per step, whereas implicit schemes require repeated field evaluations within a nonlinear solve.Newton–Krylov solvers may also require field evaluations for Jacobian–vector products.Source: §2.2 · p. 10

2.3. An abstract particle–field method

• Particle–field methods share an interface: particles generate fields and receive fields, and a field solver maps the particle measure to fields at particle positions.The methods differ in their discrete kernels and evaluation procedures, though each approximates the same continuum operator.Source: §2.3 · pp. 10, 12, 13
• PIC represents fields on a mesh, PIF uses a truncated Fourier basis and nonuniform transforms, and Barnes–Hut and FMM use hierarchical tree interactions without a mesh.Barnes–Hut combines near-field particle interactions with far-field multipole approximations; FMM uses multipole-to-local translations between separated cells.Source: §2.3 · pp. 11–12

3.1. Learning the particle–field–particle map

• NEOPIC replaces the field solver with a neural operator that maps the particle charge measure to fields at particle positions, while leaving particle evolution and time integration in place.The operator can be used with existing particle codes without changing their remaining components.Source: §3.1 · pp. 13–14
• Learning only the instantaneous field map leaves the particle update to the time integrator, which preserves charge and permits structure-preserving or implicit schemes.The field operator does not depend on time or time step, although its particle-configuration inputs may move away from the training distribution.Source: §3.1 · p. 14
• The learned operator acts on the empirical particle measure itself, rather than treating particle positions as samples of an underlying input function.It must apply to measures with varying particle counts and approximate the continuum field operator.Source: §3.1 · p. 15

3.2. Training on several particle methods

• A single operator can train on particle configurations and fields from multiple teachers and mini-apps, whose data may be generated with the method best suited to each scenario.The learned operator can be evaluated with particle counts different from those used in training, and inherits the teacher’s error as well as its own learning error.Source: §3.2 · p. 15

3.3. Requirements on the learned operator

• The requirements exclude grid deposition and fixed sensor locations; NUFNO instead uses global Fourier layers with nonuniform transforms between particle features and modes.For a fixed number of modes, NUFFT-based transforms cost O(Np), and the layers are smooth in particle positions; the framework can also use any operator satisfying the stated requirements.Source: §3.3 · p. 16
• Figure 2’s NUFNO processes a normalized particle measure through pointwise lifting, Fourier layers, and pointwise projection, without using a mesh.Each Fourier layer transforms particle features to Fourier modes, applies spectral weights, transforms back to particles, and updates pointwise.Source: Figure 2 · p. 17

4.2. Particles as a measure

• The operator is defined independently of particle count and applies to both unit-charge empirical measures and continuous charge densities.The network evaluates its field at any location, and normalized particle weights make it applicable across particle counts.Source: §4.2 · p. 19, §4.3 · p. 20

4.3. The operator

• NEOPIC applies a lifting, Fourier layers, and a projection to a probability measure, with pointwise maps and spectral weights defining its network.In this work, the lifting is a linear layer on normalized positions, the skip maps are kernel-size-one convolutions, and the projection is a two-layer MLP.Source: §4.3 · p. 19

4.4. Discrete evaluation

• The discrete evaluation exactly computes the operator on the input particle measure at every particle count and is invariant to particle permutations, splitting, or merging at coincident positions.These invariances follow because outputs depend on particle data only through the measure.Source: §4.4 · p. 21

4.5. Representation equivalence

• With NUDFTs and distinct particle positions, NEOPIC is a representation-equivalent neural operator (ReNO) with zero aliasing error [50].The theorem establishes this for each layer and the composed network.Source: §4.5 · p. 22
• Across particle samples from one density, output differences are unbiased sampling errors of order N_p^-1/2, unlike grid-FNO aliasing, which does not decrease with particle count.Particle sampling fluctuations behave as for the continuous operator applied to the empirical measure and are bounded independently of integrand smoothness.Source: §4.5 · pp. 23–24
• With NUFFTs, NEOPIC is an ε-ReNO whose error is proportional to the transform tolerance uniformly over particle discretizations, independent of particle count, positions, and weights.NUFFT spreading onto an oversampled grid introduces kernel aliasing, which the kernel width and upsampling factor control.Source: §4.5 · pp. 24–25
• The representation-error bound does not grow with particle count, and NUFFT aliasing is controlled by tolerance rather than fixed by a grid.The particle-noise contribution reaches the field through the continuous operator applied to the particle measure; its coefficient fluctuations are controlled by Lemma 3 and, in training, Section 6.Source: §4.5 · p. 26
• Representation equivalence establishes consistency of the operator and its discrete evaluation, not trained-operator approximation accuracy, teacher-bias control, or scheme stability.The passage says these other issues are treated in Section 5.Source: §4.5 · p. 26
• NEOPIC removes aliasing as one cause of numerical heating, but does not rule out every heating mechanism or guarantee energy conservation or stability.Scheme stability is governed by the learned operator’s Lipschitz constant Lθ; energy errors in the mini-apps do not grow in extrapolation windows.Source: §4.5 · p. 26

5.1. Setting and notation

• The analysis concerns periodic Vlasov–Poisson systems in dimensions one through three and uses the learned measure operator, whose particle evaluation matches its empirical-measure evaluation.It states the error analysis for kick–drift–kick leapfrog; other integrators are discussed separately.Source: §5.1 · pp. 27–28

5.2. Properties of the learned field operator

• The learned operator’s sampling error depends on architecture but not particle number or grid resolution, unlike PIC’s sampling error.The analysis assumes iid equal-weight initial samples, regular exact solutions, a Lipschitz architecture, and bounded model error.Source: §5.2 · p. 29
• NEOPIC can represent PIF with the same Fourier modes and approximate continuous field operators uniformly over probability measures, but these approximation results do not guarantee successful training or control Lθ.For discontinuous Barnes–Hut targets, small expected error can still exist under the data distribution.Source: §5.2 · pp. 29–30
• The trained model’s error against the exact field includes both training error and teacher discretization bias; for models trained on multiple teachers, the bias is the largest among the teachers on the data seen.For PIC with shape functions of order α, the teacher discretization error is O(h^α).Source: §5.2 · pp. 30–31

5.3. Convergence of trajectories and fields

• Trajectory and field errors are bounded by contributions from sampling, time stepping, and model error, using independent exact-field reference particles without requiring numerical particles to be independent.The reference particles remain iid samples of the evolving distribution, enabling sampling bounds at each time.Source: §5 · p. 27, §5.3 · pp. 31–32

5.4. Convergence of the distribution function

• For bounded Lipschitz observables, particle errors combine trajectory error with Monte Carlo sampling, while the expected Wasserstein-1 error adds a dimension-dependent sampling rate.The rate is Np^-1/(2d) for d ≥2 and Np^-1/2 log(1 + Np) for d = 1; distribution-function convergence is no slower than field convergence, with a caveat that particles are not independent samples of the particle-dependent transported solution.Source: §5.4 · p. 33
• Model error is amplified by self-consistent dynamics, and the error floor persists as particle number grows and time step shrinks; stability depends on the learned operator’s Lipschitz constant.The growth rate depends on Lθ, so controlling it can improve stability; the model-error floor includes teacher bias.Source: §5.4 · p. 34
• The bounds retain Theorem 3 and Corollary 1 rates when the reference norm used for normalization is bounded away from zero.The experiments confirm the bounds’ error structure and rates.Source: §5.4 · p. 35
• Measured against the teacher, teacher-forced error in the reference norm is a lower estimate of the L∞ learning error at the given particle number.Source: §5.4 · p. 35
• The explicit constants are pessimistic: worst-case Gronwall growth misses phase-mixing saturation in Landau damping, while filamentation makes Kf grow over time.For distribution-function convergence, Theorem 4 statements (i)–(ii) are more natural measures.Source: §5.4 · p. 35

• Because labels use the same noisy particles as inputs, sampling noise regularizes the mismatch between the network and teacher responses rather than driving the network response toward zero.If the teacher is realizable, the noise does not bias the minimizer.Source: §6 · p. 37
• For density perturbations, small Q forces the network and teacher to respond alike, while finite network dimension can leave irreducible teacher fluctuation variance outside its response subspace.When Dθ is below the teacher-response rank, Q0 is positive.Source: §6 · pp. 37–38
• Training with more particles lowers noise-free-input error but increases sensitivity; expected loss is minimized by matching training and inference particle counts, or using their harmonic mean across inference counts.For inference over all counts at least Nmin, the worst-case loss is minimized by training with Nmin particles.Source: §6 · p. 38
• Training noise regularizes residual sensitivity in the teacher-relative stability exponent, but Q controls an average over fluctuation directions rather than a supremum.The exponent depends on teacher stability and residual Lipschitzness, not the network's overall Lipschitz constant.Source: §6 · p. 39
• Choose the smallest Ntr that keeps relevant modes above noise over the training horizon, and avoid exceeding the smallest or harmonic-mean particle count planned for inference.The expansion assumes SNRk(t) = √Ntr |ρ̂k(t)| ≫ 1 in modes carrying the field.Source: §6 · p. 40

7.3. Total cost

• T(Np) = O(Np) for both NUFFT-based and direct evaluation, with coefficients independent of Np at fixed architecture parameters and tolerance.The mode count is a property of the learned operator, not the particle count; a fixed per-step cost dominates only for small Np.Source: §7.3 · pp. 41–42
• For 3D with ε = 10^-6, NUFFT uses w^3 = 343 versus M = 32768 modes, reducing arithmetic per particle by roughly two orders of magnitude.Direct evaluation can still be faster for small Np, small M, or in 2D because it avoids the fixed fine-grid FFT cost.Source: §7.3 · p. 42
• The transform-cost ratio of NEOPIC to PIF reflects the learned operator’s additional depth, width, and nonlinearity, and can be reduced by choosing smaller L or dv.Measured ratios are smaller than the formula’s estimate because channel batches share sorting, kernel evaluation, and coordinate-memory traffic.Source: §7.3 · p. 42

7.4. Memory requirements

• The dimension-wise factorization cuts 3D transform memory by a factor of 16 when dv = 2m = 32, without changing the O(NpMdv) operation count.In 2D, dv(2m) = M, so factorization brings no memory reduction and the explicit matrix is used.Source: §7.4 · p. 43
• Peak inference memory is O(Np) for both evaluation forms, with NUFFT memory per particle independent of the number of modes.For 3D with m = 16 and Np = 10^6, dimension-wise factorization reduces direct-evaluation storage from about 260 GB to about 17 GB, still over an order of magnitude above NUFFT.Source: §7.4 · pp. 43–44

8.1. 2D-2V weak Landau damping

• On the Np = Ntr diagonal, 2D weak-Landau teacher-forced field error fell from 2.3 × 10−3 to 7.3 × 10−4; inference below training particle counts was substantially worse.The error matrix was asymmetric: inference above the training particle count was benign, and error decreased with Np until reaching a model-dependent floor. Training-data sampling noise regularized models against noisier inputs, at the cost of a higher error floor on clean inputs.Source: §8.1 · pp. 46–47, Figure 3 · p. 47
• Against the CIC teacher, 2D weak-Landau in-loop diagonal field error fell from 6.7 × 10−3 to 1.5 × 10−3, while high-fidelity-reference errors were indistinguishable across models.The in-loop error was 3.0 times the teacher-forced error at 5 × 10^4 particles and 2.0 times at 8 × 10^5; against the high-fidelity reference, all three errors decreased approximately as Np^-0.51, with the teacher’s discretization error dominating at the largest Np tested.Source: §8.1 · pp. 47–48, Figure 4 · p. 48
• A weak-Landau model generalized to nearby wave numbers and some other mini-apps, but failed on strong Landau damping outside its training input regime.At k = 0.25 and k = 0.125, field errors rose to 4.5 × 10−3 and 8.7 × 10−3, respectively. It remained accurate during bump-on-tail’s linear phase and until two-stream saturation, after which error jumped; strong-Landau position error was 4.6 × 10−2.Source: §8.1 · p. 49

8.2. 2D-2V multi-app model

• In the 2D multi-app model, field errors grew beyond training, most sharply for saturated two-stream, where field and position errors rose by factors of 12 and 39.Within the training horizon, mean field errors were 4 × 10−3 for both Landau cases, 1.5 × 10−2 for two-stream, 8 × 10−3 for bump-on-tail, and 0.11 for cyclotron. The two-stream divergence reflects trajectory sensitivity after particles become trapped in nonlinear potential wells.Source: §8.2 · p. 50, Figure 7 · p. 51
• In 2D multi-app tests, NEOPIC closely matched field-energy behavior across periodic cases, including instability growth rates and saturation peaks.For weak Landau damping, fitted damping rates were −0.290 (PIF) and −0.293 (NEOPIC), versus −0.307 from linear theory; two-stream growth rates were 0.614 versus 0.613, and bump-on-tail rates were 0.297 versus 0.301. The first saturation peaks occurred at t = 28.70 versus 28.75 and t = 36.35 versus 36.30, with amplitudes differing by 4% and 1%.Source: §8.2 · p. 51, Figure 8 · p. 52
• In periodic 2D mini-apps, NEOPIC’s energy error was comparable to its teacher’s and did not grow during extrapolation; mean-force removal conserved momentum to single-precision round-off.For weak Landau damping, energy errors were 1.5 × 10−5 versus 2.5 × 10−5; for two-stream, 1.1 × 10−2 versus 9.8 × 10−3. Momentum conservation was at 10−9–10−6 relative.Source: §8.2 · pp. 51, 53, Figure 9 · p. 53
• For the open-boundary cyclotron case, NEOPIC reproduced spiral-arm formation through training, but later field error rose to about 20% and the arms over-wound in the shear layer.The field error was 1–4% through training versus 0.5% for Barnes–Hut; at t = 11.98, arm winding was 14% greater overall, while winding for 0.14 ≤ r ≤ 0.29 agreed within 1%.Source: §8.2 · pp. 53–54, Figure 12 · p. 55

8.3. 3D-3V multi-app models

• In 3D, both multi-app models reproduced several field-energy trends, with two-stream growth rates within 2% of PIF but later and higher saturation peaks.Weak-Landau damping rates were −0.191 (Barnes–Hut), −0.190 (four-app) and −0.179 (five-app). Two-stream rates were 0.547 (PIF), 0.539 (four-app) and 0.537 (five-app); the saturation peak was t = 16.4 and Wz = 509 for PIF versus t = 17.7 and Wz = 594 and t = 17.3 and Wz = 642 for the models. Bump-on-tail’s linear phase was below the noise level at this particle count; strong-Landau field energies followed the reference through the nonlinear phase.Source: §8.3 · pp. 56–57
• Adding the Penning-trap teacher makes the response consistent across mini-apps and time, whereas averaging or switching between teachers can reduce accuracy or cause temporal inconsistency.The five-app model has better energy conservation in seven of eight cases, most clearly for bump-on-tail; constructing complementary training sets with smooth transitions is left for future work.Source: §8.3 · pp. 57–58
• The five-app model stays between Barnes–Hut and PIF fluctuation levels, while the four-app model switches from Barnes–Hut-like to PIF-like behavior in bump-on-tail.In weak Landau damping, PIF fluctuations are about six times Barnes–Hut’s (1.08 versus 0.17), and the five-app model is at 0.31. In bump-on-tail, the five-app model’s non-growing components remain at 25–33, while the four-app model changes from 17 to 65 and reaches 47–64 in extrapolation.Source: §8.3 · p. 57, Figure 14 · p. 59
• Within training, instability energy errors are 8–11 times PIF for the four-app model and 2–5 times PIF for the five-app model; the five-app error is smaller except in two-stream extrapolation.The two-stream extrapolation errors are 5.1 × 10^-3 for the five-app model and 3.5 × 10^-3 for the four-app model. NEOPIC conserves momentum to 10^-5 relative in these periodic mini-apps, unlike Barnes–Hut.Source: §8.3 · p. 58, Table 3 · p. 59
• In the Penning trap, NEOPIC’s total energy gains up to 6.6% during initial relaxation, then remains bounded with a maximum 5.0% change in extrapolation.At t = 60, self-field energy is 1.04 × 10^4 for NEOPIC versus 1.23 × 10^4 for PIC, while external potential energy is 4.7 × 10^3 versus 2.5 × 10^3.Source: Figure 15 · p. 60, §8.3 · p. 60

8.4. Computational cost

• NEOPIC is not competitive with optimized particle methods at equal particle counts, but its cost is independent of the teacher and particle distribution and depends mainly on architecture.Per particle, it is about 120–200 times more expensive than optimized FFT-based PIC and about six times more expensive than Barnes–Hut at 1.6 × 10^7 particles; its cost grows linearly with Ldv and hardly depends on mode count. The comparison with Barnes–Hut depends on opening angle and is meaningful only at matched accuracy.Source: §8.4 · p. 61
• NEOPIC’s NUFFT evaluation scales linearly: fitted slopes are 79 ns per particle in 2D and 300 ns in 3D, with intercepts of 64 ms and 49 ms, respectively.Direct evaluation costs 155 ns per particle in 2D and NUDFT costs 1.4 µs in 3D; direct is faster up to about 7 × 10^5 particles in 2D, while NUFFT is faster beyond about 4 × 10^4 in 3D.Source: §8.4 · p. 60, Figure 16 · p. 61

9. Conclusions

• NEOPIC's error separates into Monte Carlo sampling, time-stepping, and model errors, matching standard Monte Carlo particle methods up to model error.The cyclotron time-step study finds the same second-order time-stepping error for NEOPIC and PIF, while model error is independent of time step.Source: §9 · pp. 62–63
• Training-data particle noise acts as an implicit regularizer, guiding the model's response to density fluctuations toward its teacher and informing training particle counts.With more diverse training data, multi-teacher models behave as a consistent discretization rather than switching between teachers.Source: §9 · pp. 62–63
• The cost per time step is O(N_p), but the current single-GPU implementation is not competitive with optimized particle codes at equal particle counts.Cost per particle scales linearly with layer count times channel width and depends little on the number of Fourier modes.Source: §9 · pp. 62–63
• Against high-fidelity references, all models converge as N_p^-1/2 through 32 times the smallest training particle number, where teacher discretization error dominates.The teacher-forced error levels off at a knee two to four times the training particle count for models trained with fewer particles; training near the smallest expected simulation size is suggested.Source: §9 · p. 62
• A single 2D model trained on five mini-apps from three methods reproduces teacher damping and growth rates within about 1%, conserves momentum to round-off, and remains stable to twice the training horizon.On periodic mini-apps, energy conservation matches the teachers; the cyclotron field error grows beyond the training horizon.Source: §9 · p. 62
• The multi-app results motivate foundation models pre-trained across methods and applications, with few application-specific samples for fine-tuning.The framework is also proposed for extensions beyond electrostatics, including electromagnetic and collisional models.Source: §9 · p. 63

Appendix A.1. Min–max normalization

• The normalization box is set by the extreme particles, so its position dependence is Lipschitz in maximum displacement rather than average displacement and is not weakly continuous in the measure.For compactly supported densities the box converges to the support as particle count grows; for unbounded densities its extent grows slowly with particle count.Source: Appendix A.1 · p. 64
• In open domains, normalizing by the bounding box makes the input invariant to translations and componentwise rescalings of the particle cloud.The box itself depends on particle count, so the same physical density can map to slightly different normalized densities as that count changes.Source: Appendix A.1 · p. 64
• The Fourier transforms treat the normalized box as periodic, imposing artificial periodicity across an open-boundary cloud; the error analysis assumes fixed normalization.For min–max normalization, that analysis applies to configurations with nondegenerate bounding boxes and does not cover the box's particle-count dependence.Source: Appendix A.1 · p. 64

Appendix A.2. Grid aliasing and particle sampling

• The proposition compares uniform-grid coefficients after deposition with coefficients computed from particles to quantify grid aliasing versus particle sampling.It considers retained Fourier modes and a uniform grid with at least twice as many points per dimension as retained modes.Source: Appendix A.2 · pp. 64–65
• Grid deposition folds high frequencies coherently onto retained modes through an alias lattice, creating deterministic error that does not decrease with particle count.In grid-based FNOs, aliasing recurs in each layer because activations are not band-limited.Source: Appendix A.2 · p. 65
• Particle quadrature contributes high-frequency effects through sampling fluctuations, yielding unbiased error that decreases as N_p^-1/2 and is bounded independently of integrand smoothness or bandwidth.This Monte Carlo behavior applies at each network layer, with propagation through layers handled by Lemma 3.Source: Appendix A.2 · p. 65
• Unlike deposition-based grid FNOs, NEOPIC's NUDFTs neither alias particle noise into retained modes nor fold activation-generated harmonics back into them.For cloud-in-cell deposition, alias weights are small well below the grid Nyquist mode but become order one near it.Source: Appendix A.2 · p. 65

Appendix A.3. Translation equivariance

• Within the normalized box, a position-independent lifting makes the operator translation equivariant for every shift on the torus.The property holds exactly at every particle number, unlike grid-based layers, which are at best equivariant under grid-spaced shifts.Source: Appendix A.3 · p. 66
• The linear lifting of positions used in this work breaks translation equivariance within the normalized box, as a property of the learned operator rather than its discretization.In open domains, translating the whole cloud leaves normalized positions and therefore the field at particles unchanged.Source: Appendix A.3 · p. 66

Appendix B. Additional remarks on the error analysis

• NEOPIC’s sampling error scales as N_p^-1/2 in every dimension with a resolution-independent constant, while PIC also has a spatial-discretization error requiring Debye-length resolution.The PIC sampling rate depends on the grid; NEOPIC’s corresponding contribution is part of ε_teacher and does not constrain inference particle resolution.Source: Appendix B · p. 66
• For stable explicit order-m integrators, the time-stepping term becomes ∆t^m; implicit midpoint is covered only when each step is solved below local-error tolerance.The bound retains a step-size restriction and does not capture implicit stability at large steps; implicit solves require multiple field evaluations.Source: Appendix B · p. 66
• A periodic lifting preserves periodicity but cannot represent non-periodic behavior; open-boundary problems need non-periodic lifting information or padding beyond the particle cloud.The stated lifting choices do not affect the representation and error theorems.Source: Appendix B · p. 67
• For Dirichlet or Neumann walls, the vacuum field from imposed potentials can be solved once and added, leaving the network to learn the self-consistent field with homogeneous boundary conditions.The kernel is no longer translation invariant; position-dependent lifting or sine/cosine expansions can represent the boundary setting.Source: Appendix B · p. 67
• The training-noise analysis is leading-order in 1/N, assumes population squared loss and a reached minimizer, and warns that correlated snapshots do not provide independent validation data.Teacher runs with multiple seeds reduce frozen-noise fitting; mixed-particle-number training needs data at each number, and subsampling labels must be recomputed.Source: Appendix B · pp. 68–69

Appendix C. Diagnostics of the cyclotron mini-app

• The discrete canonical angular momentum is conserved exactly when the total particle torque vanishes; for 2D pairwise Coulomb forces, the constant virial makes L_can exactly conserved.For other fields, L_can changes through non-central-field torques and changes in the virial; the latter are smaller here and carry a ∆t^2 factor.Source: Appendix C · p. 69
• At ω_c∆t = 3, the Boris energy diagnostic measures gyration energy with a factor of about 0.31, explaining the Barnes–Hut energy oscillation as an integrator-and-diagnostic effect.The effect disappears when the gyration is resolved.Source: Appendix C · p. 69
• At t = 11.98, NEOPIC/Barnes–Hut winding ratios approach 1 at larger radii, while the excess winding is localized to the shear layer just outside the core.The ratio is 1.14 for 0.09 ≤ r ≤ 0.29 and 0.995 for 0.14 ≤ r ≤ 0.29.Source: Appendix C · p. 70

Appendix D. Cyclotron mini-app: convergence in the time step

• For ∆t ≤ 0.005, NEOPIC and PIF guiding-centre errors converge at local orders 1.7–2.0, with fitted order 1.75 for ∆t ≤ 0.0025.Their curves lie on top of each other, and the main error floor at smaller steps is NEOPIC’s model error of 3.5 × 10^-3 at t = 3.Source: Appendix D · p. 71
• At the training timestep, the time-stepping error is 9.9 × 10^-3, about three times the model error; reducing the timestep 32-fold does not increase model error.For ∆t ≤ 0.01, the model error is independent of timestep, and the error against the finest PIF run levels off for ∆t ≲ 0.005.Source: Appendix D · p. 71
• For ∆t ≤ 0.0025, NEOPIC’s relative energy and canonical-angular-momentum changes remain 3.2–3.4 × 10^-2 and 4.3–4.5 × 10^-2, respectively, independent of timestep.PIF’s changes fall below 10^-3 for energy and 10^-4 for L_can as ∆t decreases.Source: Appendix D · p. 71

Appendix E. Implicit time integration

• NEOPIC’s differentiable field map provides exact Jacobians by automatic differentiation, enabling combination with implicit integrators that require field-position derivatives.The passage contrasts this with implicit PIC methods that typically use Jacobian-free Newton–Krylov methods and finite-difference Jacobian-vector products.Source: Appendix E · p. 71
• With implicit midpoint, NEOPIC remains stable through ∆t = 2.5 and to twice the training horizon, with maximum relative energy change at most 1.2 × 10^-3 through ∆t = 2.At ∆t = 2.5 the maximum relative energy change is 1.1 × 10^-2; the ∆t = 2.5 run does not reproduce the damping rate because the plasma oscillation and resonant particles are unresolved.Source: Appendix E · p. 73, Figure E.18 · p. 74
• Leapfrog NEOPIC reproduces leapfrog PIF’s behavior at every tested timestep, with energy errors agreeing within 6%; both lose stability at ∆t ≥ 2.Beyond the leapfrog stability limit, field energy rises to about 7 and total energy changes by a factor of six to seven.Source: Appendix E · p. 73, Figure E.18 · p. 74
• At ∆t = 1, an implicit step costs 0.64 s versus 2.8 ms for leapfrog NEOPIC, and efficient implementation and higher-dimensional preconditioners remain future work.The experiment demonstrates compatibility with implicit integrators, not cost competitiveness.Source: Appendix E · p. 73

Appendix F. Cost and architecture

• The study timed 18 subnetworks of the trained 3D five-app model for 20 steps of weak Landau damping at 10^6 and 4 × 10^6 particles.Subnetwork sizes varied in layers, channels, and modes; the full-model timings reproduced those reported in Figure 16.Source: Appendix F · p. 75
• Across 18 architectures, cost per particle follows 2Ldv within 3%; time tracks Ldv, not parameter count, while halving m changes time by at most 4%.Parameters vary 256-fold while time varies 6.8-fold; the smallest subnetwork’s cost is comparable to Barnes–Hut and high-accuracy PIF, but it is untrained and its accuracy may be lower than the full model’s.Source: Appendix F · p. 75

Availability

• The implementation, experiment configurations, analysis scripts, figures, and tagged code versions are available in the cited repositories.Training datasets, trained models, and reference outputs are available from the corresponding author upon reasonable request.Source: Availability · p. 76

Declaration of generative AI and AI-assisted technologies in the manuscript preparation process

• The authors report using AI tools for manuscript and code-related work, while stating that they designed and ran the experiments and checked the AI-assisted scripts.They reviewed and edited the content and take full responsibility for the published article.Source: Declaration of generative AI and… · p. 76