Low-Latency Activation-Regularized Sparse Neural Operators with Distillation Assistance Towards Real-Time Edge-Deployable Virtual Sensing

# Low-Latency Activation-Regularized Sparse Neural Operators with Distillation Assistance Towards Real-Time Edge-Deployable Virtual Sensing
Source: [https://arxiv.org/html/2608.23987](https://arxiv.org/html/2608.23987)
William HowesAffiliation:Grainger College of Engineering, Nuclear, Plasma & Radiological Engineering Department, University of Illinois Urbana\-Champaign, Urbana, IL, USAFarid AhmedAffiliation:Grainger College of Engineering, Nuclear, Plasma & Radiological Engineering Department, University of Illinois Urbana\-Champaign, Urbana, IL, USASyed Bahauddin Alam\*Affiliation:Grainger College of Engineering, Nuclear, Plasma & Radiological Engineering Department, University of Illinois Urbana\-Champaign, Urbana, IL, USAAffiliation:National Center for Supercomputing Applications, Urbana, IL, USA

###### Abstract

Virtual sensing enables digital twins and safety\-critical systems to reconstruct and forecast spatial\-temporal physics in real time\. However, conventional computational and data\-driven methods often face challenges in generalization, latency, and energy efficiency for edge deployment\. Neural operators offer a promising alternative but remain reliant on power\-intensive hardware\. Spiking neurons and neuromorphic computing can improve efficiency, yet surrogate\-gradient training and multi\-step spiking introduce convergence and latency challenges\. We propose the Sparse\-Activation\-ReLU \(SAR\) layer, a single\-step alternative that promotes activation sparsity without surrogate\-gradient training while remaining compatible with event\-based computing\. Within a trunk\-based NOMAD architecture, SAR achieves over a fivefold improvement in the combined Latency\-Error\-Energy \(LEE\) metric compared with Variable Spiking Neuron \(VSN\) and Leaky Integrate\-and\-Fire \(LIF\) implementations\. We further analyze spiking entropy and feature usage and introduce synthetic knowledge distillation, reducing the LEE score by more than twofold\. Finally, we improve VSN through a ReLU\-based spiking loss and graph\-neighbor thresholding\. On the Heat Exchanger dataset, these approaches reduce L2 error by more than twofold and nearly sevenfold, respectively, while reducing spiking and spatial aggregation\. Overall, the work presented is a step towards energy\-efficient virtual sensing by providing an alternative framework that can be positioned towards neuromorphic or other edge device integration that can be a gold standard to compare latency, energy, and error performance for future efficient designs that are sparsity or brain\-inspired spiking based\.

## 1Introduction

With the rapid advancement of artificial intelligence, the unprecedented computational demands of modern datacenter infrastructure have significantly increased global electricity consumption, renewing interest in nuclear energy as a reliable, carbon\-free, and scalable power source capable of supporting next\-generation computing technologies\[[1](https://arxiv.org/html/2608.23987#bib.bib1),[2](https://arxiv.org/html/2608.23987#bib.bib2)\]\. In particular, Small Modular Reactors \(SMRs\) and microreactors provide flexible, modular deployment strategies that enable efficient integration of nuclear power into distributed applications, including military installations, remote communities, industrial facilities, and other high\-power environments where conventional large\-scale nuclear reactors are impractical\[[3](https://arxiv.org/html/2608.23987#bib.bib3)\]\. As these reactor systems continue to decrease in physical scale while increasing in deployment flexibility, there is a corresponding need to improve monitoring, control, and operational awareness to reduce outage time, enhance safety margins, and improve overall economic viability\. Virtual sensing, specifically the real\-time spatial\-temporal reconstruction of governing physical fields through computational methods, offers a pathway toward achieving these goals by providing insight into system states that are either impossible or prohibitively expensive to measure directly\. By mapping sparse physical boundary measurements, such as inlet flow sensors, to dense spatial\-temporal fields including transient fluid flow, temperature, and pressure distributions, virtual sensing enables enhanced material degradation assessment, anomaly detection, predictive maintenance, and operational decision\-making throughout the lifetime of a nuclear system\. Existing approaches are primarily based on traditional numerical methods, such as the Finite Element Method \(FEM\), or machine learning frameworks including Physics\-Informed Neural Networks \(PINNs\)\[[4](https://arxiv.org/html/2608.23987#bib.bib4)\]\. While these methodologies have demonstrated success across a variety of scientific applications, they often either struggle to satisfy the strict real\-time requirements of online monitoring due to their computational complexity, limited ability to generalize to unseen operating conditions, or strict dependence on governing equations that may be unavailable or difficult to formulate for experimental systems\[[4](https://arxiv.org/html/2608.23987#bib.bib4),[5](https://arxiv.org/html/2608.23987#bib.bib5),[6](https://arxiv.org/html/2608.23987#bib.bib6)\]\. To overcome these limitations, neural operators have emerged as a powerful class of data\-driven architectures capable of learning mappings between infinite\-dimensional function spaces without requiring explicit equation integration, thereby enabling rapid inference for previously unseen boundary and initial conditions without retraining\[[7](https://arxiv.org/html/2608.23987#bib.bib7)\]\. Among the earliest and most influential neural operator architectures are DeepONet\[[8](https://arxiv.org/html/2608.23987#bib.bib8)\], which independently processes boundary conditions and evaluation coordinates through separate branch and trunk networks, and the Fourier Neural Operator \(FNO\)\[[9](https://arxiv.org/html/2608.23987#bib.bib9)\], which performs iterative spectral convolutions over structured computational grids\. More recent neural operator designs have extended these capabilities to increasingly complex multiphysics systems and irregular computational meshes\. Despite these significant advances in reconstruction accuracy and generalization, comparatively little attention has been devoted to computational efficiency\. Most existing neural operator studies assume the availability of datacenter\-class accelerators, such as NVIDIA A100 and H200 GPUs, whereas many practical virtual sensing applications require deployment on resource\-constrained edge hardware with strict latency, memory, and power budgets\[[10](https://arxiv.org/html/2608.23987#bib.bib10),[11](https://arxiv.org/html/2608.23987#bib.bib11)\]\. Consequently, the development of energy\-efficient neural operators capable of maintaining high reconstruction fidelity while substantially reducing computational cost remains an important and largely unexplored research direction\.

Spiking Neural Networks \(SNNs\)\[[12](https://arxiv.org/html/2608.23987#bib.bib12),[13](https://arxiv.org/html/2608.23987#bib.bib13),[14](https://arxiv.org/html/2608.23987#bib.bib14),[15](https://arxiv.org/html/2608.23987#bib.bib15)\]have emerged as one of the most promising paradigms for energy\-efficient artificial intelligence due to their inherently sparse and event\-driven mode of computation\. Inspired by biological neural systems, SNNs temporally integrate input and communicate information only when neurons emit spikes, allowing computation to occur selectively rather than through the dense activation patterns characteristic of conventional artificial neural networks\. This event\-driven spatial\-temporal behavior significantly reduces unnecessary computation by naturally exploiting the varying importance of individual neurons during inference\. Among the available neuron models, the Leaky Integrate\-and\-Fire \(LIF\) neuron has become one of the most widely adopted because it provides an effective balance between biological realism and computational efficiency\[[16](https://arxiv.org/html/2608.23987#bib.bib16)\]\. However, conventional LIF neurons exhibit limited performance on regression\-based tasks, including spatial\-temporal field reconstruction mainly due to precision loss, motivating the development of the Variable Spiking Neuron \(VSN\), which handles directly encoded continuous\-valued inputs, avoiding precision loss, while producing graded spike outputs on top of the memory/integration dynamics of the previous LIF neuron that preserve the computational advantages of event\-driven communication\[[17](https://arxiv.org/html/2608.23987#bib.bib17)\]\. Current neural operator implementations incorporating VSNs are almost exclusively trained using surrogate gradients, introducing an inherent mismatch between the forward and backward passes during optimization that frequently leads to unstable training dynamics and degraded reconstruction performance\[[18](https://arxiv.org/html/2608.23987#bib.bib18)\]\. An alternative strategy commonly explored within the neuromorphic community is ANN\-to\-SNN conversion, although these methods sometimes require long spike time windows and depend upon spike\-rate encoding and binary spike communication to accurately represent continuous information\[[19](https://arxiv.org/html/2608.23987#bib.bib19),[20](https://arxiv.org/html/2608.23987#bib.bib20)\]\. With the introduction of graded spike outputs in the VSN, which is shown to improve regression performance in neural operator applications, communication shifts away from rate\-based coding and back toward the direct transmission of continuous\-valued information, making it difficult to map the standard ANN\-SNN conversion formulation to the VSN framework\. Although VSN\-based neural operators have demonstrated encouraging performance for virtual sensing applications\[[21](https://arxiv.org/html/2608.23987#bib.bib21),[17](https://arxiv.org/html/2608.23987#bib.bib17),[22](https://arxiv.org/html/2608.23987#bib.bib22)\], they remain fundamentally constrained by surrogate\-gradient optimization\. Increasing the number of spike timesteps can improve reconstruction accuracy, as demonstrated later in this work, but does so at the expense of increased inference latency, greater computational cost, higher overall spiking activity, and consequently increased energy consumption\.

For this paper, we chose to stay within the variable/continuous communication framework to allow for ideal regression performance\. In addition, we focus on low\-latency performance\. The VSN has had initial success in terms of performance and potential efficiency and has even been demonstrated to perform within a single spike step\[[17](https://arxiv.org/html/2608.23987#bib.bib17),[21](https://arxiv.org/html/2608.23987#bib.bib21)\]but this formulation completely removes the intended memory dynamics\. With an interest to avoid the unstable surrogate training and stay within the one\-step realm and remove the unnecessary memory integration, keeping the low\-latency performance, we present a sparse\-activation layer within the neural operator framework that allows variable communication similar to VSN and leverages the natural sparsity induced by the ReLU activation function to achieve improved control over the energy\-accuracy tradeoff and latency performance\. Essentially, our proposed Sparse\-Activation\-ReLU \(SAR\) layer applies activation\-based sparsity regularization to the sparse output of the ReLU activation function, enabling explicit control over activation sparsity and consequently reducing neuron communication activity below that achieved by existing VSN\-based approaches at low\-latency while keeping variable communication key for regression neural operator performance and simultaneously maintaining superior reconstruction accuracy through surrogate\-free optimization\. The SAR layer framework can easily be mapped to an ANN\-to\-neuromorphic conversion formulation but with its removal of memory dynamics, our proposed work can also be directed towards any hardware/software that can leverage activation sparsity especially with the difficulty of neuromorphic implementation for densely connected neural networks\. Another issue, relatively unexplored within spiking operator research which has focused on attaching spiking layers to model architectures, is the feasibility of neuromorphic integration or general edge deployment\. To address this concern, we explore a synthetic neural operator distillation framework for edge\-deployable virtual sensing that centers around improving the performance of weaker, neuromorphic\-friendly or other general hardware portable and activation\-sparse operator architectures with synthetic knowledge distillation from sophisticated operator models that are less natural for edge\-deployment implementation, providing not only sparse communication and energy\-efficiency but low memory/latency requirements\. Knowledge distillation has had initial introduction within the neural operator field \(\[[23](https://arxiv.org/html/2608.23987#bib.bib23),[24](https://arxiv.org/html/2608.23987#bib.bib24)\]\), but we introduce the technique, specifically synthetic generation, with the unique perspective of neural operator edge deployment within neuromorphic/sparsity\-exploiting hardware, utilizing sophisticated, less device\-friendly architectures as teacher models\.

Although the SAR layer is not intended to replace conventional spiking dynamics, which may ultimately provide richer feature representations and improved modeling of transient physical phenomena once more effective optimization techniques become available, it serves as a practical low\-latency alternative for deployment in strict real\-time environments while also establishing a strong benchmark against which future spiking neural operator developments can be evaluated\. With this philosophy in mind and the interest of preserving more traditional spiking dynamics or flexibility in the variable spiking signal that is removed in the SAR layer, we also explore activation\-ReLU\-based regularization in the context of the Variable Spiking Neuron\. We explore two opportunities to reduce the latency and improve the performance of the VSN despite its surrogate gradient\. We first replace the original percentage based loss term with the activation\-based output of the neuron memory before spiking minus the threshold which is fed through a ReLU function, providing a more enriched loss term that explores the strength of neuron spiking event for improved performance without the need for longer spike trains\. In addition, we explore a ReLU\-based thresholding technique within the gating mechanism of the spatial aggregation layer for the Variable Spiking Graph Neural Operator \(VS\-GNO\)\[[21](https://arxiv.org/html/2608.23987#bib.bib21)\]that performs online filtering of a node’s neighbors, essentially removing unnecessary edges from the overall aggregation and improving computational speed and latency\. Although the original mechanism adapts the graph and forces gating weights towards zero it ultimately includes all neighbors in the overall computation\. As a result, we attempt neighbor thresholding with a separate activation\-based loss term to allow for a controllable accuracy\-edge count tradeoff, resulting in improved spatial aggregation efficiency\.

The contributions of this work are summarized as follows:

1. 1\.Surrogate\-free variable\-signal neural operators towards neuromorphic/sparse sensing\.We introduce the Sparse\-Activation\-ReLU \(SAR\) framework, which provides sparse, event\-based activation outputs for energy efficiency, rivaling traditional spiking models\. The SAR layer exists with a direct ANN\-to\-neuromorphic conversion methodology that eliminates surrogate\-gradient optimization while enabling explicit control over the energy\-accuracy tradeoff for neural operator\-based virtual sensing\. With the absence of memory dynamics, SAR can also be applicable towards non\-neuromorphic hardware/software that leverage the activation sparsity and reduce neuron\-neuron communication\.
2. 2\.Comprehensive sparsity analysis and regularization term comparison\.We investigate two activation regularization techniques, theL1L\_\{1\}norm and Hoyer norm, within sparse neural operators and compare their performance against surrogate\-gradient\-based VSN and LIF implementations across multiple spike timestep configurations as well introduce a novel entropy\-based characterization of sparse feature distributions and propose a unified evaluation metric that jointly balances reconstruction accuracy, inference latency, and energy efficiency\.
3. 3\.Synthetic knowledge distillation for assisting regression performance\.We present an initial synthetic distillation neural operator framework in which a sophisticated graph\-based teacher neural operator transfers knowledge to an activation\-regularized student model that is more natural to integrate on neuromorphic hardware or other edge devices that leverage sparsity\. Although simplified in the present work, this framework establishes a foundation for future synthetic data generation techniques capable of reducing dependence on computationally expensive high\-fidelity simulations\.
4. 4\.Activation regularization for traditional variable spiking neurons\.We extend the proposed activation\-based sparsity analysis back to variable\-spiking neural operators to improve performance without long spike trains and computational latency by replacing the original percentage\-based VSN regularization objective with a SAR\-inspired Hoyer\-based activation loss, demonstrating improved optimization behavior while highlighting the remaining limitations imposed by surrogate\-gradient training\. We also attempt activation ReLU\-based thresholding with VS\-GNO’s gating mechanism in order to remove unnecessary neighbors and improve the computational efficiencies of spatial layer in the variable spiking graph operator\.

Collectively, these contributions advance the development of energy\-efficient, edge\-deployable neural operators for scientific machine learning\. Beyond introducing a surrogate\-free reduction in activation communication that can map to neuromorphic, event\-based hardware or other sparsity\-leveraging device, this work provides new insight into activation\-based sparsity as an effective regularization mechanism, establishes novel analysis tools for evaluating sparse neural representations and performance, demonstrates the potential of synthetic knowledge distillation for reducing data requirements with neural operators while bridging the gap between spiking algorithms and realistic hardware integration, and investigates the applicability of activation\-based regularization within existing spiking neural operator architectures\. Taken together, these developments represent a meaningful step toward practical, low\-power virtual sensing systems capable of performing accurate real\-time spatial\-temporal reconstruction while reducing both computational cost and dependence on large\-scale training datasets\.

## 2Methods

### 2\.1Variable Spiking Neural Operators

To improve the regression capabilities of conventional LIF\-based spiking neural operators, previous work introduced the Variable Spiking Neuron \(VSN\), a modification of the standard LIF neuron that preserves event\-driven communication while transmitting continuous\-valued signals more suitable for scientific machine learning and operator regression tasks\[[17](https://arxiv.org/html/2608.23987#bib.bib17)\]\. Rather than transmitting only a binary spike, the VSN gates a continuous signal using the binary spiking event, allowing information to propagate only when the neuron fires while maintaining higher representational fidelity than traditional LIF neurons\. This formulation has previously been integrated into neural operator architectures, resulting in Variable Spiking Neural Operators capable of trading predictive accuracy for communication efficiency through sparse neuron activity\[[17](https://arxiv.org/html/2608.23987#bib.bib17),[21](https://arxiv.org/html/2608.23987#bib.bib21)\]\.

Given an input feature vectorz\(t\)z^\{\(t\)\}at spike time steptt, the membrane potential evolves according to

M\(t\)=β​M\(t−1\)\+z\(t\)y~\(t\)=\{1;M\(t\)≥Θ0;M\(t\)<Θ\}if​y~\(t\),M\(t\)←0y\(t\)=σ\(z\(t\)y~\(t\)\),whereσ\(0\)=0,\\begin\{gathered\}M^\{\(t\)\}=\\beta M^\{\(t\-1\)\}\+z^\{\(t\)\}\\\\ \\tilde\{y\}^\{\(t\)\}=\\left\\\{\\begin\{array\}\[\]\{lr\}1;\\quad M^\{\(t\)\}\\geq\\Theta\\\\ 0;\\quad M^\{\(t\)\}<\\Theta\\end\{array\}\\right\\\}\\quad\\text\{if \}\\tilde\{y\}^\{\(t\)\},M^\{\(t\)\}\\leftarrow 0\\\\ y^\{\(t\)\}=\\sigma\(z^\{\(t\)\}\\tilde\{y\}^\{\(t\)\}\),\\quad\\text\{where \}\\sigma\(0\)=0,\\end\{gathered\}\(1\)
whereβ\\betadenotes the membrane leakage parameter,Θ\\Thetarepresents the neuron threshold,y~\(t\)\\tilde\{y\}^\{\(t\)\}is the binary spike event\. When the membrane potential exceeds the threshold, a spike is generated, the membrane is reset, and the neuron transmits the continuous\-valued signaly\(t\)=σ⁡\(z\(t\)\)y^\{\(t\)\}=\\sigma\(z^\{\(t\)\}\)\. Otherwise, no communication occurs\. In contrast to conventional LIF neurons, this continuous output formulation can significantly improve regression performance while preserving sparse communication behavior required for neuromorphic implementations\[[17](https://arxiv.org/html/2608.23987#bib.bib17)\]\.

To regulate communication activity throughout the network, previous Variable Spiking Neural Operators introduce a spike\-percentage regularization term during training\. Rather than directly penalizing activation values, this objective minimizes the average fraction of observed spike events to total possible spikes across all neuron layers\. For a neural operator containingLLspiking layers, the overall optimization objective is written as

ℒ=α⋅ℒL​2\+γ⋅∑l=1Lℒspike\(l\),\\mathcal\{L\}=\\alpha\\cdot\\mathcal\{L\}\_\{L2\}\+\\gamma\\cdot\\sum\_\{l=1\}^\{L\}\\mathcal\{L\}\_\{\\mathrm\{spike\}\}^\{\(l\)\},\(2\)
whereℒL​2\\mathcal\{L\}\_\{L2\}denotes the reconstruction error,ℒspike\(l\)\\mathcal\{L\}\_\{\\mathrm\{spike\}\}^\{\(l\)\}is the average spike fraction measured at spiking layerll\(how many spikes occurred over how many spikes could have occurred\) over all spike time steps \(STS\) and the input batch,α\\alphacontrols the reconstruction objective, andγ\\gammadetermines the strength of the sparsity fraction\-based regularization\.ℒspike\(l\)\\mathcal\{L\}\_\{\\mathrm\{spike\}\}^\{\(l\)\}is essentially calculated by counting the binary output,y~\(t\)\\tilde\{y\}^\{\(t\)\}, at time steptt\. Increasingγ\\gammaencourages lower neuron firing rates, reducing communication and expected energy consumption while generally degrading reconstruction accuracy\. This formulation therefore provides a controllable accuracy\-efficiency tradeoff that has formed the basis of previous Variable Spiking Neural Operator architectures\.

![Refer to caption](https://arxiv.org/html/2608.23987v1/Figure1.jpg)Figure 1:Sparse\-Activation\-ReLU Model\(a\)Visual demonstration of the Variable Spiking Neuron\[[17](https://arxiv.org/html/2608.23987#bib.bib17)\]operating after a linear fully\-connected layer\. Synaptic weights and unique bias terms accumulate the input for each neuron at spike time stepttwhich is aggregated into the neuron’s leaky memory, mirroring LIF dynamics\. VSN replaces the binary spike output with an activation function operating on the neuron input\.\(b\)Sparse\-Activation\-ReLU for a fully\-connected linear layer\. We apply a ReLU activation function to the linear layer output which only allows positive signals subsequently operated on by a chosen activation function\. Such forward pass communication allows for surrogate\-free ANN\-to\-neuromorphic conversion that allows for variable signals \(unlike traditional binary spike rate to activation matching\), included in the original VSN framework\.\(c\)Sparse\-Activation\-ReLU for a generic input \(spectral/spatial convolution or norm\)\. We apply a ReLU activation function to the input subtracted by an optional threshold utilized for better spiking control\. Similar to the linear layers, the ReLU only allows positive signals subsequently operated on by a chosen activation function\.
### 2\.2Sparse\-Activation\-ReLU Alternative

The Variable Spiking Neuron \(VSN\) was originally introduced as a sparse communication mechanism with variable communication that regulates information flow through thresholded spiking behavior that is controlled by a spiking percentage loss term\. The VSN allows for improved regression performance compared to the binary LIF neuron\[[17](https://arxiv.org/html/2608.23987#bib.bib17)\]\. Previous attempts to implement VSN into operator architectures typically explore error and spiking behavior for a single spike time step \(1 STS\)\[[21](https://arxiv.org/html/2608.23987#bib.bib21),[17](https://arxiv.org/html/2608.23987#bib.bib17)\]\. Increasing spiking time steps adds increased memory and computational burden to training and risks increasing inference latency as well as energy consumption which can go against the real\-time, edge constraints that dominate virtual sensing\[[17](https://arxiv.org/html/2608.23987#bib.bib17)\]\. As a result, it seems prudent, when exploring regression\-based spiking neural operators for real\-time virtual sensing, to reduce spiking time to one step\. However, as a consequence, when operating in a single spike time step \(1 STS\), the temporal dynamics that traditionally motivate spiking neural networks no longer exist\. Since only a single spike evaluation is performed, there is no meaningful membrane memory to maintain and the leakage dynamics become irrelevant\. Consequently, the neuron functions primarily as a filtering operation, selectively allowing information to propagate when an activation exceeds a learned threshold\. This transition in framework leaves the VSN design as underutilized\. While a 1 STS VSN has shown promising effectiveness and improvement over LIF\-based neural operators\[[17](https://arxiv.org/html/2608.23987#bib.bib17)\], training the VSN still requires surrogate\-gradient approximations to back\-propagate through the non\-differentiable spike function\. This introduces a mismatch between the forward and backward passes and can complicate optimization\[[18](https://arxiv.org/html/2608.23987#bib.bib18)\]\. Traditional ANN\-to\-SNN methods tend to have long latencies and even alternatives that shorten latency requirements\[[20](https://arxiv.org/html/2608.23987#bib.bib20)\]are still based on binary communication, and for this paper we avoid such binary signaling and chose graded/variable spiking due to the previous shown improved regression performance that avoids precision loss\. Motivated by the observation that the 1 STS VSN acts primarily as a sparsity\-inducing filter, we replace the spiking neuron with a differentiable ReLU\-based activation sparsity layer that can be trained directly using conventional gradient descent and runs in only a single step with variable communication\. This layer does not rely on percentages but the direct sparsity\-norm calculation of output activations of the ReLU layer\.

Given an input feature vector𝐳\\mathbf\{z\}which is typically the result of a particular computational layer within the neural operator, the proposed Sparse\-Activation\-ReLU \(SAR\) Layer is simply defined as

𝐲=σ⁡\(ReLU⁡\(𝐳−𝝉\)\),where​σ​\(0\)=0\.\\mathbf\{y\}=\\sigma\(\\mathrm\{ReLU\}\(\\mathbf\{z\}\-\\boldsymbol\{\\tau\}\)\),\\quad\\text\{where \}\\sigma\(0\)=0\.\(3\)
where𝝉\\boldsymbol\{\\tau\}is an optional learnable subtraction parameter that can provide enriched learning dynamics towards sparse activations especially when𝐳\\mathbf\{z\}is not from a linear layer with a built in bias parameter\. For any linear layer input,𝐳=W​𝐱\+𝐛\\mathbf\{z\}=W\\mathbf\{x\}\+\\mathbf\{b\}, the included bias term acts as our thresholding, meaning we do not include𝝉\\boldsymbol\{\\tau\}in our calculation, but when applying SAR to spectral\-convolution, spatial aggregation layers, or after normalization layers without a clear bias term, we utilize𝝉\\boldsymbol\{\\tau\}as a learnable filtering term analogous to the thresholding behavior of the original VSN to allow for better control over spiking and efficiency\. The result of the ReLU function can then be fed into a subsequent activation functionσ\\sigmasuch thatσ⁡\(0\)=0\\sigma\(0\)=0, allowing for more enhanced nonlinear computations if desired \(e\.g\. softmax in transformer\) that does not add to the efficiency of the spiking framework of VSN\. This is where the variable communication seen in the VSN is partially preserved differing from traditional SNN methods that rely on binary communication\. The resulting variable activations𝐲\\mathbf\{y\}, which preserve the variable communication framework of VSN for improved performance, contain a level of sparsity \(zero\-count\) due to the definition of the ReLU function and are then passed to the subsequent network layers\. Regularization of the sparsity\-level of𝐲\\mathbf\{y\}can help achieve a controllable and desired sparsity level\.

Previous sparsity inducing methodologies similarly exploit the ReLU function but typically rely on fixing𝝉\\boldsymbol\{\\tau\}as a predefined value greater than 0, such as shifted ReLU\[[25](https://arxiv.org/html/2608.23987#bib.bib25)\], or manipulating the ReLU dynamics to admit only values above a threshold, post\-training, such as FATReLU\[[26](https://arxiv.org/html/2608.23987#bib.bib26)\]\. Both methods treat thresholding within the ReLU dynamics as a fixed hyperparameter while SAR aligns closer with the VSN by optimizing the threshold parameter\. These methods are also applied towards standard LIF/binary communication and have not been approached towards variable communication such as SAR\. Techniques such as shifted ReLU or FATReLU can be applied post training to further induce sparsity but risk severe degradation\.

Following training, the SAR layer can be incorporated into the neuromorphic systems with the following procedure: define synapse weights with the multiplicative weights utilized in computational layers within neural operators, combine the existing bias or the optional subtraction parameters if needed \(𝐛/𝝉\\mathbf\{b\}/\\boldsymbol\{\\tau\}\), and simply set a zero\-threshold spiking neuron that corresponds to the ReLU function designed to operate in one step \(no stored memory/leakage\)\. The resulting signal is simply a positive filter with a value, for a specific neuron/dimensionii, that is either zero \(no signal\) orσ⁡\(zi−τi\)/σ⁡\(zi\)\\sigma\(z\_\{i\}\-\\tau\_\{i\}\)/\\sigma\(z\_\{i\}\)if greater than zero, essentially meaning𝐲\\mathbf\{y\}from Equation[3](https://arxiv.org/html/2608.23987#S2.E3)is the signal sent aligning naturally with the originalσ\\sigma\-ReLU\-based ANN computation\. Figure[1](https://arxiv.org/html/2608.23987#S2.F1)visualizes the neuromorphic implementation of the SAR layer as well as the Variable Spiking Neuron\. It is important to note that because the ReLU is converted to the neuromorphic thresholding, ourσ\\sigmafunction does not add additional computation compared to the VSN preserving the original efficiency characteristics\. This provides a simple ANN\-to\-neuromorphic conversion mechanism without requiring surrogate\-gradient optimization during training\. Within this framework, activation sparsity becomes directly analogous to spiking activity, allowing standard sparsification techniques to regulate the communication behavior of the network\. In addition, the SAR layer operates only in a single step similar to the previously trained 1 STS VSN models\. As a result, SAR\-based neural operators provide a low latency, efficient alternative that can provide better power consumption and reconstruction performance due to the lack of surrogate gradient training\. The summarized dynamics of the SAR spiking layer for neuron/dimensioniiis below whereyi=0y\_\{i\}=0corresponds to no signal sent and no further downstream computation:

yi=σ⁡\(zi−τi\)ifzi−τi\>0otherwiseyi=0\.\\begin\{gathered\}y\_\{i\}=\\sigma\(z\_\{i\}\-\\tau\_\{i\}\)\\quad\\text\{if\}\\quad z\_\{i\}\-\\tau\_\{i\}\>0\\quad\\text\{otherwise\}\\quad y\_\{i\}=0\.\\end\{gathered\}\(4\)
Alternatively, hardware that supports sparse activations and reduced neuron\-neuron communication without memory dynamics unlike the existing neuromorphic devices \(e\.g\. Loihi 2\[[27](https://arxiv.org/html/2608.23987#bib.bib27)\]\) are another alternative for the SAR implementation and potentially a preferable one since neuromorphic hardware does not naturally match with the dense connectivity and high memory requirements of neural operators\.

A potential limitation of this formulation arises from its reliance on the ReLU activation which allows for the zero threshold neuromorphic conversion\. In the original VSN formulation, the spike gate is applied to the incoming signal with an optimized threshold that is isolated from the signal sent and can potentially allow negative values\. The resulting signal can also be operated on any subsequent activation function \(that is zero at zero input\) with a positive or negative value depending on the neuron input and threshold\. Our sparsity formulation includes a subsequent general activation function, but the output of Eq\.[3](https://arxiv.org/html/2608.23987#S2.E3)is strictly non\-negative, restricting all subsequent activations to positive\-valued inputs which could limit the expressive capabilities of the ReLU\-sparsity layer in terms of optimization\. If the thresholds employed by the original VSN are predominantly positive, this behavior may resemble the original VSN formulation, however, this is just speculation and may not hold at all and represents a potential limitation of the proposed approach\.

Consequently, future work should investigate alternative sparsity\-inducing layers capable of preserving both positive and negative signal information that also avoiding surrogate\-gradient optimization\. Such approaches may provide a closer approximation to the original VSN formulation while retaining the optimization benefits of conventional artificial neural networks\. Despite this limitation, the results presented in this work demonstrate improved or similar predictive performance relative to the original VSN framework across multiple activation functions, suggesting that the practical impact of this restriction, relative to the gradient mismatch issue, remains unclear and warrants further investigation\.

To regulate activation sparsity, previous methods explore initialization\[[28](https://arxiv.org/html/2608.23987#bib.bib28)\]as well as the previously mentioned fixed thresholding with shifted ReLU\[[25](https://arxiv.org/html/2608.23987#bib.bib25)\]and FATReLU\[[26](https://arxiv.org/html/2608.23987#bib.bib26)\], but to allow guaranteed and optimal control over the sparsification/efficiency level without hyperparameter tuning, regularization\-based training in the form of a sparsity loss term, similar to the VSN loss objection, is most effective\. We replace the average spike\-percentage regularization term utilized by the original VSN framework with norm\-based activation regularization applied directly to the output of Eq\.[3](https://arxiv.org/html/2608.23987#S2.E3)\. When applied post\-ReLU, the norm\-based terms will attempt to force vector elements towards zero which improves sparsity\. Two activation terms, previous utilized for activation sparsity\[[26](https://arxiv.org/html/2608.23987#bib.bib26)\], are explored:L1L\_\{1\}and Hoyer norm\. For a feature vector𝐯⁡\(𝐱i\)∈ℝd\\mathbf\{v\}\(\\mathbf\{x\}\_\{i\}\)\\in\\mathbb\{R\}^\{d\}corresponding to node or evaluation point𝐱i\\mathbf\{x\}\_\{i\}, theL1L\_\{1\}activation sparsity term is defined as

ℒL​1\(i\)=\|𝐯⁡\(𝐱i\)\|1=∑j=1d\|vj​\(𝐱i\)\|,\\mathcal\{L\}\_\{L1\}^\{\(i\)\}=\|\\mathbf\{v\}\(\\mathbf\{x\}\_\{i\}\)\|\_\{1\}=\\sum\_\{j=1\}^\{d\}\|v\_\{j\}\(\\mathbf\{x\}\_\{i\}\)\|,\(5\)
which encourages sparse activations by penalizing activation magnitude\. We additionally investigate the Hoyer sparsity measure,

ℒH\(i\)=\|𝐯⁡\(𝐱i\)\|12\|𝐯⁡\(𝐱i\)\|22\+ϵ=\(∑j=1d\|vj​\(𝐱i\)\|\)2∑j=1dvj​\(𝐱i\)2\+ϵ,\\mathcal\{L\}\_\{H\}^\{\(i\)\}=\\frac\{\|\\mathbf\{v\}\(\\mathbf\{x\}\_\{i\}\)\|\_\{1\}^\{2\}\}\{\|\\mathbf\{v\}\(\\mathbf\{x\}\_\{i\}\)\|\_\{2\}^\{2\}\+\\epsilon\}=\\frac\{\\left\(\\sum\_\{j=1\}^\{d\}\|v\_\{j\}\(\\mathbf\{x\}\_\{i\}\)\|\\right\)^\{2\}\}\{\\sum\_\{j=1\}^\{d\}v\_\{j\}\(\\mathbf\{x\}\_\{i\}\)^\{2\}\+\\epsilon\},\(6\)
whereϵ\\epsilonis a small constant introduced for numerical stability\. With the loss term defined for a single node𝐱i\\mathbf\{x\}\_\{i\}, we can subsequently average or sum over the provided grid to provide the sparsity loss contribution for the specific layer of compuation\. Both regularization strategies have previously been explored for activation sparsification\[[26](https://arxiv.org/html/2608.23987#bib.bib26)\]\. However, the Hoyer formulation provides improved control over sparsity due to its scale\-invariant nature, whereas theL1L\_\{1\}norm primarily suppresses activation magnitude across the network\. In the following section, we present the complete ANN\-to\-neuromorphic training methodology\.

![Refer to caption](https://arxiv.org/html/2608.23987v1/Figure2.jpg)Figure 2:Sparse\-Activation\-ReLU Nonlinear Manifold Decoder for Operator LearningArchitecture of SAR\-NOMAD specifically for the 2D Heat Exchanger\. In green, exists every SAR layer which naturally replaces the ReLU layers\. Since every computation before SAR is a simple linear layer, we do not include a threshold parameter which is redundant in the presence of the bias term\. As shown, the full\-field output provides theL2L\_\{2\}error for the total objective function while the individual sparse activation outputs of the ReLU \(from all SAR layers in the branch, trunk and combined networks not fully depicted in the figure\) with the Hoyer/L1L\_\{1\}loss term provide the control on energy efficiency\.![Refer to caption](https://arxiv.org/html/2608.23987v1/Figure3.jpg)Figure 3:Sparse\-Activation\-ReLU Graph Neural OperatorArchitecture of SAR\-GNO for the 2D Heat Exchanger\. In green, exists every SAR layer included in the neural operator which does not utilize the optional threshold parameter\. The green SAR layers replace either identity mappings or ReLU layers\. Components in dark purple represent SAR with the optional threshold parameters𝝉\\boldsymbol\{\\tau\}included for improved spiking control since they follow a normalization layer and non\-linear computational layers\. These threshold layers replace an identity mapping \(spatial\) and GeLU \(spectral\)\. SAR\-GNO takes in boundary input with the input embedding mappingMMand combines it with the geometry coordinates to produce an input for the latent projection mappingPP\. Subsequent spectral\-spatial blocks \(10 layers total\) provide global and local analysis that is combined through a collaboration layerff\. Final a downlift layerQQprovides the final full\-field reconstructed multi\-physics output\.
### 2\.3ANN\-to\-Neuromorphic Conversion for SAR Neural Operators

To evaluate the proposed SAR framework, we integrate the layer and activation sparsity regularization into both the Nonlinear Manifold Decoder for Operator Learning \(NOMAD\)\[[29](https://arxiv.org/html/2608.23987#bib.bib29)\]and Virtual Irregular Real\-Time Sparse Operator \(VIRSO\)\[[6](https://arxiv.org/html/2608.23987#bib.bib6)\]architectures, resulting in SAR\-NOMAD \(Figure[2](https://arxiv.org/html/2608.23987#S2.F2)\) and SAR\-GNO \(Figure[3](https://arxiv.org/html/2608.23987#S2.F3)\)\. The trunk\-branch structure of NOMAD naturally utilizes ReLU activations throughout the network, providing a straightforward transition to the proposed SAR implementation shown in Figure[2](https://arxiv.org/html/2608.23987#S2.F2)\. We additionally investigate VIRSO, which employed GeLU activations within the architecture\. To maintain as much consistency with the original VIRSO design, Sparse\-Activation\-ReLU layers are inserted at locations where VSNs were previously utilized shown in Figure[3](https://arxiv.org/html/2608.23987#S2.F3)in green, followed by a GeLU activation to preserve the original activation behavior as closely as possible\.

For conventional linear layers, which comprise the entirety of the NOMAD architecture, the SAR formulation is implemented as just a ReLU with the bias parameters providing a threshold\-like dynamic\. Within VS\-GNO, certain layers do not contain bias parameters, specifically spectral/spatial convolutions followed by normalization layers\. In these cases, optional subtraction parameter introduced in Eq\.[3](https://arxiv.org/html/2608.23987#S2.E3)is utilized to provide additional flexibility and improve training dynamics\.

Activation sparsity regularization is applied to the outputs of each SAR layer\. Specifically, the activations immediately following the ReLU operation are utilized when computing the sparsity penalty\. For a network containingLLSAR layers, the sparsity contribution from each layer is computed independently and then summed to obtain a global measure of network sparsity\. This produces a direct analogue to the energy\-accuracy tradeoff provided by the VSN framework, where increased sparsity corresponds to reduced neuron activity and, consequently, reduced communication within the network\.

The final training objective combines prediction accuracy and activation sparsity into a single loss function depicted in Figure[2](https://arxiv.org/html/2608.23987#S2.F2),

ℒ=α⋅ℒL​2\+γ⋅∑l=1Lℒactivation\(l\),\\mathcal\{L\}=\\alpha\\cdot\\mathcal\{L\}\_\{L2\}\+\\gamma\\cdot\\sum\_\{l=1\}^\{L\}\\mathcal\{L\}\_\{\\mathrm\{activation\}\}^\{\(l\)\},\(7\)
whereℒL​2\\mathcal\{L\}\_\{L2\}denotes the reconstruction error,ℒactivation\(l\)\\mathcal\{L\}\_\{\\mathrm\{activation\}\}^\{\(l\)\}represents either theL1L\_\{1\}or Hoyer activation sparsity regularization term computed at layerll\(this could be a layer in a branch network from a trunk\-branch model or a spectral block from a graph operator\) and averaged over the input batch,α\\alphacontrols the emphasis placed on predictive accuracy, andγ\\gammadetermines the strength of sparsity regularization\. Theℒactivation\(l\)\\mathcal\{L\}\_\{\\mathrm\{activation\}\}^\{\(l\)\}term, which is defined over the features for each grid point, might also be averaged over the evaluation domain as well as the batch, especially in the trunk and combined networks for NOMAD\. Increasingγ\\gammaencourages greater activation sparsity and therefore lower neuron activity throughout the network\. The sparsity hyperparameterγ\\gammacan be made unique for specific sparsity layers which is explored below\.

Training is performed entirely within the ANN framework using the proposed SAR layers and sparsity regularization terms\. Following training, the learned model can be transferred directly to a neuromorphic implementation operating with a threshold of zero\. Under this formulation, neuron spiking activity becomes directly equivalent to the sparsity pattern learned during ANN optimization, allowing the resulting neuromorphic system to reproduce the same communication behavior observed during training\. With the single\-step framework, SAR\-based layers not only reduce the latency of neuromorphic\-based inference, but the computational overhead during training, providing an overall faster alternative to large STS VSN\-based neural operators\.

![Refer to caption](https://arxiv.org/html/2608.23987v1/Figure4.jpg)Figure 4:Synthetic DistillationWe show the synthetic distillation framework utilized for neuromorphic virtual sensing\. A graph\-based VIRSO model, not native for neuromorphic hardware due to high connectivity and difficult integration, generates synthetic Heat Exchanger examples by randomly sampling input parameters for the Heat Exchanger dataset\. These synthetic examples are compared against SAR\-NOMAD’s predictions, a more neuromorphic friendly model, allowing for improvedL2L\_\{2\}error performance in SAR\-NOMAD, ideally keeping efficiency the same\.
### 2\.4Model Synthetic Distillation

Increasing sparsity within the proposed framework introduces a natural tradeoff with reconstruction accuracy, as higher sparsity restricts the amount of information that can be transmitted through the network\. This effect is particularly pronounced in regimes with limited training data, which is a common constraint in virtual sensing applications where data generation via finite element methods or similar high\-fidelity solvers is computationally expensive and time\-consuming\. Under such conditions, the reduced availability of training samples may prevent the model from learning an efficient communication strategy that is robust under high sparsity constraints\. Also, sophisticated model designs, such as FNO\[[9](https://arxiv.org/html/2608.23987#bib.bib9)\], Wavelet Neural Operator \(WNO\)\[[30](https://arxiv.org/html/2608.23987#bib.bib30)\], and VIRSO\[[6](https://arxiv.org/html/2608.23987#bib.bib6)\]that can better handle limited data, might present challenges when implemented within existing neuromorphic hardware \(e\.g\. Loihi 2\[[27](https://arxiv.org/html/2608.23987#bib.bib27)\]\) or other edge\-deployable hardware that handle our SAR framework but exist with strict memory constraints, so we are further limited with less sophisticated, more neuromorphic/edge\-friendly designs such as trunk\-branch networks which present with simple FCN layers and independent grid point evaluation\.

This behavior motivates the following unique neural operator framework: a non\-neuromorphic\-native or non\-edge\-friendly neural operator first learns an accurate approximation of the nonlinear operator mapping and subsequently transfers this knowledge to a more edge\-deployable neural operator architecture, thereby improving its performance beyond what can be achieved through direct training alone\.

To address this limitation, we employ knowledge distillation, a widely used technique in machine learning for transferring knowledge from a large teacher model to a more compact student model\[[31](https://arxiv.org/html/2608.23987#bib.bib31)\]which has seen initial introduction within neural operators\[[23](https://arxiv.org/html/2608.23987#bib.bib23),[24](https://arxiv.org/html/2608.23987#bib.bib24)\]\. In particular, we investigate synthetic distillation\[[31](https://arxiv.org/html/2608.23987#bib.bib31),[32](https://arxiv.org/html/2608.23987#bib.bib32)\], where additional training data is generated using a pretrained teacher model\. In this setting, a large and complex neural operator that is not directly suitable for sparse or neuromorphic deployment is used to produce synthetic field outputs corresponding to inputs sampled from a known distribution\.

Given that the input distribution for the target application is assumed to be known, new training samples are generated by drawing inputs from this distribution and evaluating them using a high\-capacity teacher model\. In this work, the teacher model is the graph\-based neural operator VIRSO without sparsity regularization\. VIRSO’s architecture renders it difficult to implement within edge\-deployable hardware, but it can be used to efficiently generate corresponding solution fields for generated synthetic inputs, which are then incorporated into the training dataset of the sparsified neural operator that better aligns with hardware integration\.

For this study, VIRSO\-generated samples are used to augment the training data for the SAR\-NOMAD implementation depicted in Figure[4](https://arxiv.org/html/2608.23987#S2.F4)where the main parameters that define the Heat Exchanger input is randomly generated and fed through both the teacher \(VIRSO\) and student \(SAR\-NOMAD\) models defining a training loss to improve SAR\-NOMAD’s performance\. While this procedure introduces additional computational overhead during training, it enables improved coverage of the input space and provides richer supervision for learning under sparsity constraints with no overhead during inference\. Further work is required to reduce the cost of synthetic data generation and training or to explore alternative distillation strategies; however, the focus of this paper is to demonstrate the transition from a large, complex, and non\-neuromorphic\-friendly model to a sparse architecture that is compatible with neuromorphic or other edge\-device deployment\. This framework is intended to improve the ability of the sparse model to approximate the underlying nonlinear mapping between input conditions and output fields under limited\-data, high\-sparsity regimes\.

### 2\.5Bringing Activation\-Regularized ReLU to Variable Spiking Operators

#### Activation\-ReLU Loss Term:

Although the SAR framework is motivated by the observation that many regression\-based neural operators with a VSN\-type neuron operate in the single spike time step \(1 STS\) regime to reduce latency, where the temporal dynamics of the VSN become largely unnecessary, this does not imply that the original VSN formulation cannot benefit from multiple spike time steps\. On the contrary, allowing neurons to accumulate information over time introduces temporal dynamics that can improve representational capacity\. In particular, for transient prediction problems, the persistent membrane memory of the VSN may better integrate temporal information and improve long\-term forecasting performance\. Moreover, VSN based operators provide flexibility in its signals, allowing negative output\. As a result, it would be also prudent to explore techniques to improve the performance VSN\-based operators without the need for high latency, allowing for memory dynamics and better signal expressiveness\. Motivated by these observations and beliefs, we additionally investigate improvements to the training methodology of Variable Spiking Neural Operators while preserving their underlying neuron dynamics\.

Rather than proposing an alternative to surrogate\-gradient optimization for neurons with persistent membrane memory, which remains an important direction for future work, we instead reconsider the sparsity objective employed during training\. Specifically, we hypothesize that the original spike\-percentage regularization term can be replaced with a more informative ReLU\-defined, activation\-based objective term similar to the SAR layer that takes difference between the current memory \(before any reset at time steptt\) similar to the input𝐳\\mathbf\{z\}from SAR and the trainable threshold𝚯\\boldsymbol\{\\Theta\}\(similar the parameter𝝉\\boldsymbol\{\\tau\}from SAR\)\. Instead of operating on binary spike events, we define a continuous pre\-spike activation corresponding to the membrane potential above the firing threshold,

𝐚\(t\)=ReLU​\(β​𝐌\(t−1\)\+𝐳\(t\)−𝚯\)∈ℝd,\\mathbf\{a\}^\{\(t\)\}=\\text\{ReLU\}\\left\(\\beta\\mathbf\{M\}^\{\(t\-1\)\}\+\\mathbf\{z\}^\{\(t\)\}\-\\boldsymbol\{\\Theta\}\\right\)\\in\\mathbb\{R\}^\{d\},\(8\)
whereβ​𝐌\(t−1\)\+𝐳\(t\)\\beta\\mathbf\{M\}^\{\(t\-1\)\}\+\\mathbf\{z\}^\{\(t\)\}is the membrane potential at timettimmediately before thresholding and any potential reset\. If a neuroniidoes not fire,ai\(t\)=𝟎a\_\{i\}^\{\(t\)\}=\\mathbf\{0\}\. Otherwise,ai\(t\)a\_\{i\}^\{\(t\)\}measures the amount by which the membrane potential exceeds the firing threshold before the membrane is reset\. This continuous activation therefore contains substantially more information than the corresponding binary spike event\.

The Hoyer sparsity measure is then applied directly to the activations defined in Eq\.[8](https://arxiv.org/html/2608.23987#S2.E8)\. Since zero\-valued activations correspond to neurons whose membrane potentials never exceed the firing threshold, encouraging sparsity in these activations naturally promotes reduced spiking activity\. Unlike the original spike\-percentage regularization, which is computed solely from binary firing events, the proposed objective exploits the continuous activation values above threshold, providing richer optimization information while maintaining the same sparsity objective\. During training, the original spike\-percentage regularization term in Eq\.[2](https://arxiv.org/html/2608.23987#S2.E2)is therefore replaced by the activation\-based sparsity objective following the same formulation as Eq\.[7](https://arxiv.org/html/2608.23987#S2.E7), with the Hoyer loss computed from the activations in Eq\.[8](https://arxiv.org/html/2608.23987#S2.E8)and averaged over the spike time steps\. The effectiveness of this activation\-based regularization strategy is evaluated in Section[3\.5](https://arxiv.org/html/2608.23987#S3.SS5)\.

#### ReLU Activations Towards Neighbor Thresholding in Spatial Graph Aggregation:

In addition to replacing activation\-regularized ReLU layers with the VSN to improve the training of low\-latency models, we extend the same thresholding concept to the graph spatial gating mechanism within variable spiking graph operators\. Within VIRSO\[[6](https://arxiv.org/html/2608.23987#bib.bib6)\], each edge \(neighbor\) connection is weighted through a gating function,γ⁡\(u,v\)\\gamma\(u,v\), which is computed as a function of the source nodeuu, destination nodevv, the associated edge attributes, and their Lipschitz positional encoding\. The learned gate is intended to adaptively construct the graph by assigning larger weights to more informative neighbors while suppressing less relevant connections\.

For all previous experiments, the graph connectivity was constructed using a fixed neighborhood size ofk=30k=30\. Consequently, every node aggregates information from 30 neighboring nodes regardless of their relative importance\. Although this dense aggregation improves representational capacity, it also increases computational cost and inference latency, particularly for future neuromorphic implementations where every neighbor corresponds to an additional synaptic event\. To alleviate this burden, we introduce a thresholded graph gating mechanism that selectively removes low\-importance neighbors\.

Specifically, we define a threshold parameter,τ\\tau, which is constrained to lie within the interval\[0,1\]\[0,1\]\. The threshold is subtracted from the learned gate value, and the resulting quantity is passed through a ReLU activation,

γ^​\(u,v\)=ReLU⁡\(γ⁡\(u,v\)−τ\)τ←clip⁡\(τ,0,1\)\.\\hat\{\\gamma\}\(u,v\)=\\mathrm\{ReLU\}\\left\(\\gamma\(u,v\)\-\\tau\\right\)\\quad\\tau\\leftarrow\\mathrm\{clip\}\(\\tau,0,1\)\.\(9\)
Only neighbors with gate values exceeding the threshold contribute to the spatial aggregation since a gate weight of zero represents no addition to the output\. The resulting thresholded gated aggregation output for the next convolution layerℓ\\ellcan therefore be expressed as

𝐡u\(ℓ\)=∑v∈𝒩⁡\(u\)γ^​\(u,v\)​𝐖𝐡v\(ℓ−1\),\\mathbf\{h\}\_\{u\}^\{\(\\ell\)\}=\\sum\_\{v\\in\\mathcal\{N\}\(u\)\}\\hat\{\\gamma\}\(u,v\)\\mathbf\{W\}\\mathbf\{h\}\_\{v\}^\{\(\\ell\-1\)\},\(10\)
where𝐖\\mathbf\{W\}denotes the linear weights applied to the spatial aggregation input before summation\.

During neuromorphic deployment, the gating network itself is removed, and the learned thresholded gate values are treated as fixed synaptic weights under the assumption of a static computational graph\. With our gating threshold technique, the spatial aggregation computation needed during neuromorphic deployed would then be significantly reduced since we only use a fraction of the original gating weight count before thresholding and forcing particular edge weights to zero\. If the underlying graph topology changes, the gating module must be recomputed; however, the proposed thresholding formulation should remain applicable which is a potential area of further investigation\.

We investigate two thresholding strategies\. The first employs a single trainable threshold shared across all graph edges\. The second predicts an edge\-dependent threshold by introducing an additional output from the gating network, allowing each edge to learn its own adaptive threshold value\.

To explicitly control the number of active neighbors, we introduce an additional sparsity loss weighted by the hyperparameterϵ\\epsilon\. This loss is defined as the average activation of the thresholded gate values over all samples, nodes, and neighboring edges,

ℒedge=1N​∑i=1N1\|𝒱\|​∑u∈𝒱∑v∈𝒩⁡\(u\)γ^i​\(u,v\),\\mathcal\{L\}\_\{\\mathrm\{edge\}\}=\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}\\frac\{1\}\{\|\\mathcal\{V\}\|\}\\sum\_\{u\\in\\mathcal\{V\}\}\\sum\_\{v\\in\\mathcal\{N\}\(u\)\}\\hat\{\\gamma\}\_\{i\}\(u,v\),\(11\)
whereNNis the batch size and𝒱\\mathcal\{V\}denotes the set of graph nodes\. Minimizing this objective encourages the network to reduce the number of active neighbors while preserving only the most informative graph connections\.

The edge sparsity loss is incorporated into the overall optimization objective with theℒL​2\\mathcal\{L\}\_\{L2\}reconstruction error and spiking percentage errorℒspike\\mathcal\{L\}\_\{\\text\{spike\}\}as

ℒtotal=α⋅ℒL​2\+γ⋅ℒspike\+ϵ⋅ℒedge\.\\mathcal\{L\}\_\{\\mathrm\{total\}\}=\\alpha\\cdot\\mathcal\{L\}\_\{L2\}\+\\gamma\\cdot\\mathcal\{L\}\_\{\\text\{spike\}\}\+\\epsilon\\cdot\\mathcal\{L\}\_\{\\mathrm\{edge\}\}\.\(12\)
The weighting parameterϵ\\epsilonprovides explicit control over the trade\-off between reconstruction accuracy and the number of active graph edges, enabling the model to learn computationally efficient graph representations while maintaining predictive performance\. This technique, which we demonstrate below, outperforms training on a graph with smaller neighbor/k value most likely because it leaves the spectral layer untouched by the reduced edge count and allows the option to focus on farther neighbors than only the k closes nodes\. Instead of generating a separate graph with smaller nodes, we reduce memory requirements and provide flexibility in neighbor choices\. Theoretically, this concept can applied to the original VIRSO graph network and any SAR\-based graph models which we have seen success in compared to lower k graph alternatives, but we chose to focus on the VSN implementation because of the emphasis of spiking in this paper and the increased usage of the spatial by the VSN\-based graph model over the SAR alternative which forces the spatial layer to minimal communication\. We found that SAR\-based models are best utilized in the spectral\-only form, completely dropping the unnecessary spatial block\. In all fairness, this is also true for the VSN implementation for our chosen benchmark which shows better regression performance with the spectral layer, but we present the option for spatial aggregation with improved computation in case local calibration is needed for future applications\.

### 2\.6Problem Formulation

Let𝒴⊂ℝd\\mathcal\{Y\}\\subset\\mathbb\{R\}^\{d\}denote thedd\-dimensional spatial domain of the physical system for which a solution field is sought\. The objective of neural operators is to learn the underlying nonlinear mapping between input boundary conditions and the corresponding physical solution fields through the following operator formulation:

𝒢:𝒰→𝒮,𝒢⁡\(𝐮\)​\(𝐱\)=𝐬⁡\(𝐱\),𝐱∈𝒴,\\mathcal\{G\}:\\mathcal\{U\}\\rightarrow\\mathcal\{S\},\\qquad\\mathcal\{G\}\(\\mathbf\{u\}\)\(\\mathbf\{x\}\)=\\mathbf\{s\}\(\\mathbf\{x\}\),\\qquad\\mathbf\{x\}\\in\\mathcal\{Y\},\(13\)
where𝒢\\mathcal\{G\}denotes a nonlinear operator acting between functional spaces\. The input is given by𝐮=\[u1,…,ub\]∈𝒰=∏i=1bℱ\\mathbf\{u\}=\[u\_\{1\},\\ldots,u\_\{b\}\]\\in\\mathcal\{U\}=\\prod\_\{i=1\}^\{b\}\\mathcal\{F\}, consisting ofbbpotentially multi\-modal input components\. The input space is defined such that eachuiu\_\{i\}may correspond either to a scalar quantity inℝ\\mathbb\{R\}or a functional input belonging to a space such asL2​\(D′\)L^\{2\}\(D^\{\\prime\}\)defined over an alternative domain\. The resulting output field is represented by𝐬⁡\(𝐱\)∈ℝk\\mathbf\{s\}\(\\mathbf\{x\}\)\\in\\mathbb\{R\}^\{k\}, describingkkphysical quantities at location𝐱\\mathbf\{x\}within the domain𝒴\\mathcal\{Y\}\. Consequently, the output space is given by𝒮⊂L2​\(𝒴,ℝk\)\\mathcal\{S\}\\subset L^\{2\}\(\\mathcal\{Y\};\\mathbb\{R\}^\{k\}\)\.

Kernel\-based neural operators, such as VIRSO, approximate the mapping𝒢\\mathcal\{G\}through a sequence of nonlinear integral transformations inspired by Green’s function formulations\[[7](https://arxiv.org/html/2608.23987#bib.bib7)\]\. The resulting iterative update can be expressed as

𝐯ℓ\+1​\(𝐱\)=σ⁡\(W​𝐯ℓ​\(𝐱\)\+∫𝒴𝒦ϕ​\(𝐱,𝐳\)​𝐯ℓ​\(𝐳\)​𝑑𝐳\),\\mathbf\{v\}\_\{\\ell\+1\}\(\\mathbf\{x\}\)=\\sigma\\left\(W\\mathbf\{v\}\_\{\\ell\}\(\\mathbf\{x\}\)\+\\int\_\{\\mathcal\{Y\}\}\\mathcal\{K\}\_\{\\phi\}\(\\mathbf\{x\},\\mathbf\{z\}\)\\mathbf\{v\}\_\{\\ell\}\(\\mathbf\{z\}\)d\\mathbf\{z\}\\right\),\(14\)
where𝐯ℓ\\mathbf\{v\}\_\{\\ell\}and𝐯ℓ\+1\\mathbf\{v\}\_\{\\ell\+1\}denote latent function representations evaluated at spatial location𝐱∈𝒴\\mathbf\{x\}\\in\\mathcal\{Y\}for layer indicesℓ=0,…,L\\ell=0,\\ldots,L\. The nonlinear activation functionσ\\sigmamay correspond to common choices such as ReLU or Sigmoid\. The termW​𝐯ℓ​\(𝐱\)W\\mathbf\{v\}\_\{\\ell\}\(\\mathbf\{x\}\)provides a learnable residual mapping, while𝒦ϕ\\mathcal\{K\}\_\{\\phi\}denotes a parameterized kernel responsible for propagating information across the spatial domain through the integral operator\. Collectively, these components enable the network to iteratively approximate the target nonlinear operator over the geometry of interest\.

VIRSO utilizes the Graph Fourier Transform to lift the input𝐯ℓ​\(𝐱\)\\mathbf\{v\}\_\{\\ell\}\(\\mathbf\{x\}\)to the spectral domain to allow for the multiplication of a directly parameterized kernel that is optimized to approximate the convolution integral\. VIRSO also utilizes a gated point\-wise spatial aggregation as a local approximation of the integral\. Combined with a projection mappingff, the spectral and spatial approximations provide both global and local analysis of the provided input while avoiding the previous scalability concerns\. With VIRSO, the boundary input is lifted, with a mappingMM, to a latent dimension and then copied to each evaluation coordinate\. Then, the entire grid with feature vectors equal to coordinates and latent condition representations is fed through a projection layerPP, to further lift the input into the defined intermediate dimension, and the subsequent convolution\-integral layers\. A final downlift layerQQmaps the grid’s intermediate feature vectors to the final output channels\. This encoder\-decoder structure processes the entire grid and evaluation points at once\.

Trunk\-branch architectures such as NOMAD provide a different approximation to Equation[13](https://arxiv.org/html/2608.23987#S2.E13)\. Multiple branch encoder networks \(typically fully\-connected\) are utilized to lift the boundary inputs \(profile and inlet values\) to individual embeddingsbi​\(ui\)∈ℝpb^\{i\}\(u\_\{i\}\)\\in\\mathbb\{R\}^\{p\}for boundary inputii\. In addition, a trunk network is utilized to process the input queries \(where we evaluate our solution\)\. It learns a set ofppbasis functions,t⁡\(𝐲\)∈ℝpt\(\\mathbf\{y\}\)\\in\\mathbb\{R\}^\{p\}, that are optimized for the problem and evaluation coordinate\. One or more evaluation points can be provided to the trunk branch, allowing for more natural handling of unseen evaluation coordinates\. Typically, models such as the Multi\-Input Operator Network \(MIONet\)\[[33](https://arxiv.org/html/2608.23987#bib.bib33)\]compute an element\-wise dot product between branches,b⁡\(u1,…,uk\)=b1​\(u1\)⊙…⊙bk​\(uk\)b\(u^\{1\},\.\.\.,u^\{k\}\)=b^\{1\}\(u^\{1\}\)\\odot\.\.\.\\odot b^\{k\}\(u^\{k\}\)for k boundary inputs, and a subsequent dot product with the trunk basis output,t⁡\(𝐲\)∈ℝpt\(\\mathbf\{y\}\)\\in\\mathbb\{R\}^\{p\}, to finalize an output solution for the given trunk evaluation points\. NOMAD recognizes that solutions might live on a low\-dimensional nonlinear submanifold\. This prompts NOMAD to map the individual branch embeddings and trunk basis \(usually concatenated\) to the final solution for each evaluation point with a nonlinear network \(typically fully\-connected\) called a "combined" network\.

### 2\.72D Heat Exchanger and Training Details

To evaluate the proposed sparsity framework, we consider the steady\-state flow analysis of a two\-dimensional cross section extracted from a three\-dimensional heat exchanger featuring enhanced heat transfer through a highly irregular dimpled surface and wavy tape insert geometry\[[34](https://arxiv.org/html/2608.23987#bib.bib34)\]\. The benchmark follows the boundary\-to\-field reconstruction framework defined in Equation[13](https://arxiv.org/html/2608.23987#S2.E13), where pressure and velocity fields are reconstructed at 3,977 evaluation nodes from only two scalar inlet conditions and a discretized heat flux profile\.

The heat exchanger benchmark represents a challenging evaluation with a high degree of geometric irregularity, complex vortex physics, and a four\-component coupled output field\. The geometry consists of a dimpled cylindrical channel with a wavy tape insert \(Figure[5](https://arxiv.org/html/2608.23987#S3.F5)a\), whose complex surface topology removes some rotational and reflective symmetries that could be exploited\. The wavy insert generates large recirculation regions and secondary vortex structures that significantly increase the complexity of the underlying flow physics and challenge operator learning architectures that lack sufficient local resolution capability\[[34](https://arxiv.org/html/2608.23987#bib.bib34)\]\.

The input consists of two scalar inlet conditions shown in Figure[5](https://arxiv.org/html/2608.23987#S3.F5), temperatureTi​nT\_\{in\}and axial velocityvi​nv\_\{in\}, together with a 100\-point discretization of an axial wall heat flux profile defined by the following equationA​sin⁡\(π​x/H\)A\\sin\(\\pi x/H\), whereHHis the height of the original 3D heat exchanger andAAis a provided amplitude input\. The output consists of four physical fields evaluated at 3,977 spatial locations on a two\-dimensional axial cross section𝒴h​x\\mathcal\{Y\}\_\{hx\}: pressurep⁡\(z,y\)p\(z,y\)and the three velocity componentsux​\(z,y\)u\_\{x\}\(z,y\),uy​\(z,y\)u\_\{y\}\(z,y\), anduz​\(z,y\)u\_\{z\}\(z,y\), resulting in the sparse\-to\-dense reconstruction ratio of3,977×4102≈156:1\\frac\{3,977\\times 4\}\{102\}\\approx 156:1\. The corresponding operator mapping can be expressed as

𝒢h​x:\(ℝ\)2×L2​\(ℝ\)→\(L2​\(𝒴h​x\)\)4,𝒴h​x⊂ℝ2\.\\mathcal\{G\}\_\{hx\}:\(\\mathbb\{R\}\)^\{2\}\\times L^\{2\}\(\\mathbb\{R\}\)\\rightarrow\(L^\{2\}\(\\mathcal\{Y\}\_\{hx\}\)\)^\{4\},\\qquad\\mathcal\{Y\}\_\{hx\}\\subset\\mathbb\{R\}^\{2\}\.\(15\)
In addition to the velocity components, we compute the resulting velocity magnitude from our predicted values and provide the loss in magnitude reconstruction as a regularization term during training and as a fifth performance metric during inference among the other four output channels\. The velocity magnitude regularization term was fixed atλm​a​g=0\.1\\lambda\_\{mag\}=0\.1for all experiments during training\. The term is defined by the following equation:

emag=λm​a​g​‖u^x2\+u^y2\+u^z2−u2‖2‖u2‖2,e\_\{\\text\{mag\}\}=\\lambda\_\{mag\}\\frac\{\\left\\\|\\hat\{u\}\_\{x\}^\{2\}\+\\hat\{u\}\_\{y\}^\{2\}\+\\hat\{u\}\_\{z\}^\{2\}\-u^\{2\}\\right\\\|\_\{2\}\}\{\\\|u^\{2\}\\\|\_\{2\}\},\(16\)
Training, validation, and test data were generated using ANSYS Fluent\[[35](https://arxiv.org/html/2608.23987#bib.bib35)\], resulting in 988 training examples, 248 validation examples, and 310 testing examples\. Due to the computational cost associated with generating high\-fidelity computational fluid dynamics solutions for complex geometries, this benchmark additionally serves as a representative virtual sensing application where limited training data motivates the use of synthetic distillation techniques\. We utilize a VIRSO model with 6 layers, width of 64, and 64 spectral modes as our teacher model to generate 1000, 2000, 4000, and 8000 synthetic examples that are added to the original 988 examples in our training dataset\. The same validation and test sets for the original performance are utilized for performance evaluation of our synthetic distillation framework\. We were able to generate synthetic inputs for our teacher model utilizing the data distributions employed for our ANSYS generated data, which are described below:

A∼𝒰⁡\(540,660\)​\[kW/m2\],Ti​n∼𝒰⁡\(536\.4,655\.6\)​\[K\],vi​n∼𝒰⁡\(4\.05,4\.95\)​\[m/s\]\.A\\sim\\mathcal\{U\}\(540,660\)\\,\[\\mathrm\{kW/m^\{2\}\}\],\\quad T\_\{in\}\\sim\\mathcal\{U\}\(536\.4,655\.6\)\\,\[\\mathrm\{K\}\],\\quad v\_\{in\}\\sim\\mathcal\{U\}\(4\.05,4\.95\)\\,\[\\mathrm\{m/s\}\]\.\(17\)
Two SAR\-based neural operator architectures were evaluated within this work: SAR\-NOMAD and SAR\-GNO shown in Figures[2](https://arxiv.org/html/2608.23987#S2.F2)and[3](https://arxiv.org/html/2608.23987#S2.F3)\. SAR\-NOMAD follows the original trunk\-branch formulation consisting of two branch networks, a trunk network, and a final combination network\. Each branch network and the combined branch contains four FCN layers while the trunk network contains three\. SAR layers are inserted after every linear layer except the final output layer in the combined network\. A hidden dimension of 256 was utilized for all sub\-networks, representing an input of3∗2563\*256for the combined network\. This architecture provides a natural integration of the proposed sparsity framework due to its existing reliance on ReLU activations throughout the network\. The optional subtraction parameter𝝉\\boldsymbol\{\\tau\}was not utilized since the bias parameters of the linear layers provide a natural threshold\-like dynamic after the application of linear weights\.

For SAR\-GNO, both a full and spectral\-only configuration are investigated\. Each model contains ten operator layers with a latent width of 64 and 100 spectral modes\. The spectral\-only model removes the spatial interaction block while retaining the spectral operator components\. SAR layers are incorporated throughout the architecture, including the input embedding layers, projection layers, downlift layers, spectral operator blocks, and spatial operator blocks\. Within the spectral and spatial layers which are unlike the straightforward linear layers and have normalization blocks, we utilized the optional subtraction parameter𝝉\\boldsymbol\{\\tau\}to enhance filtering by mimicking a threshold like behavior\. To preserve behavior consistent with the original VIRSO architecture, each sparsity layer is implemented as a ReLU filtering operation followed by a GeLU activation\.

For the sparsity loss term, we add the batch\-wise averaged activation loss for every sparsity layer in SAR\-NOMAD and SAR\-GNO, multiply ourγ\\gammahyperparameter, and then add this term to the L2 loss term to define our entire objective function\.

All sparsity activation regularization strategies discussed in Sections[2\.2](https://arxiv.org/html/2608.23987#S2.SS2)and[2\.3](https://arxiv.org/html/2608.23987#S2.SS3)are applied identically across both architectures\. For each SAR layer, activations immediately following the ReLU operation are collected and utilized in the sparsity regularization term\. The resulting sparsity penalties are summed across all layers and incorporated into the final loss function\. Training is performed entirely within the ANN framework\. As discussed before, conversion to a neuromorphic implementation is simple when operating with a zero\-threshold spiking mechanism\. This formulation allows us to compare directly the efficiency performance of SAR layers compared to traditional spiking layers such as VSN and LIF since the sparsity percentage \(essentially binary count of which dimension is zero or nonzero\) is identical in comparison to the spiking percentage\.

SAR\-NOMAD and SAR\-GNO were trained using the Adam optimizer with an initial learning rate of10−310^\{\-3\}and a batch size of 16\. Learning rate scheduling was performed using a step decay strategy with a step size of 40 epochs and a decay factor of 0\.5\. Early stopping was employed with a patience of 40 epochs with all results presented representing the best relative L2 error with the validation dataset\. It should be noted that all graph\-based models \(VIRSO, VS\-GNO, SAR\-GNO\), unless specified, operated on a graph generated from the Heat Exchanger grid using the KNN algorithm with a k value of3030nearest neighbors\. Moreover, the graph Laplacian utilized for the spectral block was chosen to be distance\-weighted\.

SAR\-NOMAD was trained for a maximum of 500 epochs with a weight decay coefficient of10−510^\{\-5\}while SAR\-GNO was trained for a maximum of 200 epochs with a weight decay coefficient of10−310^\{\-3\}\. Unless otherwise stated, all remaining training procedures, optimization settings, and evaluation methodologies were identical between SAR\-NOMAD and SAR\-GNO\. In this paper, we set theL2L\_\{2\}norm hyperparmaeterα\\alphato 1 in Equation[7](https://arxiv.org/html/2608.23987#S2.E7)while we vary theγ\\gammaparameter to explore different levels of sparsity and accuracy and demonstrate the controllability of our framework\. We explore both the and Hoyer loss terms and their effectiveness towards accurate and energy\-efficient neural operators\. Success is not defined by a specificL2L\_\{2\}error threshold but providing better accuracy at similar or improved efficiency compared to the VSN and demonstrating success in the controlling the accuracy\-energy tradeoff without severe reconstruction error degradation\. As mentioned, the relativeL2L\_\{2\}error, defined below withnninterior nodes and output channeloo,

erel=‖𝐬^o−𝐬o‖2‖𝐬o‖2=∑p=1n\(s^o​\(𝐱p\)−so​\(𝐱p\)\)2∑p=1nso​\(𝐱p\)2,e\_\{\\text\{rel\}\}=\\frac\{\\\|\\hat\{\\mathbf\{s\}\}\_\{o\}\-\\mathbf\{s\}\_\{o\}\\\|\_\{2\}\}\{\\\|\\mathbf\{s\}\_\{o\}\\\|\_\{2\}\}=\\frac\{\\sqrt\{\\sum\_\{p=1\}^\{n\}\\left\(\\hat\{s\}\_\{o\}\(\\mathbf\{x\}\_\{p\}\)\-s\_\{o\}\(\\mathbf\{x\}\_\{p\}\)\\right\)^\{2\}\}\}\{\\sqrt\{\\sum\_\{p=1\}^\{n\}s\_\{o\}\(\\mathbf\{x\}\_\{p\}\)^\{2\}\}\},\(18\)
was utilized as the accuracy\-based loss term for training and for performance evaluation during testing\. All training and evaluation were performed \(with spiking simulated by Snntorch\[[36](https://arxiv.org/html/2608.23987#bib.bib36)\]\) using an NVIDIA GH200 hardware provided through the DeltaAI cluster and an NVIDIA H200/A100 with the Delta cluster at the National Center for Supercomputing Applications \(NCSA\)\[[37](https://arxiv.org/html/2608.23987#bib.bib37)\]\. The resulting experiments provide a challenging assessment of SAR neural operators on a highly irregular, multi\-output reconstruction problem characterized by extreme boundary\-to\-field mapping requirements and complex nonlinear flow physics\. Lastly, as previously mentioned, we explore the performance of SAR\-NOMAD and VS\-NOMAD under synthetic distillation in Section[3\.4](https://arxiv.org/html/2608.23987#S3.SS4)and Figure[4](https://arxiv.org/html/2608.23987#S2.F4)\. Using trained weights from the VIRSO architecture previously optimized under the same training hyperparameters \(6 spectral\-spatial layers, 64 modes, and a width of 64\), we synthetically generated 2D Heat Exchanger samples then used them, alongside the original data used to train VIRSO, to improve the convergence of SAR\-NOMAD with a gamma value equal to0\.0010\.001,0\.0050\.005, and0\.010\.01\.

In addition to SAR\-based neural operators, we also present results from a variable spiking NOMAD \(VS\-NOMAD\) and variable spiking GNO \(VS\-GNO\)\[[21](https://arxiv.org/html/2608.23987#bib.bib21)\]that implement the VSN to provide as a spiking benchmark and demonstrate the enhanced training of the activation sparsity layer over the variable spiking implementation\. The VS\-GNO results utilized the same previously explored graph architecture\[[21](https://arxiv.org/html/2608.23987#bib.bib21)\]and training hyperparameters as its SAR\-GNO counterpart\. The VS\-NOMAD was also implemented with the same parameters as its SAR counterpart\. For activation functions, the VS\-NOMAD uses ReLU while the VS\-GNO utilizes GeLU layers\. The VSN was implemented with unique thresholds and leakage parameters along the hidden feature dimension \(256 for VS\-NOMAD and 100 for VS\-GNO\)\. We explored two surrogate gradient functions implemented by Snntorch\[[36](https://arxiv.org/html/2608.23987#bib.bib36)\]: fast sigmoid and arctangent\. We also explored three different spiking time steps for the VS\-NOMAD: 1, 10, and 20 STS to provide analysis of the spiking temporal dynamics and two different time steps \(1 and 10 STS\) for the spectral\-only VS\-GNO \(full version could not fit on the H200’s available memory\)\. The loss function utilized for the VSN\-based neural operators is regularized by the spiking percentages which is averaged over each layer within a sub\-component \(branch, spectral\-block, etc\.\) and then summed and multiplied by the spiking hyperparameterγ\\gammawhich provides the same energy\-accuracy tradeoff as the sparsity regularization term\. We also further explore in Section[3\.5](https://arxiv.org/html/2608.23987#S3.SS5)replacing the percentage\-based regularization term for variable spiking with an activation\-inspired Hoyer loss term utilizing a ReLU function applied to the difference between neuron memory and threshold at each spike pass\. We compare the resulting reconstruction error and spiking percentages with the original VSN framework for various different gamma values\. Moreover, we utilize the presented full VS\-GNO implementations for the 2D Heat Exchanger and experiment with our gated\-aggregation threshold technique in Section[2\.5](https://arxiv.org/html/2608.23987#S2.SS5)\. We utilize the sameγ\\gammachoises for the 1 STS full VS\-GNO as well as two values for our edge\-count hyperparameterϵ\\epsilon:0\.010\.01and0\.0010\.001\. We utilized the two presented gate threshold definitions and also present VS\-GNO results with KNN graphs having 5, 10, and 15 neighbors with a distance\-weighted graph Laplacian to demonstrate the superior performance and flexibility of our presented filter\.

Lastly, we also include the Leaky\-Integrate\-Fire neuron \(LIF\) as another surrogate gradient based alternative that essentially replaces the VSN \(and the activation function it utilizes\), defining the LIF\-NOMAD and LIF\-GNO neural operators\. We explore direct and rate based encoding for our boundary and grid input with the latter generated by Snntorch’sspikegen\.ratefunction\[[36](https://arxiv.org/html/2608.23987#bib.bib36)\]\. For the LIF\-NOMAD, the direct encoding results utilized the sameγ\\gammaparameters as VSN with 1, 10, and 20 STS while the rate encoding results were performed with 10 and 20 STS\. LIF\-GNO \(spectral and full\) explored 1 and 10 STS for direct encoding and 10 STS for rate encoding withγ\\gammavalue of00and0\.50\.5\. We only trained the spectral LIF\-GNO for 10 STS due to memory constraints of the H200\. For LIF\-NOMAD and LIF\-GNO, the arctangent surrogate gradient function was utilized\. In addition, all spiking percentages presented for the LIF\-based models are the result of the original calculated percentage divided by1\.71\.7in order to match the energy efficiency of the VSN which, due to its variable communication, does consume more energy based on a SpiNNaker2 analysis\[[17](https://arxiv.org/html/2608.23987#bib.bib17)\]\. This correction is admittedly crude and might present an exaggeration of the LIF’s energy reduction, but without physical hardware results it is difficult to implement\. Further integration on physical neuromorphic hardware itself is needed to truly test the energy efficiency performance of SAR, VSN, or LIF based models\.

### 2\.8Lid Driven Cavity and Training Details

To further test the performance of the activation\-based sparsity regularization, the Lid\-Driven Cavity \(LDC\) benchmark was utilized\. Unlike conventional field\-to\-field benchmarks, the LDC problem requires reconstructing high\-dimensional spatial fields from a low\-dimensional temporal forcing signal while operating on a structured computational mesh\.

The computational domain consists of a two\-dimensional square cavity, depicted in Figure[7](https://arxiv.org/html/2608.23987#S3.F7)a, in which fluid motion is generated by a moving upper boundary\. Rather than prescribing a constant lid velocity, the network receives a time\-varying lid velocity profile,V⁡\(t\)V\(t\), sampled over 90 discrete time steps, as the input shown in Figure[7](https://arxiv.org/html/2608.23987#S3.F7)a,b\. The target output consists of three coupled flow quantities evaluated over 4,225 interior spatial nodes: the pressure fieldp⁡\(x,y\)p\(x,y\), the velocity magnitude\|𝐯⁡\(x,y\)\|\|\\mathbf\{v\}\(x,y\)\|, and the turbulent kinetic energy fieldk⁡\(x,y\)k\(x,y\)\. This corresponds to reconstructing 12,675 output quantities from 90 input values, resulting in an approximate reconstruction ratio of141:1141\{:\}1\.

The underlying fluid dynamics are governed by the incompressible Reynolds Averaged Navier Stokes \(RANS\) equations together with a standardkk\-ε\\varepsilonturbulence closure,

∇⋅𝒗=0,\\nabla\\cdot\\boldsymbol\{v\}=0,\(19\)∂𝒗∂t\+\(𝒗⋅∇\)𝒗=−∇p\+∇⋅\[\(ν\+νt\)\(∇𝒗\+∇𝒗⊤\)\],\\frac\{\\partial\\boldsymbol\{v\}\}\{\\partial t\}\+\(\\boldsymbol\{v\}\\cdot\\nabla\)\\boldsymbol\{v\}=\-\\nabla p\+\\nabla\\cdot\\left\[\\left\(\\nu\+\\nu\_\{t\}\\right\)\\left\(\\nabla\\boldsymbol\{v\}\+\\nabla\\boldsymbol\{v\}^\{\\top\}\\right\)\\right\],\(20\)∂k∂t\+𝒗⋅∇k=Pk−ε\+∇⋅\[\(ν\+νtσk\)∇k\]\.\\frac\{\\partial k\}\{\\partial t\}\+\\boldsymbol\{v\}\\cdot\\nabla k=P\_\{k\}\-\\varepsilon\+\\nabla\\cdot\\left\[\\left\(\\nu\+\\frac\{\\nu\_\{t\}\}\{\\sigma\_\{k\}\}\\right\)\\nabla k\\right\]\.\(21\)
The objective of the neural operator is to approximate the nonlinear operator

𝒢l​d​c:L2​\(ℝ≥0\)→\(L2​\(𝒴l​d​c\)\)3,𝒴l​d​c⊂ℝ2\.\\mathcal\{G\}\_\{ldc\}:L^\{2\}\(\\mathbb\{R\}\_\{\\geq 0\}\)\\to\(L^\{2\}\(\\mathcal\{Y\}\_\{ldc\}\)\)^\{3\},\\quad\\mathcal\{Y\}\_\{ldc\}\\subset\\mathbb\{R\}^\{2\}\.\(22\)
Unless otherwise stated, all models are trained using a batch size of 16 for a maximum of 500 epochs with the Adam optimizer\. The initial learning rate is set to5×10−45\\times 10^\{\-4\}and is reduced by a factor of 0\.5 every 20 epochs using a StepLR scheduler\. A weight decay of10−310^\{\-3\}is employed for regularization, and all inputs and outputs are normalized using range normalization\. Early stopping is performed using a validation patience of 20 epochs to mitigate overfitting with all results presented representing the best relative L2 error with the validation dataset\. All training was performed on the Delta supercomputer at NCSA, utilizing NVIDIA A100 and H200 GPUs\.

To evaluate the proposed sparsity framework, experiments are performed using both the proposed SAR\-NOMAD architecture and the original VS\-NOMAD architecture\. Since the preceding Heat Exchanger experiments demonstrate \(shown below\) that the Hoyer activation sparsity regularization consistently provides a superior accuracy\-sparsity tradeoff compared to theL1L\_\{1\}formulation, only the Hoyer loss term is considered for SAR\-NOMAD throughout the remainder of this work\. For comparison, VS\-NOMAD is evaluated using both the fast sigmoid and arctangent surrogate gradient functions under the 1, 10, 20, and 30 spike step operating regimes\. We also implement LIF\-NOMAD with direct and rate encoding for the same time steps \(1, 10, 20, 30 for direct and 10, 20, 30 for rate\)\. Similar to the Heat Exchanger, we divided all calculated percentages from the LIF models by1\.71\.7to ideally match the efficiency of variable based alternatives\.

To ensure a fair comparison, all models employ an identical architecture consisting of a hidden width of 256 with three trunk FCN layers, four branch FCN layers, and four combined FCN layers\. Since the Lid\-Driven Cavity benchmark is parameterized solely by the temporal lid velocity profile, only a single branch network is required to encode the input forcing signal\. Consistent with the 2D Heat Exchanger experiments, the reconstruction weighting parameter is fixed atα=1\\alpha=1for all models, while the sparsity weighting coefficientγ\\gammais varied to investigate the resulting accuracy\-efficiency tradeoff\. For SAR\-NOMAD,γ\\gammacontrols the Hoyer activation sparsity regularization, whereas for VS\-NOMAD it controls the spike\-percentage regularization objective\. Increasingγ\\gammatherefore encourages reduced communication activity in both formulations, enabling a direct comparison of reconstruction accuracy as network sparsity increases\.

## 3Results

![Refer to caption](https://arxiv.org/html/2608.23987v1/Figure5.jpg)Figure 5:2D Heat Exchanger\(a\)We define the 2D Heat Exchanger geometry which is a cross\-sectional slice of a 3D grid at axial positionx=0\.789x=0\.789with a dimpled wall surface and tape insert utilized for better heat transfer, creating large vortex flow behavior and an irregular geometry with reduced symmetry\(b\)The boundary input utilized for our 2D Heat Exchanger\. It includes a wall heat profile defined along the axial direction of the original 3D geometry and two inlet quantities: temperature and axial velocity/speed\.\(c\)The 50th percentile test example for the SAR\-NOMAD model with the Hoyer loss term andγ=0\.005\\gamma=0\.005\. We see that the main physics and vortex behavior is closely captured with mostly stochastic absolute error behavior and only small absolute error hotspots located around the wall surface\.### 3\.1Heat Exchanger Performance \(NOMAD\)

Table 1:L2L\_\{2\}errors and spiking percentages for SAR\-NOMAD \(L1L\_\{1\}and Hoyer loss terms\) and VS/LIF\-NOMAD \(fast sigmoid and arctanget surrogate gradients\) with the 2D Heat Exchanger\.\*Training loss utilized 3γ\\gammafor Trunk and Comb terms andγ\\gamma/3 for Branch terms#### SAR\-NOMAD Results:

Table[1](https://arxiv.org/html/2608.23987#S3.T1)summarizes the accuracy of the different loss terms for SAR\-NOMAD \(L1L\_\{1\}and Hoyer\) together with different surrogate gradient functions for VS\-NOMAD on the 2D Heat Exchanger benchmark\. For SAR\-NOMAD the reported spiking percentage is equivalent to the spike activity obtained after neuromorphic conversion using a zero activation threshold, making the comparison to VSN/LIF spiking results easier\. As a result, we will discuss the sparsity of SAR layers in terms of spiking based off the previously discussed neuromorphic formulation\. The baseline NOMAD model achieves approximately 1% average relativeL2L\_\{2\}error across the output channels\. Despite being trained without any explicit sparsity regularization, the baseline architecture exhibits nearly 30% average spiking across the network components when converted to an event\-driven representation\. This indicates that the original ReLU\-based architecture is inherently sparse and that a substantial amount of unnecessary computation can be eliminated through neuromorphic\.

Table[1](https://arxiv.org/html/2608.23987#S3.T1)compares two activation sparsity regularization strategies, namely theL1L\_\{1\}activation loss and the Hoyer sparsity loss, over varying values ofγ\\gammawhile holdingα=1\\alpha=1\. Both approaches successfully reduce network spiking relative to the baseline model, resulting in more computationally efficient inference\. As expected, increasingγ\\gammadecreases the spiking activity while increasing the reconstruction error\. However, the manner in which this tradeoff occurs differs significantly between the two regularization strategies\.

For theL1L\_\{1\}activation loss, the reduction in spiking occurs gradually while the average relativeL2L\_\{2\}error remains below approximately 1% untilγ=0\.005\\gamma=0\.005\. At this point, the reconstruction error increases abruptly to over 7% while the average spiking percentage remains approximately 13%\. Examination of the component\-wise spiking statistics in Table[2](https://arxiv.org/html/2608.23987#S3.T2)explains this behavior\. The baseline model exhibits the largest spiking percentages within the trunk and combined sub\-networks, indicating that these components perform the majority of the informative computation\. TheL1L\_\{1\}loss primarily suppresses activity within the branch networks while leaving the trunk and combined networks largely unaffected\. Eventually, both branch networks approach zero spiking activity, effectively collapsing into constant bias terms because only the final layer bias remains active\. Consequently, the operator loses its ability to communicate boundary\-condition information through the branch pathways, resulting in the dramatic increase in reconstruction error observed atγ=0\.005\\gamma=0\.005\.

In contrast, the Hoyer sparsity loss exhibits a considerably more favorable sparsity\-accuracy tradeoff\. Asγ\\gammaincreases from10−410^\{\-4\}to0\.0050\.005, the average relative L2 error increases gradually from approximately 1% to approximately 5%, while the average spiking percentage decreases to below 5%\. Notably, the Hoyer regularization simultaneously achieves lower spiking activity and lower reconstruction error than theL1L\_\{1\}formulation atγ=0\.005\\gamma=0\.005especially\. Furthermore, the Hoyer loss substantially reduces activity within the trunk and combined sub\-networks, decreasing their spiking percentages by more than 20% across the explored regularization strengths\. This behavior explains why theL1L\_\{1\}loss maintains low reconstruction error for smallerγ\\gammavalues; although the branch sub\-networks are heavily suppressed, the majority of the information continues to be transmitted through the trunk and combined sub\-networks\. However, this also makes theL1L\_\{1\}regularization more difficult to tune since the computationally dominant sub\-networks remain relatively active\.

The superior behavior of the Hoyer loss can be understood by considering the objective that it optimizes\. Unlike theL1L\_\{1\}norm, which penalizes only the overall activation magnitude, the Hoyer loss explicitly promotes sparse activation patterns\. For example, the activation vectors\[0\.5,0\.5,0\.5,0\.5\]\[0\.5,0\.5,0\.5,0\.5\]and\[1,0,1,0\]\[1,0,1,0\]both have anL1L\_\{1\}norm equal to 2 despite the latter producing significantly fewer spikes\. Under the Hoyer metric, these vectors have values of 4 and 2, respectively, thereby favoring the genuinely sparse representation\. Although the Hoyer loss produces slightly larger reconstruction errors for some operating points, it provides significantly better control over the energy\-accuracy tradeoff\. Since the acceptable reconstruction error is ultimately application dependent, controllable sparsity/spiking dynamics are more valuable than simply minimizingL2L\_\{2\}error\. Furthermore, as demonstrated later through synthetic distillation, the Hoyer\-based models can recover additional accuracy while maintaining similarly low spiking percentages\.

Additional evidence for the superiority of the Hoyer formulation is obtained by independently weighting the activation penalties for each NOMAD component\. Rather than assigning an identical coefficient ofγ=5×10−4\\gamma=5\\times 10^\{\-4\}to every sub\-network, the trunk and combined sub\-networks were weighted using3​γ3\\gamma, while the branch sub\-networks were assignedγ/3\\gamma/3to more aggressively suppress the dominant communication pathways\. Under this modified weighting, theL1L\_\{1\}formulation reduced trunk and combined spiking by only 15% and 7%, respectively\. In contrast, the Hoyer formulation reduced trunk and combined spiking by 45% and 58%, respectively\. These results further demonstrate that the Hoyer loss directly targets sparsity rather than simply reducing activation magnitudes\.

In Figure[5](https://arxiv.org/html/2608.23987#S3.F5)c, we see the 50th percentile performance of SAR\-NOMAD with the Hoyer loss andγ=0\.005\\gamma=0\.005\. The main vortices caused by the 3D geometry’s wavy tape insert are shown in high detail\. The absolute error is mostly stochastic in nautre with only small hotspots located on the wall surface\. This visualization shows that with SAR\-based layers, we can still accurately simulate the governing physics while providing highly efficient computation\.

Table 2:Per\-layer spiking percentages for each network component for Sparse\-NOMAD and VS\-NOMAD with the 2D Heat Exchanger\.\*Training loss utilized 3γ\\gammafor Trunk and Comb terms andγ\\gamma/3 for Branch termsTable[2](https://arxiv.org/html/2608.23987#S3.T2)provides the layer\-wise spiking percentages for the SAR\-NOMAD models and the 2D Heat Exchanger\. For the trunk network and the second branch network, which process the spatial coordinates and inlet boundary conditions, the largest spiking activity generally occurs within the first and final hidden layers\. The elevated activity in the first layer is expected because this layer lifts a low\-dimensional input into the higher\-dimensional latent representation, requiring denser communication to preserve representational capacity\. Increased activity in the final layer, which is also observed within the combined network, may indicate progressively stronger feature selection as the network approaches the output representation\. Although this interpretation remains speculative, it suggests that different layers may exhibit distinct functional roles within the sparse operator and warrants further investigation\.

More importantly, the layer\-wise analysis highlights the tendency of theL1L\_\{1\}regularization to completely suppress the branch sub\-networks\. Atγ=0\.005\\gamma=0\.005, both branches contain layers with zero spiking activity, causing each branch to collapse into a constant bias term independent of the input boundary conditions\. In contrast, the Hoyer regularization maintains nonzero activity throughout the branch sub\-networks\. The minimum observed spiking percentage of approximately 0\.39% corresponds to roughly one active neuron out of 256, allowing boundary\-condition information to continue propagating through the network while substantially reducing computational cost\. This preservation of sparse yet meaningful communication provides another explanation for the improved accuracy obtained using the Hoyer formulation\.

Table 3:Average fractional feature\-spiking entropy for each evaluation node for each network component of Sparse\-NOMAD and VS\-NOMAD with the 2D Heat Exchanger\.\*Training loss utilized 3γ\\gammafor Trunk and Comb terms andγ\\gamma/3 for Branch termsFinally, Table[3](https://arxiv.org/html/2608.23987#S3.T3)presents an entropy\-based analysis of feature utilization during spiking inference on the 2D Heat Exchanger\. For each layer, spike counts are accumulated across the entire test dataset to produce a feature\-wise spike count vector where for each feature we have the number of spikes that this feature neuron emits\. Normalizing by the total number of spikes along the feature vector yields a probability distribution,

piℓ=siℓ∑j=1dsjℓ,p\_\{i\}^\{\\ell\}=\\frac\{s\_\{i\}^\{\\ell\}\}\{\\sum\_\{j=1\}^\{d\}s\_\{j\}^\{\\ell\}\},\(23\)
wheresis\_\{i\}denotes the total spike count of featureiiin layerℓ\\ellandddis the feature dimension for layerℓ\\ell\. The corresponding spiking entropy for layerℓ\\ellis

Hℓ=−∑i=1dpiℓlog\(piℓ\),H^\{\\ell\}=\-\\sum\_\{i=1\}^\{d\}p\_\{i\}^\{\\ell\}\\log\(p\_\{i\}^\{\\ell\}\),\(24\)
which is normalized by the maximum possible entropy,

Hnormℓ=Hℓlog⁡\(d\)\.H^\{\\ell\}\_\{\\mathrm\{norm\}\}=\\frac\{H^\{\\ell\}\}\{\\log\(d\)\}\.\(25\)
For the trunk and combined sub\-networks, the normalized entropyHℓH^\{\\ell\}for a feature vector in layerℓ\\ellis also averaged over all node points evaluated in that layer \(3,977 for the 2D Heat Exchanger\)\. A normalized entropy approaching one indicates that communication is distributed across nearly the entire feature space, whereas values approaching zero indicate that the network communicates through only a small subset of intermediate latent features\.

![Refer to caption](https://arxiv.org/html/2608.23987v1/Figure6.jpg)Figure 6:Spiking Activity ResultsWe depict the probability/frequency of different feature dimensions spiking for the SAR\-NOMAD model with the Hoyer loss andγ=0\.005\\gamma=0\.005\. For the branch layers, the sparsity regularization forces only a small subset, as a low as one dimension, to fire and communication information, indicating that the model is forced to collapse it feature dimensions in order to provide low spiking\. The trunk and combined subnetworks have a lot more activity across the entire dimension\.The entropy analysis reveals a consistent trend across all SAR\-NOMAD models\. Asγ\\gammaincreases and the overall spiking percentage decreases, the normalized entropy also decreases, indicating that the model increasingly relies on a smaller subset of latent features for communication\. This behavior is particularly evident within the first branch network\. UnderL1L\_\{1\}regularization, the entropy rapidly decreases to zero as the branch sub\-networks completely collapse, particularly atγ=0\.005\\gamma=0\.005, where both branches lose all expressive capacity\. Conversely, the Hoyer formulation generally maintains nonzero entropy within the branch sub\-networks, with occasional zero entries only atγ=0\.005\\gamma=0\.005, corresponding to communication through a single active feature\. This behavior again aligns with the improved reconstruction accuracy observed for the Hoyer models\. The entropy of the Hoyer results still show a collapse of layer dimensions to a small subset with the branch networks further indicated by the spiking activity results in Figure[6](https://arxiv.org/html/2608.23987#S3.F6)for the Hoyer run withγ=0\.005\\gamma=0\.005\.

The entropy analysis presents both opportunities and challenges\. Low entropy suggests considerable potential for structured network pruning since many latent features, and as a result synapse weights, are rarely utilized\. In fact, for the branch networks in the Hoyer results could easily be eliminated to the single feature and single vector of synapse weights utilized based off our test run, saving a lot of memory\. However, it also indicates that the model may be minimizing spike activity by reducing the effective latent dimensionality rather than learning richer threshold\-based dynamics distributed across the full feature space\. Although preliminary experiments explored entropy as an additional regularization objective, limited improvements were observed, indicating that further investigation is necessary\.

#### VS/LIF\-NOMAD Results:

In addition to SAR\-NOMAD, we investigate the VS\-NOMAD and LIF\-NOMAD architectures using the spiking percentage loss introduced in\[[17](https://arxiv.org/html/2608.23987#bib.bib17)\]\. Two surrogate gradient functions, fast sigmoid and arctangent, are evaluated together for the VSN implementation with spike time step \(STS\) values of 1, 10, and 20 and regularization coefficients ofγ=0\\gamma=0and10−610^\{\-6\}to assess the effect of introducing temporal dynamics into VS\-NOMAD\. For the LIF\-NOMAD model, we explore the arctangent function with direct encoding for 1, 10, and 20 STS and rate encoding for 10 and 20 STS\.

As expected, increasingγ\\gammareduces the average spiking percentage \(except for 20 STS LIF\-NOMAD with direct encoding\) while increasing the reconstruction error, as summarized in Table[1](https://arxiv.org/html/2608.23987#S3.T1)\. A clear distinction emerges between the surrogate gradient functions for the VSN implementation\. Fast sigmoid produces lower average spiking percentages, ranging from approximately 6\-12%, but suffers from substantially larger mean reconstruction errors ranging from 20\-46%\. Conversely, arctangent produces significantly lower reconstruction errors between approximately 8\-35% while exhibiting higher spiking percentages between approximately 7\.5\-15%\. Although fast sigmoid achieves greater computational efficiency, its severe degradation in reconstruction accuracy makes arctangent the preferable surrogate function, particularly since it still achieves substantially lower spiking activity than the baseline NOMAD model\. Comparing VS\-NOMAD and LIF\-NOMAD, we find that the rate\-encoded implementations with reduced precision for input and neuron\-neuron communication for LIF\-NOMAD resulted in severe error degradation compared to the far better reconstruction performance with the VS\-NOMAD with meanL2L\_\{2\}errors over 67% even though the equivalent spiking percentages are reduced\. For example, VS\-NOMAD with arctangent for 20 STS andγ=10−6\\gamma=10^\{\-6\}does have around 1\.7 times more spiking that the most efficient rate\-encoded LIF implementation but over 5 times improvedL2L\_\{2\}error\. Contrastingly, the results with the LIF\-NOMAD model utilizing direct encoding indicate a different story due to the improved precision of the model input\. Although VS\-NOMAD withγ=0\\gamma=0for 10 and 20 STS resulted with lowerL2L\_\{2\}error, the overall performance of the LIF\-NOMAD with direct encoding is at least similar to arctangent VS\-NOMAD with similar error and similar or lower spiking percentage\. Directly encoded LIF\-NOMAD presents no clear performance distinction from VS\-NOMAD unlike its rate\-encoded counterpart which utilizes reduced precision in its input\.

Increasing the number of spike time steps generally improves reconstruction accuracy \(not for LIF\-NOMAD, directly encoded for 10 to 20 STS\), suggesting that temporal dynamics enrich the latent representation and improve the reconstruction of the underlying physics\. However, the average spiking percentage also increases because additional spike time steps allow greater membrane potential accumulation and consequently more spike events\. With inference on an NVIDIA A100, we find that 1, 10, 20, and 100 STS result in an average latency of 7, 57, 114, and 556 milliseconds \(approximately linear scaling\)\. With applications requiring real time performance and potentially multiple auto\-regressive calculations for long\-term forecasting, a low STS with lower latency might be desirable in terms of edge\-deployment and integration of energy\-efficient neural operators\. In addition, for the 2D Heat Exchanger, we were unable to train an STS of 30 due to memory constraints for a single H200, further demonstrating the difficulty of dealing with high spike step dynamics and the benefits of reducing spiking to the single step regime\.

Despite the performance of LIF and VSN implementations for NOMAD, the Hoyer/L1L\_\{1\}\-based SAR\-NOMAD provides significantly improved reconstruction and efficiency results over the surrogate gradient alternatives\. The best\-reconstruction\-performing VS/LIF configuration employs the arctangent surrogate and VSN with ten spike time steps andγ=0\\gamma=0, producing approximately 8% average relativeL2L\_\{2\}error with roughly 15% average spiking\. The lowest spiking VS/LIF result is the 10 STS LIF\-NOMAD with rate encoding which presented with 73\.39% averageL2L\_\{2\}error and 6\.41% spiking\. In comparison, the Hoyer\-based SAR\-NOMAD model atγ=0\.005\\gamma=0\.005achieves approximately 5\.41% average relativeL2L\_\{2\}error while requiring only approximately 4\.87% average spiking and only a single spike step\. Not only does the SAR based model perform with better power consumption and reconstruction error than all VS/LIF\-NOMAD results, it also significantly reduces the latency that might be required by VS/LIF\-NOMAD to achieve similar performance\. This performance gap is likely attributable to the surrogate gradient mismatch that commonly limits the optimization of spiking neural operators\. In conclusion, SAR\-NOMAD provides an efficient alternative that avoids the long spike trains and surrogate gradient training while providing desired energy consumption performance\. SAR\-NOMAD is guaranteed to be the closest to the desired speed for real\-time monitoring with only a single spike step, making it an ideal candidate for efficient virtual sensing and a benchmark towards proposed spiking models in the future\.

An additional observation is that VS\-NOMAD primarily suppresses activity within the first branch network and the combined network rather than within the second branch network like SAR\-NOMAD and LIF\-NOMAD\. Since the combined network integrates information from all preceding sub\-networks before projecting to the nonlinear output manifold, excessive suppression of this component may explain the larger degradation in reconstruction accuracy compared to SAR\-NOMAD and LIF\-NOMAD with direct encoding\.

Table[2](https://arxiv.org/html/2608.23987#S3.T2)further illustrates the layer\-wise spiking activity for VS\-NOMAD using the arctangent surrogate gradient, the best performing function\. Similar to SAR\-NOMAD, the highest activity occurs within the first hidden layer of both branch networks and the trunk network due to the lifting of low\-dimensional inputs into the latent space\. However, unlike the Hoyer\-based SAR\-NOMAD models, the VS\-NOMAD architecture permits complete collapse of the first branch network for the one\-STS configuration withγ=10−6\\gamma=10^\{\-6\}, where zero spiking is observed within the first hidden layer for fast sigmoid and arctangent\. Consequently, this branch reduces to a constant term independent of the heat profile input\. The Hoyer\-based SAR\-NOMAD models, which have no gradient mismatch avoid this collapse completely\.

Finally, Table[3](https://arxiv.org/html/2608.23987#S3.T3)presents the same entropy\-based analysis for arctangent VS\-NOMAD\. Consistent with the SAR\-NOMAD observations, reduced spiking activity corresponds to lower normalized entropy, indicating that the network communicates through a progressively smaller subset of latent features as spike activity is suppressed\. Although entropy for VS\-NOMAD is higher than the SAR models, the spiking percentages of the variable spiking results is higher\. All indications show that with lower spiking both VS\-NOMAD and SAR\-NOMAD tend to reduce entropy and reduce the scale of the feature dimensions utilized for inference\.

Overall, the results for the 2D Heat Exchanger strongly demonstrate the potential for SAR\-based spiking as an alternative for surrogate gradient trained VSN/LIF operators\. Not only does it provide improved accuracy and spiking performance, the ReLU layer also removes the need for high STS and keeps latency down closer to the desired real\-time properties for virtual sensing applications\. At the very least, as improvements in training for traditional spiking neural operators occur, the SAR layer can provide a gold\-standard benchmark in terms of efficiency, latency, and accuracy reconstruction\.

### 3\.2Heat Exchanger Performance \(VIRSO\)

Table 4:L2L\_\{2\}errors and spiking percentages between the spectral\-only SAR\-GNO and spectral\-only VS/LIF\-GNO with the 2D Heat Exchanger\[[21](https://arxiv.org/html/2608.23987#bib.bib21)\]\. \(Thresh\.\) corresponds to the inclusion of the threshold parameter for SAR and an additional GeLU layer\. \(No Thresh\.\) corresponds to removing the threshold and keeping the GeLU layer while \(Thresh\. No GeLU\) represents the removal of the GeLU function while keeping the threshold\.Mean Spiking Percentages \(%\)Spiking ModelMean Rel\.L2L\_\{2\}\(%\)MPSpectralSpatialfQMeanArctan\. LIF\-GNOγ=0\\gamma=0DE 1 STS10\.9325\.5131\.8721\.52\-\-14\.8523\.44Arctan\. LIF\-GNOγ=0\.5\\gamma=0\.5DE 1 STS10\.850\.1\.173\.37\-\-9\.193\.43Arctan\. LIF\-GNOγ=0\.5\\gamma=0\.5DE 10 STS9\.710\.1\.173\.80\-\-8\.423\.35Arctan\. LIF\-GNOγ=0\\gamma=0RE 10 STS21\.007\.816\.4418\.47\-\-11\.3511\.01Arctan\. LIF\-GNOγ=0\.5\\gamma=0\.5RE 10 STS16\.390\.0\.0\.94\-\-6\.541\.87Fast Sig\. VS\-GNOγ=0\\gamma=01 STS0\.4010\.5314\.6425\.69\-\-15\.8416\.68Fast Sig\. VS\-GNOγ=0\.5\\gamma=0\.51 STS0\.490\.782\.1720\.06\-\-8\.817\.96Fast Sig\. VS\-GNOγ=0\.5\\gamma=0\.510 STS0\.452\.341\.7020\.73\-\-9\.478\.56Arctan\. VS\-GNOγ=0\\gamma=01 STS0\.368\.649\.4127\.42\-\-13\.0314\.62Arctan\. VS\-GNOγ=0\.5\\gamma=0\.51 STS0\.470\.393\.7513\.10\-\-4\.625\.47Arctan\. VS\-GNOγ=0\.5\\gamma=0\.510 STS0\.481\.563\.3013\.98\-\-4\.675\.88SAR\-GNOγ=0\.001\\gamma=0\.001\(Thresh\.\)0\.453\.958\.8511\.04\-\-24\.5412\.10SAR\-GNOγ=0\.003\\gamma=0\.003\(Thresh\.\)0\.540\.391\.576\.67\-\-16\.246\.22SAR\-GNOγ=0\.001\\gamma=0\.001\(No Thresh\.\)0\.564\.187\.4516\.50\-\-24\.5713\.18SAR\-GNOγ=0\.003\\gamma=0\.003\(No Thresh\.\)6\.320\.0\.8\.42\-\-15\.726\.03SAR\-GNOγ=0\.001\\gamma=0\.001\(Thresh\. No GeLU\)0\.413\.1711\.6213\.53\-\-30\.2014\.63SAR\-GNOγ=0\.003\\gamma=0\.003\(Thresh\. No GeLU\)0\.462\.554\.569\.51\-\-24\.2110\.21

Table 5:L2L\_\{2\}errors and spiking percentages between the full SAR\-GNO architecture and the Full VS/LIF\-GNO architecture with the 2D Heat Exchanger\[[21](https://arxiv.org/html/2608.23987#bib.bib21)\]\. \(Thresh\.\) corresponds to the inclusion of the threshold parameter for SAR and an additional GeLU layer\. \(No Thresh\.\) corresponds to removing the threshold and keeping the GeLU layer while \(Thresh\. No GeLU\) represents the removal of the GeLU function while keeping the threshold\.Following the success of SAR layer within NOMAD, we extend the approach to the graph\-based VIRSO architecture by inserting SAR layers followed by GeLU activations, consistent with the original VIRSO design, resulting in the proposed SAR\-GNO architecture\. We evaluate SAR\-GNO on the 2D Heat Exchanger benchmark and directly compare its performance against the previously proposed VS\-GNO using the VSN formulation\[[21](https://arxiv.org/html/2608.23987#bib.bib21)\]and LIF\-GNO which utilizes the traditional LIF neuron\. Consistent with the VS\-GNO study, both the spectral\-only and full VIRSO architectures are investigated for VSN, LIF, and SAR implementations\. Because SAR\-NOMAD existed with already built in bias terms within its linear layers, we chose not to include the optional threshold𝝉\\boldsymbol\{\\tau\}\. SAR\-GNO has non\-linear layers with normalization blocks such as the output of the spectral block and the spatial aggregation\. With those layers, we decided to add the threshold which we test its performance in this section\.

Tables[4](https://arxiv.org/html/2608.23987#S3.T4)and[5](https://arxiv.org/html/2608.23987#S3.T5)summarize the average relative L2 error across all output channels together with the average spiking percentages for each component of the original VIRSO architecture for the spectral\-only and full version of LIF/VS\-GNO and SAR\-GNO\. For the spectral\-only architecture \(Table[4](https://arxiv.org/html/2608.23987#S3.T4)\), the previously proposed VS\-GNO achieved strong performance compared to LIF\-GNO, obtaining average relativeL2L\_\{2\}errors below 0\.5% while reducing the average spiking percentage to approximately 5\.5% for the arctangent implementation and 8% for fast sigmoid\. For the graph architecture, the VSN demonstrates far improved performance over the LIF implementation unlike the VS\-NOMAD results\. The LIF\-GNO did achieve lower spiking with the 1 STS directly encoded model withγ=0\.5\\gamma=0\.5having 3\.43% spiking and the 10 STS,γ=0\.5\\gamma=0\.5rate encoded model having 1\.87% spiking, but the error is over 20 times higher, representing severe performance degradation\. The VSN, with its variable communication most likely naturally aligns with the hybrid SNN/ANN structure of VS\-GNO, provides far better reconstruction performance while also still significantly reducing spiking albeit potentially not as low as the LIF implementations although the 1\.7 correction factor is not as accurate for a hybrid architecture unlike the purely spiking NOMAD model since there exists standard downstream computations that don’t have reduced energy between VSN or LIF since they would be identical\. The concept of the 1\.7 different originated from analysis of how much energy a LIF/VSN neuron needs for computation\[[17](https://arxiv.org/html/2608.23987#bib.bib17)\]which is the most relevant for a completely spiking model such as VS/LIF\-NOMAD\. Further integration is needed to understand the true difference in energy consumption\.

For the full\-version of VS\-GNO, we see even more improved performance over the LIF\-GNO implementation which has, at best, twice as less spiking but over 20 times worse reconstruction performance, strongly indicating the superior communication enrichment from the Variable Spiking Neuron\.

With the SAR layer, we see similar or slightly improved performance over the VSN implementation\. For the spectral\-only architecture, atγ=0\.003\\gamma=0\.003, SAR\-GNO with a threshold parameter𝝉\\boldsymbol\{\\tau\}and added GeLU layer achieves an average relativeL2L\_\{2\}error below 0\.6% while maintaining an average spiking percentage around 6\.22%\. This spiking is lower than the fast sigmoid VSN implementations but the SAR\-GNO’s reconstruction performance is slightly worse \(0\.05% higher than the lowest spiking fast sigmoid VSN result\), concluding with an overall similar performance between SAR\-GNO and VS\-GNO\. For the arctangent implementation of VS\-GNO, we have slightly betterL2L\_\{2\}\(at most 0\.18%\) and spiking \(at most 0\.75%\) than the SAR\-GNO implementation\. In addition, for both fast sigmoid and arctangent, the 10 STS implementation for VS\-GNO resulted in similarL2L\_\{2\}error with worse spiking percentages indicating that the 1 STS VS\-GNO architecture with no memory temporal dynamics was the preferred option\. Comparing both a step SAR layer and the 1 STS VS\-GNO, the results potentially indicate that the spectral\-only VSN implementation was able to avoid surrogate gradient instability and provided slightly better communication enrichment with the ability to permit negative signals unlike the SAR layer and its dependence on the ReLU function\. Despite the comparison, the SAR implementation with the spectral\-only architecture is essentially the same, and the reliance on ReLU only provides a slight difference in performance\. Although further investigation into sparsity\-inducing layers that allow negative signals might be warranted, the magnitude of the improved communication’s impact is not clear, and in the context of surrogate gradient functions, might be surpassed by the training instability\.

For the full graph architecture, the SAR\-GNO provides slightly improved results over VS\-GNO\. For the fast arctangent implementation of VS\-GNO, as gamma increases the spiking is reduced but the mean relativeL2L\_\{2\}error degrades to over 5%\. Meanwhile, SAR provides a reconstruction performance of at most 4% with spiking below 3% and far below the arctangent spiking performance\. In contrast, the fast sigmoid VSN implementation, somewhat surprisingly, provides lowL2L\_\{2\}error with a spiking performance as low as 9\.82% with 0\.83% error, avoiding the same degradation that the arctangent implementation experienced\. While SAR\-GNO provides spiking below 3%, the resulting L2 is slightly higher reaching at most 4%\. Despite this, we see that atγ=0\.01\\gamma=0\.01with the threshold and attached GeLU layer, the full SAR\-GNO architecture with threshold parameters achieves an average relativeL2L\_\{2\}error of only 1\.17% with spiking 6\.49% resulting in a combined Error\-Energy score \(described in more detail in the next section\) of 7\.59 which is slightly lower than the Error\-Energy score of 8\.15 for the fast sigmoid VS\-GNO\. Based on these results, SAR\-GNO does seem to provide slightly improved performance over VS\-GNO, although the overall reconstruction performance remains relatively close between the two approaches\. More importantly, these results demonstrate that the reliance on ReLU within the SAR layer does not severely degrade performance compared to the GeLU activation used with VS\-GNO\. In contrast, the use of surrogate gradients in VS\-GNO can introduce additional instability and training degradation, as demonstrated by the performance degradation observed with the arctangent implementation\. Thus, SAR provides a competitive alternative that avoids surrogate\-gradient training while maintaining comparable or slightly improved reconstruction performance\.

We also examine the impact of the threshold parameter for the full and spectral\-only version of SAR\-GNO as well as the inclusion of the subsequent GeLU layer in Tables[4](https://arxiv.org/html/2608.23987#S3.T4)and[5](https://arxiv.org/html/2608.23987#S3.T5)\. For the spectral\-only model, removing the threshold for theγ=0\.003\\gamma=0\.003implementation significantly worsensL2L\_\{2\}error with similar spiking performance\. We can see that the removal of the threshold parameter allows the spectral spiking to increase forcing the model to collapse the input embedding spiking \(M\) to zero which hurts performance\. Adding a threshold allows the spectral SAR\-GNO to better regulate spiking for the spectral layer and avoid the collapse for the input embedding, resulting in an improved optimization and efficiency performance\. This is further emphasized with the performance of spectral\-only SAR\-GNO atγ=0\.001\\gamma=0\.001which displays a less dramatic difference between no threshold and a threshold but still exhibits slightly better spiking and accuracy when the threshold is included because it is able to significantly reduce the spectral layer percentage due to better spiking control\.

For both gamma values explored with the full graph architecture, adding the threshold slightly worsens relativeL2L\_\{2\}error compared to the non\-threshold version, but can also reduce spiking overall\. We see that adding threshold brings down the spectral layer spiking from 9\.64% to 4\.73% forγ=0\.01\\gamma=0\.01and 3\.21% to 1\.59% forγ=0\.05\\gamma=0\.05\. The spatial layer does observe less dramatic behavior and worse performance with the threshold, but these results still indicate how the optional threshold𝝉\\boldsymbol\{\\tau\}can be used to improve spiking control for certain layers and potentially reduce overall spiking\.

In addition, we explore ablation analysis with the GeLU activation function\. We remove the GeLU layer while keeping the threshold parameters for both spectral\-only and full architectures\. For the first model version, we find that adding the GeLU layer slightly weakens the L2 error performance \(around 0\.1% increase forγ=0\.003\\gamma=0\.003\) but the spiking percentage shows a significant decrease\. This result indicates that for the spectral\-only model the addition of the GeLU layer does allow for noticeable improvements in neuron communication and training, resulting in almost similar reconstruction with lower spiking\. For the full\-version, the performance of the GeLU layer demonstrates similar characteristics\. When adding the GeLU layer forγ=0\.05\\gamma=0\.05, the L2 error increases by around 0\.82% while the spiking decreases by around 1\.34%\. Forγ=0\.01\\gamma=0\.01, the L2 error between GeLU and no\-GeLU differs by only 0\.22% while the spiking with the GeLU layer was lower by 1\.72%, indicating that although the GeLU layer can introduce a small increase in reconstruction error, it consistently reduces the required spiking activity\. This suggests that the GeLU layer potentially improves the efficiency of neuron communication while maintaining comparable reconstruction performance, further supporting its use after the ReLU thresholding in the SAR layer depicted in Equation[3](https://arxiv.org/html/2608.23987#S2.E3)\. Despite these results, the conclusion is not definitive and further benchmarking and analysis is required\.

Overall, SAR\-GNO further reinforces the potential advantages of sparse activation regularization with ReLU over surrogate\-gradient\-based variable spiking formulations\. By directly regularizing the activation outputs, SAR\-GNO effectively suppresses unnecessary computation while preserving information flow throughout the network\. Although the SAR layer’s performance is not as strong compared to VS\-GNO as seen with the NOMAD implementation, our presented sparsity layer provides similar to slightly better reconstruction performance and efficiency, indicating the success of SAR within the context of non\-linear layers/architectures\. The inclusion of the threshold and subsequent GeLU does show noticeable improvement in spiking performance with only slightly worse or similarL2L\_\{2\}error, demonstrating their potential effectiveness when included in the SAR framework introduced in Equation[3](https://arxiv.org/html/2608.23987#S2.E3)\. Overall, the reliance of ReLU has not demonstrated severe performance degradation within the context of SAR\-GNO\.

Finally, it is worth reemphasizing that implementing 10 STS did not drastically change performance while worsening computational efficiency compared to its 1 STS counterpart\. Unlike the improvements observed for VS\-NOMAD, the introduction of additional temporal dynamics provided little benefit for the graph neural operator architecture, suggesting that temporal spiking dynamics are less effective within the VIRSO framework\.

![Refer to caption](https://arxiv.org/html/2608.23987v1/Figure7.jpg)Figure 7:2D Lid\-Driven Cavity\(a\)We define the boundary conditions for the Lid\-Driven Cavity geometry which shows zero velocity on all faces except for a temporal signal at the top surface\.\(b\)We display example time\-varying velocity profiles for the top surface which is defined for 90 time steps\.\(c\)The 50th percentile test example for the SAR\-NOMAD model with the Hoyer loss term andγ=0\.05\\gamma=0\.05\. We see that the main vortex and boundary behavior is closely captured as well as the corner eddies with mostly stochastic absolute error behavior\. We do see small hotspots of absolute error located at the top surface and corners indicating some difficulty in capturing surface physics with a regular discretized grid\.
### 3\.3Lid Driven Cavity \(NOMAD\)

To further demonstrate the improved convergence behavior of SAR\-based operators relative to surrogate gradient\-based VSN and LIF neurons, both as a potential neuromorphic alternative mechanism and as a benchmark for future development, we evaluated both approaches on the LDC dataset\. Consistent with the 2D Heat Exchanger study, Hoyer\-based SAR\-NOMAD models were trained using anα\\alphavalue of 1 while varying the sparsity coefficientγ\\gamma\. In parallel, VS\-NOMAD models employing both the fast sigmoid and arctangent surrogate gradients were trained across multiple values ofγ\\gammaand STS\. We also utilized LIF\-NOMAD implementations with the arctangent gradient function utilizing direct encoding with 1, 10, 20, and 30 STS and rate encoding with 10, 20, 30 STS\. A summary of the resulting performance is provided in Table[6](https://arxiv.org/html/2608.23987#S3.T6)\.

As expected, the baseline model without explicit spiking regularization already exhibits a degree of inherent sparsity across its sub\-networks \(almost 40% average spiking\)\. The VS\-NOMAD models reduce the average spiking activity to as low as 14\.78% using the arctangent surrogate gradient and 15\.02% using the fast sigmoid surrogate gradient, corresponding to a reduction of more than half of the baseline spiking activity\. However, this reduction is accompanied by a substantial degradation in reconstruction accuracy, with the lowest observed mean relativeL2L\_\{2\}error remaining above 13%\. The spiking behavior of VS\-NOMAD is similar to what is observed for the 2D Heat Exchanger dataset\. For the arctangent surrogate gradient, increasingγ\\gammagenerally produces lower spiking activity\. The fast sigmoid also shows this behavior except for the implementation with 10 STS which exhibits increased spiking atγ=10−5\\gamma=10^\{\-5\}compared to lower gamma\. Contrastingly, increasing STS does not result in the exact same behavior with the 2D Heat Exchanger\. From 20 to 30 STS for the Fast Sigmoid, the L2 error increases from 20\.44% to 31\.15% forγ=10−5\\gamma=10^\{\-5\}and 12\.86% to 14\.83% forγ=0\\gamma=0while increasing or keeping the same spiking percentage\. Fast Sigmoid also exhibits similar behavior at different spike time transitions\. Although the arctangent function is more aligned with our expectations, the VS\-NOMAD results indicate that the spiking dynamic behavior of VS\-NOMAD is not guaranteed to improve performance with higher STS potentially due to its high dependence on the surrogate gradient formulation\. In addition, we were unable to train an STS of 40 due to memory constraints with a single H200 further indicating the difficulties of training with high spike time neural operators for virtual sensing\.

Suprisingly, the LIF\-NOMAD results for direct encoding slightly outperformed the VS\-NOMAD implementation\. Similar to VSN, higher STS improvedL2L\_\{2\}performance while generally increasing spiking slighlty but the increase in gamma produced less predictable results with 1 STS and 20 STS seeing improvedL2L\_\{2\}error with higher gamma\. Either way, with 30 STS andγ=0\\gamma=0, the directly\-encoded LIF\-NOMAD had anL2L\_\{2\}error of 8\.76% with spiking around 10\.32% which is improved in both reconstruction and effeciency over all VSN implementations\. Of course the rate\-encoded, with far reduced precision in its input, presents sever error degradation similar to the Heat Exchanger\. Overall, the LIF\-NOMAD with direct\-encoding provides a more competitive output than the VS\-NOMAD unlike the VS/LIF\-GNO results in the previous section which might have benefited with the variable communication more due to the SNN\-ANN hybrid architecture design\.

Despite the previous results, the SAR layer with its variable communication outperforms the VSN and LIF implementations\. Similar to the Heat Exchanger application, SAR\-NOMAD exhibits a strong, consistent, and predictable sparsity\-accuracy tradeoff\. Asγ\\gammaincreases from 0\.001 to 0\.05, the average spiking activity decreases monotonically from approximately 19% to 4%, while the mean relativeL2L\_\{2\}error increases gradually from approximately 1% to 10%\. Although reconstruction error increases with stronger sparsity regularization, SAR\-NOMAD models withγ\>0\.005\\gamma\>0\.005outperform every evaluated VS\-NOMAD configuration in both reconstruction accuracy and spiking efficiency\. In comparsion with LIF\-NOMAD, the LIF neuron implemention with direct encoding does present slightly lower L2 error at 30 STS than theγ=0\.05\\gamma=0\.05SAR result but at the cost of higher spiking and significantly higher latency\. Overally, SAR\-NOMAD provides improved performance over all LIF results \(even outperforming the spiking percentages of the rate\-encoded model\) in terms of error, spiking, and latency\. Similar to the results obtained for the 2D Heat Exchanger, these findings further establish SAR\-based spiking as a low\-latency, high\-accuracy alternative to VSN\-based networks while simultaneously providing a practical gold\-standard benchmark for future spiking neural operator development due to its single spike\-step operation\. In Figure[7](https://arxiv.org/html/2608.23987#S3.F7)c, we see the 50th percentile performance of SAR\-NOMAD with the Hoyer loss term atγ=0\.05\\gamma=0\.05\. The main vortex and corner behaviors are captured accurately, further indicating that SAR is capable of providing accurate reconstruction with high energy\-efficiency\.

Because SAR\-NOMAD and VS/LIF\-NOMAD share the same underlying network architecture and differ only in their spiking formulation, their performance can be directly compared using a simple Latency\-Error\-Energy \(LEE\) metric similar to the energy\-delay product in concept but including reconstruction performance: an important pillar of efficient virtual sensing\. The proposed LEE score is defined as the product of the latency measured in spike time steps or STS, the mean relativeL2L\_\{2\}reconstruction error \(%\), and the average spiking percentage across all sub\-networks \(%\), with lower values indicating superior overall performance under the assumption that latency, accuracy, and energy efficiency are equally weighted\. To further isolate the tradeoffs between these quantities, three additional metrics are considered: the Latency\-Error \(LErr\), Error\-Energy \(EE\), and Latency\-Energy \(LEn\) scores\. While latency and energy consumption could ultimately be expressed in physical units following hardware implementation, the identical network architecture and implementation of both SAR\-NOMAD and VS/LIF\-NOMAD allow STS and average spiking percentage to serve as appropriate proxies for latency and energy consumption in the present software\-based evaluation\. Moreover, although equal weighting is assumed in these metrics, application\-specific requirements may prioritize one objective over another\. For example, safety\-critical applications may favor lower reconstruction error, whereas energy\-constrained deployments may prioritize reduced energy consumption\.

Using these metrics, Tables[7](https://arxiv.org/html/2608.23987#S3.T7)and[8](https://arxiv.org/html/2608.23987#S3.T8)clearly demonstrate that the selected SAR\-NOMAD configuration outperforms every evaluated VS/LIF\-NOMAD model across all four performance metrics\. For the comprehensive LEE score, the selected SAR\-NOMAD model achieves a score of 26\.35 for the 2D Heat Exchanger, representing at least a fivefold improvement over the best\-performing VS/LIF\-NOMAD configuration\. Similarly, for the Lid Driven Cavity dataset, SAR\-NOMAD achieves an LEE score of 37\.96, corresponding to an improvement of also more than fivefold compared to the best VS/LIF\-NOMAD model\. SAR\-NOMAD also consistently achieves superior performance across the individual LErr, EE, and LEn metrics, demonstrating improvements not only in reconstruction accuracy but also in computational efficiency while maintaining the low\-latency advantages of single\-step spiking\. Given that real\-time prediction, and ultimately faster\-than\-real\-time prediction, is essential for providing timely information to system operators, latency represents a critical performance metric\. Consequently, the combination of low latency, high reconstruction accuracy, and improved energy efficiency makes SAR\-NOMAD and the SAR layer in general a highly effective candidate for efficient virtual sensing and field reconstruction\.

Table 6:L2L\_\{2\}errors and spiking percentages for Hoyer\-based SAR\-NOMAD and VS/LIF\-NOMAD \(Fast Sigmoid and Arctanget surrogate gradients\) with the Lid Driven Cavity\.Mean Relative L2 Errors \(%\)Mean Spiking Percentages \(%\)Spiking ModelppvvkkMeanTrunkBranchCombMeanBase Model0\.891\.691\.051\.2145\.3531\.4137\.9238\.23SAR Hoyerγ=0\.001\\gamma=0\.0011\.031\.931\.101\.3526\.217\.2424\.9519\.47SAR Hoyerγ=0\.005\\gamma=0\.0051\.693\.401\.812\.3016\.593\.8316\.7912\.40SAR Hoyerγ=0\.01\\gamma=0\.012\.885\.362\.883\.7113\.032\.7412\.289\.35SAR Hoyerγ=0\.05\\gamma=0\.055\.1916\.108\.159\.815\.900\.744\.983\.87VS Fast Sig\. 1 STSγ=0\\gamma=06\.6514\.0222\.5614\.4131\.4535\.7212\.6326\.60VS Fast Sig\. 1 STSγ=10−5\\gamma=10^\{\-5\}30\.0245\.1642\.8239\.3319\.5514\.982\.6712\.40VS Fast Sig\. 10 STSγ=0\\gamma=012\.2520\.6625\.7919\.5728\.9310\.229\.2816\.14VS Fast Sig\. 10 STSγ=10−5\\gamma=10^\{\-5\}12\.2727\.8423\.9121\.3430\.0819\.547\.1218\.91VS Fast Sig\. 20 STSγ=0\\gamma=06\.3615\.0517\.1712\.8633\.2016\.3210\.5820\.03VS Fast Sig\. 20 STSγ=10−5\\gamma=10^\{\-5\}13\.2123\.1324\.9920\.4430\.4222\.587\.9720\.32VS Fast Sig\. 30 STSγ=0\\gamma=07\.9516\.5719\.9714\.8342\.4027\.3317\.2629\.00VS Fast Sig\. 30 STSγ=10−5\\gamma=10^\{\-5\}15\.5437\.2840\.6431\.1529\.9212\.149\.2917\.12VS Arctan\. 1 STSγ=0\\gamma=010\.0511\.5221\.2814\.2833\.8918\.087\.9619\.98VS Arctan\. 1 STSγ=10−5\\gamma=10^\{\-5\}19\.9233\.4745\.7533\.0518\.6323\.242\.4814\.78VS Arctan\. 10 STSγ=0\\gamma=05\.9217\.2923\.6415\.6235\.5731\.6014\.8527\.34VS Arctan\. 10 STSγ=10−5\\gamma=10^\{\-5\}10\.0132\.4627\.7123\.3928\.4916\.049\.0217\.85VS Arctan\. 20 STSγ=0\\gamma=06\.5613\.0221\.8113\.8035\.1926\.8611\.2424\.43VS Arctan\. 20 STSγ=10−5\\gamma=10^\{\-5\}11\.2821\.5924\.7419\.2129\.3518\.298\.9518\.86VS Arctan\. 30 STSγ=0\\gamma=04\.3213\.5022\.5213\.4535\.5327\.4612\.8125\.27VS Arctan\. 30 STSγ=10−5\\gamma=10^\{\-5\}7\.5315\.9323\.2915\.5829\.4120\.908\.4719\.59LIF Arctan\. 1 STS DEγ=0\\gamma=027\.2036\.9333\.1532\.423\.904\.8210\.576\.43LIF Arctan\. 1 STS DEγ=10−5\\gamma=10^\{\-5\}26\.2435\.1429\.4030\.263\.898\.0010\.297\.39LIF Arctan\. 10 STS DEγ=0\\gamma=09\.0013\.3810\.8911\.098\.778\.5511\.269\.53LIF Arctan\. 10 STS DEγ=10−5\\gamma=10^\{\-5\}9\.1512\.4612\.2811\.308\.8110\.9011\.4610\.39LIF Arctan\. 20 STS DEγ=0\\gamma=07\.9611\.3910\.8910\.089\.898\.8511\.3310\.02LIF Arctan\. 20 STS DEγ=10−5\\gamma=10^\{\-5\}7\.4211\.6910\.089\.739\.599\.1610\.769\.84LIF Arctan\. 30 STS DEγ=0\\gamma=08\.249\.808\.258\.769\.849\.4111\.7110\.32LIF Arctan\. 30 STS DEγ=10−5\\gamma=10^\{\-5\}8\.0510\.469\.029\.189\.899\.3411\.2310\.15LIF Arctan\. 10 STS REγ=0\\gamma=030\.8356\.3263\.9750\.375\.125\.234\.935\.09LIF Arctan\. 10 STS REγ=10−5\\gamma=10^\{\-5\}29\.9955\.8663\.7049\.855\.034\.434\.744\.73LIF Arctan\. 20 STS REγ=0\\gamma=026\.2054\.1862\.0847\.497\.796\.124\.276\.06LIF Arctan\. 20 STS REγ=10−5\\gamma=10^\{\-5\}23\.1050\.2660\.7944\.717\.015\.984\.535\.84LIF Arctan\. 30 STS REγ=0\\gamma=021\.1047\.9057\.1042\.047\.896\.024\.516\.14LIF Arctan\. 30 STS REγ=10−5\\gamma=10^\{\-5\}20\.7648\.1057\.8242\.237\.685\.414\.515\.87

Table 7:LEE, LErr, EE, and LEn scores for the Hoyer\-Based SAR\-NOMAD withγ=0\.005\\gamma=0\.005and various VS/LIF\-NOMAD results for the 2D Heat Exchanger\.Table 8:LEE, LErr, EE, and LEn scores for the Hoyer\-Based SAR\-NOMAD withγ=0\.005\\gamma=0\.005and various VS\-NOMAD results for the Lid Driven Cavity\.
### 3\.4Synthetic Distillation Performance

As discussed in Section[2\.4](https://arxiv.org/html/2608.23987#S2.SS4), within virtual sensing, and particularly spatial\-temporal field reconstruction, the availability of training data is frequently limited for several reasons\. For experimentally acquired datasets, measurements may be restricted due to security considerations or may be prohibitively expensive to obtain\. Furthermore, when the desired reconstruction cannot be directly measured, high\-fidelity computational methods, such as the Finite Element Method \(FEM\), are required to generate training data, resulting in substantial computational cost\. Consequently, neural operators for virtual sensing are often trained using relatively small datasets\.

Advanced neural operator architectures, such as VIRSO and Geo\-FNO, can effectively approximate nonlinear operator mappings under limited data conditions due to their strong inductive biases, sophisticated spectral representations, and ability to operate on irregular geometries\. However, these architectures require convolution over the entire computational grid, making them less suitable for deployment on edge devices and neuromorphic hardware\. In contrast, NOMAD and other trunk\-branch neural operator architectures consist primarily of fully connected network \(FCN\) layers and independently evaluate each spatial location\. While this structure is well suited for neuromorphic implementation, it lacks the enriched feature extraction and communication dynamics present in convolution\-based architectures\.

This behavior is reflected in the performance difference between VS\-GNO and VS\-NOMAD on the 2D Heat Exchanger dataset, despite both models employing the same Variable Spiking Neuron \(VSN\)\. Atγ=0\.5\\gamma=0\.5, VS\-GNO achieves both a lower mean relativeL2L\_\{2\}error and lower average spiking activity than the majority of VS\-NOMAD configurations\. These observations bring us towards synthetic distillation, as discussed in Section[2\.4](https://arxiv.org/html/2608.23987#S2.SS4)\. The large graph\-based VIRSO model provides new data examples, based on the distributions in Equation[17](https://arxiv.org/html/2608.23987#S2.E17)and the model’s outputs, for SAR\-NOMAD which results in better optimization and convergence while ideally preserving the energy\-efficiency of the model architecture\.

To establish the feasibility of this framework, the synthetic input\-output pairs from VIRSO are then combined with the original training dataset and used to train a SAR\-NOMAD model with the Hoyer sparsity loss andγ=0\.001\\gamma=0\.001,0\.0050\.005, and0\.010\.01\. Tables[9](https://arxiv.org/html/2608.23987#S3.T9),[10](https://arxiv.org/html/2608.23987#S3.T10), and[11](https://arxiv.org/html/2608.23987#S3.T11)summarize the resulting performance using the same test dataset employed throughout the previous experiments\.

For all gamma examples, as we add more synthetic examples, the mean relativeL2L\_\{2\}error decreases with either similar or improved spiking efficiency\. With the addition of 8,000 synthetic training examples, the mean relativeL2L\_\{2\}error is reduced by approximately a factor of two for bothγ=0\.001\\gamma=0\.001and0\.0050\.005, while the average spiking activity decreases by approximately 1% forγ=0\.001\\gamma=0\.001and0\.6%0\.6\\%forγ=0\.005\\gamma=0\.005\. These results indicate that additional knowledge generated by the teacher model enables SAR\-NOMAD to maintain its energy\-efficient sparse representations while learning significantly more accurate field reconstructions\. Consequently, the overall Latency\-Error\-Energy \(LEE\) score for theγ=0\.005\\gamma=0\.005result is reduced by more than a factor of two to 11\.34\.

Forγ=0\.01\\gamma=0\.01in Table[11](https://arxiv.org/html/2608.23987#S3.T11), the baseline model trained without synthetic data exhibits substantial performance degradation, with the relativeL2L\_\{2\}error increasing to nearly19%19\\%while reducing the average spiking activity to below2%2\\%\. Introducing 1,000 synthetic training samples dramatically improves reconstruction accuracy, reducing theL2L\_\{2\}error to below5%5\\%while increasing the average spiking activity by only1\.3%1\.3\\%\. As additional synthetic samples are incorporated, the model continues to follow the same accuracy trend, with 8,000 synthetic samples further reducing the meanL2L\_\{2\}error by approximately a factor of two relative to the 1,000\-sample case while requiring only a slight increase in spiking activity\. Although the 8,000\-sample model exhibits higher average spiking than the originalγ=0\.01\\gamma=0\.01baseline, the substantial improvement in reconstruction accuracy demonstrates the effectiveness of our initial synthetic distillation approach for SAR\-NOMAD\. In contrast, applying synthetic distillation to VS\-NOMAD using the 10 STS,γ=10−6\\gamma=10^\{\-6\}, arctangent configuration yielded minimal improvements in eitherL2L\_\{2\}error or average spiking percentage\. These results suggest that extending synthetic distillation to variable\-spiking architectures is considerably more challenging, potentially due to a mismatch between the surrogate gradient used during optimization and the behavior of the synthesized training data\.

However, generating 8,000 synthetic samples increases the effective training dataset to nearly nine times its original size\. Although this straightforward synthetic data generation strategy substantially reduces the dependence on expensive experimental measurements or computationally intensive high\-fidelity simulations, it still introduces a considerable computational burden during training because a large quantity of synthetic data is required when random input sampling is employed\. Future work should therefore investigate more sophisticated knowledge distillation strategies that transfer information more efficiently while requiring fewer synthetic examples\. Nevertheless, the present results demonstrate the viability of a distillation framework specifically tailored for spiking neural operators, wherein a teacher model that is not inherently suitable for neuromorphic deployment learns from limited ground\-truth data and subsequently generates synthetic training data that enables a neuromorphic\-friendly student model to achieve substantially better performance than training on the original dataset alone\.

Table 9:L2L\_\{2\}errors and spiking percentages for SAR\-NOMAD and Hoyer loss term atγ=0\.001\\gamma=0\.001for different synthetically generated training example counts with the 2D Heat Exchanger\.Mean Relative L2 Errors \(%\)Mean Spiking Percentages \(%\)Spiking Modelppvzv\_\{z\}vyv\_\{y\}vxv\_\{x\}vmagv\_\{\\mathrm\{mag\}\}MeanTrunkB1B2CombMeanSAR Hoyer \(Originalγ=0\.001\\gamma=0\.001\)2\.144\.423\.552\.392\.963\.0917\.811\.862\.0714\.098\.96SAR Hoyer \(Samples = 1000\)1\.834\.082\.962\.482\.842\.8417\.551\.952\.2213\.458\.79SAR Hoyer \(Samples = 2000\)1\.582\.782\.412\.122\.652\.3117\.281\.561\.7113\.368\.48SAR Hoyer \(Samples = 4000\)1\.432\.451\.901\.462\.141\.8716\.920\.971\.9113\.138\.23SAR Hoyer \(Samples = 8000\)1\.171\.701\.561\.241\.891\.5115\.911\.071\.2312\.687\.72

Table 10:L2L\_\{2\}errors and spiking percentages for SAR\-NOMAD and Hoyer loss term atγ=0\.005\\gamma=0\.005for different synthetically generated training example counts with the 2D Heat Exchanger\.Mean Relative L2 Errors \(%\)Mean Spiking Percentages \(%\)Spiking Modelppvzv\_\{z\}vyv\_\{y\}vxv\_\{x\}vmagv\_\{\\mathrm\{mag\}\}MeanTrunkB1B2CombMeanSAR Hoyer \(Originalγ=0\.005\\gamma=0\.005\)4\.116\.575\.824\.845\.725\.418\.910\.600\.839\.154\.87SAR Hoyer \(Samples = 1000\)3\.115\.634\.843\.514\.424\.308\.550\.390\.758\.984\.67SAR Hoyer \(Samples = 2000\)2\.924\.343\.762\.903\.563\.508\.770\.490\.718\.674\.66SAR Hoyer \(Samples = 4000\)2\.613\.753\.392\.703\.403\.178\.230\.390\.688\.884\.55SAR Hoyer \(Samples = 8000\)2\.203\.062\.922\.193\.022\.687\.970\.390\.587\.994\.23

Table 11:L2L\_\{2\}errors and spiking percentages for SAR\-NOMAD and Hoyer loss term atγ=0\.01\\gamma=0\.01for different synthetically generated training example counts with the 2D Heat Exchanger\.
### 3\.5Sparse Activation ReLU towards VS\-GNO

#### Improving The VSN Loss Term

In order to preserve the dynamic spiking behavior of the VSN, enabling richer feature representations and improved modeling of transient physics, we investigate an alternative regularization strategy rather than applying SAR\-based neural operators directly\. Specifically, we introduce an activation\-based regularization loss inspired by SAR, described in Section[2\.5](https://arxiv.org/html/2608.23987#S2.SS5)and Equation[8](https://arxiv.org/html/2608.23987#S2.E8)\. Unlike the original percentage\-based spiking loss, this formulation employs the Hoyer loss to regularize the activation strength of the spiking neurons, providing a more informative measure of spiking behavior\. To evaluate the effectiveness of this approach, we compare the original VS\-NOMAD with Activation Regularized VS\-NOMAD \(ARVS\-NOMAD\) on both the 2D Heat Exchanger and Lid\-Driven Cavity datasets using 1 and 10 STS configurations and an arctangent surrogate gradient across a range ofγ\\gammavalues, as summarized in Tables[12](https://arxiv.org/html/2608.23987#S3.T12)and[13](https://arxiv.org/html/2608.23987#S3.T13)\.

For the 1 STS experiments on the 2D Heat Exchanger, ARVS\-NOMAD consistently achieves substantially lower mean relativeL2L\_\{2\}error than the original VS\-NOMAD\. The highest\-performing ARVS\-NOMAD configuration reduces the meanL2L\_\{2\}error to below9%9\\%, whereas all VS\-NOMAD configurations remain above30%30\\%\. Furthermore, the spiking activity of ARVS\-NOMAD is comparable to or lower than that of VS\-NOMAD, with theγ=5⋅10−7\\gamma=5\\cdot 10^\{\-7\}configuration achieving both lower spiking and lower reconstruction error than every VS\-NOMAD result\. Similar improvements are observed for the 10 STS experiments, where ARVS\-NOMAD consistently produces lower meanL2L\_\{2\}errors than the original VS\-NOMAD\. Although the average spiking activity is generally slightly higher, a direct comparison of identicalγ\\gammavalues is not entirely appropriate because the percentage\-based and activation\-based regularization terms differ in scale and magnitude\. Nevertheless, ARVS\-NOMAD is able to reduce the average spiking activity to below6\.52%6\.52\\%while achieving anL2L\_\{2\}error significantly lower than the16\.05%16\.05\\%obtained by the lowest\-spiking VS\-NOMAD configuration\. Interestingly, the reconstruction error temporarily improves asγ\\gammaincreases from5⋅10−85\\cdot 10^\{\-8\}to10−610^\{\-6\}, suggesting that a moderate amount of activation regularization improves optimization by encouraging more efficient utilization of the learned feature representations\.

The Lid\-Driven Cavity results exhibit similar trends\. Across both the 1 and 10 STS configurations, ARVS\-NOMAD consistently achieves mean relativeL2L\_\{2\}errors of approximately11%11\\%, whereas the corresponding VS\-NOMAD models remain above23%23\\%\. For the 1 STS configuration, ARVS\-NOMAD also reduces the average spiking activity, providing improvements in both reconstruction accuracy and energy efficiency\. For the 10 STS configuration withγ=10−4\\gamma=10^\{\-4\}, ARVS\-NOMAD achieves a substantial reduction inL2L\_\{2\}error compared to VS\-NOMAD while maintaining nearly identical spiking activity\.

The results on both the Heat Exchanger and Lid\-Driven Cavity datasets demonstrate the improved learning dynamics enabled by the proposed activation\-based regularization\. By leveraging the superior sparsity\-promoting properties of the Hoyer loss, as demonstrated in our comparison between theL1L\_\{1\}and Hoyer regularizers, ARVS\-NOMAD regularizes the magnitude of neuron activations rather than treating spikes as purely binary events, as in the original percentage\-based loss\. Specifically, the proposed regularization operates on the difference between the pre\-reset membrane potential and the firing threshold, thereby encoding richer information about spike intensity\. This formulation preferentially penalizes neurons with consistently large membrane potentials/low thresholds that are likely to spike continuously, encouraging more adaptive threshold behavior and improved feature utilization\. Consequently, this targeted regularization produces both lower reconstruction errors and reduced spiking activity\.

Despite these improvements, ARVS\-NOMAD still underperforms the SAR\-NOMAD architectures because the underlying optimization remains dependent on surrogate gradients, resulting in an inherent gradient mismatch\. Consequently, these results primarily demonstrate the effectiveness of combining the SAR\-inspired activation encoding with Hoyer regularization for developing energy\-efficient neural operators, rather than resolving the optimization challenges associated with variable\-spiking neural operators\. Future work should therefore investigate alternative training algorithms that eliminate or mitigate the reliance on surrogate gradients\.

Table 12:L2L\_\{2\}errors and spiking percentages for Activation\-Regularized VS\-NOMAD and Percentage\-Regularized VS\-NOMAD with the 2D Heat Exchanger trained using the Arctangent surrogate gradient function\.Table 13:L2L\_\{2\}errors and spiking percentages for Activation\-Regularized VS\-NOMAD and Percentage\-Regularized VS\-NOMAD with the Lid Driven Cavity trained using the Arctangent surrogate gradient function\.
#### Reducing Spatial Aggregation

Table 14:Average error and spiking with the layer\-wise neighbor count and average count for the VS\-GNO model on the Heat Exchanger with the arctangent surrogate function utilizing both spatial thresholding gating mechanisms: shared threshold \(Exp\. 1\) and edge\-dependent threshold \(Exp\.2\)\.To demonstrate the potential of ReLU thresholding to reduce spatial aggregation costs and, consequently, the latency of VS\-GNO implementations, Table[14](https://arxiv.org/html/2608.23987#S3.T14)presents the full VS\-GNO results using the arctangent surrogate gradient withγ=0\\gamma=0andγ=0\.5\\gamma=0\.5for the two proposed thresholding techniques: a single trainable threshold value \(Exp\. 1\) and an edge\-dependent trainable threshold value \(Exp\. 2\)\. We utilize twoϵ\\epsilonvalues \(0\.01 and 0\.001\) to demonstrate different levels of emphasis on retaining neighbor connections\. For comparison, we also report the full VS\-GNO with a reduced number of nearest neighbors within the KNN graph construction utilized for these results\. Forγ=0\\gamma=0, all thresholding experiments outperform the reduced\-neighbor configurations in terms ofL2L\_\{2\}error while maintaining slightly higher or similar spiking percentages and similar or reduced neighbor counts\. Forϵ=0\.01\\epsilon=0\.01andγ=0\\gamma=0, with the constant trainable threshold \(Exp\. 1\), we achieve an almost sevenfold lowerL2L\_\{2\}error than the best\-performing reduced\-neighbor configuration based on KNN construction, while also obtaining lower spiking and an average neighbor count of less than 5\. Examining the layer\-wise neighbor statistics shows that the thresholding technique effectively removes neighbors in specific layers while maintaining aggregation in others, essentially finding a favorable tradeoff between reconstruction error and the overall need for spatial aggregation\. With the edge\-dependent thresholding atγ=0\\gamma=0, we observe cases where the neighbor count is reduced to zero, effectively reducing the model to a spectral\-only version\. Despite this, the resulting models maintain errors below 1% with similar spiking percentages\. Forγ=0\.5\\gamma=0\.5, the edge\-dependent thresholding \(Exp\. 2\) demonstrates the better overall performance, providing greater control over which neighbors are blocked based on their spatial positioning and edge weights\. In particular, theϵ=0\.001\\epsilon=0\.001result forγ=0\.5\\gamma=0\.5in Exp\. 2 achieves an averageL2L\_\{2\}error of approximately 1%, a spiking percentage below 9%, and an average neighbor count lower than all reduced\-neighbor alternatives forγ=0\.5\\gamma=0\.5\. These results suggest that trainable edge thresholding is preferable over simply reducing neighbor count which not only affects the amount of spatial aggregation but also limits the model’s ability to selectively determine which neighbors should be removed, while simultaneously affecting the spectral component\. Experimenting with post\-training thresholding, which essentially entailed finetuning of a constant threshold, resulted in poor performance which makes the proposed trainable ReLU thresholding approach a more favorable alternative for maintaining lowL2L\_\{2\}error while reducing spatial aggregation\. We also experimented the thresholding concept with the full SAR\-GNO model but the results demonstrated more error degradation\. Since the SAR implementations were shown, in Table[5](https://arxiv.org/html/2608.23987#S3.T5), to have far less spiking \(at least 10 times less compared to the arctangent VS\-GNO results\) for the spatial blocks, the inclusion of the neighbor ReLU\-based loss term could have thrown off ideal convergence for the SAR implementation especially since it has been shown with SAR\-NOMAD that lower spiking leads to collapse in the feature dimension\. This, of course, is purely speculative and further analysis into the spiking behavior of SAR\-GNO and the impact of the neighbor thresholding is required\. Nevertheless, if spatial aggregation is required for local spatial calibration, the proposed method provides a potential mechanism that does require further exploration to determine its true efficacy\. At least in the context of the VS\-GNO algorithm, which does rely more on the spatial aggregation in terms of spiking amount, our proposed neighbor "filter" has promising performance towards optimizing the tradeoff between accuracy and the overall number of neighbor connections\. Further research is also required to fully understand the impact of the proposed thresholding techniques and the differences between the two thresholding approaches\. Also, further applications are required to determine whether spatial aggregation is necessary at all, given that the spectral\-only layer achieves betterL2L\_\{2\}error and that lower neighbor counts correspond to Exp\. 1 withγ=0\\gamma=0and Exp\. 2 withγ=0\.5\\gamma=0\.5\.

## 4Conclusions and Further Work

In this paper, we addressed several key limitations of previously presented energy\-efficient solutions, mainly spiking neural operators and the Variable Spiking Neuron\. Owing to surrogate gradient mismatch, the reconstruction performance of single\-spike\-timestep \(1 STS\) model tends to degrade in specific circumstances, often requiring additional spike timesteps to recover accuracy\. Increasing the number of spike timesteps, however, introduces undesirable latency and significantly increases the computational burden during training\. We demonstrate this limitation through memory constraints encountered on a single H200 GPU when training 30 STS and 40 STS models for the 2D Heat Exchanger and Lid\-Driven Cavity use cases\. To overcome these challenges, we move away from the traditional spiking\-based models proposed for energy efficiency and propose an activation\-sparsity framework\. We present a low\-latency, single\-step framework based on a Sparse\-Activation\-ReLU \(SAR\) layer that avoids any surrogate gradient mismatch entirely and allows for variable communication similar to the VSN while removing the temporal integration that can risk increased latency, operating in a one\-step forward pass\. The SAR layer provides a alternative solution to the typical spiking implementations\. Although other hardware implementations that take advantage of activation sparsity, our energy\-efficient models can fit within an ANN\-to\-neuromorphic framework, utilizing a zero\-thresholding event policy, that is unlike traditional ANN\-SNN algorithms which rely on rate coding and binary communication and risk precision loss\. Although SAR is somewhat more restrictive than standard VSNs, as it only permits positive activations beyond the zero threshold whereas VSNs can utilize both positive and negative thresholds and signals, it consistently outperformed or performed similarly to VSN\-based and LIF\-based implementations for both the trunk\-branch NOMAD architecture and the graph\-based VIRSO architecture \(which utilized GeLU layers\) on the presented benchmark problems\.

To quantify the overall tradeoff between reconstruction accuracy, latency, and energy efficiency, we introduced the Latency\-Energy\-Error \(LEE\) score, which equally weights latency \(number of spike time steps or STS\), reconstruction error \(mean relativeL2L\_\{2\}error \(%\)\), and energy consumption \(average spiking \(%\)\)\. Using this metric, we found that SAR\-NOMAD improved upon VS/LIF\-NOMAD by at least a factor of five on the 2D Heat Exchanger and the Lid\-Driven Cavity benchmark\.

We also introduced analysis based on spiking entropy to quantify how much of a model’s feature representation is utilized across the test dataset\. Evaluating the Heat Exchanger benchmark, we observed that both SAR\-NOMAD and VS\-NOMAD exhibit decreasing feature entropy as the target spiking percentage is reduced\. For sufficiently small sparsity coefficients, such asγ=0\.005\\gamma=0\.005for SAR\-NOMAD, the entropy approaches zero, indicating that nearly all information is concentrated within a single feature dimension\. These results suggest that, in order to satisfy increasingly restrictive sparsity constraints, both SAR and VSN models progressively collapse their latent feature representations to achieve lower activation or spiking rates\.

Building upon these findings, we proposed a synthetic distillation framework for neuromorphic/edge device virtual sensing\. Synthetic input\-output samples were generated using the graph\-based VIRSO model and combined with the original ground\-truth dataset to train SAR\-NOMAD on the 2D Heat Exchanger benchmark\. This approach reduced the mean relativeL2L\_\{2\}error by approximately a factor of two using 8,000 synthetic training samples while maintaining nearly identical or lower average spiking percentages\. Consequently, the corresponding LEE score was reduced by approximately a factor of two, further improving SAR\-NOMAD’s overall performance and increasing its advantage over VS\-NOMAD to more than an order of magnitude on the 2D Heat Exchanger benchmark\.

Finally, we investigated methods to improve the performance of the VSN implementation which provides potentially beneficially temporal integration, especially in the context of transient dynamics, as well as improved signaling \(no ReLU reliance\)\. Rather than replacing the neuron itself, we substituted the original spike\-percentage regularization term with a SAR\-inspired sparse activation loss\. Although these results were not as pronounced as those achieved with SAR\-NOMAD, the modified VS\-NOMAD models consistently achieved lower reconstruction errors while maintaining comparable or lower spiking percentages on both the 2D Heat Exchanger and Lid\-Driven Cavity benchmarks, especially in the 1 STS domain\. For the Heat Exchanger with 1 STS andγ=5⋅10−7\\gamma=5\\cdot 10^\{\-7\}, the ReLU activation loss implementation provided anL2L\_\{2\}error at least 2 times lower than the presented VSN results with a fraction\-based spiking loss term while also having an overall lower spiking percentage\. These findings suggest that sparse activation regularization provides richer information encoding than the original binary spike\-percentage objective while preserving the temporal processing capabilities unique to spiking neural networks\. We also explored ReLU\-based thresholding within the context of gated neighborhood aggregation of the VS\-GNO algorithm for the Heat Exchanger, which showed stronger dependence on the spatial blocks, indicated by higher spiking percentages\. Trainable thresholds, either shared or edge\-dependent, were able to reduce overall neighbor count with mostly improved L2 error and spiking over lower neighbor constructed graphs or manual post\-training thresholding\. Forγ=0\\gamma=0, the single threshold \(Exp\. 1\) result withϵ=0\.01\\epsilon=0\.01provided almost 7 times lower L2 error than the lower k\-value results for the same gamma with also lower spiking and overall average neighbor count\. Although the true efficacy of the spatial layer is in doubt, we present a methodology to mitigate the neighbor connectivity while still providing accurate reconstruction and low spiking, overall improving the integration of VS\-GNO into neuromorphic/edge device hardware and its overall latency\. Further research should delve into transient or generally varying graph structures instead of the single, stationary graph topology which might represent only a subset of virtual sensing applications\.

As discussed throughout this work, SAR\-based neural operators are not intended to replace conventional spiking neural operators, but rather to serve as a complementary low\-latency sparsity\-based alternative and a strong performance benchmark while improvements in spiking dynamics and optimization continue to emerge\. Initial research should explore modifications to SAR that allow the expressiveness seen with the VSN\. Currently, the ReLU function filters out negative inputs, and although a subsequent activation functionσ\\sigmaexists that can be important for sequential/transformer models, using tanh/softmax, the layer has less expression than the VSN which permits negative signals\. Despite this limitation, SAR still performs similar or slightly outperforms VSN with the GeLU activation in VS\-GNO, indicating that gradient mismatch poses a larger disadvantage than layer expressiveness\. Moving beyond SAR, future research should also focus on overcoming the limitations of surrogate\-gradient training and developing optimization techniques that preserve the inherent advantages of spiking computation\. Ideally, neuron membrane integration, temporal memory, and dynamic spike generation should provide richer feature representations and more effective processing of transient physical phenomena\. Spiking neural networks have consistently demonstrated strong performance on time\-series problems, making it natural to investigate whether aligning neuron memory dynamics with the temporal behavior of transient physics can further improve virtual sensing performance\. In particular, an important research question is whether membrane potentials should be reset after each physical timestep or propagated continuously throughout a transient simulation to better capture long\-term temporal dependencies\.

ANN\-to\-SNN conversion remains one of the primary surrogate\-free training methodologies and has demonstrated considerable success for LIF\-based neural networks\. However, existing conversion techniques are fundamentally based on rate\-coded communication and therefore do not naturally extend to Variable Spiking Neurons \(VSNs\), whose graded spike outputs encode continuous\-valued information rather than firing rates and binary communication\. The proposed SAR neural operator framework provides a alternative solution that is based on sparsity events not standard spiking dynamics but can still enable direct ANN\-to\-neuromorphic conversion for single\-spike\-timestep using zero thresholding while avoiding surrogate\-gradient optimization and providing variable communication vital for performance on regression based spatial\-temporal reconstruction\. Nevertheless, this simplification sacrifices the temporal memory dynamics that make VSNs particularly attractive for multi\-timestep inference\. Consequently, extending surrogate\-free optimization techniques such as SAR\-NOMAD to incorporate membrane memory and temporal spiking dynamics \(and exist not as a simple sparsity filter\) represents a promising direction for future investigation\. More broadly, this challenge motivates exploration of fundamentally different optimization paradigms that eliminate gradient mismatch altogether, including gradient\-free optimization techniques, reinforcement learning, and emerging quantum computing approaches for large\-scale search\-based training\[[38](https://arxiv.org/html/2608.23987#bib.bib38),[39](https://arxiv.org/html/2608.23987#bib.bib39)\]\.

Beyond training methodology, our entropy\-based analysis of sparse neural representations and VSN implementations suggests additional opportunities for improving overall event\-based model efficiency and expressiveness\. In particular, the branch layers of the Hoyer\-regularized SAR\-NOMAD model exhibited entropy values approaching zero in several intermediate layers, indicating that the learned representations collapsed onto only one or a few dominant feature dimensions\. While this behavior promotes sparsity, it also suggests that the available feature space is not being fully utilized\. These observations motivate two complementary research directions\. The first is model pruning\[[40](https://arxiv.org/html/2608.23987#bib.bib40)\], where inactive neurons, unused synapses connections, or even entire model components can be removed to reduce both parameter count and memory requirements\. Such compression is especially important for deployment on highly resource\-constrained edge platforms, such as neuromorphic devices \(Loihi 2\[[27](https://arxiv.org/html/2608.23987#bib.bib27)\]\), where memory capacity often represents a more significant limitation than computational throughput\. The second direction involves developing new regularization strategies that explicitly consider feature entropy in addition to activation sparsity\. Encouraging a more balanced utilization of the latent feature space while maintaining equivalent spiking activity could improve reconstruction accuracy by preventing feature collapse and increasing the representational capacity of intermediate neural embeddings\.

In parallel with algorithmic advances, continued work toward realistic hardware deployment remains essential\. Implementing SAR\-based neural operators on edge device hardware that exploits the sparsity, which we have discussed could include neuromorphic chips, will not only demonstrate and answer the question of the practical feasibility of neural operators for energy\-efficient virtual sensing but also quantify the latency, throughput, and energy\-efficiency improvements predicted by the proposed framework\. Hardware deployment will additionally highlight the importance of the synthetic knowledge distillation strategy introduced in this work\. As neural operators continue to increase in architectural complexity, direct deployment onto current edge device hardware becomes increasingly difficult due to memory, connectivity, and architectural constraints\. Knowledge transfer from large, expressive teacher models to compact trunk\-branch student architectures therefore represents a practical pathway for translating state\-of\-the\-art neural operators into deployable edge systems\.

Further development of the proposed synthetic distillation framework is also warranted\. Although synthetic training data substantially improved reconstruction performance on the limited Heat Exchanger dataset, particularly for SAR\-NOMAD, these gains were accompanied by a significant increase in training cost, especially when utilizing 8,000 synthetic samples\. Future work should therefore investigate more efficient synthetic data generation strategies capable of achieving comparable performance improvements with significantly fewer generated examples\. Likewise, performing knowledge distillation directly on the original training dataset rather than relying exclusively on synthetic samples may further reduce computational overhead while preserving the benefits of teacher\-student learning\. Finally, the current framework assumes that the input parameter distributions are known a priori, enabling straightforward sampling of synthetic training examples\. Extending this methodology to more general functional inputs with unknown or highly complex distributions remains an important open problem, particularly for experimental measurements and real\-world datasets where governing parameterizations are unavailable\.

Overall, this work represents an important step toward practical, energy\-efficient, real\-time virtual sensing through the development of sparse neural operators\. By introducing surrogate\-free activation\-sparsity regularization, analyzing sparse feature representations through entropy\-based metrics, investigating improved optimization strategies for spiking neural operators, and proposing a synthetic knowledge distillation framework for reducing data requirements, this work advances several key challenges associated with deploying neural operators on low\-power edge hardware\. Collectively, these contributions help bridge the gap between increasingly sophisticated scientific machine learning architectures and practical hardware implementations capable of accurate, low\-latency spatial\-temporal physics reconstruction\.

## 5Acknowledgments

This work was made possible by support from the National Center for Supercomputing Applications \(NCSA\) and the U\.S\. Department of Energy Office of Nuclear Energy, specifically the University Nuclear Leadership Program’s Graduate Fellowship\. LLMs were utilized solely for language modification and structuring\.

## References

- \[1\]Electric Power Research Institute \(EPRI\)\.Powering intelligence: Analyzing artificial intelligence and data center energy consumption, May 2024\.
- \[2\]Sophia Chen\.Data centres will use twice as much energy by 2030 driven by ai\.[https://www\.nature\.com/articles/d41586\-025\-01113\-z](https://www.nature.com/articles/d41586-025-01113-z), 2025\.Nature News, 10 April 2025\.
- \[3\]Raffaella Testoni et al\.Review of nuclear microreactors: Status, potentialities and challenges\.Progress in Nuclear Energy, 138:103822, 2021\.
- \[4\]M\.Z\. Naser\.Fundamental flaws of physics\-informed neural networks and explainability methods in engineering systems\.Computers & Industrial Engineering, 212:111704, 2026\.
- \[5\]Kazuma Kobayashi, Farid Ahmed, and Syed Bahauddin Alam\.Virtual sensing to enable real\-time monitoring of inaccessible locations & unmeasurable parameters, 2024\.
- \[6\]William Howes, Jason Yoo, Kazuma Kobayashi, Subhankar Sarkar, Farid Ahmed, Souvik Chakraborty, and Syed Bahauddin Alam\.Real\-time sensing of inaccessible physical fields via an edge\-deployable hardware\-portable graph neural operator, 2026\.
- \[7\]Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar\.Neural operator: Learning maps between function spaces with applications to pdes\.Journal of Machine Learning Research, 24\(89\):1–97, 2023\.
- \[8\]Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis\.Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators\.Nat\. Mach\. Intell\., 3\(3\):218–229, March 2021\.
- \[9\]Zongyi Li, Nikola B\. Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew M\. Stuart, and Anima Anandkumar\.Fourier neural operator for parametric partial differential equations\.CoRR, abs/2010\.08895, 2020\.
- \[10\]Haozhou Zhai, Jinwei Du, Yuhui Ai, and Tianjiang Hu\.Edge deployment of deep networks for visual detection: A review\.IEEE Sensors Journal, 25\(11\):18662–18683, 2025\.
- \[11\]Chiara Contoli, Lorenzo Calisti, Giacomo Di Fabrizio, Nicholas Kania, Alessandro Bogliolo, and Emanuele Lattanzi\.A power\-aware vision\-based virtual sensor for real\-time edge computing\.J\. Real\-Time Image Process\., 21\(4\), May 2024\.
- \[12\]SAMANWOY GHOSH\-DASTIDAR and HOJJAT ADELI\.Spiking neural networks\.International Journal of Neural Systems, 19\(04\):295–308, 2009\.PMID: 19731402\.
- \[13\]Jiadong Wu, Yinan Wang, Zhiwei Li, Lun Lu, and Qingjiang Li\.A review of computing with spiking neural networks\.Computers, Materials and Continua, 78\(3\):2909–2939, 2024\.
- \[14\]Kashu Yamazaki, Viet\-Khoa Vo\-Ho, Darshan Bulsara, and Ngan Le\.Spiking neural networks and their applications: A review\.Brain Sciences, 12\(7\), 2022\.
- \[15\]Parijata Majumdar\.Spiking neural networks: a comprehensive review of diverse applications, research progress, challenges and future research directions\.Evol\. Syst\., 16\(4\), December 2025\.
- \[16\]Jason K\. Eshraghian, Max Ward, Emre O\. Neftci, Xinxin Wang, Gregor Lenz, Girish Dwivedi, Mohammed Bennamoun, Doo Seok Jeong, and Wei D\. Lu\.Training spiking neural networks using lessons from deep learning\.Proceedings of the IEEE, 111\(9\):1016–1054, 2023\.
- \[17\]Shailesh Garg and Souvik Chakraborty\.Neuroscience inspired neural operator for partial differential equations\.Journal of Computational Physics, 515:113266, 2024\.
- \[18\]Kungjui Hou, Kunlun Wu, and Yongcheng Zhou\.Adaptive and lightweight surrogate gradients: enhancing training efficiency of spiking neural networks\.Frontiers in Neuroscience, Volume 20 \- 2026, 2026\.
- \[19\]Qian Zhang, Chenxi Wu, Adar Kahana, Youngeun Kim, Yuhang Li, George Em Karniadakis, and Priyadarshini Panda\.Artificial to spiking neural networks conversion for scientific machine learning, 2023\.
- \[20\]Haiyan Jiang, Srinivas Anumasa, Giulia De Masi, Huan Xiong, and Bin Gu\.A unified optimization framework of ANN\-SNN conversion: Towards optimal mapping from activation values to firing rates\.In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett, editors,Proceedings of the 40th International Conference on Machine Learning, volume 202 ofProceedings of Machine Learning Research, pages 14945–14974\. PMLR, 23–29 Jul 2023\.
- \[21\]William Howes, Farid Ahmed, Kazuma Kobayashi, Souvik Chakraborty, and Syed Bahauddin Alam\.Neuroscience inspired graph operators towards edge\-deployable virtual sensing for irregular geometries, 2026\.
- \[22\]Isha Jain, Shailesh Garg, Shaurya Shriyam, and Souvik Chakraborty\.Hybrid variable spiking graph neural networks for energy\-efficient scientific machine learning\.Journal of the Mechanics and Physics of Solids, 200:106152, 2025\.
- \[23\]Han Wan, Rui Zhang, and Hao Sun\.Spectral\-inspired operator learning with limited data and unknown physics, 2026\.
- \[24\]Lan Chen, Houhui Yi, Zhiyang Zhang, Muwei Liu, and Wenjun Liu\.Fit\-deepokan: enhancing neural operator learning of soliton dynamics via FitNets distillation\.Nonlinear Dyn\., 114\(4\), February 2026\.
- \[25\]Ilan Price, Nicholas Daultry Ball, Samuel C\. H\. Lam, Adam C\. Jones, and Jared Tanner\.Deep neural network initialization with sparsity inducing activations, 2024\.
- \[26\]Mark Kurtz, Justin Kopinsky, Rati Gelashvili, Alexander Matveev, John Carr, Michael Goin, William Leiserson, Sage Moore, Nir Shavit, and Dan Alistarh\.Inducing and exploiting activation sparsity for fast inference on deep neural networks\.In Hal Daumé III and Aarti Singh, editors,Proceedings of the 37th International Conference on Machine Learning, volume 119 ofProceedings of Machine Learning Research, pages 5533–5543\. PMLR, 13–18 Jul 2020\.
- \[27\]Garrick Orchard, E\. Paxon Frady, Daniel Ben Dayan Rubin, Sophia Sanborn, Sumit Bam Shrestha, Friedrich T\. Sommer, and Mike Davies\.Efficient neuromorphic signal processing with loihi 2\.In2021 IEEE Workshop on Signal Processing Systems \(SiPS\), pages 254–259, 2021\.
- \[28\]Hongru Yang, Ziyu Jiang, Ruizhe Zhang, Yingbin Liang, and Zhangyang Wang\.Neural networks with sparse activation induced by large bias: Tighter analysis with bias\-generalized ntk, 2024\.
- \[29\]Jacob H\. Seidman, Georgios Kissas, Paris Perdikaris, and George J\. Pappas\.Nomad: Nonlinear manifold decoders for operator learning, 2022\.
- \[30\]Tapas Tripura and Souvik Chakraborty\.Wavelet neural operator: a neural operator for parametric partial differential equations, 2022\.
- \[31\]Myeonginn Kang and Seokho Kang\.Knowledge distillation with insufficient training data for regression\.Engineering Applications of Artificial Intelligence, 132:108001, 2024\.
- \[32\]Tianxun Zhou and Keng\-Hwee Chiam\.Synthetic data generation method for data\-free knowledge distillation in regression neural networks\.Expert Systems with Applications, 227:120327, October 2023\.
- \[33\]Pengzhan Jin, Shuai Meng, and Lu Lu\.Mionet: Learning multiple\-input operators via tensor product, 2022\.
- \[34\]Farid Ahmed, Md\. Hasan Nasrullah, Istiak Ahmad, Kazuma Kobayashi, and Syed Bahauddin Alam\.Enhancing thermo\-hydraulic performance in dimpled channels with wavy tape inserts for heat pipe & heat exchanger design with complex energy systems\.Case Studies in Thermal Engineering, 60:104583, 2024\.
- \[35\]ANSYS Inc\.Ansys fluent user’s guide, 2024\.
- \[36\]Jason K\. Eshraghian\.snntorch documentation, 2021\.
- \[37\]National Center for Supercomputing Applications \(NCSA\)\.Delta user documentation, 2025\.Accessed: 2025\-11\-10\.
- \[38\]Sneha Aenugu\.Training spiking neural networks using reinforcement learning, 2020\.
- \[39\]Yanhu Chen, Cen Wang, Hongxiang Guo, Xiong Gao, and Jian Wu\.Accelerating spiking neural networks using quantum algorithm with high success probability and high calculation accuracy\.Neurocomputing, 493:435–444, 2022\.
- \[40\]Hongrong Cheng, Miao Zhang, and Javen Qinfeng Shi\.A survey on deep neural network pruning\-taxonomy, comparison, analysis, and recommendations, 2024\.