Neural Particle Automata

• Neural Particle Automata are a particle-based extension of neural cellular automata featuring continuous positions and internal states for dynamic, local interactions.
• They use differentiable Smoothed Particle Hydrodynamics (SPH) operators to compute local densities, smoothing, and gradients, enabling asynchronous state and position updates.
• The framework leverages CUDA acceleration and stochastic update rules to efficiently simulate large-scale, self-organizing behaviors in morphogenesis, texture synthesis, and classification tasks.

Neural Particle Automata (NPA) are a Lagrangian extension of Neural Cellular Automata (NCA), generalizing the concept of neural automata from static lattice-based cells to dynamic systems of particles with continuous positions and internal states. NPA replaces the Eulerian paradigm—where state is tied to fixed pixels or voxels—with a particle-based approach wherein each particle is individually simulated and interacts locally via neural rules. This framework enables explicit individuation, heterogeneous dynamics across particles, and focused computation in regions of activity. A differentiable Smoothed Particle Hydrodynamics (SPH) formulation, with CUDA-accelerated kernels, allows NPA to scale efficiently to large particle ensembles and supports end-to-end learning of self-organizing behaviors across morphogenesis, texture synthesis, and classification tasks (Kim et al., 22 Jan 2026).

1. Formal Model Specification

NPA operates on a set of $N$ particles indexed by $i=1,\dots,N$ in $D$-dimensional space. Each particle is defined by a continuous position $x_i \in \mathbb{R}^D$ and an internal state $S_i \in \mathbb{R}^C$ representing hidden channels. The evolution of the system is governed by a shared, learnable local update rule $f_\theta$, typically implemented as a two-layer MLP. At each discrete timestep $t$:

1. For each particle $i$, a local perception vector $Z_i^t$ is constructed from its $\epsilon$-neighborhood $N(i)$:
   • $S_i^t$: current state
   • $\bar{S}_i^t$: density-normalized, locally smoothed state
   • $\nabla S_i^t$: state gradient
   • $\nabla \rho_i^t$: density gradient
   • Expressed as $$Z_i^t = [ S_i^t ; \bar S_i^t ; \nabla S_i^t ; \nabla\rho_i^t ]$$
2. The neural rule computes $\Delta y_i^t = f_\theta(Z_i^t)$, yielding updates for state (static) or both position and state (dynamic tasks).
3. A stochastic update is applied per particle using $\delta_i^t \sim \mathrm{Bernoulli}(p)$, typically $p = 0.5$ for asynchrony:
   • $S_i^{t+1} = S_i^t + \delta_i^t \cdot \Delta S_i^t$
   • $x_i^{t+1} = x_i^t + \delta_i^t \cdot \Delta x_i^t$

This model permits both heterogeneity and locality in updates, supporting individuation and dynamic adaptation even in unstructured particle ensembles.

2. Differentiable SPH Perception Operators

To enable fully differentiable, locality-aware perception and interaction, NPA utilizes SPH for estimating densities, smoothing states, and calculating gradients:

• Density Estimation:

$$\rho_i = \sum_{j \in N(i)} m_j W(\|r_{ji}\|; \epsilon)$$

where $r_{ji} = x_j - x_i$, kernel $W$ is Poly6 in 2D, and $m_i$ is usually normalized $(m_i = 1/N)$.

• State Smoothing:

$$\bar{S}_i = \sum_{j \in N(i)} \frac{m_j}{\rho_j} S_j W(\|r_{ji}\|; \epsilon)$$

• Gradient Estimation (Difference Form):

$$\nabla_0 S_i = \sum_{j \in N(i)} \frac{m_j}{\rho_j} (S_j - S_i) \otimes W^\nabla(r_{ji}; \epsilon)$$

• Gradient Correction (Moment Matrix):

$$M_i = \sum_j \frac{m_j}{\rho_j} r_{ji} W^\nabla(r_{ji}; \epsilon)^\top$$

$\nabla_1 S_i = M_i^{-1} \nabla_0 S_i$

In practice, this is implemented as: $$\nabla S_i = \nabla_0 S_i + \mathrm{detach}(\nabla_1 S_i - \nabla_0 S_i)$$ to ensure numerical stability in backpropagation.

• Density Gradient:

$$\nabla \rho_i = \sum_{j \in N(i)} m_j W^\nabla(r_{ji}; \epsilon)$$

SPH Kernel Functions

Kernel	Formula	Notes
Poly6 $W(r;\epsilon)$	$\frac{4}{\pi \epsilon^8}(\epsilon^2 - r^2)^3$	Support: $0 \leq \|r\| \leq \epsilon$
Spiky Gradient $W^\nabla(r;\epsilon)$	$$\frac{10}{\pi \epsilon^5}(\epsilon - \|r\|)^2 (r/\|r\|)$$	Support: $0 < \|r\| \leq \epsilon$

A uniform hash-grid bins particles by position for efficient neighbor querying; two CUDA kernel strategies are implemented:

• Particle-centric: One thread per particle, scanning adjacent grid cells.
• Grid-centric: One block per cell, staging per-cell particles in shared memory for coalesced access.

