Iterative Refinement Neural Operators are Learned Fixed-Point Solvers: A Principled Approach to Spectral Bias Mitigation

# Iterative Refinement Neural Operators are Learned Fixed-Point Solvers: A Principled Approach to Spectral Bias Mitigation
Source: [https://arxiv.org/html/2605.24041](https://arxiv.org/html/2605.24041)
Shuyuan ShangCUHK ShenzhenXiaopeng WangDartmouth CollegePu RenLawrence Berkeley National LabYaoqing YangDartmouth College

###### Abstract

Neural operators serve as fast, data\-driven surrogates for scientific modeling but typically rely on a monolithic, single\-pass inference procedure that struggles to resolve high\-frequency details, a limitation known as spectral bias\. We introduce the Iterative Refinement Neural Operator \(IRNO\), which augments pre\-trained operators with a learned refinement module iteratively applied via fixed\-point iteration\. IRNO decomposes the prediction into a coarse initialization followed by successive residual corrections, paralleling classical numerical solvers\. Under local assumptions, we establish contraction of the induced operator, ensuring convergence to a unique fixed point\. To explicitly target high\-frequency errors, we propose a progressive spectral loss that adaptively increases penalty on high\-frequency components over refinement steps during training\. Across physical systems, IRNO consistently lowers error, with up to 56\.05% improvement on turbulent flow\. On Active Matter, spectral analysis reveals that, relative to base operator, the normalized error ratios decrease to 27\.72–36\.10% in low\-, 5\.07–6\.68% in mid\-, and 1\.48–2\.04% in high\-frequencies, remaining stable beyond the trained iteration count\. Code is available at[https://github\.com/xiaotianliu\-dartmouth/Iterative\_Refinement\_Neural\_Operator](https://github.com/xiaotianliu-dartmouth/Iterative_Refinement_Neural_Operator)

## 1Introduction

![Refer to caption](https://arxiv.org/html/2605.24041v1/x1.png)

\(a\) Ground Truth![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/teaser_fno.png)

\(b\) Base Operator \(FNO\)
![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/teaser_fig_2_step.png)Iteration 2
![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/teaser_fig_4_step.png)Iteration 4
![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/teaser_fig_6_step.png)Iteration 6
Training Cutoff:K=6K=6

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/teaser_fig_8_step.png)Iteration 8
![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/teaser_fig_10_step.png)Iteration 10
![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/teaser_fig_12_step.png)Iteration 12

\(c\) Iterative Refinement Process \(k=1→12k=1\\to 12\)

Figure 1:Iterative Refinement on ERA5 16×\\timesSuper\-Resolution\.The base operator \(FNO\) captures large\-scale atmospheric structure but under\-resolves fine\-scale details\. Successive refinement steps \(k=1→12k=1\\rightarrow 12\) from IRNO progressively reduce high\-frequency error while preserving global coherence\. Horizontal dashed line indicates the training cutoff atK=6K=6; iterations beyond this point demonstrate stable extrapolation and continue to recover fine\-scale details\. All panels share a common color scale\. Shown field: kinetic energy\.Neural operators have emerged as an effective approach for learning mappings between function spaces, enabling fast surrogate models for parametric Partial Differential Equations \(PDEs\) and complex physical systems\. Architectures such as the Fourier Neural Operator \(FNO\)Liet al\.\([2021](https://arxiv.org/html/2605.24041#bib.bib3)\)and DeepONetLuet al\.\([2021](https://arxiv.org/html/2605.24041#bib.bib17)\)achieve high accuracy across a range of scientific tasks while maintaining tractable inference computation costs\. In this work, we introduce theIterative Refinement Neural Operator \(IRNO\), a learned refinement operator that augments a pre\-trained neural operator with iterative application of a shared\-weight update as a test\-time inference loop\. Instead of producing a single solution estimate, IRNO defines a dynamical process in function space that progressively corrects residual error\. IRNO enables consistent error reduction across refinement steps, remains stable when evaluated for iteration counts beyond those used during training, and improves spectral precision without retraining the base model\.

Neural operator architectures are typically trained to approximate the solution operator through a single forward evaluation\. Improvements in accuracy are therefore primarily achieved through*training\-time scaling*, including increased model capacity, higher\-resolution data, or larger training sets\. Unlike monolithic operators, IRNO decouples accuracy improvement from retraining through*test\-time iteration*\. Furthermore, the function\-space formulation allows IRNO trained with one base operator to refine predictions from different operators\.

From the perspective of numerical analysis, this refinement process defines a learned fixed\-point iteration in function space\. IRNO applies a sequence of residual corrections that progressively reduce the remaining error, aligning the method with classical residual\-based solvers\(Almgrenet al\.,[2013](https://arxiv.org/html/2605.24041#bib.bib23)\)while preserving the expressiveness and efficiency of modern neural operators\.

We provide a theoretical analysis showing that IRNO converges as a contraction mapping in function space in Section[3](https://arxiv.org/html/2605.24041#S3)\. The analysis models each refinement update as a locally affine map in a neighborhood of the solution manifold\. We establish conditions for monotonic convergence toward a unique fixed point and show presence of a residual floor when bias is present\. These results provide a formal interpretation of the convergence and saturation behavior when extrapolated beyond training iterations\.

Across systems and tasks, IRNO yields consistent error reduction across refinement steps, strong improvements in mid\-to\-high frequency errors, stable extrapolation to iteration counts up to twice those used during training, transferability across base operators, and Pareto\-dominant accuracy–compute trade\-offs relative to capacity\-matched monolithic baselines\. Figure[1](https://arxiv.org/html/2605.24041#S1.F1)illustrates the refinement process on ERA5\(Hersbachet al\.,[2020](https://arxiv.org/html/2605.24041#bib.bib22); Renet al\.,[2025](https://arxiv.org/html/2605.24041#bib.bib20)\), where successive iterations recover fine\-scale structure while preserving large\-scale coherence\.

Our contributions are summarized as follows:

1. 1\.Iterative Refinement Operator\.We introduce a modular, architecture\-agnostic iterative refinement operator that enhances predictions of pre\-trained neural operators, enabling multi\-step error reduction at inference time\.
2. 2\.Contraction\-Based Analysis\.We model refinement as a fixed\-point iteration in function space and derive conditions for monotonic convergence, stable extrapolation over iteration, and a bias\-controlled residual error floor, which we empirically validate by correlating the fixed\-point bias with the observed error plateau in Figure[3](https://arxiv.org/html/2605.24041#S3.F3)\.
3. 3\.Consistent Error Reduction, Spectral Bias Mitigation, and Cross\-Operator Transfer\.We conduct a systematic empirical study on efficacy of IRNO across multiple physical systems in Section[4](https://arxiv.org/html/2605.24041#S4)\. IRNO achieves consistent error reductions across benchmarks \(Section[4\.2](https://arxiv.org/html/2605.24041#S4.SS2)\), demonstrates effective spectral bias mitigation \(Section[4\.3](https://arxiv.org/html/2605.24041#S4.SS3)\), exhibits cross\-operator transferability between distinct neural operator architectures \(Section[4\.5](https://arxiv.org/html/2605.24041#S4.SS5)\), and extends effectively to unstructured mesh settings \(Section[4\.6](https://arxiv.org/html/2605.24041#S4.SS6)\)\. Cost–performance analysis \(Section[4\.8](https://arxiv.org/html/2605.24041#S4.SS8)\) shows that IRNO traces a Pareto\-dominant frontier relative to capacity\-matched baselines\.

Detailed discussions of related work are available in Appendix[B](https://arxiv.org/html/2605.24041#A2)\.

## 2Methodology

In this section, we formalize IRNO as an inference\-time fixed\-point scheme in which a pre\-trained base operator provides an initial estimate that is iteratively refined by residual correction\. Let𝒳\\mathcal\{X\}andℋ\\mathcal\{H\}be Banach spaces representing the input \(e\.g\., initial conditions\) and solution function spaces, respectively\. Our goal is to learn the solution operator𝒢:𝒳→ℋ\\mathcal\{G\}:\\mathcal\{X\}\\rightarrow\\mathcal\{H\}from a finite dataset𝒟=\{\(xi,yi\)\}i=1N\\mathcal\{D\}=\\\{\(x\_\{i\},y\_\{i\}\)\\\}\_\{i=1\}^\{N\}, whereyi=𝒢​\(xi\)y\_\{i\}=\\mathcal\{G\}\(x\_\{i\}\)\.

### 2\.1Iterative Refinement Neural Operator \(IRNO\)

#### 2\.1\.1Inference\-time Refinement Scheme

The inference process is a two\-stage dynamical system:

1. 1\.Initialization:A pre\-trained base operatorTbase:𝒳→ℋT\_\{\\text\{base\}\}:\\mathcal\{X\}\\rightarrow\\mathcal\{H\}\(e\.g\., FNO\) produces an initial coarse ansatzh0=Tbase​\(x\)h\_\{0\}=T\_\{\\text\{base\}\}\(x\)\. This captures the dominant low\-frequency structure but may lack local fidelity\.
2. 2\.Iterative Correction:A learned refinement operatorΦθ:𝒳×ℋ→ℋ\\Phi\_\{\\theta\}:\\mathcal\{X\}\\times\\mathcal\{H\}\\rightarrow\\mathcal\{H\}estimates the local residual\. The solution state is updated via the fixed\-point iteration: hk\+1=hk\+α⋅Φθ​\(x,hk\),k=0,…,K−1,h\_\{k\+1\}=h\_\{k\}\+\\alpha\\cdot\\Phi\_\{\\theta\}\(x,h\_\{k\}\),\\quad k=0,\\dots,K\-1,whereα∈\(0,1\]\\alpha\\in\(0,1\]is a step size of choice, controlling the trade\-off between convergence speed and stability\.

This formulation enablesΦθ\\Phi\_\{\\theta\}to refine predictions from different base operatorsTbaseT\_\{\\text\{base\}\}\(Section[4\.5](https://arxiv.org/html/2605.24041#S4.SS5)\) and allows us to establish IRNO as a contraction mapping in function space \(Section[3](https://arxiv.org/html/2605.24041#S3)\)\. The comparison between the workflows of a single\-pass and the iterative scheme is shown in Figure[2](https://arxiv.org/html/2605.24041#S2.F2)\.

##### Architectural Choice\.

The refinement operatorΦθ\\Phi\_\{\\theta\}is designed to represent a stable update rule that maps the current iteratehkh\_\{k\}to a corrective residual, rather than directly approximating the solution operator itself\. Architectural considerations are therefore guided by numerical properties of the induced inference dynamics: \(i\)*smoothness*to ensure stable iteration, \(ii\)*multi\-scale expressiveness*to capture spectral corrections, and \(iii\)*computational efficiency*through weight sharing across iterations\.

These requirements can be satisfied by various multi\-scale encoder\-decoder architectures with skip connections\. In our experiments \(Section[4](https://arxiv.org/html/2605.24041#S4)\), we instantiateΦθ\\Phi\_\{\\theta\}as a lightweight U\-NetRonnebergeret al\.\([2015](https://arxiv.org/html/2605.24041#bib.bib45)\), though the framework is architecture\-agnostic and compatible with any operator backbone that meets these design criteria\.

### 2\.2Training the Refinement Operator

We trainΦθ\\Phi\_\{\\theta\}using a composite objective that enforces trajectory control, frequency recovery, and convergence stability\.

#### 2\.2\.1Multi\-step Supervision \(Trajectory Control\)

To ensure the iterative process creates a stable dynamical trajectory toward the solution, we impose supervision on the output of every refinement step\. We minimize theL2L^\{2\}error norm averaged over the trajectory:

ℒspatial=1K​∑k=1K‖hk−y‖2\.\\mathcal\{L\}\_\{\\text\{spatial\}\}=\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}\\\|h\_\{k\}\-y\\\|^\{2\}\.This deep supervision prevents the network from learning unstable intermediate updates that might drift off the solution manifold before being corrected\. Appendix[E](https://arxiv.org/html/2605.24041#A5)demonstrates that without trajectory control, the learned operator fails to satisfy the strong monotonicity condition underlying the contraction theory \(Section[3](https://arxiv.org/html/2605.24041#S3)\)\.

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/flow_chart.png)Figure 2:Overview of the Iterative Refinement Neural Operator \(IRNO\)\.\(a\) Standard neural operators approximate the solution in a monolithic, single forward evaluation, often losing fine\-scale details\. \(b\) IRNO reformulates inference as a dynamic process\. A Base Operator provides a coarse initialization, which is then iteratively corrected by a shared\-weight refinement operatorΦθ\\Phi\_\{\\theta\}\. At each stepkk, the network concatenates the original inputxxwith the current estimatehkh\_\{k\}to predict a residual update, progressively resolving fine\-scale local physics\.
#### 2\.2\.2Progressive Spectral Loss \(Frequency Recovery\)

StandardL2L^\{2\}loss in the physical space often fails to capture high\-frequency details\. We apply a progressive spectral loss that dynamically pivots from coarse to fine scales as iterations proceed\. Early refinement steps focus on correcting coarse structure, while later steps progressively recover local features, implementing multi\-resolution numerical refinement\.

Leth^k=ℱ​\(hk\)\\hat\{h\}\_\{k\}=\\mathcal\{F\}\(h\_\{k\}\)be the Fourier transform of the estimate\. The spectral loss at stepkkis defined as

ℒspectral\(k\)=1H​W​1ρ¯k​∑ωρ​\(ω,λk\)⋅\|\|h^k​\(ω\)\|−\|y^​\(ω\)\|\|2,\\mathcal\{L\}\_\{\\text\{spectral\}\}^\{\(k\)\}=\\frac\{1\}\{HW\}\\frac\{1\}\{\\bar\{\\rho\}\_\{k\}\}\\sum\_\{\\omega\}\\rho\(\\omega,\\lambda\_\{k\}\)\\cdot\\left\|\|\\hat\{h\}\_\{k\}\(\\omega\)\|\-\|\\hat\{y\}\(\\omega\)\|\\right\|^\{2\},where the weighting functionρ\\rhopenalizes high frequencies exponentially asρ​\(ω,λk\)=1\+\(\|ω\|/\|ω\|nyq\)λk,\\rho\(\\omega,\\lambda\_\{k\}\)=1\+\\left\(\|\\omega\|/\|\\omega\|\_\{\\text\{nyq\}\}\\right\)^\{\\lambda\_\{k\}\},andρ¯k\\bar\{\\rho\}\_\{k\}is the averaged weights\. The exponentλk\\lambda\_\{k\}increases linearly fromλstart\\lambda\_\{\\text\{start\}\}toλend\\lambda\_\{\\text\{end\}\}across iterations\. The total spectral loss is the average over the trajectory,

ℒspectral=1K​∑kℒspectral\(k\)\.\\mathcal\{L\}\_\{\\text\{spectral\}\}=\\frac\{1\}\{K\}\\sum\_\{k\}\\mathcal\{L\}\_\{\\text\{spectral\}\}^\{\(k\)\}\.

#### 2\.2\.3Fixed\-point Regularization \(Convergence Stability\)

Finally, we explicitly enforce the stability condition of a fixed\-point solver\. If the input is the exact solutionyy, the residual correction should be zero\. We add a regularization term

ℒfp=‖Φθ​\(x,y\)‖2\.\\mathcal\{L\}\_\{\\text\{fp\}\}=\\\|\\Phi\_\{\\theta\}\(x,y\)\\\|^\{2\}\.This ensures thatyyis a fixed point of the learned dynamics, preventing the model from pushing the state away from the true solution once converged\. The fixed\-point regularizationℒfp\\mathcal\{L\}\_\{\\text\{fp\}\}directly minimizes the bias term in Theorem[3\.1](https://arxiv.org/html/2605.24041#S3.Thmtheorem1), thereby reducing the error floor from Corollary[3\.3](https://arxiv.org/html/2605.24041#S3.Thmtheorem3)\.

##### Total Objective\.

The full loss function is a weighted sum:

ℒtotal=ℒspatial\+βspectral​ℒspectral\+βfp​ℒfp,\\mathcal\{L\}\_\{\\text\{total\}\}=\\mathcal\{L\}\_\{\\text\{spatial\}\}\+\\beta\_\{\\text\{spectral\}\}\\mathcal\{L\}\_\{\\text\{spectral\}\}\+\\beta\_\{\\text\{fp\}\}\\mathcal\{L\}\_\{\\text\{fp\}\},whereβspectral\\beta\_\{\\text\{spectral\}\}andβfp\\beta\_\{\\text\{fp\}\}are hyperparameters balancing spectral corrections, convergence stability, and spatial fidelity\.

Importantly,Φθ\\Phi\_\{\\theta\}learns the iteration\-invariant*update dynamics*rather than the solution mapping itself and can be applied fork\>Kk\>Ksteps at inference time\.

## 3Theoretical Analysis

We analyze the convergence behavior of the proposed iterative refinement operatorΦ​\(x,⋅\)\\Phi\(x,\\cdot\)\. Letek=y−hke\_\{k\}=y\-h\_\{k\}denote the residual error at iterationkk\. The analysis characterizes how‖ek‖\\\|e\_\{k\}\\\|evolves under mild regularity assumptions on the learned refinement operatorΦ​\(x,⋅\)\\Phi\(x,\\cdot\)and the pre\-trained base operatorTbaseT\_\{\\text\{base\}\}\.

### 3\.1Theoretical Setup

Formally, we assume the following properties hold in a local neighborhood of the true solutionyy:

1. 1\.Local affine approximation:In a neighborhoodBδ​\(y\)=\{hi∈ℋ:‖ei‖<δ\}B\_\{\\delta\}\(y\)=\\\{h\_\{i\}\\in\\mathcal\{H\}:\\\|e\_\{i\}\\\|<\\delta\\\}of the true solution, we assume the refinement operatorΦ​\(x,⋅\)\\Phi\(x,\\cdot\)admits an affine linearization of the form Φ​\(x,h\)=b​\(x\)\+A​\(x,h\)​e\+R​\(x,h\),\\Phi\(x,h\)=b\(x\)\+A\(x,h\)e\+R\(x,h\),wheree=y−he=y\-h,A​\(x,h\):ℋ→ℋA\(x,h\):\\mathcal\{H\}\\to\\mathcal\{H\}is a bounded linear operator,b​\(x\)b\(x\)is the bias term, and the remainderR​\(x,h\)R\(x,h\)satisfies‖R​\(x,h\)‖≤L2​‖e‖2\\\|R\(x,h\)\\\|\\leq\\tfrac\{L\}\{2\}\\\|e\\\|^\{2\}, for someL\>0L\>0\. For notational simplicity, we suppress thexxdependence and writeb​\(x\)b\(x\)asbb\. By evaluating the affine approximationΦ​\(x,h\)\\Phi\(x,h\)ath=yh=y, we have the bias termb=Φ​\(x,y\)b=\\Phi\(x,y\)\. Ideally, a perfectly learned refinement operator would satisfy‖b‖=0\\\|b\\\|=0\. This decomposition is analogous to a first\-order Taylor expansion ofΦ​\(x,⋅\)\\Phi\(x,\\cdot\)around the true solutionyy, withA​\(x,h\)A\(x,h\)being a linear operator acting on residualeeandR​\(x,h\)R\(x,h\)capturing higher\-order remainder terms\.
2. 2\.Lipschitz continuity of the operator\-valued map and strong monotonicity:The mappingh↦A​\(x,h\)h\\mapsto A\(x,h\)should vary smoothly in a neighborhood of the fixed point\. Assume there existsμ\>0\\mu\>0such that for allh∈Bδ​\(y\)h\\in B\_\{\\delta\}\(y\), ‖A​\(x,h\)−A​\(x,y\)‖op≤μ​‖e‖\\\|A\(x,h\)\-A\(x,y\)\\\|\_\{\\text\{op\}\}\\leq\\mu\\\|e\\\|whereA​\(x,y\)≡A​\(x\)A\(x,y\)\\equiv A\(x\)\. Furthermore, we assume the linearization at the solution,A​\(x,y\)A\(x,y\), is locallyboundedandstrongly monotone, i\.e\., there exist constants0<m≤M<∞0<m\\leq M<\\inftysuch that ⟨A​\(x,y\)​e,e⟩≥m​‖e‖2,‖A​\(x,y\)‖op≤M\.\\langle A\(x,y\)e,e\\rangle\\geq m\\\|e\\\|^\{2\},\\quad\\\|A\(x,y\)\\\|\_\{\\text\{op\}\}\\leq M\.Under this condition, choosing0<α<2​mM20<\\alpha<\\frac\{2m\}\{M^\{2\}\}guarantees‖I−α​A​\(x,y\)‖op=q<1\\\|I\-\\alpha A\(x,y\)\\\|\_\{\\text\{op\}\}=q<1, since for any unit vectoree, ‖\(I−α​A\)​e‖2≤1−2​α​m\+α2​M2<1\.\\\|\(I\-\\alpha A\)e\\\|^\{2\}\\leq 1\-2\\alpha m\+\\alpha^\{2\}M^\{2\}<1\. The Lipschitz continuity assumption onh↦A​\(x,h\)h\\mapsto A\(x,h\)is standard in convergence analysisWerner and Hofmann \([2019](https://arxiv.org/html/2605.24041#bib.bib46)\); Rastogiet al\.\([2020](https://arxiv.org/html/2605.24041#bib.bib47)\); Kovachkiet al\.\([2023a](https://arxiv.org/html/2605.24041#bib.bib4)\), while the strong monotonicity condition is empirically validated in Appendix[E](https://arxiv.org/html/2605.24041#A5)and ensures the spectral radius conditionq<1q<1is achievable by explicit choices ofα\\alpha\.
3. 3\.Initialization quality and invariant\-ball:The base operatorTbaseT\_\{\\text\{base\}\}provides a sufficiently accurate initial ansatz such that ‖e0‖<min⁡\{δ,1−q2​c\},c=α​\(L2\+μ\),\\\|e\_\{0\}\\\|<\\min\\left\\\{\\delta,\\frac\{1\-q\}\{2c\}\\right\\\},\\quad c=\\alpha\\left\(\\frac\{L\}\{2\}\+\\mu\\right\),whereδ\\deltais the radius in Assumption[1](https://arxiv.org/html/2605.24041#S3.I1.i1)\. We further assume there existsrrsuch that‖e0‖≤r<min⁡\{δ,1−q2​c\}\\\|e\_\{0\}\\\|\\leq r<\\min\\left\\\{\\delta,\\frac\{1\-q\}\{2c\}\\right\\\}and α​‖b‖≤r​\(1−q−c​r\),\\alpha\\\|b\\\|\\leq r\(1\-q\-cr\),which is a small\-bias condition ensuring the iterates remain inBr​\(y\)⊂Bδ​\(y\)B\_\{r\}\(y\)\\subset B\_\{\\delta\}\(y\)\. This bias magnitude is directly minimized by the fixed\-point regularizationℒfp\\mathcal\{L\}\_\{\\mathrm\{fp\}\}in Section[2\.2\.3](https://arxiv.org/html/2605.24041#S2.SS2.SSS3)\. A sufficient condition and derivation are given in Appendix[F](https://arxiv.org/html/2605.24041#A6)\.

Intuitively, Assumption[1](https://arxiv.org/html/2605.24041#S3.I1.i1)reflects the fact thatΦ\\Phiis trained by deep supervision to approximate the residuale=y−he=y\-h, and its local linearization is encouraged to produce residual\-correcting updates near the solution manifold\. Assumption[2](https://arxiv.org/html/2605.24041#S3.I1.i2)captures architectural smoothness induced by convolutional layers inΦ\\Phi, which empirically limit abrupt changes in the refinement dynamics\. The strong monotonicity condition further ensures that a valid step sizeα\\alphaexists, and is empirically validated in Appendix[E](https://arxiv.org/html/2605.24041#A5)\. Finally, Assumption[3](https://arxiv.org/html/2605.24041#S3.I1.i3)formalizes the role ofTbaseT\_\{\\text\{base\}\}as providing a reasonable initialization within the basin of attraction, and requires the learned bias‖b‖\\\|b\\\|to remain small enough to keep iterates from escaping the local neighborhood, a condition directly enforced by the fixed\-point regularizationℒfp\\mathcal\{L\}\_\{\\text\{fp\}\}in Section[2\.2\.3](https://arxiv.org/html/2605.24041#S2.SS2.SSS3)\.

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/AM_bias_error.png)\(a\) Active Matter \(FNO\)
![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/TR_bias_error.png)\(b\) TR\-2D \(TFNO\)

Figure 3:Empirical validation of bias–error floor relationship\.Scatter plots of minimum errormink⁡‖ek‖\\min\_\{k\}\\\|e\_\{k\}\\\|as a function of bias magnitude‖b‖=‖Φθ​\(x,y\)‖\\\|b\\\|=\\\|\\Phi\_\{\\theta\}\(x,y\)\\\|overk=24k=24refinement steps\.\(a\)Active Matter: Pearsonr=0\.933r=0\.933\(p≪10−10p\\ll 10^\{\-10\}\)\.\(b\)TR\-2D: Pearsonr=0\.949r=0\.949\(p≪10−10p\\ll 10^\{\-10\}\)\. Both systems exhibit a strong linear dependence between the asymptotic error floor and the bias, consistent with Corollary[3\.3](https://arxiv.org/html/2605.24041#S3.Thmtheorem3)\.
### 3\.2Main Results

###### Theorem 3\.1\(Quadratic–Linear Convergence\)\.

Under the above assumptions, the iteration satisfies

‖ek\+1‖≤q​‖ek‖\+c​‖ek‖2\+α​‖b‖\.\\\|e\_\{k\+1\}\\\|\\leq q\\\|e\_\{k\}\\\|\+c\\\|e\_\{k\}\\\|^\{2\}\+\\alpha\\\|b\\\|\.Whenb=0b=0, the scheme is locally contractive:

‖ek\+1‖≤\(q\+c​‖ek‖\)​‖ek‖,with​q\+c​‖ek‖<1\.\\\|e\_\{k\+1\}\\\|\\leq\(q\+c\\\|e\_\{k\}\\\|\)\\\|e\_\{k\}\\\|,\\quad\\text\{with \}q\+c\\\|e\_\{k\}\\\|<1\.Thus, the errors decrease monotonically andhk→yh\_\{k\}\\to y\.

Proof of Theorem[3\.1](https://arxiv.org/html/2605.24041#S3.Thmtheorem1)is available in Appendix[C](https://arxiv.org/html/2605.24041#A3)\. The result in Theorem[3\.1](https://arxiv.org/html/2605.24041#S3.Thmtheorem1)characterizes IRNO as a learned iterative solver whose convergence behavior is governed by the contraction factorqq, smoothness of its update dynamics, and the bias, providing a theoretical basis for the monotonic error decay and stable extrapolation observed empirically in Section[4](https://arxiv.org/html/2605.24041#S4)\. Detailed interpretations of the error recursion in Theorem[3\.1](https://arxiv.org/html/2605.24041#S3.Thmtheorem1)are provided in Appendix[C](https://arxiv.org/html/2605.24041#A3)\.

###### Corollary 3\.2\(Geometric Convergence and Iteration Complexity\)\.

In the linear\-dominant regime withb=0b=0, the error decays geometrically:

‖ek‖≲qk​‖e0‖\.\\\|e\_\{k\}\\\|\\lesssim q^\{k\}\\\|e\_\{0\}\\\|\.To achieve‖ek‖≤ε\\\|e\_\{k\}\\\|\\leq\\varepsilon, it suffices to perform

k=O​\(log⁡\(‖e0‖/ε\)log⁡\(1/q\)\)k=O\\\!\\left\(\\frac\{\\log\(\\\|e\_\{0\}\\\|/\\varepsilon\)\}\{\\log\(1/q\)\}\\right\)iterations of refinement\.

Corollary[3\.2](https://arxiv.org/html/2605.24041#S3.Thmtheorem2)assumesb=0b=0, allowing arbitrarily small error with sufficient iterations\. In practice, however,Φθ\\Phi\_\{\\theta\}may exhibit non\-zero biasb≠0b\\neq 0, introducing a limiting error floor\. The following corollary quantifies the upper bound on the bias\-induced error floor\.

###### Corollary 3\.3\(Convergence with Bias\)\.

Ifb=Φ​\(x,y\)≠0b=\\Phi\(x,y\)\\neq 0and‖b‖\\\|b\\\|is sufficiently small, then there exists a unique fixed pointh∗h^\{\*\}satisfyingh∗=h∗\+α​Φ​\(x,h∗\)h^\{\*\}=h^\{\*\}\+\\alpha\\Phi\(x,h^\{\*\}\), and the iteration converges linearly toh∗h^\{\*\}for any initialization within a neighborhood ofyy\. The limiting error satisfies

‖e∗‖≤α​‖b‖1−q\+O​\(‖b‖2\)\.\\\|e^\{\*\}\\\|\\leq\\frac\{\\alpha\\\|b\\\|\}\{1\-q\}\+O\(\\\|b\\\|^\{2\}\)\.

Proofs of Corollaries[3\.2](https://arxiv.org/html/2605.24041#S3.Thmtheorem2)and[3\.3](https://arxiv.org/html/2605.24041#S3.Thmtheorem3)are available in Appendix[D](https://arxiv.org/html/2605.24041#A4)\. Under mild assumptions, the refinement mapT​\(h\)=h\+α​Φ​\(x,h\)T\(h\)=h\+\\alpha\\Phi\(x,h\)is locally contractive\. The contraction factorqqgoverns the convergence speed and is directly affected by the choice of step sizeα\\alpha\(Section[4\.7\.3](https://arxiv.org/html/2605.24041#S4.SS7.SSS3)\)\. In presence of non\-zero bias, convergence remains guaranteed with a limiting error floor proportional to‖b‖\\\|b\\\|\.

Empirical validation in Figure[3](https://arxiv.org/html/2605.24041#S3.F3)supports the conclusions of Corollary[3\.3](https://arxiv.org/html/2605.24041#S3.Thmtheorem3)across two physical systems \(Section[1](https://arxiv.org/html/2605.24041#S4.T1)\)\. For both Active Matter and TR\-2D, the minimum attainable error over the refinement trajectory exhibits a strong linear dependence on the bias magnitude, with Pearson correlations exceeding0\.930\.93in both cases \(p≪10−10p\\ll 10^\{\-10\}\)\. Least\-squares fits indicate that variation in the asymptotic error floor is largely explained by the magnitude of the fixed\-point residualbb, motivating fixed\-point regularization to directly reduce‖b‖\\\|b\\\|\(Section[2\.2\.3](https://arxiv.org/html/2605.24041#S2.SS2.SSS3)\)\.

## 4Experiments

We evaluate IRNO from the perspective of a learned iterative solver rather than a single\-pass predictor\. Experiments are designed to assess six core areas: \(i\) global convergence behavior across distinct physical systems and extrapolation beyond the training iteration cutoff, \(ii\) spectral dynamics of refinement and preferential reduction of mid\-to\-high frequency spectral error, \(iii\) transferability across base operators and operator architectures, \(iv\) the role of progressive spectral supervision and step sizeα\\alphathrough targeted ablations, \(v\) the cost–performance trade\-off relative to capacity scaling, and \(vi\) generalization to graph\-based operators on irregular meshes\.

### 4\.1Experimental Setup

Table 1:Iterative error reduction across physical systems\.Comparison of single\-pass base operator performance and IRNO after refinement, with FNO evaluated atK=6K=6and TFNO/WDSR atK=4K=4\. Metrics reported are VRMSE for TR\-2D and Active Matter, and ACC and RFNE for ERA5 with standard errors<0\.001<0\.001for all entries\.↑\\uparrow/↓\\downarrowindicate better direction\.##### Physical Systems\.

We evaluate IRNO on four scientific benchmarks: \(i\) Turbulent Radiative Layer \(TR\-2D\) and \(ii\) Active Matter from the WellOhanaet al\.\([2024](https://arxiv.org/html/2605.24041#bib.bib44)\); Morelet al\.\([2025](https://arxiv.org/html/2605.24041#bib.bib53)\); Wuet al\.\([2025](https://arxiv.org/html/2605.24041#bib.bib51)\); Holzschuhet al\.\([2025](https://arxiv.org/html/2605.24041#bib.bib50)\), \(iii\) ERA5 global weather16×16\\timessuper\-resolution from SuperBench\(Renet al\.,[2025](https://arxiv.org/html/2605.24041#bib.bib20); Chenet al\.,[2024](https://arxiv.org/html/2605.24041#bib.bib49); Hassanet al\.,[2023](https://arxiv.org/html/2605.24041#bib.bib54)\), and \(iv\) CE\-Gauss\(Mousaviet al\.,[2025](https://arxiv.org/html/2605.24041#bib.bib59)\), an irregular unstructured\-mesh benchmark\. These datasets are widely adopted for assessing operator stability, multi\-scale error dynamics, and high\-resolution reconstruction\. Detailed dataset description can be found in Appendix[G\.1](https://arxiv.org/html/2605.24041#A7.SS1)\.

We report variance\-scaled RMSE \(VRMSE\) for TR\-2D and Active Matter, and anomaly correlation coefficient \(ACC\) and relative Frobenius norm error \(RFNE\) for ERA5\. Definitions of the metrics are detailed in Appendix[G\.5\.1](https://arxiv.org/html/2605.24041#A7.SS5.SSS1)\.

Base Operators and Baselines\.For TR\-2D and Active Matter, we use the FNO\(Liet al\.,[2021](https://arxiv.org/html/2605.24041#bib.bib3)\)and the Tucker\-Factorized FNO \(TFNO\)\(Kossaifiet al\.,[2025](https://arxiv.org/html/2605.24041#bib.bib42)\)as base operators\. For ERA5, we evaluate both FNO and the Wide\-Activated Deep Super\-Resolution network \(WDSR\)Yuet al\.\([2018](https://arxiv.org/html/2605.24041#bib.bib25)\)\.

Table 2:ERA516×16\\timessuper\-resolution against spectral baselines\.Comparison of IRNO against state\-of\-the\-art spectral methods on ERA5, with IRNO evaluated using WDSR as the base operator\. Metrics reported are ACC and RFNE\.↑\\uparrow/↓\\downarrowindicate better direction\.We compare IRNO against three baselines\.\(i\) Standard Base Operator:single\-pass inference using the pre\-trained backbone without refinement\.\(ii\) Capacity\-Matched Residual Model:single\-shot residual correction networks with parameter count and FLOPs exceed those of IRNO executed forKKrefinement steps\. This baseline controls for improvements due to increased model capacity rather than iterative refinement\.\(iii\) State\-of\-the\-Art Models for Spectral Bias Mitigation:Hierarchical Neural Operator Transformer \(HiNOTE\)\(Luoet al\.,[2024](https://arxiv.org/html/2605.24041#bib.bib57)\)and High Frequency Scaling \(HFS\)\(Khodakaramiet al\.,[2026](https://arxiv.org/html/2605.24041#bib.bib32)\), compared on ERA5 in Table[2](https://arxiv.org/html/2605.24041#S4.T2)\.

Training and Inference\.The refinement operator is trained with a horizon ofK=4,6K=4,6steps using composite objective defined in Section[2\.2](https://arxiv.org/html/2605.24041#S2.SS2)\. The base operatorTbaseT\_\{\\text\{base\}\}is pre\-trained and frozen\. Models trained withK=4K=4are evaluated fork∈\[0,8\]k\\in\[0,8\], and models trained withK=6K=6are evaluated fork∈\[0,12\]k\\in\[0,12\]\. Detailed experiments setup is available in Appendix[G](https://arxiv.org/html/2605.24041#A7)\.

### 4\.2Global Convergence Behavior Across Physical Regimes

We evaluate whether IRNO performs stable and monotonic error reduction across distinct physical systems, consistent with the contractive dynamics predicted by Theorem[3\.1](https://arxiv.org/html/2605.24041#S3.Thmtheorem1)\.

Quantitative Performance\.Table[1](https://arxiv.org/html/2605.24041#S4.T1)reports aggregate performance metrics for TR\-2D, Active Matter, and ERA5\. Across all benchmarks and base architectures, IRNO consistently improves upon base operator inference, achieving 45–56% reductions on TR\-2D, 51–80% on Active Matter in VRMSE and 34% in RFNE while increasing 19% ACC on ERA5\.

Comparison with State\-of\-the\-Art Models for Spectral Bias Mitigation\.Table[2](https://arxiv.org/html/2605.24041#S4.T2)compares IRNO against recent spectral methods on ERA5 16×\\timessuper\-resolution\. IRNO \(WDSR\) achieves the best ACC \(0\.910\) and RFNE \(0\.195\), outperforming both High\-Frequency Scaling \(HFS\)\(Khodakaramiet al\.,[2026](https://arxiv.org/html/2605.24041#bib.bib32)\)and the hierarchical attention operator HiNOTE\(Luoet al\.,[2024](https://arxiv.org/html/2605.24041#bib.bib57)\)\. Furthermore, IRNO is complementary to these architectural spectral methods\. On Active Matter, combining IRNO with HFS further reduces VRMSE from 0\.0631 to 0\.0486, suggesting the iterative refinement mechanism compounds gains from frequency\-aware architectures\. Table[17](https://arxiv.org/html/2605.24041#A8.T17)and Figure[19](https://arxiv.org/html/2605.24041#A8.F19)in the appendix provide a detailed step\-by\-step breakdown, showing VRMSE decreasing from 0\.0631 \(HFS base\) to 0\.0486 atk=6k=6and remaining stable atk=8k=8\(0\.04870\.0487\), with IRNO achieving lower spectral error energy across the full radial frequency range\.

Error Trajectories vs\. Iteration\.Figure[4](https://arxiv.org/html/2605.24041#S4.F4)shows that error decreases monotonically during early iterations and reaches a stable plateau\. IRNO remains stable beyond the training horizon over the evaluated extrapolation range \(2×2\\times\), while longer\-horizon behavior \(8×8\\times\) depends on the step size; smaller or scheduled step sizes ensure stability\. The effect of step sizeα\\alphaon long\-horizon stability, evaluated up tok=48k=48, is analyzed in Appendix[H](https://arxiv.org/html/2605.24041#A8)\(Figure[18](https://arxiv.org/html/2605.24041#A8.F18)\)\.

### 4\.3Spectral Dynamics of Iterative Refinement

![Refer to caption](https://arxiv.org/html/2605.24041v1/x2.png)\(a\)TR\-2D\(VRMSE↓\\downarrow\)
![Refer to caption](https://arxiv.org/html/2605.24041v1/x3.png)\(b\)ACTIVE\(VRMSE↓\\downarrow\)
![Refer to caption](https://arxiv.org/html/2605.24041v1/x4.png)\(c\)ERA5\(RFNE↓\\downarrow\)
![Refer to caption](https://arxiv.org/html/2605.24041v1/x5.png)\(d\)ERA5\(ACC↑\\uparrow\)

Figure 4:Convergence behavior across physical systems\.VRMSE for TR\-2D and Active Matter, and RFNE and ACC for ERA5, plotted as a function of refinement stepkk\. Vertical dashed lines indicate the training cutoff \(FNO atK=6K=6and TFNO/WDSR atK=4K=4\)\. Error metrics decrease or reach a stable plateau fork\>Kk\>K, while ACC increases and stabilizes\.To characterize the spectral dynamics of IRNO, we analyze how error evolves in the frequency domain across refinement steps at both the dataset and instance levels\.

#### 4\.3\.1Dataset\- and Instance\-Level Spectral Error Distribution

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/average_spectral_error.png)\(a\)Dataset\-level median normalized spectral error ratiosE~\(k\)​\(ω\)\\tilde\{E\}^\{\(k\)\}\(\\omega\)\.
![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/per_sample.png)\(b\)Instance\-level spectral MSE trajectories\.

Figure 5:Spectral Error Evolution under iterative refinement \(Active Matter, FNO\)\.\(a\)Median normalized spectral error ratiosE~\(k\)​\(ω\)\\tilde\{E\}^\{\(k\)\}\(\\omega\)across the test set, with shaded interquartile ranges \(25–75%\)\. IRNO exhibits consistent attenuation of mid\-to\-high frequency error with increasing refinement steps, with stable behavior near the Nyquist limitω=128\\omega=128, indicated by a vertical dashed line\.\(b\)Spectral Mean Squared Error \(MSE\) for a representative test instance fork∈\[0,12\]k\\in\[0,12\]refinement steps\. Each curve shows the radial\-spectral error of the base operator and successive IRNO refinements, illustrating monotonic reduction of spectral error and stable behavior beyond the training cutoff\. Across both panels, mid\-to\-high frequency error decreases with refinement, with the largest relative reductions observed near the Nyquist limitω=128\\omega=128\(vertical dashed line\)\.For each refinement stepkkand radial frequencyω\\omega, we compute the normalized spectral error ratio

E~\(k\)​\(ω\)=E\(k\)​\(ω\)E\(0\)​\(ω\)\+ε,\\tilde\{E\}^\{\(k\)\}\(\\omega\)=\\frac\{E^\{\(k\)\}\(\\omega\)\}\{E^\{\(0\)\}\(\\omega\)\+\\varepsilon\},whereE\(0\)​\(ω\)E^\{\(0\)\}\(\\omega\)denotes the spectral error of the base operator andε=10−10\\varepsilon=10^\{\-10\}is a small constant for numerical stability\. We report the median and interquartile range \(25–75%\) across the test set\.

Figure[5](https://arxiv.org/html/2605.24041#S4.F5)\(a\) shows that refinement reduces error primarily in mid\-to\-high frequency regimes where base operators exhibit strongest spectral bias\. We further quantify this observation in Section[4\.4](https://arxiv.org/html/2605.24041#S4.SS4)\. Figure[5](https://arxiv.org/html/2605.24041#S4.F5)\(b\) confirms monotonic spectral error reduction per\-sample with largest improvements nearω=128\\omega=128\(Nyquist limit\)\.

![Refer to caption](https://arxiv.org/html/2605.24041v1/x6.png)\(a\)Low\-frequency band
![Refer to caption](https://arxiv.org/html/2605.24041v1/x7.png)\(b\)Mid\-frequency band
![Refer to caption](https://arxiv.org/html/2605.24041v1/x8.png)\(c\)High\-frequency band
![Refer to caption](https://arxiv.org/html/2605.24041v1/x9.png)

Figure 6:Frequency\-band error ratios across refinement steps \(Active Matter, FNO\)\.Normalized error ratiosRk\(B\)R\_\{k\}^\{\(B\)\}for each tensor field, computed separately over \(a\) low\-, \(b\) mid\-, and \(c\) high\-frequency bands as a function of refinement stepkk\. The vertical dashed line marks the training cutoff atK=6K=6\. Lower values indicate stronger error reduction relative to the base operator\.

### 4\.4Frequency\-Band Convergence

To summarize frequency\-dependent convergence, we partition the spectrum into low\-, mid\-, and high\-frequency bands using resolution\-invariant percentile cutoffs \(bottom 33%, middle 33%, and top 33%, respectively\)\. For each bandBBand refinement stepkk, we compute the normalized error ratio aggregated over frequencies in the band and across the test set\.

LetΩB\\Omega\_\{B\}denote the set of frequencies in bandBB,i∈𝒟t​e​s​ti\\in\\mathcal\{D\}\_\{test\}index test samples, andEi\(k\)​\(ω\)E\_\{i\}^\{\(k\)\}\(\\omega\)denote the spectral error of sampleiiat frequencyω\\omegaand refinement stepkk\. We define

Rk\(B\)=mediani∈𝒟t​e​s​t​\(∑ω∈ΩBEi\(k\)​\(ω\)∑ω∈ΩBEi\(0\)​\(ω\)\+ε\)\.R\_\{k\}^\{\(B\)\}=\\mathrm\{median\}\_\{i\\in\\mathcal\{D\}\_\{test\}\}\\left\(\\frac\{\\sum\_\{\\omega\\in\\Omega\_\{B\}\}E\_\{i\}^\{\(k\)\}\(\\omega\)\}\{\\sum\_\{\\omega\\in\\Omega\_\{B\}\}E\_\{i\}^\{\(0\)\}\(\\omega\)\+\\varepsilon\}\\right\)\.
Figure[6](https://arxiv.org/html/2605.24041#S4.F6)reportsRk\(B\)R\_\{k\}^\{\(B\)\}across refinement steps for orientation and strain\-rate tensor fields\. Table[3](https://arxiv.org/html/2605.24041#S4.T3)summarizes averaged values at the initial, cutoff, and final iterations\. Across tensor fields, IRNO achieves 98–99% error reduction in high\-frequency bands, 93–96% in mid\-frequency bands, and 62–73% in low\-frequency bands at the final refinement step\.

Table 3:Frequency\-band error ratios across refinement steps \(Active Matter, FNO\)\.Normalized error ratiosRk\(B\)R\_\{k\}^\{\(B\)\}for low\-, mid\-, and high\-frequency bands, reported for orientation and strain tensor fields at the initial step \(k=1k=1\), training cutoff \(k=6k=6\), and final step \(k=12k=12\)\. Lower values indicate stronger error reduction\.
### 4\.5Transferability Across Base Operators

The function\-space formulation of IRNO in Section[2\.1](https://arxiv.org/html/2605.24041#S2.SS1)suggests that refinement dynamics depend primarily on the local residual geometry rather than the specific architecture producing the initialization\. Under Theorem[3\.1](https://arxiv.org/html/2605.24041#S3.Thmtheorem1), convergence is guaranteed whenever the initialization lies within the basin of attraction \(Assumption[3](https://arxiv.org/html/2605.24041#S3.I1.i3)\), and the operator is well\-conditioned\.

As shown in Table[4](https://arxiv.org/html/2605.24041#S4.T4), IRNO trained on one base operator successfully transfers without retraining to refine predictions from different operators, achieving consistent error reduction across operators and tasks\. Notably, IRNO trained on low\-performing base operators often outperforms same\-operator configurations when transferred to high\-performing base operators\. For example, IRNOTFNO\{\}\_\{\\text\{TFNO\}\}improves FNO’s initial prediction VRMSE by 58\.53% on TR\-2D, 13\.21 percentage points higher than IRNOFNO\{\}\_\{\\text\{FNO\}\}’s improvement\. This suggests that less accurate operators generate larger, more diverse residual structures during training, forcing IRNO to learn more robust error\-correction strategies that generalize effectively to the smaller, more structured residuals of higher\-performing operators\.

Table 4:Transferability across base operators\.IRNO trained on one base operator successfully transfers to refine predictions from different base operators\. Subscripts indicate the operator IRNO was trained with \(e\.g\., IRNOTFNO\{\}\_\{\\text\{TFNO\}\}was trained with TFNO as the base\)\. Standard errors are<0\.001<0\.001for all entries\.
### 4\.6Generalization to Graph\-Based Operators and Irregular Meshes

To evaluate whether IRNO extends beyond structured\-grid operators, we apply it to RIGNO\(Mousaviet al\.,[2025](https://arxiv.org/html/2605.24041#bib.bib59)\), a graph\-based neural operator that operates on irregular unstructured meshes\. We use the CE\-Gauss benchmark, an irregular unstructured mesh dataset with 16,384 nodes and 4 physical variables, and evaluate autoregressive rollout over 7 timesteps withK=4K=4refinement steps andα=0\.3\\alpha=0\.3\. Dataset details are provided in Appendix[G\.1](https://arxiv.org/html/2605.24041#A7.SS1)\.

Table 5:Autoregressive rollout on irregular meshes\.IRNO \(K=4K=4\) with RIGNO as the base operator on the CE\-Gauss unstructured mesh benchmark, evaluated over a 7\-step autoregressive rollout\. Metric reported is relativeL2L^\{2\}error \(%\) against ground truth\. Lower is better\.As shown in Table[5](https://arxiv.org/html/2605.24041#S4.T5), IRNO reducesL2L^\{2\}error at every timestep, with improvements compounding from 12\.5% att=1t=1to 21\.3% att=7t=7\. This suggests that early refinement suppresses error accumulation in autoregressive rollout\. Implementation details are provided in Appendix[G](https://arxiv.org/html/2605.24041#A7)\.

### 4\.7Ablation Study

#### 4\.7\.1Progressive vs\. Fixed Spectral Supervision

We ablate the proposed progressive spectral loss to evaluate the role of frequency curriculum in stable multiscale refinement\. Table[6](https://arxiv.org/html/2605.24041#S4.T6)compares models trained with a linearly increasing spectral exponentλk:1\.0→2\.0\\lambda\_\{k\}:1\.0\\rightarrow 2\.0against fixed\-weight baselines with constantλ∈\{1\.0,1\.25,1\.75,2\.0\}\\lambda\\in\\\{1\.0,1\.25,1\.75,2\.0\\\}\.

The progressive spectral loss schedule achieves VRMSE of 0\.039 versus 0\.051\-0\.070 for fixedλ\\lambda\(23\.5\-44\.3% reduction\), with particularly strong high\-frequency improvements \(normalized error 0\.24 vs 0\.60\-0\.88\)\.

Table 6:Ablation of progressive spectral loss \(Active Matter, TFNO\)\.Comparison of a progressive schedule \(λk:1→2\\lambda\_\{k\}:1\\rightarrow 2\) against fixed spectral weights, with TFNO \(4 steps\) as base operator\. We report VRMSE and frequency\-band normalized error ratios for low\-, mid\-, and high\-frequency bands\. Lower values indicate stronger error reduction\.
#### 4\.7\.2Robustness Across Refinement Architectures and Normalization

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/ablation_alpha.png)Figure 7:Ablation: Step sizeα\\alphasensitivity \(TR\-2D, TFNO\)\.Convergence behavior forα∈\{0\.1,0\.2,0\.3,0\.4,0\.5,0\.6\}\\alpha\\in\\\{0\.1,0\.2,0\.3,0\.4,0\.5,0\.6\\\}, trained withK=4K=4steps \(dashed line\)\. Small step sizes converge slowly; moderate step sizes \(α∈\[0\.2,0\.4\]\\alpha\\in\[0\.2,0\.4\]\) achieve optimal balance;α=0\.6\\alpha=0\.6diverges rapidly beyond the training horizon, violating the contraction conditionq=‖I−α​A​\(x\)‖op<1q=\\\|I\-\\alpha A\(x\)\\\|\_\{\\text\{op\}\}<1from Theorem[3\.1](https://arxiv.org/html/2605.24041#S3.Thmtheorem1)\.IRNO’s iterative mechanism is robust across refinement backbone choices, with ResNet, ConvNext, and FNO backbones all achieving\>71%\>71\\%VRMSE reduction on Active Matter \(TFNO base,K=4K=4\); see Appendix[H](https://arxiv.org/html/2605.24041#A8)\(Table[18](https://arxiv.org/html/2605.24041#A8.T18)\)\. Normalization choice does not qualitatively affect IRNO’s behavior\. BatchNorm, LayerNorm, and GroupNorm all yield consistent error reduction on TR\-2D \(TFNO base\); see Appendix[H](https://arxiv.org/html/2605.24041#A8)\(Table[19](https://arxiv.org/html/2605.24041#A8.T19)\)\.

#### 4\.7\.3Varying Step Sizeα\\alpha

Figure[7](https://arxiv.org/html/2605.24041#S4.F7)shows how step sizeα\\alphaaffects convergence on TR\-2D with TFNO\-based IRNO trained atk=4k=4\. Small value \(α=0\.1\\alpha=0\.1\) converges slowly but stably, while moderate values \(α=0\.2\\alpha=0\.2\) achieve optimal convergence speed and stable extrapolation beyond the training horizon\. Moderately large step size \(0\.3≤α≤0\.50\.3\\leq\\alpha\\leq 0\.5\) initially decreases error but diverges fork\>6k\>6, and large value \(α=0\.6\\alpha=0\.6\) diverges within training horizon with VRMSE increasing from 0\.12 to 0\.24\. This instability occurs whenα\\alphamoves outside the local contraction regime depicted in the theory section[3](https://arxiv.org/html/2605.24041#S3)\. We useα∈\{0\.2,0\.25\}\\alpha\\in\\\{0\.2,0\.25\\\}for our experiments, balancing convergence speed and stability\.

### 4\.8Cost–Performance Pareto Frontier

Finally, we evaluate whether iterative refinement provides a more efficient cost\-performance trade\-off than increasing model capacity\. Figure[8](https://arxiv.org/html/2605.24041#S4.F8)compares IRNO against a capacity\-matched monolithic residual model in terms of ACC and RFNE versus FLOPs and memory consumption\.

Across all metrics, IRNO traces a Pareto\-dominant frontier\. At 1100 GFLOPs, IRNO attains ACC = 0\.84, while a 15×\\timesU\-Net achieves ACC = 0\.79 at comparable computational cost\. Similarly, IRNO reduces RFNE to 0\.28 at 1000 GFLOPs, whereas monolithic baselines remain above 0\.295 in the same cost range\. These results indicate that performance gains are driven by the refinement mechanism itself rather than by increased model capacity\.

We additionally compare against F\-Adapter\(Zhanget al\.,[2025](https://arxiv.org/html/2605.24041#bib.bib58)\), a parameter\-efficient spectral fine\-tuning method\. F\-Adapter achieves a2\.31%2\.31\\%VRMSE reduction at low computational overhead, with gains concentrated in mid\- and high\-frequency bands; IRNO trades additional training cost for substantially larger gains \(50\.73%50\.73\\%\), representing complementary design regimes\. A full spectral comparison is provided in Appendix[H](https://arxiv.org/html/2605.24041#A8)\(Table[20](https://arxiv.org/html/2605.24041#A8.T20), Figure[20](https://arxiv.org/html/2605.24041#A8.F20)\)\.

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/acc-flops.png)\(a\)ACC vs\. FLOPs
![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/acc-memory.png)\(b\)ACC vs\. Memory
![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/rfne-flops.png)\(c\)RFNE vs\. FLOPs
![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/rfne-memory.png)\(d\)RFNE vs\. Memory

Figure 8:Cost–performance trade\-offs under iterative refinement \(ERA5, FNO\)\.ACC and RFNE plotted against computational cost \(FLOPs\) and memory usage during inference\. Each marker corresponds to a refinement step of the iterative U\-Net, with step indices annotated\. Baselines \(FNO, SRCNN, and 15×\\timesU\-Net\) are shown for reference\. IRNO achieves improved ACC or lower error at comparable or lower computational cost\. See Appendix[H](https://arxiv.org/html/2605.24041#A8)for detailed comparisons and wall\-clock times\.

## 5Conclusion

We propose IRNO, a learned iterative refinement operator that progressively reduces prediction error across physical systems\. We provide a contraction\-based theoretical analysis that interprets refinement as a fixed\-point iteration in function space, establishing conditions under which monotonic convergence and stability beyond the training horizon are expected\. Experimental results validate that IRNO achieves consistent spectral error reduction, stable extrapolation to iteration counts beyond those seen in training, and improved cost–performance trade\-offs relative to capacity\-matched baselines\.

## 6Acknowledgement

We thank our colleagues and funding agencies\. This work is supported by the DARPA AIQ program, the U\.S\. Department of Energy under Award Number DE\-SC0025584, and Dartmouth College\.

## References

- A\. S\. Almgren, J\. B\. Bell, M\. J\. Lijewski, Z\. Lukić, and E\. Van Andel \(2013\)Nyx: a massively parallel amr code for computational cosmology\.The Astrophysical Journal765\(1\),pp\. 39\.External Links:ISSN 1538\-4357,[Link](http://dx.doi.org/10.1088/0004-637X/765/1/39),[Document](https://dx.doi.org/10.1088/0004-637x/765/1/39)Cited by:[§1](https://arxiv.org/html/2605.24041#S1.p3.1)\.
- S\. Bai, J\. Z\. Kolter, and V\. Koltun \(2019\)Deep equilibrium models\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Cited by:[§B\.2](https://arxiv.org/html/2605.24041#A2.SS2.p2.1)\.
- R\. T\. Q\. Chen, Y\. Rubanova, J\. Bettencourt, and D\. Duvenaud \(2018\)Neural ordinary differential equations\.Advances in Neural Information Processing Systems\.Cited by:[§B\.2](https://arxiv.org/html/2605.24041#A2.SS2.p2.1)\.
- W\. Chen, J\. Song, P\. Ren, S\. Subramanian, D\. Morozov, and M\. W\. Mahoney \(2024\)Data\-efficient operator learning via unsupervised pretraining and in\-context learning\.InAdvances in Neural Information Processing Systems,A\. Globerson, L\. Mackey, D\. Belgrave, A\. Fan, U\. Paquet, J\. Tomczak, and C\. Zhang \(Eds\.\),Vol\.37,pp\. 6213–6245\.External Links:[Document](https://dx.doi.org/10.52202/079017-0201)Cited by:[§4\.1](https://arxiv.org/html/2605.24041#S4.SS1.SSS0.Px1.p1.1)\.
- D\. B\. Fielding, E\. C\. Ostriker, G\. L\. Bryan, and A\. S\. Jermyn \(2020\)Multiphase gas and the fractal nature of radiative turbulent mixing layers\.The Astrophysical Journal Letters894\(2\),pp\. L24\.Cited by:[§G\.1\.1](https://arxiv.org/html/2605.24041#A7.SS1.SSS1.Px1.p1.1)\.
- E\. Haber and L\. Ruthotto \(2017\)Stable architectures for deep neural networks\.Inverse Problems34\(1\),pp\. 014004\.External Links:[Document](https://dx.doi.org/10.1088/1361-6420/aa9a90),[Link](https://doi.org/10.1088/1361-6420/aa9a90)Cited by:[§B\.2](https://arxiv.org/html/2605.24041#A2.SS2.p1.1)\.
- S\. M\. S\. Hassan, A\. Feeney, A\. Dhruv, J\. Kim, Y\. Suh, J\. Ryu, Y\. Won, and A\. Chandramowlishwaran \(2023\)BubbleML: a multiphase multiphysics dataset and benchmarks for machine learning\.InAdvances in Neural Information Processing Systems,A\. Oh, T\. Naumann, A\. Globerson, K\. Saenko, M\. Hardt, and S\. Levine \(Eds\.\),Vol\.36,pp\. 418–449\.Cited by:[§4\.1](https://arxiv.org/html/2605.24041#S4.SS1.SSS0.Px1.p1.1)\.
- K\. He, X\. Zhang, S\. Ren, and J\. Sun \(2016\)Deep Residual Learning for Image Recognition\.InProceedings of 2016 IEEE Conference on Computer Vision and Pattern Recognition,CVPR ’16,pp\. 770–778\.External Links:[Document](https://dx.doi.org/10.1109/CVPR.2016.90),ISSN 1063\-6919,[Link](http://ieeexplore.ieee.org/document/7780459)Cited by:[§B\.2](https://arxiv.org/html/2605.24041#A2.SS2.p1.1)\.
- H\. Hersbach, B\. Bell, P\. Berrisford, S\. Hirahara, A\. Horányi, J\. Muñoz\-Sabater, J\. Nicolas, C\. Peubey, R\. Radu, D\. Schepers,et al\.\(2020\)The era5 global reanalysis\.Quarterly journal of the royal meteorological society146\(730\),pp\. 1999–2049\.Cited by:[§G\.1\.3](https://arxiv.org/html/2605.24041#A7.SS1.SSS3.p1.1),[§1](https://arxiv.org/html/2605.24041#S1.p5.1)\.
- B\. Holzschuh, Q\. Liu, G\. Kohl, and N\. Thuerey \(2025\)PDE\-transformer: efficient and versatile transformers for physics simulations\.InForty\-second International Conference on Machine Learning,External Links:[Link](https://openreview.net/forum?id=3BaJMRaPSx)Cited by:[§4\.1](https://arxiv.org/html/2605.24041#S4.SS1.SSS0.Px1.p1.1)\.
- S\. Khodakarami, V\. Oommen, A\. Bora, and G\. E\. Karniadakis \(2026\)Mitigating spectral bias in neural operators via high\-frequency scaling for physical systems\.Neural Networks193,pp\. 108027\.External Links:ISSN 0893\-6080,[Document](https://dx.doi.org/https%3A//doi.org/10.1016/j.neunet.2025.108027),[Link](https://www.sciencedirect.com/science/article/pii/S0893608025009074)Cited by:[§B\.1](https://arxiv.org/html/2605.24041#A2.SS1.p1.1),[§G\.3](https://arxiv.org/html/2605.24041#A7.SS3.SSS0.Px3.p1.1),[§4\.1](https://arxiv.org/html/2605.24041#S4.SS1.SSS0.Px1.p4.1),[§4\.2](https://arxiv.org/html/2605.24041#S4.SS2.p3.4),[Table 2](https://arxiv.org/html/2605.24041#S4.T2.8.3.1.1)\.
- J\. Kossaifi, N\. Kovachki, Z\. Li, D\. Pitt, M\. Liu\-Schiaffini, V\. Duruisseaux, R\. J\. George, B\. Bonev, K\. Azizzadenesheli, J\. Berner, and A\. Anandkumar \(2025\)A library for learning neural operators\.arXiv preprint arXiv:2412\.10354\.Cited by:[§4\.1](https://arxiv.org/html/2605.24041#S4.SS1.SSS0.Px1.p3.1)\.
- N\. B\. Kovachki, Z\. Li, B\. Liu, K\. Azizzadenesheli, K\. Bhattacharya, A\. M\. Stuart, and A\. Anandkumar \(2023a\)Neural operator: learning maps between function spaces\.Journal of Machine Learning Research24\(89\),pp\. 1–97\.External Links:[Link](https://www.jmlr.org/papers/v24/21-1524.html)Cited by:[item 2](https://arxiv.org/html/2605.24041#S3.I1.i2.p3.3)\.
- N\. Kovachki, Z\. Li, B\. Liu, K\. Azizzadenesheli, K\. Bhattacharya, A\. Stuart, and A\. Anandkumar \(2023b\)Neural operator: learning maps between function spaces\.Journal of Machine Learning Research24\(89\),pp\. 1–97\.Cited by:[§G\.2\.1](https://arxiv.org/html/2605.24041#A7.SS2.SSS1.Px2.p1.1)\.
- Z\. Li, N\. Kovachki, K\. Azizzadenesheli, B\. Liu, K\. Bhattacharya, A\. Stuart, and A\. Anandkumar \(2021\)Fourier neural operator for parametric partial differential equations\.InInternational Conference on Learning Representations \(ICLR\),Note:See also arXiv:2010\.08895 \(2020\)External Links:[Link](https://openreview.net/forum?id=c8P9NQVtmnO)Cited by:[§B\.2](https://arxiv.org/html/2605.24041#A2.SS2.p1.1),[§G\.2\.1](https://arxiv.org/html/2605.24041#A7.SS2.SSS1.Px1.p1.1),[§1](https://arxiv.org/html/2605.24041#S1.p1.1),[§4\.1](https://arxiv.org/html/2605.24041#S4.SS1.SSS0.Px1.p3.1)\.
- X\. Liu, B\. Xu, S\. Cao, and L\. Zhang \(2024\)Mitigating spectral bias for the multiscale operator learning\.Journal of Computational Physics506,pp\. 112944\.External Links:ISSN 0021\-9991,[Document](https://dx.doi.org/https%3A//doi.org/10.1016/j.jcp.2024.112944),[Link](https://www.sciencedirect.com/science/article/pii/S0021999124001931)Cited by:[§B\.1](https://arxiv.org/html/2605.24041#A2.SS1.p1.1)\.
- L\. Lu, P\. Jin, G\. Pang, Z\. Zhang, and G\. E\. Karniadakis \(2021\)Learning nonlinear operators via deeponet based on the universal approximation theorem of operators\.Nature Machine Intelligence3\(3\),pp\. 218–229\.External Links:[Document](https://dx.doi.org/10.1038/s42256-021-00302-5)Cited by:[§B\.2](https://arxiv.org/html/2605.24041#A2.SS2.p1.1),[§1](https://arxiv.org/html/2605.24041#S1.p1.1)\.
- X\. Luo, X\. Qian, and B\. Yoon \(2024\)Hierarchical neural operator transformer with learnable frequency\-aware loss prior for arbitrary\-scale super\-resolution\.arXiv preprint arXiv:2405\.12202\.Cited by:[§G\.3](https://arxiv.org/html/2605.24041#A7.SS3.SSS0.Px3.p1.1),[§4\.1](https://arxiv.org/html/2605.24041#S4.SS1.SSS0.Px1.p4.1),[§4\.2](https://arxiv.org/html/2605.24041#S4.SS2.p3.4),[Table 2](https://arxiv.org/html/2605.24041#S4.T2.8.4.2.1)\.
- R\. Morel, J\. Han, and E\. Oyallon \(2025\)DISCO: learning to DISCover an evolution operator for multi\-physics\-agnostic prediction\.InForty\-second International Conference on Machine Learning,External Links:[Link](https://openreview.net/forum?id=6EZ3MDDf6p)Cited by:[§4\.1](https://arxiv.org/html/2605.24041#S4.SS1.SSS0.Px1.p1.1)\.
- S\. Mousavi, S\. Wen, L\. Lingsch, M\. Herde, B\. Raonić, and S\. Mishra \(2025\)RIGNO: a graph\-based framework for robust and accurate operator learning for pdes on arbitrary domains\.InAdvances in Neural Information Processing Systems,Vol\.38\.Cited by:[§G\.1\.4](https://arxiv.org/html/2605.24041#A7.SS1.SSS4.Px1.p1.1),[§4\.1](https://arxiv.org/html/2605.24041#S4.SS1.SSS0.Px1.p1.1),[§4\.6](https://arxiv.org/html/2605.24041#S4.SS6.p1.2)\.
- R\. Ohana, M\. McCabe, L\. Meyer, R\. Morel, F\. Agocs, M\. Beneitez, M\. Berger, B\. Burkhart, S\. Dalziel, D\. Fielding,et al\.\(2024\)The well: a large\-scale collection of diverse physics simulations for machine learning\.Advances in Neural Information Processing Systems37,pp\. 44989–45037\.Cited by:[§G\.5\.1](https://arxiv.org/html/2605.24041#A7.SS5.SSS1.Px1.p1.3),[§4\.1](https://arxiv.org/html/2605.24041#S4.SS1.SSS0.Px1.p1.1)\.
- V\. Oommen, A\. Bora, Z\. Zhang, and G\. E\. Karniadakis \(2025\)Integrating neural operators with diffusion models improves spectral representation in turbulence modelling\.Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences481\(2309\),pp\. 20240819\.External Links:ISSN 1364\-5021,[Document](https://dx.doi.org/10.1098/rspa.2024.0819),[Link](https://doi.org/10.1098/rspa.2024.0819),https://royalsocietypublishing\.org/rspa/article\-pdf/doi/10\.1098/rspa\.2024\.0819/513944/rspa\.2024\.0819\.pdfCited by:[§B\.1](https://arxiv.org/html/2605.24041#A2.SS1.p2.1)\.
- S\. Qin, F\. Lyu, W\. Peng, D\. Geng, J\. Wang, X\. Tang, S\. Leroyer, N\. Gao, X\. Liu, and L\. L\. Wang \(2024\)Toward a better understanding of fourier neural operators from a spectral perspective\.External Links:2404\.07200,[Link](https://arxiv.org/abs/2404.07200)Cited by:[§B\.1](https://arxiv.org/html/2605.24041#A2.SS1.p1.1)\.
- N\. Rahaman, A\. Baratin, D\. Arpit, F\. Draxler, M\. Lin, F\. A\. Hamprecht, Y\. Bengio, and A\. Courville \(2019\)On the spectral bias of neural networks\.InInternational Conference on Machine Learning \(ICML\),Proceedings of Machine Learning Research, Vol\.97,pp\. 5301–5310\.External Links:[Link](https://proceedings.mlr.press/v97/rahaman19a.html)Cited by:[§B\.1](https://arxiv.org/html/2605.24041#A2.SS1.p1.1)\.
- A\. Rastogi, G\. Blanchard, and P\. Mathé \(2020\)Convergence analysis of Tikhonov regularization for non\-linear statistical inverse problems\.Electronic Journal of Statistics14\(2\),pp\. 2798 – 2841\.External Links:[Document](https://dx.doi.org/10.1214/20-EJS1735),[Link](https://doi.org/10.1214/20-EJS1735)Cited by:[item 2](https://arxiv.org/html/2605.24041#S3.I1.i2.p3.3)\.
- P\. Ren, N\. B\. Erichson, J\. Guo, S\. Subramanian, O\. San, Z\. Lukic, and M\. W\. Mahoney \(2025\)SuperBench: a super\-resolution benchmark dataset for scientific machine learning\.Journal of Data\-centric Machine Learning Research\.Note:Dataset Certification, Reproducibility CertificationExternal Links:[Link](https://openreview.net/forum?id=OJ6zUcWldW)Cited by:[§G\.1\.3](https://arxiv.org/html/2605.24041#A7.SS1.SSS3.p1.1),[§1](https://arxiv.org/html/2605.24041#S1.p5.1),[§4\.1](https://arxiv.org/html/2605.24041#S4.SS1.SSS0.Px1.p1.1)\.
- O\. Ronneberger, P\. Fischer, and T\. Brox \(2015\)U\-net: convolutional networks for biomedical image segmentation\.InMedical Image Computing and Computer\-Assisted Intervention – MICCAI 2015,N\. Navab, J\. Hornegger, W\. M\. Wells, and A\. F\. Frangi \(Eds\.\),Cham,pp\. 234–241\.External Links:ISBN 978\-3\-319\-24574\-4Cited by:[§2\.1\.1](https://arxiv.org/html/2605.24041#S2.SS1.SSS1.Px1.p2.1)\.
- W\. Wang and F\. Zhang \(2025\)Breaking scale anchoring: frequency representation learning for accurate high\-resolution inference from low\-resolution training\.External Links:2512\.05132,[Link](https://arxiv.org/abs/2512.05132)Cited by:[§B\.1](https://arxiv.org/html/2605.24041#A2.SS1.p2.1)\.
- F\. Werner and B\. Hofmann \(2019\)Convergence analysis of \(statistical\) inverse problems under conditional stability estimates\.Inverse Problems36\(1\),pp\. 015004\.External Links:[Document](https://dx.doi.org/10.1088/1361-6420/ab4cd7),[Link](https://doi.org/10.1088/1361-6420/ab4cd7)Cited by:[item 2](https://arxiv.org/html/2605.24041#S3.I1.i2.p3.3)\.
- R\. Wienands and W\. Joppich \(2005\)Practical fourier analysis for multigrid methods\.Numerical Insights,Taylor & Francis\.External Links:ISBN 9781584884927,LCCN 2004055113,[Link](https://books.google.com/books?id=IOSux5GxacsC)Cited by:[§B\.2](https://arxiv.org/html/2605.24041#A2.SS2.p1.1)\.
- Z\. Wu, S\. Zhang, S\. He, S\. Wang, M\. Zhu, A\. Jiao, L\. Lu, and D\. van Dijk \(2025\)COAST: intelligent time\-adaptive neural operators\.In2nd AI for Math Workshop @ ICML 2025,External Links:[Link](https://openreview.net/forum?id=cz99ALd3yv)Cited by:[§4\.1](https://arxiv.org/html/2605.24041#S4.SS1.SSS0.Px1.p1.1)\.
- J\. Yu, Y\. Fan, J\. Yang, N\. Xu, Z\. Wang, X\. Wang, and T\. Huang \(2018\)Wide activation for efficient and accurate image super\-resolution\.arXiv preprint arXiv:1808\.08718\.Cited by:[§G\.2\.1](https://arxiv.org/html/2605.24041#A7.SS2.SSS1.Px3.p1.1),[§4\.1](https://arxiv.org/html/2605.24041#S4.SS1.SSS0.Px1.p3.1)\.
- H\. Zhang, C\. Kang, Y\. Wang, and D\. Zou \(2025\)F\-adapter: frequency\-adaptive parameter\-efficient fine\-tuning in scientific machine learning\.arXiv preprint arXiv:2509\.23173\.Cited by:[§B\.1](https://arxiv.org/html/2605.24041#A2.SS1.p3.1),[Appendix H](https://arxiv.org/html/2605.24041#A8.SS0.SSS0.Px10.p1.2),[§4\.8](https://arxiv.org/html/2605.24041#S4.SS8.p3.2)\.

## Appendix ALimitations and Future Work

The proposed Iterative Refinement Neural Operator \(IRNO\) demonstrates stable convergence, reduced spectral error, and improved cost–performance trade\-offs across multiple physical systems\. Several open questions remain that point toward productive directions for future research\.

Beyond Local Neighborhoods\.The present analysis establishes convergence under smoothness and spectral radius assumptions in a neighborhood of the true solution, depending on the quality of the initialization from the base operator\. Extending these results to broader regions of the solution space remains an open problem\. Tools from monotone operator theory, non\-expansive mappings, or Lyapunov stability analysis may provide a path toward characterizing global behavior and basin structure for learned refinement dynamics\.

Adaptive Control of Refinement Dynamics\.IRNO uses a fixed step size chosen empirically to satisfy a contraction condition derived from the local Jacobian\. This choice offers a balance between stability and convergence speed, though it does not incorporate information from the refinement trajectory itself\. A natural criterion for adaptive stopping is when‖Φθ​\(x,hk\)‖\\\|\\Phi\_\{\\theta\}\(x,h\_\{k\}\)\\\|falls below a threshold tied to the bias level, which avoids over\-refinement once the error floor has been reached\. Future work may explore adaptive step\-size policies or learned controllers that respond to local error geometry or spectral content, thereby drawing connections to classical line\-search and trust\-region methods\.

Extensions to Stochastic and Inverse Settings\.This work focuses on deterministic forward prediction tasks\. Extending iterative refinement to stochastic operators, uncertainty\-aware prediction, or inverse problems, such as data assimilation or parameter identification, represents a promising direction\. In these settings, refinement dynamics could be coupled with probabilistic inference mechanisms, potentially bridging learned solvers with Bayesian filtering or variational frameworks\.

Cross\-Resolution and Distribution Generalization\.Current evaluation uses matched train and test resolutions, and all benchmarks draw from the same data distribution\. Whether the learned update rule transfers to finer grids than those seen during training, or remains accurate under distributional shift \(e\.g\., out\-of\-distribution forcing terms or boundary conditions\), is an open question\. Studying resolution generalization and distributional robustness is a key direction for establishing the broader applicability of IRNO in scientific workflows\.

## Appendix BRelated Work

### B\.1Spectral Bias and Multi\-Scale Operator Learning

A fundamental challenge in training neural operators for SciML is the phenomenon of spectral bias, where neural networks prioritize learning low\-frequency components of the target function while struggling to capture high\-frequency detailsRahamanet al\.\[[2019](https://arxiv.org/html/2605.24041#bib.bib6)\]\. Early efforts to address this limitation in the context of operator learning focused on architectural innovations that explicitly route information across scales\. For instance, the Hierarchical Attention Neural Operator \(HANO\)Liuet al\.\[[2024](https://arxiv.org/html/2605.24041#bib.bib33)\]was proposed to enable nested feature computation that better captures multiscale solution spaces, by employing a scale\-adaptive interaction range and self\-attention mechanisms over a hierarchy of levels\. More recent approaches have sought to mitigate spectral bias through frequency\-specific boosting and scaling mechanisms\. SpecBoost was introduced as an ensemble framework designed to enhance high\-frequency capture in FNOsQinet al\.\[[2024](https://arxiv.org/html/2605.24041#bib.bib34)\]\. Their spectral analysis revealed that FNOs exhibit a distinct parameterization bias that favors dominant low frequencies\. To overcome this, SpecBoost trains a secondary residual operator specifically tasked with learning the high\-frequency residuals left by the primary model\. Similarly, High\-Frequency Scaling \(HFS\) modulates the latent space of convolutional\-based neural operators to amplify high\-frequency modesKhodakaramiet al\.\[[2026](https://arxiv.org/html/2605.24041#bib.bib32)\]\. Unlike Fourier\-based interventions, HFS operates directly in the latent representation, avoiding the computational overhead of the FFT while successfully recovering fine\-scale features in complex multiphase flow and turbulence problems\.

The integration of generative models has emerged as a powerful strategy for recovering high\-frequency spectral content, particularly in turbulent and multiscale systems\. Conditioning diffusion models on the output of neural operators enables the generation of fine\-scale fluctuations that are typically smoothed out by deterministic operator predictions, leading to improved alignment with the true energy spectrumOommenet al\.\[[2025](https://arxiv.org/html/2605.24041#bib.bib31)\]\. From a resolution\-generalization perspective, scale anchoring has been identified as a fundamental limitation in which models trained on low\-resolution data fail to extrapolate to higher Nyquist frequencies, and Frequency Representation Learning has been proposed as a means to align spectral features across resolutions and reduce this dependencyWang and Zhang \[[2025](https://arxiv.org/html/2605.24041#bib.bib30)\]\. Our work complements these approaches but differs fundamentally in its inference mechanism\. While these methods improve spectral coverage through enhanced architectures or training objectives, they retain a*single\-pass*inference structure where all frequency components must be resolved simultaneously\. In contrast, IRNO reframes operator inference as an*iterative refinement process*, progressively correcting residual errors across scales through learned fixed\-point iterations\. This enables systematic error reduction at inference time without increasing the base model’s capacity or retraining\.

A concurrent line of work, F\-Adapter\[Zhanget al\.,[2025](https://arxiv.org/html/2605.24041#bib.bib58)\]\(NeurIPS 2025\), addresses spectral bias through parameter\-efficient adaptation of Fourier layers using low\-rank decompositions\. Its theoretical grounding concerns the approximation capacity of LoRA\-style adapters in the Fourier domain\. IRNO’s theoretical grounding is distinct, concerning the dynamics of fixed\-point iteration and establishing contraction conditions for monotonic convergence\. Empirically, F\-Adapter targets minimal\-parameter adaptation \(2\.31% VRMSE reduction on Active Matter\) while IRNO trades additional training cost for substantially larger gains \(50\.73%\), representing complementary design regimes rather than competing approaches\. See Appendix[H](https://arxiv.org/html/2605.24041#A8)for a direct quantitative comparison\.

### B\.2Numerical Principles as Inductive Bias

Deep learning architectures often exhibit structural similarity with classical numerical methods\. Residual networks, for example, are mathematically equivalent to the forward Euler discretization of continuous\-time dynamical systemsHeet al\.\[[2016](https://arxiv.org/html/2605.24041#bib.bib38)\]\. Similarly, encoder\-decoder structures with skip connections, such as U\-Nets, mirror the restriction and prolongation operations in multigrid methods, enabling efficient error reduction across scalesHaber and Ruthotto \[[2017](https://arxiv.org/html/2605.24041#bib.bib39)\], Wienands and Joppich \[[2005](https://arxiv.org/html/2605.24041#bib.bib41)\]\. Neural operators such as FNOLiet al\.\[[2021](https://arxiv.org/html/2605.24041#bib.bib3)\]and DeepONetLuet al\.\[[2021](https://arxiv.org/html/2605.24041#bib.bib17)\]parameterize mappings in spectral or basis\-function spaces, consistent with classical spectral/pseudo\-spectral methods\. In most cases, these models are trained as*direct solvers*that approximate the solution operator through a single forward evaluation\.

A closer conceptual link to our framework arises in implicit layers, Deep Equilibrium Models \(DEQs\)Baiet al\.\[[2019](https://arxiv.org/html/2605.24041#bib.bib35)\], and Neural Ordinary Differential Equations \(Neural ODEs\)Chenet al\.\[[2018](https://arxiv.org/html/2605.24041#bib.bib36)\], which explicitly cast network inference as the solution of a fixed\-point or continuous\-time dynamical system\. These approaches define the*network itself*as an equilibrium or flow, with iterative solvers embedded into training and inference\. IRNO differs by preserving a standard, explicit base operator and introducing a learned*refinement dynamics*applied at inference time\. Opposed to learning an implicit representation, IRNO learns a residual update rule that iteratively corrects the output of a pre\-trained operator, retaining architectural modularity while inheriting the convergence and stability properties of classical fixed\-point iterations\.

## Appendix CProof of Theorem[3\.1](https://arxiv.org/html/2605.24041#S3.Thmtheorem1)\(Quadratic–Linear Convergence\)

We sketch the main steps establishing the error recursion and local contractivity\.

Letek=y−hke\_\{k\}=y\-h\_\{k\}\. From the iteration

hk\+1=hk\+α​Φ​\(x,hk\),h\_\{k\+1\}=h\_\{k\}\+\\alpha\\Phi\(x,h\_\{k\}\),and the local affine decomposition in Assumption[1](https://arxiv.org/html/2605.24041#S3.I1.i1),

Φ​\(x,hk\)=A​\(x,hk\)​ek\+b\+R​\(x,hk\),\\Phi\(x,h\_\{k\}\)=A\(x,h\_\{k\}\)e\_\{k\}\+b\+R\(x,h\_\{k\}\),we obtain

ek\+1=ek−α​A​\(x,hk\)​ek−α​b−α​R​\(x,hk\)\.e\_\{k\+1\}=e\_\{k\}\-\\alpha A\(x,h\_\{k\}\)e\_\{k\}\-\\alpha b\-\\alpha R\(x,h\_\{k\}\)\.DecomposingA​\(x,hk\)=A​\(x\)\+\(A​\(x,hk\)−A​\(x\)\)A\(x,h\_\{k\}\)=A\(x\)\+\\big\(A\(x,h\_\{k\}\)\-A\(x\)\\big\)yields

ek\+1=\(I−α​A​\(x\)\)​ek−α​\(A​\(x,hk\)−A​\(x\)\)​ek−α​b−α​R​\(x,hk\)\.e\_\{k\+1\}=\(I\-\\alpha A\(x\)\)e\_\{k\}\-\\alpha\\big\(A\(x,h\_\{k\}\)\-A\(x\)\\big\)e\_\{k\}\-\\alpha b\-\\alpha R\(x,h\_\{k\}\)\.
Taking norms and applying the Jacobian stability and remainder bounds from Assumptions[1](https://arxiv.org/html/2605.24041#S3.I1.i1)and[2](https://arxiv.org/html/2605.24041#S3.I1.i2), we obtain

‖ek\+1‖≤‖I−α​A​\(x\)‖op​‖ek‖\+α​μ​‖ek‖2\+α​L2​‖ek‖2\+α​‖b‖\.\\\|e\_\{k\+1\}\\\|\\leq\\\|I\-\\alpha A\(x\)\\\|\_\{\\mathrm\{op\}\}\\\|e\_\{k\}\\\|\+\\alpha\\mu\\\|e\_\{k\}\\\|^\{2\}\+\\alpha\\tfrac\{L\}\{2\}\\\|e\_\{k\}\\\|^\{2\}\+\\alpha\\\|b\\\|\.Definingq=‖I−α​A​\(x\)‖opq=\\\|I\-\\alpha A\(x\)\\\|\_\{\\mathrm\{op\}\}andc=α​\(L2\+μ\)c=\\alpha\\left\(\\tfrac\{L\}\{2\}\+\\mu\\right\)gives the claimed bound

‖ek\+1‖≤q​‖ek‖\+c​‖ek‖2\+α​‖b‖\.\\\|e\_\{k\+1\}\\\|\\leq q\\\|e\_\{k\}\\\|\+c\\\|e\_\{k\}\\\|^\{2\}\+\\alpha\\\|b\\\|\.
Whenb=0b=0and‖e0‖<1−q2​c\\\|e\_\{0\}\\\|<\\tfrac\{1\-q\}\{2c\}, we have

q\+c​‖ek‖<1,q\+c\\\|e\_\{k\}\\\|<1,which implies local contractivity\. By induction, the error decreases monotonically and the iteration converges toyy\. The resulting behavior is linear for small errors, with a quadratic correction term governing the transient regime\.

□\\square

The error recursion in Theorem[3\.1](https://arxiv.org/html/2605.24041#S3.Thmtheorem1)decomposes the refinement dynamics into three contributions with clear numerical interpretations\. The linear termq​‖ek‖q\\\|e\_\{k\}\\\|governs the convergence rate and is determined by the choice of step sizeα\\alphaand the spectral properties of the local Jacobian throughq=‖I−α​A​\(x\)‖opq=\\\|I\-\\alpha A\(x\)\\\|\_\{\\mathrm\{op\}\}\. The quadratic termc​‖ek‖2c\\\|e\_\{k\}\\\|^\{2\}captures higher\-order deviations from the local linearization and dominates in the transient regime, explaining the accelerated error reduction observed in early refinement steps\. Finally, the bias termα​‖b‖\\alpha\\\|b\\\|sets a limiting error floor in the presence of model mismatch, leading to convergence toward a neighborhood of the true solution as formalized in Corollary[3\.3](https://arxiv.org/html/2605.24041#S3.Thmtheorem3)\.

## Appendix DProofs of Corollaries

### Corollary[3\.2](https://arxiv.org/html/2605.24041#S3.Thmtheorem2)Geometric Convergence and Iteration Complexity

From Theorem[3\.1](https://arxiv.org/html/2605.24041#S3.Thmtheorem1)withb=0b=0, the error satisfies

‖ek\+1‖≤q​‖ek‖\+c​‖ek‖2\.\\\|e\_\{k\+1\}\\\|\\leq q\\\|e\_\{k\}\\\|\+c\\\|e\_\{k\}\\\|^\{2\}\.As‖ek‖→0\\\|e\_\{k\}\\\|\\to 0, the quadratic term becomes negligible\. For anyε\>0\\varepsilon\>0, there existsk0k\_\{0\}such that for allk≥k0k\\geq k\_\{0\},

c​‖ek‖<ε,c\\\|e\_\{k\}\\\|<\\varepsilon,and hence

‖ek\+1‖≤\(q\+ε\)​‖ek‖\.\\\|e\_\{k\+1\}\\\|\\leq\(q\+\\varepsilon\)\\\|e\_\{k\}\\\|\.Iterating this inequality yields

‖ek‖≤C​\(q\+ε\)k​‖e0‖,\\\|e\_\{k\}\\\|\\leq C\(q\+\\varepsilon\)^\{k\}\\\|e\_\{0\}\\\|,for some constantC\>0C\>0\. Sinceε\\varepsilonis arbitrary, this implies the asymptotic geometric rate‖ek‖≲qk​‖e0‖\\\|e\_\{k\}\\\|\\lesssim q^\{k\}\\\|e\_\{0\}\\\|\.

To achieve‖ek‖≤ε′\\\|e\_\{k\}\\\|\\leq\\varepsilon^\{\\prime\}, it suffices to choosekksuch that

qk​‖e0‖≤ε′,q^\{k\}\\\|e\_\{0\}\\\|\\leq\\varepsilon^\{\\prime\},which gives

k=O​\(log⁡\(‖e0‖/ε′\)log⁡\(1/q\)\)\.k=O\\\!\\left\(\\frac\{\\log\(\\\|e\_\{0\}\\\|/\\varepsilon^\{\\prime\}\)\}\{\\log\(1/q\)\}\\right\)\.□\\square

### Corollary[3\.3](https://arxiv.org/html/2605.24041#S3.Thmtheorem3)Convergence with Bias

Part 1: Existence and Uniqueness of Fixed Point, and Convergence\.

Define the operatorT:Bδ​\(y\)→ℋT:B\_\{\\delta\}\(y\)\\to\\mathcal\{H\}by

T​\(h\)=h\+α​Φ​\(x,h\)\.T\(h\)=h\+\\alpha\\Phi\(x,h\)\.We now verify explicitly thatTTmaps a ball aroundyyinto itself and is a contraction, so the Banach Fixed Point Theorem applies\. The derivation also reveals how the admissible radius depends onqqand the biasbb\.

Let a positiver<δr<\\deltabe chosen later and considerh∈B¯r​\(y\)h\\in\\overline\{B\}\_\{r\}\(y\)\. Using Assumption[1](https://arxiv.org/html/2605.24041#S3.I1.i1)gives

‖T​\(h\)−y‖\\displaystyle\\\|T\(h\)\-y\\\|=‖h−y\+α​\[b\+A​\(x,h\)​\(y−h\)\+R​\(x,h\)\]‖\\displaystyle=\\\|h\-y\+\\alpha\[b\+A\(x,h\)\(y\-h\)\+R\(x,h\)\]\\\|=‖\(I−α​A​\(x,h\)\)​\(h−y\)\+α​b\+α​R​\(x,h\)‖\.\\displaystyle=\\\|\(I\-\\alpha A\(x,h\)\)\(h\-y\)\+\\alpha b\+\\alpha R\(x,h\)\\\|\.Applying Assumptions[2](https://arxiv.org/html/2605.24041#S3.I1.i2)and triangular inequality, we have

‖T​\(h\)−y‖\\displaystyle\\\|T\(h\)\-y\\\|=‖\(I−α​A​\(x,h\)\)​\(h−y\)\+α​b\+α​R​\(x,h\)‖\\displaystyle=\\\|\(I\-\\alpha A\(x,h\)\)\(h\-y\)\+\\alpha b\+\\alpha R\(x,h\)\\\|≤‖I−α​A​\(x,h\)‖op​‖h−y‖\+α​‖b‖\+α​‖R​\(x,h\)‖\\displaystyle\\leq\\\|I\-\\alpha A\(x,h\)\\\|\_\{\\text\{op\}\}\\,\\\|h\-y\\\|\+\\alpha\\\|b\\\|\+\\alpha\\\|R\(x,h\)\\\|≤\(‖I−α​A​\(x\)‖op\+α​‖A​\(x,h\)−A​\(x\)‖op\)​‖h−y‖\+α​‖b‖\+α​‖R​\(x,h\)‖\\displaystyle\\leq\\big\(\\\|I\-\\alpha A\(x\)\\\|\_\{\\text\{op\}\}\+\\alpha\\\|A\(x,h\)\-A\(x\)\\\|\_\{\\text\{op\}\}\\big\)\\\|h\-y\\\|\+\\alpha\\\|b\\\|\+\\alpha\\\|R\(x,h\)\\\|≤\(q\+α​μ​r\)​r\+α​‖b‖\+α​L2​r2\\displaystyle\\leq\\big\(q\+\\alpha\\mu r\\big\)r\+\\alpha\\\|b\\\|\+\\alpha\\frac\{L\}\{2\}r^\{2\}=q​r\+α​\(μ\+L2\)​r2\+α​‖b‖\.\\displaystyle=qr\+\\alpha\\\!\\left\(\\mu\+\\frac\{L\}\{2\}\\right\)r^\{2\}\+\\alpha\\\|b\\\|\.Thus,TTmapsB¯r​\(y\)\\overline\{B\}\_\{r\}\(y\)into itself providedrrsatisfies

r≥q​r\+α​\(μ\+L2\)​r2\+α​‖b‖\.r\\;\\geq\\;qr\+\\alpha\\\!\\left\(\\mu\+\\tfrac\{L\}\{2\}\\right\)r^\{2\}\+\\alpha\\\|b\\\|\.This quadratic inequality inrradmits a positive solution wheneverq<1q<1and‖b‖\\\|b\\\|is sufficiently small to make the discriminant positive\. In particular,

r−,r\+=1−q±\(1−q\)2−4​α2​\(μ\+L2\)​‖b‖2​α​\(μ\+L2\),r\_\{\-\},r\_\{\+\}=\\frac\{1\-q\\pm\\sqrt\{\(1\-q\)^\{2\}\-4\\alpha^\{2\}\\\!\\left\(\\mu\+\\tfrac\{L\}\{2\}\\right\)\\\|b\\\|\}\}\{2\\alpha\\\!\\left\(\\mu\+\\tfrac\{L\}\{2\}\\right\)\},and operatorTTis self\-mapping ifr∈\[r−,r\+\]∪\(0,δ\)r\\in\[r\_\{\-\},r\_\{\+\}\]\\cup\(0,\\delta\)\.

For small enough‖b‖\\\|b\\\|, letε:=4​α2​\(μ\+L2\)​‖b‖\(1−q\)2\\varepsilon:=\\frac\{4\\alpha^\{2\}\\\!\\left\(\\mu\+\\tfrac\{L\}\{2\}\\right\)\\\|b\\\|\}\{\\\!\\left\(1\-q\\right\)^\{2\}\}\. Then

r−=\(1−q\)−\(1−q\)​1−ε2​α​\(μ\+L2\)=\(1−q\)​\(1−1\+12​ε\+O​\(ε2\)\)2​α​\(μ\+L2\)=α​‖b‖1−q\+O​\(‖b‖2\),r\_\{\-\}=\\frac\{\(1\-q\)\-\(1\-q\)\\sqrt\{1\-\\varepsilon\}\}\{2\\alpha\\\!\\left\(\\mu\+\\tfrac\{L\}\{2\}\\right\)\}=\\frac\{\(1\-q\)\\left\(1\-1\+\\frac\{1\}\{2\}\\varepsilon\+O\(\\varepsilon^\{2\}\)\\right\)\}\{2\\alpha\\\!\\left\(\\mu\+\\tfrac\{L\}\{2\}\\right\)\}=\\frac\{\\alpha\\\|b\\\|\}\{1\-q\}\+O\(\\\|b\\\|^\{2\}\),by taking the first\-order Taylor expansion of1−ε\\sqrt\{1\-\\varepsilon\}\.

Forh1,h2∈B¯r​\(y\)h\_\{1\},h\_\{2\}\\in\\bar\{B\}\_\{r\}\(y\), we bound the contraction factor ofTTusing the Fréchet differentiability ofΦ\\Phiand Assumption[1](https://arxiv.org/html/2605.24041#S3.I1.i1)\-[2](https://arxiv.org/html/2605.24041#S3.I1.i2)\.

By the fundamental theorem of calculus for Fréchet derivatives,

Φ​\(x,h1\)−Φ​\(x,h2\)=∫01Dh​Φ​\(x,h2\+t​\(h1−h2\)\)​\(h1−h2\)​𝑑t\.\\Phi\(x,h\_\{1\}\)\-\\Phi\(x,h\_\{2\}\)=\\int\_\{0\}^\{1\}D\_\{h\}\\Phi\\bigl\(x,h\_\{2\}\+t\(h\_\{1\}\-h\_\{2\}\)\\bigr\)\(h\_\{1\}\-h\_\{2\}\)\\,dt\.
Under Assumption[1](https://arxiv.org/html/2605.24041#S3.I1.i1), the derivativeDh​Φ​\(x,h\)D\_\{h\}\\Phi\(x,h\)exists and satisfies

Dh​Φ​\(x,h\)=−A​\(x,h\)\+E​\(x,h\),D\_\{h\}\\Phi\(x,h\)=\-A\(x,h\)\+E\(x,h\),where the remainder termE​\(x,h\)E\(x,h\)arises from the quadratic remainderRRand fulfills‖E​\(x,h\)‖op≤C​‖y−h‖\\\|E\(x,h\)\\\|\_\{\\mathrm\{op\}\}\\leq C\\\|y\-h\\\|for someC\>0C\>0\(sinceRRis quadratic in‖e‖\\\|e\\\|\)\.

Hence,

‖T​\(h1\)−T​\(h2\)‖\\displaystyle\\\|T\(h\_\{1\}\)\-T\(h\_\{2\}\)\\\|=‖\(h1−h2\)\+α​∫01Dh​Φ​\(x,h2\+t​\(h1−h2\)\)​\(h1−h2\)​𝑑t‖\\displaystyle=\\biggl\\\|\(h\_\{1\}\-h\_\{2\}\)\+\\alpha\\int\_\{0\}^\{1\}D\_\{h\}\\Phi\\bigl\(x,h\_\{2\}\+t\(h\_\{1\}\-h\_\{2\}\)\\bigr\)\(h\_\{1\}\-h\_\{2\}\)\\,dt\\biggr\\\|=‖∫01\[I\+α​Dh​Φ​\(x,h2\+t​\(h1−h2\)\)\]​\(h1−h2\)​𝑑t‖\\displaystyle=\\biggl\\\|\\int\_\{0\}^\{1\}\\Bigl\[I\+\\alpha D\_\{h\}\\Phi\\bigl\(x,h\_\{2\}\+t\(h\_\{1\}\-h\_\{2\}\)\\bigr\)\\Bigr\]\(h\_\{1\}\-h\_\{2\}\)\\,dt\\biggr\\\|≤∫01‖I\+α​Dh​Φ​\(x,h2\+t​\(h1−h2\)\)‖op​𝑑t⋅‖h1−h2‖\.\\displaystyle\\leq\\int\_\{0\}^\{1\}\\Bigl\\\|I\+\\alpha D\_\{h\}\\Phi\\bigl\(x,h\_\{2\}\+t\(h\_\{1\}\-h\_\{2\}\)\\bigr\)\\Bigr\\\|\_\{\\mathrm\{op\}\}\\,dt\\;\\cdot\\;\\\|h\_\{1\}\-h\_\{2\}\\\|\.
Now, for anyξ∈B¯r​\(y\)\\xi\\in\\bar\{B\}\_\{r\}\(y\),

‖I\+α​Dh​Φ​\(x,ξ\)‖op\\displaystyle\\bigl\\\|I\+\\alpha D\_\{h\}\\Phi\(x,\\xi\)\\bigr\\\|\_\{\\mathrm\{op\}\}=‖I−α​A​\(x,ξ\)\+α​E​\(x,ξ\)‖op\\displaystyle=\\bigl\\\|I\-\\alpha A\(x,\\xi\)\+\\alpha E\(x,\\xi\)\\bigr\\\|\_\{\\mathrm\{op\}\}≤‖I−α​A​\(x,ξ\)‖op\+α​‖E​\(x,ξ\)‖op\.\\displaystyle\\leq\\\|I\-\\alpha A\(x,\\xi\)\\\|\_\{\\mathrm\{op\}\}\+\\alpha\\\|E\(x,\\xi\)\\\|\_\{\\mathrm\{op\}\}\.
Using Assumption[2](https://arxiv.org/html/2605.24041#S3.I1.i2)and the bound onEE,

‖I−α​A​\(x,ξ\)‖op\\displaystyle\\\|I\-\\alpha A\(x,\\xi\)\\\|\_\{\\mathrm\{op\}\}≤‖I−α​A​\(x,y\)‖op\+α​‖A​\(x,ξ\)−A​\(x,y\)‖op\\displaystyle\\leq\\\|I\-\\alpha A\(x,y\)\\\|\_\{\\mathrm\{op\}\}\+\\alpha\\\|A\(x,\\xi\)\-A\(x,y\)\\\|\_\{\\mathrm\{op\}\}≤q\+α​μ​‖y−ξ‖,\\displaystyle\\leq q\+\\alpha\\mu\\\|y\-\\xi\\\|,and‖E​\(x,ξ\)‖op≤C​‖y−ξ‖\\\|E\(x,\\xi\)\\\|\_\{\\mathrm\{op\}\}\\leq C\\\|y\-\\xi\\\|\. Since‖y−ξ‖≤r\\\|y\-\\xi\\\|\\leq r, we obtain

‖I\+α​Dh​Φ​\(x,ξ\)‖op≤q\+α​\(μ\+C\)​r\.\\bigl\\\|I\+\\alpha D\_\{h\}\\Phi\(x,\\xi\)\\bigr\\\|\_\{\\mathrm\{op\}\}\\leq q\+\\alpha\(\\mu\+C\)r\.
Choosingrrsufficiently small so that

ρ:=q\+α​\(μ\+C\)​r<1,\\rho:=q\+\\alpha\(\\mu\+C\)r<1,we have for allh1,h2∈B¯r​\(y\)h\_\{1\},h\_\{2\}\\in\\bar\{B\}\_\{r\}\(y\)

‖T​\(h1\)−T​\(h2\)‖≤ρ​‖h1−h2‖,\\\|T\(h\_\{1\}\)\-T\(h\_\{2\}\)\\\|\\leq\\rho\\\|h\_\{1\}\-h\_\{2\}\\\|,which establishes thatTTis a contraction onB¯r​\(y\)\\bar\{B\}\_\{r\}\(y\)\.

Conclusion\.Therefore,

TTis a contraction mapping

B¯r​\(y\)\\overline\{B\}\_\{r\}\(y\)into itself\. By theBanach Fixed Point Theorem,

TTadmits a unique fixed point

h∗h^\{\*\}in

B¯r​\(y\)\\overline\{B\}\_\{r\}\(y\), and the iteration

hk\+1=T​\(hk\)h\_\{k\+1\}=T\(h\_\{k\}\)converges linearly to

h∗h^\{\*\}for any initial

h0∈Bδ​\(y\)h\_\{0\}\\in B\_\{\\delta\}\(y\)\.

Part 2: Limiting Error Bound\.

At the fixed pointh∗h^\{\*\}, one can trivially conclude that

0=α​Φ​\(x,h∗\)=α​\(A​\(x,h∗\)​e∗\+b\+R​\(x,h∗\)\)\.0=\\alpha\\Phi\(x,h^\{\*\}\)=\\alpha\\left\(A\(x,h^\{\*\}\)e^\{\*\}\+b\+R\(x,h^\{\*\}\)\\right\)\.Rearranging gives

α​A​\(x,h∗\)​e∗=−α​b−α​R​\(x,h∗\)\.\\alpha A\(x,h^\{\*\}\)e^\{\*\}=\-\\alpha b\-\\alpha R\(x,h^\{\*\}\)\.Taking norms on both sides and by Assumption[1](https://arxiv.org/html/2605.24041#S3.I1.i1),

‖α​A​\(x,h∗\)​e∗‖≤α​‖b‖\+α​L2​‖e∗‖2\.\\\|\\alpha A\(x,h^\{\*\}\)e^\{\*\}\\\|\\leq\\alpha\\\|b\\\|\+\\frac\{\\alpha L\}\{2\}\\\|e^\{\*\}\\\|^\{2\}\.Since

‖I−α​A​\(x,h∗\)‖op≤‖I−α​A​\(x\)‖op\+α​‖A​\(x,h∗\)−A​\(x\)‖op≤q\+α​μ​‖e∗‖,\\\|I\-\\alpha A\(x,h^\{\*\}\)\\\|\_\{\\text\{op\}\}\\leq\\\|I\-\\alpha A\(x\)\\\|\_\{\\text\{op\}\}\+\\alpha\\\|A\(x,h^\{\*\}\)\-A\(x\)\\\|\_\{\\text\{op\}\}\\leq q\+\\alpha\\mu\\\|e^\{\*\}\\\|,ifq\+α​μ​‖e∗‖<1q\+\\alpha\\mu\\\|e^\{\*\}\\\|<1, by the Neumann series, we haveα​A​\(x,h∗\)\\alpha A\(x,h^\{\*\}\)invertible and with

σmin​\(α​A​\(x,h∗\)\)≥1−‖I−α​A​\(x,h∗\)‖op≥1−q−α​μ​‖e∗‖\.\\sigma\_\{\\text\{min\}\}\\left\(\\alpha A\(x,h^\{\*\}\)\\right\)\\geq 1\-\\\|I\-\\alpha A\(x,h^\{\*\}\)\\\|\_\{\\text\{op\}\}\\geq 1\-q\-\\alpha\\mu\\\|e^\{\*\}\\\|\.Plugging back, we obtain

‖e∗‖≤α​‖b‖\+α​L2​‖e∗‖2σmin​\(α​A​\(x,h∗\)\)≤a​‖b‖\+α​L2​‖e∗‖21−q−α​μ​‖e∗‖\.\\\|e^\{\*\}\\\|\\leq\\frac\{\\alpha\\\|b\\\|\+\\tfrac\{\\alpha L\}\{2\}\\\|e^\{\*\}\\\|^\{2\}\}\{\\sigma\_\{\\text\{min\}\}\\left\(\\alpha A\(x,h^\{\*\}\)\\right\)\}\\leq\\frac\{a\\\|b\\\|\+\\tfrac\{\\alpha L\}\{2\}\\\|e^\{\*\}\\\|^\{2\}\}\{1\-q\-\\alpha\\mu\\\|e^\{\*\}\\\|\}\.This is the same quadratic inequality we solved in proving self\-mapping, and with first\-order Taylor expansion, we have

‖e∗‖≤α​‖b‖1−q\+O​\(‖b‖2\)\.\\\|e^\{\*\}\\\|\\leq\\frac\{\\alpha\\\|b\\\|\}\{1\-q\}\+O\(\\\|b\\\|^\{2\}\)\.
□\\square

## Appendix EEmpirical Validation of Strong Monotonicity

Assumption[2](https://arxiv.org/html/2605.24041#S3.I1.i2)requires the linearization at the solution,A​\(x,y\)A\(x,y\), to be bounded and strongly monotone\. Concretely, we assume there exist constants0<m≤M<∞0<m\\leq M<\\inftysuch that

⟨A​\(x,y\)​e,e⟩≥m​‖e‖2,‖A​\(x,y\)‖op≤M,\\langle A\(x,y\)e,e\\rangle\\geq m\\\|e\\\|^\{2\},\\qquad\\\|A\(x,y\)\\\|\_\{\\mathrm\{op\}\}\\leq M,for alle∈ℋe\\in\\mathcal\{H\}\. A natural question is whether this local assumption is reasonable for a trained refinement network\. Since strong monotonicity is not automatic for a learned neural network, we examine it empirically in a controlled 1\-D synthetic setting and on the full U\-Net models used in our main experiments\.

##### 1\-D synthetic experiment\.

We consider the periodic problem

\(I−ε​L\)​y=tanh⁡\(x\),\(I\-\\varepsilon L\)y=\\tanh\(x\),withn=128n=128andε=0\.3\\varepsilon=0\.3, whereLLis the discrete Laplacian with periodic boundary conditions and the solution is computed exactly by matrix inversion\. We train a two\-hidden\-layer MLP refinement operatorΦθ\\Phi\_\{\\theta\}with no architectural constraint enforcing monotonicity of the Jacobian\. We compare two variants:Full IRNOandNo deep supervision\.

For each test pair\(x,y\)\(x,y\), we compute the exact Jacobian at the solution,

A​\(x,y\)=−Dh​Φθ​\(x,y\),A\(x,y\)=\-D\_\{h\}\\Phi\_\{\\theta\}\(x,y\),using automatic differentiation, and measure

m=λmin​\(A\+A⊤2\),M=‖A‖op=σmax​\(A\)\.m=\\lambda\_\{\\min\}\\\!\\left\(\\frac\{A\+A^\{\\top\}\}\{2\}\\right\),\\qquad M=\\\|A\\\|\_\{\\mathrm\{op\}\}=\\sigma\_\{\\max\}\(A\)\.Since⟨A​e,e⟩=⟨A\+A⊤2​e,e⟩\\langle Ae,e\\rangle=\\left\\langle\\frac\{A\+A^\{\\top\}\}\{2\}e,e\\right\\ranglefor allee, the quantitymmis the tightest strong\-monotonicity constant\.

Table[7](https://arxiv.org/html/2605.24041#A5.T7)summarizes the results\. Full IRNO yields uniformly positivemmand tightly concentrated operator norms across all test samples\. In contrast, removing deep supervision leads to predominantly negativemmvalues and substantially larger, more variable operator norms, indicating that the bounded\-and\-strongly\-monotone regime generally does not hold in that setting\.

Table 7:Empirical monotonicity and boundedness of the learned Jacobian at the true solution in the 1\-D synthetic experiment\.![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/strongmono_m_small.png)

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/strongmono_M.png)

Figure 9:Distribution of the estimated strong\-monotonicity constantm=λmin​\(\(A\+A⊤\)/2\)m=\\lambda\_\{\\min\}\(\(A\+A^\{\\top\}\)/2\)\(left\) and operator normM=σmax​\(A\)M=\\sigma\_\{\\max\}\(A\)\(right\) in the 1\-D synthetic experiment, comparing Full IRNO and training without deep supervision\. Full IRNO consistently yields positivemmand tightly boundedMM, whereas removing deep supervision produces predominantly negativemmvalues and much larger variability inMM\. The dashed line in the left panel marks the boundarym=0m=0\. Both y\-axes are log\-scaled\.Figure[9](https://arxiv.org/html/2605.24041#A5.F9)visualizes the distributions ofmmandMM\. The contrast between the two training procedures suggests that the proposed training objective, especially deep supervision, strongly promotes a locally monotone and well\-conditioned refinement map near the solution\.

Table 8:Estimated strong monotonicity of the learned U\-Net refinement operator on the main experimental benchmarks\.
##### U\-Net models on TR\-2D and Active Matter\.

We further evaluate the same criterion on the U\-Net refinement models used in the main paper\. Because exact Jacobian computation is expensive at this scale, we estimate the spectrum using power iteration\. As shown in Table[8](https://arxiv.org/html/2605.24041#A5.T8), the estimated monotonicity constant is positive on essentially all test samples, yielding100%100\\%on TR\-2D and99\.5%99\.5\\%on Active Matter\. These results support the relevance of the strong monotonicity assumption in realistic settings\.

Overall, these experiments show that the bounded\-and\-strongly\-monotone regime required by Assumption[2](https://arxiv.org/html/2605.24041#S3.I1.i2)emerges empirically under the proposed training procedure\.

## Appendix FInvariant\-ball Condition

We derive conditions under which the ballBr​\(y\)B\_\{r\}\(y\)is forward invariant under the iterationT​\(h\)=h\+α​Φ​\(x,h\)T\(h\)=h\+\\alpha\\Phi\(x,h\)\. Leth∈Br​\(y\)h\\in B\_\{r\}\(y\)withr≤δr\\leq\\delta\. Using the affine decomposition from Assumption[1](https://arxiv.org/html/2605.24041#S3.I1.i1),

‖T​\(h\)−y‖\\displaystyle\\\|T\(h\)\-y\\\|=‖h−y\+α​Φ​\(x,h\)‖\\displaystyle=\\\|h\-y\+\\alpha\\Phi\(x,h\)\\\|=‖\(I−α​A​\(x,h\)\)​\(h−y\)\+α​b\+α​R​\(x,h\)‖\\displaystyle=\\\|\(I\-\\alpha A\(x,h\)\)\(h\-y\)\+\\alpha b\+\\alpha R\(x,h\)\\\|≤‖I−α​A​\(x,h\)‖op​‖e‖\+α​‖b‖\+α​‖R​\(x,h\)‖\.\\displaystyle\\leq\\\|I\-\\alpha A\(x,h\)\\\|\_\{\\text\{op\}\}\\\|e\\\|\+\\alpha\\\|b\\\|\+\\alpha\\\|R\(x,h\)\\\|\.
Applying Assumption[2](https://arxiv.org/html/2605.24041#S3.I1.i2)to bound‖I−α​A​\(x,h\)‖op\\\|I\-\\alpha A\(x,h\)\\\|\_\{\\text\{op\}\},

‖I−α​A​\(x,h\)‖op≤‖I−α​A​\(x\)‖op\+α​‖A​\(x,h\)−A​\(x\)‖op≤q\+α​μ​r,\\\|I\-\\alpha A\(x,h\)\\\|\_\{\\text\{op\}\}\\leq\\\|I\-\\alpha A\(x\)\\\|\_\{\\text\{op\}\}\+\\alpha\\\|A\(x,h\)\-A\(x\)\\\|\_\{\\text\{op\}\}\\leq q\+\\alpha\\mu r,
and applying the remainder bound‖R​\(x,h\)‖≤L2​‖e‖2≤L2​r2\\\|R\(x,h\)\\\|\\leq\\frac\{L\}\{2\}\\\|e\\\|^\{2\}\\leq\\frac\{L\}\{2\}r^\{2\}from Assumption[1](https://arxiv.org/html/2605.24041#S3.I1.i1), we obtain

‖T​\(h\)−y‖≤q​r\+α​μ​r2\+α​L2​r2\+α​‖b‖=q​r\+c​r2\+α​‖b‖,\\\|T\(h\)\-y\\\|\\leq qr\+\\alpha\\mu r^\{2\}\+\\alpha\\frac\{L\}\{2\}r^\{2\}\+\\alpha\\\|b\\\|=qr\+cr^\{2\}\+\\alpha\\\|b\\\|,
wherec=α​\(L2\+μ\)≥0c=\\alpha\\left\(\\frac\{L\}\{2\}\+\\mu\\right\)\\geq 0\. ForBr​\(y\)B\_\{r\}\(y\)to be forward invariant, we need‖T​\(h\)−y‖≤r\\\|T\(h\)\-y\\\|\\leq r, i\.e\.,

q​r\+c​r2\+α​‖b‖≤r⇔α​‖b‖≤\(1−q\)​r−c​r2\.qr\+cr^\{2\}\+\\alpha\\\|b\\\|\\leq r\\iff\\alpha\\\|b\\\|\\leq\(1\-q\)r\-cr^\{2\}\.
Case 1\.c=0c=0\.The right\-hand side is\(1−q\)​r\(1\-q\)r, giving the condition

‖b‖≤\(1−q\)​rα\.\\\|b\\\|\\leq\\frac\{\(1\-q\)r\}\{\\alpha\}\.
Case 2\.c\>0c\>0\.The function\(1−q\)​r−c​r2\(1\-q\)r\-cr^\{2\}is maximized atr∗=1−q2​cr^\{\*\}=\\frac\{1\-q\}\{2c\}, yielding maximum value\(1−q\)24​c\\frac\{\(1\-q\)^\{2\}\}\{4c\}\. Therefore, a validr∈\(0,δ\]r\\in\(0,\\delta\]exists whenever

‖b‖≤\(1−q\)24​α​c,\\\|b\\\|\\leq\\frac\{\(1\-q\)^\{2\}\}\{4\\alpha c\},in which case the two positive roots of\(1−q\)​r−c​r2=α​‖b‖\(1\-q\)r\-cr^\{2\}=\\alpha\\\|b\\\|are

r±=\(1−q\)±\(1−q\)2−4​α​c​‖b‖2​c,r\_\{\\pm\}=\\frac\{\(1\-q\)\\pm\\sqrt\{\(1\-q\)^\{2\}\-4\\alpha c\\\|b\\\|\}\}\{2c\},andBr​\(y\)B\_\{r\}\(y\)is forward invariant for anyr∈\[r−,r\+\]r\\in\[r\_\{\-\},r\_\{\+\}\]\.

## Appendix GExperimental Setup and Implementation Details

### G\.1Dataset Specifications

#### G\.1\.1Turbulent Radiative Layer\-2D \(TR\-2D\)

##### Data Generation\.

The TR\-2D dataset simulates Kelvin–Helmholtz instability using the ATHENA\+\+ hydrodynamics code\[Fieldinget al\.,[2020](https://arxiv.org/html/2605.24041#bib.bib55)\]\. Initial conditions consist of two gas layers with different temperatures and velocities, creating shear\-driven turbulence\.

##### Physical system\.

This dataset models the interaction between hot, dilute gas and cold, dense gas moving at subsonic velocities, a configuration unstable to the Kelvin\-Helmholtz instability\. The turbulent mixing populates intermediate\-temperature gas, which rapidly cools as heating and cooling become imbalanced\. This process is fundamental to understanding phase structure in the interstellar and circumgalactic medium\. The governing equations are given by

∂ρ∂t\+∇⋅\(ρ​v→\)\\displaystyle\\frac\{\\partial\\rho\}\{\\partial t\}\+\\nabla\\cdot\(\\rho\\vec\{v\}\)=0,\\displaystyle=0,∂\(ρ​v→\)∂t\+∇⋅\(ρ​v→​v→\)\\displaystyle\\frac\{\\partial\(\\rho\\vec\{v\}\)\}\{\\partial t\}\+\\nabla\\cdot\(\\rho\\vec\{v\}\\vec\{v\}\)=−∇P,\\displaystyle=\-\\nabla P,∂E∂t\+∇⋅\(\(E\+P\)​v→\)\\displaystyle\\frac\{\\partial E\}\{\\partial t\}\+\\nabla\\cdot\(\(E\+P\)\\vec\{v\}\)=−Etcool,\\displaystyle=\-\\frac\{E\}\{t\_\{\\text\{cool\}\}\},whereρ\\rhois density,v→\\vec\{v\}is the 2D velocity field,PPis pressure,EEis total energy, andtcoolt\_\{\\text\{cool\}\}is the cooling time parameter\.

##### Dataset details\.

The dataset contains 90 trajectories \(10 random initializations for each of 9 cooling time values:tcool∈\{0\.03,0\.06,0\.1,0\.18,0\.32,0\.56,1\.00,1\.78,3\.16\}t\_\{\\text\{cool\}\}\\in\\\{0\.03,0\.06,0\.1,0\.18,0\.32,0\.56,1\.00,1\.78,3\.16\\\}\)\. Each trajectory consists of 101 timesteps at resolution 384×128, with spatial domainx∈\[−0\.5,0\.5\]x\\in\[\-0\.5,0\.5\],y∈\[−1,2\]y\\in\[\-1,2\]and temporal snapshots separated byΔ​t=1\.597\\Delta t=1\.597in simulation time\. Fields include density, pressure, and velocity components\.

##### Task\.

Given 4 consecutive frames\{st−3,st−2,st−1,st\}\\\{s\_\{t\-3\},s\_\{t\-2\},s\_\{t\-1\},s\_\{t\}\\\}wherest=\(ρt,v→t,Pt\)s\_\{t\}=\(\\rho\_\{t\},\\vec\{v\}\_\{t\},P\_\{t\}\), predict the next framest\+1s\_\{t\+1\}\. This conditional prediction task tests the model’s ability to capture short\-term dynamics and turbulent evolution\.

#### G\.1\.2Active Matter \(AM\)

##### Physical system\.

This dataset simulates continuum dynamics of rod\-like active particles immersed in a Stokes fluid\. The system captures collective behavior, including energy transfer across scales, vorticity\-orientation coupling, and phase transitions from isotropic to nematic states\. The governing equations \(Eqs\. 1\-5 in the associated paper\) describe the evolution of concentration, velocity, orientation tensor, and strain\-rate tensor fields\.

##### Dataset details\.

The dataset contains 225 simulations spanning parameter spaceα∈\{−1,−2,−3,−4,−5\}\\alpha\\in\\\{\-1,\-2,\-3,\-4,\-5\\\}\(dipole strength\),β=0\.8\\beta=0\.8, andζ∈\{1,3,5,7,9,11,13,15,17\}\\zeta\\in\\\{1,3,5,7,9,11,13,15,17\\\}\(alignment strength\), with 5 trajectories per parameter set using different random initializations\. Each trajectory consists of 81 timesteps at resolution 256×256, covering temporal ranget∈\[0,20\]t\\in\[0,20\]seconds with snapshots every 0\.25 seconds\. The spatial domain isLx=Ly=10L\_\{x\}=L\_\{y\}=10with periodic boundary conditions\.

##### Task\.

Given 4 consecutive frames\{st−3,st−2,st−1,st\}\\\{s\_\{t\-3\},s\_\{t\-2\},s\_\{t\-1\},s\_\{t\}\\\}wherest=\(ρt,v→t,Dt,Ut\)s\_\{t\}=\(\\rho\_\{t\},\\vec\{v\}\_\{t\},D\_\{t\},U\_\{t\}\)includes concentration, velocity, orientation tensor, and strain\-rate tensor, predict the next framest\+1s\_\{t\+1\}\.

#### G\.1\.3ERA5 Weather Super\-Resolution

We conduct our experiments on the climate super\-resolution benchmark provided by SuperBench\[Renet al\.,[2025](https://arxiv.org/html/2605.24041#bib.bib20)\], which is constructed from the ERA5 reanalysis dataset\[Hersbachet al\.,[2020](https://arxiv.org/html/2605.24041#bib.bib22)\]\. In this work, we strictly follow the dataset construction and experimental protocol defined in SuperBench\.

##### Data Source and Task Definition\.

ERA5 is a global atmospheric reanalysis dataset produced by the European Centre for Medium\-Range Weather Forecasts \(ECMWF\) on a regular latitude–longitude grid\. Following SuperBench, we focus on three climate variables: kinetic energy, temperature, and total column water vapor\. We follow Scenario \(i\) \(the general computer vision setting\) from SuperBench, in which low\-resolution \(LR\) atmospheric states are obtained by bicubic downsampling without noise\. The task is formulated as supervised spatial super\-resolution, mapping LR climate fields of resolution45×9045\\times 90to corresponding high\-resolution \(HR\) targets of resolution720×1440720\\times 1440, resulting in a×16\\times 16super\-resolution factor along each spatial dimension\.

##### Dataset construction\.

We follow the official data construction and temporal splits provided by SuperBench\. Training data are selected from the years 2008, 2010, 2011, and 2013\. The validation sets for interpolation and extrapolation \(look\-back\) evaluation are drawn from 2012 and 2007, respectively\. The test sets for the corresponding interpolation and extrapolation tasks are from 2009 and 2014–2015, respectively\. All methods are trained and evaluated on the same splits to ensure fair comparison\. For visualizations in the main paper, we use samples from the interpolation test set \(year 2009\)\.

##### Data Preprocessing\.

Following SuperBench, we standardize each variable using statistics computed from the training split, and apply the same normalization to the validation and test sets\. No additional data augmentation or task\-specific preprocessing is performed\.

#### G\.1\.4CE\-Gauss \(Irregular Mesh\)

##### Physical system\.

CE\-Gauss is an unstructured\-mesh benchmark from the RIGNO dataset\[Mousaviet al\.,[2025](https://arxiv.org/html/2605.24041#bib.bib59)\]\. It simulates convection\-diffusion dynamics on irregular Gaussian\-distributed point clouds, capturing transport phenomena on non\-uniform spatial domains that are not representable on regular grids\. The dataset provides 4 physical variables per node\.

##### Dataset details\.

The mesh contains 16,384 nodes with irregular spatial distribution\. We follow the official train/validation/test splits provided with the RIGNO benchmark\. The evaluation task is autoregressive rollout over 7 timesteps\.

##### Task\.

Given the current state on the irregular mesh, predict the next timestep\. IRNO is applied withK=4K=4refinement steps andα=0\.3\\alpha=0\.3, using RIGNO as the base operator\. This experiment shows that the IRNO framework extends to graph\-based operators on irregular meshes\.

The refinement model for CE\-Gauss shares the same RIGNO architecture as the frozen base operator\. It takes an 8\-channel input formed by concatenating the current estimateh^k\\hat\{h\}\_\{k\}with the frozen model’s initial predictionh^0\\hat\{h\}\_\{0\}, and outputs a 4\-channel correction field\. The update is applied in normalized solution space, making the step\-sizeα\\alphadimensionless and consistent across output variables\. Training minimizes a trajectory supervision loss averaged overKKunrolled refinement steps rather than only the final\-step loss, which stabilizes the iterative dynamics\. At inference, refinement is applied at every autoregressive step as zero\-shot generalization; the model is trained only for single\-step prediction\.

### G\.2Architecture Details

#### G\.2\.1Base Operators

Table 9:Hyperparameters of FNO and TFNO across different datasets\.##### Fourier Neural Operator \(FNO\)\.

We adopt the FNO2D from the Well and SuperBench as the base operator, which implements the Fourier Neural Operator\[Liet al\.,[2021](https://arxiv.org/html/2605.24041#bib.bib3)\]\. Table[9](https://arxiv.org/html/2605.24041#A7.T9)summarizes the hyperparameters used for FNO and TFNO across all datasets\.

##### Tucker\-Factorized Fourier Neural Operator \(TFNO\)\.

To improve parameter efficiency and generalization, we also employ the TFNO\[Kovachkiet al\.,[2023b](https://arxiv.org/html/2605.24041#bib.bib56)\]\. TFNO replaces the dense spectral convolution layers in the standard FNO with Tucker\-factorized tensor contractions\. This decomposition allows for significant compression of the network weights while preserving expressivity\. As shown in Table[9](https://arxiv.org/html/2605.24041#A7.T9), the TFNO follows a similar architectural configuration to the FNO \(in terms of depth and channel width\)\.

##### Wide Activation Super\-Resolution \(WDSR\)\.

We adopt the WDSR\-A architecture\[Yuet al\.,[2018](https://arxiv.org/html/2605.24041#bib.bib25)\]as a CNN\-based base operator for scientific super\-resolution within the SuperBench framework\. The model follows a head–body–tail design and emphasizes local feature extraction through wide\-activation residual blocks\. Hyperparameter details are provided in Tables[10](https://arxiv.org/html/2605.24041#A7.T10)\.

- •Residual Blocks\.The network body consists of 18 lightweight residual blocks, each using an expansion ratio of 4 and residual scaling factor of 0\.1\. Each block contains two3×33\\times 3convolutional layers with ReLU activation and weight normalization\.
- •Upsampling and Skip Connection\.The output is produced by a tail convolution followed by pixel\-shuffle upsampling\. In parallel, a skip branch directly upsamples the input using a separate convolution and pixel\-shuffle, and the two paths are summed to form the final output\.

Table 10:Hyperparameters of the WDSR\-A model across different datasets\.

#### G\.2\.2Refinement OperatorΦθ\\Phi\_\{\\theta\}

The refinement operator is implemented as a lightweight U\-Net with the following specifications:

##### Encoder Path\.

The encoder receives the concatenation of the original inputxxand the current estimatehkh\_\{k\}, resulting in2​C2Cinput channels\. It consists of 3 levels of downsampling\. Each level contains aconvolutional blockfollowed by a2×22\\times 2max\-poolinglayer\. Each convolutional block is composed of two successive sequences of:

Conv2d​\(3×3\)→BatchNorm2d→GELU\\text\{Conv2d\}\(3\\times 3\)\\rightarrow\\text\{BatchNorm2d\}\\rightarrow\\text\{GELU\}The number of feature channels doubles at each level, starting fromCbaseC\_\{\\text\{base\}\}and reaching4​Cbase4C\_\{\\text\{base\}\}before the bottleneck\.

##### Bottleneck\.

The bottleneck bridges the encoder and decoder at the lowest spatial resolution\. It consists of a single convolutional block \(two3×33\\times 3convolutions with BatchNorm and GELU\) that processes the features in a2depth​Cbase2^\{\\text\{depth\}\}C\_\{\\text\{base\}\}\-dimensional latent space\.

##### Decoder Path\.

The decoder is symmetric to the encoder and performs 3 levels of upsampling\. At each level, the feature maps are first upsampled usingbilinear interpolation\(scale factor of 2\) followed by a3×33\\times 3convolution that halves the channel dimension\. These features are then concatenated with the corresponding skip connections from the encoder\. The merged features pass through a standard dual\-convolutional block\.

##### Output Head & Initialization\.

A final1×11\\times 1convolution maps theCbaseC\_\{\\text\{base\}\}hidden channels back to the output dimension\. To ensure the iterative process starts with stable, small\-magnitude corrections, we applyXavier uniform initializationto the weights with a gain of0\.10\.1and initialize all biases to zero\.

##### Padding\.

To respect the periodic nature of the scientific datasets used in our benchmark, all3×33\\times 3convolutional layers inΦθ\\Phi\_\{\\theta\}utilizecircular paddinginstead of zero padding, ensuring spatial continuity\.

##### Parameter Count\.

Table[11](https://arxiv.org/html/2605.24041#A7.T11)reports the parameter counts \(in millions\) of all models across different datasets\. A cross \(×\\times\) indicates that the corresponding model is not used for that dataset\. For our proposed Iterative Refinement Neural Operator \(IRNO\), we report the cumulative number of parameters required for a complete inference cycle consisting ofKKrefinement iterations\. The last three rows correspond to large monolithic baselines with different refinement models, which are introduced and discussed in the main paper[4\.7](https://arxiv.org/html/2605.24041#S4.SS7)\.

Table 11:Model parameter counts \(in millions\) across different datasets\.ModelTR\-2D\(M\)Active Matter\(M\)ERA5\(M\)FNO \(base\)19\.019\.04\.8TFNO \(base\)19\.319\.3×\\timesWDSR \(base\)×\\times×\\times1\.6RefinementΦθ\\Phi\_\{\\theta\}8\.68\.72\.1FNO \+ IRNO \(K=6K=6\)27\.627\.76\.9FNO \+ SRCNN×\\times×\\times6\.6FNO \+ UNet \(×15\\times 15\)×\\times×\\times37\.3FNO \+ FNO×\\times×\\times9\.5

### G\.3Training Configuration

Algorithm 1Training Iterative Refinement Neural Operator \(IRNO\)1:Input:Dataset

𝒟=\{\(xi,yi\)\}i=1N\\mathcal\{D\}=\\\{\(x\_\{i\},y\_\{i\}\)\\\}\_\{i=1\}^\{N\}, base operator

𝒯base\\mathcal\{T\}\_\{\\text\{base\}\}, refinement steps

KK, step size

α\\alpha, weights

βsp,βfp\\beta\_\{\\text\{sp\}\},\\beta\_\{\\text\{fp\}\}
2:Initialize:Refinement parameters

θ\\theta
3:forepoch

=1,…,E=1,\\dots,Edo

4:forbatch

\(x,y\)∈𝒟\(x,y\)\\in\\mathcal\{D\}do

5:

h0←𝒯base​\(x\)h\_\{0\}\\leftarrow\\mathcal\{T\}\_\{\\text\{base\}\}\(x\)⊳\\trianglerightBase prediction \(No gradient\)

6:

ℒaccum←0\\mathcal\{L\}\_\{\\text\{accum\}\}\\leftarrow 0
7:for

k=1k=1to

KKdo

8:

hk←hk−1\+α⋅Φθ​\(x,hk−1\)h\_\{k\}\\leftarrow h\_\{k\-1\}\+\\alpha\\cdot\\Phi\_\{\\theta\}\(x,h\_\{k\-1\}\)⊳\\trianglerightRefinement step

9:

ℒk←ℒspatial​\(hk,y\)\+βsp⋅ℒspectral​\(hk,y\)\\mathcal\{L\}\_\{k\}\\leftarrow\\mathcal\{L\}\_\{\\text\{spatial\}\}\(h\_\{k\},y\)\+\\beta\_\{\\text\{sp\}\}\\cdot\\mathcal\{L\}\_\{\\text\{spectral\}\}\(h\_\{k\},y\)
10:

ℒaccum←ℒaccum\+ℒk\\mathcal\{L\}\_\{\\text\{accum\}\}\\leftarrow\\mathcal\{L\}\_\{\\text\{accum\}\}\+\\mathcal\{L\}\_\{k\}
11:endfor

12:

ℒfp←βfp⋅‖Φθ​\(x,y\)‖2\\mathcal\{L\}\_\{\\text\{fp\}\}\\leftarrow\\beta\_\{\\text\{fp\}\}\\cdot\\\|\\Phi\_\{\\theta\}\(x,y\)\\\|^\{2\}⊳\\trianglerightFixed\-point regularization

13:

ℒtotal←1K​ℒaccum\+ℒfp\\mathcal\{L\}\_\{\\text\{total\}\}\\leftarrow\\frac\{1\}\{K\}\\mathcal\{L\}\_\{\\text\{accum\}\}\+\\mathcal\{L\}\_\{\\text\{fp\}\}
14:Update

θ\\thetaby minimizing

ℒtotal\\mathcal\{L\}\_\{\\text\{total\}\}
15:endfor

16:endfor

##### Base Operator Training\.

In our experiments, the base operators are initialized from a pretrained checkpoint\. The base operators are kept frozen during all subsequent training and evaluation stages\.

##### Refinement Operator Training\.

Algorithm[1](https://arxiv.org/html/2605.24041#alg1)provides the complete training procedure for IRNO, detailing the progressive refinement steps, combined spatial\-spectral losses, and fixed\-point regularization\.

The refinement operator is optimized using AdamW with an initial learning rate of3×10−43\\times 10^\{\-4\}and weight decay1×10−51\\times 10^\{\-5\}\. We employ a cosine learning rate scheduler with a minimum learning rate of10−610^\{\-6\}\. Training is conducted for 250 epochs with a per\-GPU batch size of 16, and gradient norms are clipped to a maximum value of 1\.0 to ensure stability\. During training, the refinement operator is unrolled forK=6K=6andK=4K=4iterative refinement steps when used with the FNO and WDSR/TFNO base models, respectively, with a fixed step sizeαr​e​f​i​n​e=0\.25\\alpha\_\{refine\}=0\.25for ERA5 dataset andαr​e​f​i​n​e=0\.2\\alpha\_\{refine\}=0\.2for TR\-2D and AM datasets\. The training objective consists of the standard reconstruction loss augmented with a spectral loss term, as described below\.

##### Baseline implementations\.

We reproduce HFS\[Khodakaramiet al\.,[2026](https://arxiv.org/html/2605.24041#bib.bib32)\]and HiNOTE\[Luoet al\.,[2024](https://arxiv.org/html/2605.24041#bib.bib57)\]from their respective official code repositories, evaluated on the ERA516×16\{\\times\}super\-resolution task under the same dataset split and evaluation protocol used for IRNO\.

##### Progressive Spectral Loss Schedule\.

We incorporate a spectral loss to encourage accurate reconstruction across different frequency bands\. Two distinct mechanisms control the spectral emphasis\. The frequency exponentλk\\lambda\_\{k\}determines the per\-step weighting of high\-frequency components, and the spectral loss weightβspectral\\beta\_\{\\text\{spectral\}\}scales the overall spectral term in the training objective\.

For the TR\-2D and Active Matter datasets, the frequency exponentλk\\lambda\_\{k\}increases linearly fromλstart=1\.0\\lambda\_\{\\text\{start\}\}=1\.0toλend=2\.0\\lambda\_\{\\text\{end\}\}=2\.0over theKKrefinement steps within each forward pass\. The spectral loss weightβspectral\\beta\_\{\\text\{spectral\}\}undergoes a linear warm\-up over the first 5 training epochs before being held fixed for the remainder of training\.

For the ERA5 dataset, we adopt the same progressive exponent schedule \(λk∈\[1\.0,2\.0\]\\lambda\_\{k\}\\in\[1\.0,2\.0\]over refinement steps\)\. Table[12](https://arxiv.org/html/2605.24041#A7.T12)shows that the progressive schedule improves ACC from0\.8750\.875to0\.8920\.892and reduces RFNE from0\.2350\.235to0\.2140\.214relative to a fixed exponent \(λk=1\.5\\lambda\_\{k\}=1\.5\), confirming that progressive spectral emphasis benefits the super\-resolution task\.

Table 12:ERA5 spectral loss schedule ablation \(FNO base,K=6K=6\)\. Progressiveλk∈\[1\.0,2\.0\]\\lambda\_\{k\}\\in\[1\.0,2\.0\]vs\. fixedλk=1\.5\\lambda\_\{k\}=1\.5\.
##### Training Compute Budget\.

Table[13](https://arxiv.org/html/2605.24041#A7.T13)reports per\-batch wall\-clock timing for the IRNO refinement operator \(U\-Net backbone\) on Active Matter at different unrolling horizonsKK\. Training overhead scales sublinearly in total wall\-clock time, withK=4K=4incurring2\.69×2\.69\\timesoverhead relative toK=1K=1andK=6K=6incurring4\.14×4\.14\\times, reflecting the amortization of optimizer and data\-loading costs across longer unrolls\.

Table 13:Per\-batch training timing \(ms\) for IRNO on Active Matter at different unrolling horizonsKK\. Measurements on a single NVIDIA RTX PRO 6000 Blackwell Max\-Q GPU\.

### G\.4Implementation Details

##### Software and Libraries\.

All experiments are implemented in Python 3\.10\.16 using PyTorch 2\.5\.1 with CUDA 12\.1\. NumPy and Matplotlib are used for numerical processing and visualization, respectively\. FLOPs are measured usingfvcore\.

##### Hardware\.

All experiments were conducted on NVIDIA GPUs\. Specifically, we used NVIDIA L40 GPUs with 48 GB memory and NVIDIA RTX PRO 6000 Blackwell GPUs with 96 GB memory\.

### G\.5Evaluation Metrics

#### G\.5\.1Primary Metrics

##### Variance\-scaled Root Mean Squared Error \(VRMSE\)\.

For TR\-2D and Active Matter, we report the variance\-scaled root mean squared error \(VRMSE\)Ohanaet al\.\[[2024](https://arxiv.org/html/2605.24041#bib.bib44)\], defined as the square root of the mean squared error normalized by the variance of the ground\-truth field\. Specifically, given a predictionuuand referencevvwith spatial meanv¯\\bar\{v\}, we compute

VRMSE​\(u,v\)=\(⟨\|u−v\|2⟩⟨\|v−v¯\|2⟩\+ε\)1/2,\\mathrm\{VRMSE\}\(u,v\)=\\left\(\\frac\{\\langle\|u\-v\|^\{2\}\\rangle\}\{\\langle\|v\-\\bar\{v\}\|^\{2\}\\rangle\+\\varepsilon\}\\right\)^\{1/2\},where⟨⋅⟩\\langle\\cdot\\rangledenotes averaging over all spatial locations and output channels, andε\\varepsilonis a small constant for numerical stability\.

##### Anomaly Correlation Coefficient \(ACC\)\.

For ERA5 dataset, we report the Anomaly Correlation Coefficient \(ACC\) to assess the spatial pattern similarity between predictions and ground truth\. Given a predicted fieldu^\\hat\{u\}and reference fielduu, we first compute their anomalies by removing the spatial mean:

u^′=u^−u^¯,u′=u−u¯,\\hat\{u\}^\{\\prime\}=\\hat\{u\}\-\\overline\{\\hat\{u\}\},\\qquad u^\{\\prime\}=u\-\\overline\{u\},where\(⋅\)¯\\overline\{\(\\cdot\)\}denotes the spatial mean\. The ACC is then defined as

ACC=⟨u^′,u′⟩⟨u^′,u^′⟩​⟨u′,u′⟩,\\mathrm\{ACC\}=\\frac\{\\langle\\hat\{u\}^\{\\prime\},u^\{\\prime\}\\rangle\}\{\\sqrt\{\\langle\\hat\{u\}^\{\\prime\},\\hat\{u\}^\{\\prime\}\\rangle\\,\\langle u^\{\\prime\},u^\{\\prime\}\\rangle\}\},where⟨⋅,⋅⟩\\langle\\cdot,\\cdot\\rangledenotes the inner product over all spatial locations and channels\. ACC measures the similarity of spatial anomaly patterns and is particularly suitable for evaluating geophysical and climate variables\.

##### Relative Frobenius Norm Error \(RFNE\)\.

We evaluate reconstruction accuracy using the Relative Frobenius Norm Error \(RFNE\), which measures the normalized discrepancy between the predicted fieldu^\\hat\{u\}and the ground\-truth fielduu\. Specifically, RFNE is defined as

RFNE=∥u^−u∥F∥u∥F,\\mathrm\{RFNE\}=\\frac\{\\lVert\\hat\{u\}\-u\\rVert\_\{F\}\}\{\\lVert u\\rVert\_\{F\}\},where∥⋅∥F\\lVert\\cdot\\rVert\_\{F\}denotes the Frobenius norm over all spatial dimensions and channels\. RFNE provides a scale\-invariant measure of global reconstruction error and is widely adopted in scientific super\-resolution benchmarks\.

### G\.6Computational Cost Analysis

##### FLOPs Calculation\.

We compute floating\-point operation counts \(FLOPs\) using PyTorch’sFlopCountAnalysiswith custom operator handlers designed to match the actual execution graph of our models\. FLOPs are estimated at the operator level using the tensor shapes observed during a forward pass\. Reported values correspond to a single forward pass of the base operator network, and cumulative FLOPs across multiple refinement steps are obtained by linear accumulation\.

##### Memory Profiling\.

We measure peak GPU memory usage during inference using PyTorch’s CUDA memory profiler\. All measurements are conducted with models in evaluation mode and undertorch\.no\_grad\(\)to exclude gradient storage\.

Peak GPU memory usage is measured during inference using PyTorch’s CUDA memory profiler\. For the base FNO model, peak memory is obtained from a single forward pass\. For iterative models with refinement steps \(IRNO\), memory usage is profiled throughout the refinement process by querying the peak allocated memory after each refinement step, and the maximum value observed across all iterations is reported as the peak memory usage for a given number of stepsKK\. Specifically, we record the maximum allocated CUDA memory viatorch\.cuda\.max\_memory\_allocated\.

### G\.7Monolithic Residual Correction Models

Table[14](https://arxiv.org/html/2605.24041#A7.T14)summarizes the architectural designs of the large monolithic refinement model baselines, performing single\-pass inference, see[4\.7](https://arxiv.org/html/2605.24041#S4.SS7)\.

Table 14:Architectural configurations of monolithic refinement model baselines used for comparison\.

## Appendix HAdditional Tables and Visualization

##### Numerical values for scatter plots\.

Figure[8](https://arxiv.org/html/2605.24041#S4.F8)in the main text presents the trade\-offs between accuracy and computational cost for all compared models\. To provide a precise quantitative reference, Table[15](https://arxiv.org/html/2605.24041#A8.T15)lists the corresponding numerical values\. For the iterative refinement baseline \(IRNO\), we report cumulative FLOPs, as well as peak memory usage across the entire inference process\.

Table 15:Detailed numerical results corresponding to the scatter plots in the main text\. The table reports ACC, RFNE, peak memory \(MB\), and cumulative FLOPs \(G\) for the base model, IRNO iterations, and the three large monolithic refinement model baselines\.
##### Inference Time per Step\.

We report the wall\-clock inference time per step for different refinement pipelines on a single GPU\. All measurements are conducted on an NVIDIA RTX PRO 6000 Blackwell Max\-Q Workstation Edition GPU\. For each method, the inference time is measured over five independent runs, and we report the mean and standard deviation\. The reported numbers correspond to the average wall\-clock time of one forward step, measured in milliseconds \(ms\)\.

Table 16:Inference time per step \(ms\) on a single NVIDIA RTX PRO 6000 Blackwell Max\-Q GPU\. Results are reported as mean±\\pmstandard deviation over five runs\.
##### Visualization for ERA5 with8×8\\timesupscale factor\.

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/f8_F1.png)

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/f8_F2.png)

Figure 10:Qualitative results on ERA5 for Kinetic Energy \(top\) and Temperature Field \(bottom\) under8×8\\timesspatial downsampling\. The model is trained with up to4 refinement steps\. We visualize the ground truth, the base FNO prediction, and the progressive refinement results after each refinement step\.Figure[10](https://arxiv.org/html/2605.24041#A8.F10)presents qualitative results on ERA5 for Kinetic Energy and Temperature Field under the8×8\\timesspatial downsampling\. The visualizations demonstrate that the proposed iterative correction framework consistently improves predictions across different physical fields and upsampling factors\.

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/acc_channel0_8.png)

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/acc_channel1_8.png)

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/acc_channel2_8.png)

Figure 11:Temporal evolution of ACC on ERA5 under8×8\\timesspatial downsampling for Field 1: Kinetic Energy, Field 2: Temperature, and Field 3: Precipitation \(left to right\)\. The model is trained with up to4 refinement steps\. We report the base FNO and predictions obtained after different numbers of iterative refinement steps\.Figure[11](https://arxiv.org/html/2605.24041#A8.F11)shows the temporal evolution of anomaly correlation coefficient \(ACC\) on ERA5 for three different physical fields under the8×8\\timesdownsampling\. Across all fields, increasing the number of refinement steps consistently improves forecasting accuracy compared to the base FNO\. These trends indicate that the iterative refinement framework generalizes well across both different fields and lower upsampling factors\.

##### Visualization for ERA5 with16×16\\timesupscale factor\.

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/f16_F1.png)Figure 12:Qualitative results on ERA5 for Temperature Field under16×16\\timesspatial downsampling\. The model is trained with up to6 refinement steps\. We visualize the ground truth, the base FNO prediction, and the progressive refinement results after each refinement step\.Figure[12](https://arxiv.org/html/2605.24041#A8.F12)presents additional qualitative results on ERA5 for the Temperature Field under the16×16\\timesspatial downsampling, complementing Figure[1](https://arxiv.org/html/2605.24041#S1.F1)in the main paper, which reports results on the Kinetic Energy Field\. Visualizations demonstrate that the proposed iterative correction framework generalizes consistently across different physical fields\.

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/acc_channel0_16.png)

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/acc_channel1_16.png)

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/acc_channel2_16.png)

Figure 13:Temporal evolution of ACC on ERA5 under16×16\\timesspatial downsampling for Field 1: Kinetic Energy, Field 2: Temperature, and Field 3: Precipitation \(left to right\)\. We report the base FNO and predictions obtained after different numbers of iterative refinement steps\.Figure[13](https://arxiv.org/html/2605.24041#A8.F13)reports the temporal evolution of anomaly correlation coefficient \(ACC\) on ERA5 for three different physical fields under the16×16\\timesdownsampling setting\. Across all fields, increasing the number of refinement steps consistently improves forecasting accuracy compared to the base FNO\. These trends are consistent across different fields, suggesting that the performance improvements are not field\-specific but arise from the proposed iterative correction framework\.

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/unet1.png)
![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/fno1.png)
![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/SRCNN1.png)

Figure 14:Qualitative comparison between Ground Truth, Base FNO, and single\-step refinement across three large monolithic baseline architectures \(Models UNet \(10×10\\times\), FNO, and SRCNN\)\.Figure[14](https://arxiv.org/html/2605.24041#A8.F14)illustrates the prediction results for three large monolithic baseline models\. These single\-step refinement \(right column\) still fail to resolve small\-scale turbulent features\. In contrast, our IRNO framework, see Figure[1](https://arxiv.org/html/2605.24041#S1.F1)is more effective at capturing fine\-grained physics\.

##### Visualization for TR\-2D and Active Matter\.

Figure[15](https://arxiv.org/html/2605.24041#A8.F15)and[16](https://arxiv.org/html/2605.24041#A8.F16)visualize predictions of FNO\-based and TFNO\-based IRNO on TR\-2D density field\. Figure[17](https://arxiv.org/html/2605.24041#A8.F17)visualizes predictions of FNO\-based and TFNO\-based IRNO on Active Matter velocity field\.

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/tr2d_fno_vis.png)Figure 15:Qualitative results on TR\-2D with FNO as base operator\. The model is trained with up to6 refinement steps\. We visualize the ground truth, the base FNO prediction, and the progressive refinement results after each refinement step\.![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/tr2d_tfno_vis.png)Figure 16:Qualitative results on TR\-2D with TFNO as base operator\. The model is trained with up to4 refinement steps\. We visualize the ground truth, the base TFNO prediction, and the progressive refinement results after each refinement step\.![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/AM_fno_vis.png)\(a\)FNO base operator \(6 refinement steps\)
![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/AM_tfno_vis.png)\(b\)TFNO base operator \(4 refinement steps\)

Figure 17:Qualitative results on Active Matter with different base operators\. We visualize the ground truth, base operator predictions, and progressive refinement results after each correction step\.
##### Long\-Horizon Extrapolation\.

Figure[18](https://arxiv.org/html/2605.24041#A8.F18)shows VRMSE as a function of refinement stepkk, evaluated beyond the training horizon ofK=6K=6up tok=48k=48on Active Matter with FNO base operator\. With step sizeα=0\.20\\alpha=0\.20, the error reaches a minimum of approximately0\.05050\.0505neark≈14k\\approx 14before diverging to0\.08060\.0806atk=48k=48\. Withα=0\.05\\alpha=0\.05, refinement remains stable throughout without divergence\. Two complementary strategies can protect against divergence\. First, step\-size scheduling reducesα\\alphaaskkincreases\. Second, adaptive stopping halts iteration when‖Φθ​\(x,hk\)‖\\\|\\Phi\_\{\\theta\}\(x,h\_\{k\}\)\\\|falls below a user\-defined threshold tied to the bias levelα​‖b‖/\(1−q\)\\alpha\\\|b\\\|/\(1\-q\)from Corollary[3\.3](https://arxiv.org/html/2605.24041#S3.Thmtheorem3)\.

![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/irno_extrapolation_alpha_comparison.png)Figure 18:VRMSE vs\. refinement stepkkup to8×8\\timesthe training horizon \(K=6K=6, dashed vertical line\) forα=0\.05\\alpha=0\.05\(stable\) andα=0\.20\\alpha=0\.20\(diverges\) with FNO\-based IRNO on Active Matter\.
##### HFS\-ResUNet and IRNO Complementarity\.

Table[17](https://arxiv.org/html/2605.24041#A8.T17)reports VRMSE at each refinement step when IRNO is applied on top of an HFS\-ResUNet base on Active Matter\. Starting from the HFS base error of0\.06310\.0631, IRNO progressively reduces error to0\.04860\.0486atk=6k=6, remaining stable atk=8k=8\(0\.04870\.0487\), demonstrating that iterative refinement compounds gains from frequency\-aware architectures\. Figure[19](https://arxiv.org/html/2605.24041#A8.F19)further shows that IRNO achieves lower spectral error energy than HFS\-ResUNet across the full radial frequency range\.

Table 17:HFS\-IRNO VRMSE per refinement step on Active Matter\.![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/spectral_error_hfs_fft_vs_irno.png)Figure 19:Spectral error energy vs\. radial frequency for HFS\-ResUNet and IRNO on Active Matter\. Shaded regions denote±1\\pm 1standard deviation\.
##### Alternative Refinement Architectures\.

Table[18](https://arxiv.org/html/2605.24041#A8.T18)reports VRMSE for IRNO with different refinement backbone architectures on Active Matter \(TFNO base,K=4K=4\)\. All architectures achieve substantial error reduction, confirming that IRNO’s gains are driven by the iterative mechanism rather than a specific backbone choice\.

Table 18:IRNO with alternative refinement architectures on Active Matter \(TFNO base,K=4K=4\)\.
##### Normalization Ablation\.

Table[19](https://arxiv.org/html/2605.24041#A8.T19)compares IRNO trained with BatchNorm, LayerNorm, and GroupNorm on TR\-2D \(TFNO base,K=4K=4\)\. All three normalization choices achieve consistent error reduction, confirming that the iterative refinement mechanism is robust to this architectural choice\.

Table 19:Normalization ablation on TR\-2D \(TFNO base\)\. VRMSE and improvement over base operator atk=4k=4andk=8k=8\.
##### F\-Adapter vs\. IRNO\.

Table[20](https://arxiv.org/html/2605.24041#A8.T20)compares F\-Adapter\[Zhanget al\.,[2025](https://arxiv.org/html/2605.24041#bib.bib58)\]and IRNO \(FNO base\) on Active Matter\. F\-Adapter targets minimal\-parameter adaptation, achieving2\.31%2\.31\\%VRMSE reduction with gains concentrated in mid\- and high\-frequency bands\. IRNO trades additional training cost for substantially larger gains \(50\.73%50\.73\\%overall\) with broad spectral improvements across all frequency bands\. Figure[20](https://arxiv.org/html/2605.24041#A8.F20)visualizes the spectral error profile of F\-Adapter and its per\-frequency improvement over the FNO baseline\.

Table 20:F\-Adapter vs\. IRNO on Active Matter \(FNO base\)\. Frequency\-band improvements are relative to the base FNO\.![Refer to caption](https://arxiv.org/html/2605.24041v1/fig/spectral_error_fadapter.png)Figure 20:Spectral error vs\. radial frequency for F\-Adapter \(left\) and per\-frequency improvement over FNO baseline \(right\) on Active Matter\.