Neighborhoods are discovered on-the-fly; memory grows $\mathcal{O}(N)$ and computational cost is $\mathcal{O}(\text{avg neighbors})$ per particle.

3. Training Regime and Loss Construction

End-to-end NPA training employs backpropagation through SPH kernels and adapts several practices from NCA:

• Persistent state pool for sampling simulation rollouts.
• Variable rollout lengths with $T \sim [T_{\min}, T_{\max}]$.
• Overflow regularization with $L_2$ penalties enforcing $S_i \in [-1,1]$.
• Vector-input scaling for stability:

$$v \leftarrow \log(1 + |v|) \cdot \frac{v}{|v| + \eta}$$

• Stop gradients through positions $x$ for perception computations.
• Displacement regularization: penalizes $\sum_t \|\Delta x_i^t\|$ for smooth motion.

Representative Tasks and Losses

Task	Particle Count ($N$)	State Channels ($C$)	Loss Function
2D Morphogenesis	4096	16	$$L_{1+2}(D, D^*) + \lambda \cdot L_{1+2}(C, C^*) \cdot \text{detach}(e^{-|D-D^*|})$$
3D Morphogenesis	16384	24	Multi-scale SSIM + isotropy regularizer ($\|A_i\|_1$)
Texture Synthesis (particles)	4096	16	VGG OT-loss (on $C$ and $D$ across channels), power law transform $\overline{D} = (D+\epsilon)^\alpha, \alpha = 0.1$
Self-classifying Point Clouds	512	16	$L_2$ on one-hot labels, result: $98.42\%$ test accuracy

Each task utilizes customized rendering—Gaussian splatting for morphogenesis and textures, GSplat-based multi-view for 3D—plus curriculum strategies (progressive difficulty via blurred/sharpened targets) and per-task regularization.

4. Observed Dynamics and Empirical Behavior

NPA demonstrates several salient behaviors across evaluated tasks:

• Robustness to Discretization: Dynamics remain stable across varied $\epsilon \in [0.05, 0.2]$ and particle number $N$, degrading only under extreme downsampling.
• Stochastic Updates: Varying Bernoulli $p$ during inference (from $0.1$ to $1.0$) preserves morphogenetic convergence.
• Self-Regeneration: Local perturbations (state zeroing, slit cuts, clumping) are autonomously repaired by continued application of local update rules.
• Heterogeneous Multi-Species Interaction: Independently trained NPA rules interact within the same simulation, exhibiting cooperation, mixing, and disruption, despite no joint training.
• Hidden State Structure and Persistent Flows: Visualization of RGB-projected hidden states reveals internally consistent, directional channel patterns ("rainbow" flows) and persistent vortex-like movement post-target attainment, which suggests emergent memory/maintenance mechanisms.
• PointMNIST Dynamics: Hidden state clusters, tracked by UMAP, form per digit class over time, indicating emergence of global semantic structure from local exchanges.

5. Essential Algorithms and Mathematical Formulations

Key pseudocode for NPA update:

1. Build hash-grid index from {x_i} 2. For each particle i (in parallel): # Gather neighborhood N(i): all j in 3^D grid cells around x_i, with ||x_j - x_i|| < ε # Compute local perception via SPH ρ_i, \bar S_i, ∇S_i, ∇ρ_i via SPH sums # Scale vectors ∇S_i ← scale(∇S_i), ∇ρ_i ← scale(∇ρ_i) # Combine perception Z_i = [S_i; \bar S_i; ∇S_i; ∇ρ_i] # Execute neural rule Δy_i = f_θ(Z_i) 3. For each i: δ_i ∼ Bernoulli(p) S_i ← S_i + δ_i·ΔS_i ; x_i ← x_i + δ_i·Δx_i

Summary mathematical constructs central to NPA:

• SPH Poly6 kernel:

$$W(r;\epsilon) = \frac{4}{\pi \epsilon^8} (\epsilon^2 - r^2)^3$$

• Spiky gradient kernel:

$$W^\nabla(r;\epsilon) = \frac{10}{\pi \epsilon^5} (\epsilon - r)^2 \frac{r}{\|r\|}$$

• Logarithmic scaling:

$$\text{scale}(v) = \log(1+|v|) \cdot \frac{v}{|v|+\eta}$$

CUDA constants:

• Each SPH kernel (forward/backward) scans $3^D$ cells per particle.
• Particle-centric (thread/particle), grid-centric (block/cell) strategies—no explicit adjacency lists, neighborhoods discovered per step, memory growth $\mathcal{O}(N)$.

6. Theoretical and Practical Significance

NPA synthesizes compactness and robustness characteristic of NCA with the flexibility and resolution of particle-based, Lagrangian models. By leveraging differentiable SPH operators and efficient CUDA implementations, NPA accommodates large-scale, self-organizing particle dynamics with local computation and persistent individuation. Its empirical properties—stability under discretization, resilience to stochastic updates, self-regenerative capacity, and emergent structure—suggest potential utility in a range of simulation, synthesis, and classification domains. The particle framework inherently avoids global synchronization, handles unstructured domains, and empowers the learning of both robust and heterogeneous behaviors from local rules (Kim et al., 22 Jan 2026).