CFO: Learning Continuous-Time PDE Dynamics via Flow-Matched Neural Operators
Authors: Xianglong Hou
Xinquan Huang
Paris Perdikaris
Venue: ICLR 2026 Poster (2026)
Keywords: Neural PDE Solver
Operator Learning
Flow Matching
Machine Learning
Long Rollouts
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: Neural operator surrogates for time-dependent partial differential equations (PDEs) conventionally employ autoregressive (AR) prediction schemes, which accumulate error over long rollouts and require uniform temporal discretization.
We introduce the Continuous Flow Operator (CFO), a framework that learns continuous-time PDE dynamics without the computational burden of standard continuous approaches, e.g., neural ODE. The key insight is repurposing flow matching to directly learn the right-hand side of PDEs without backpropagating through ODE solvers. CFO fits temporal splines to trajectory data, using finite-difference estimates of time derivatives at knots to construct probability paths whose velocities closely approximate the true PDE dynamics. A neural operator is then trained via flow matching to predict these analytic velocity fields. This approach is inherently time-resolution invariant: training accepts trajectories sampled on trajectory-specific, non-uniform time grids while inference queries solutions at any temporal resolution through ODE integration. 
Across four benchmarks (Lorenz, 1D Burgers, 2D diffusion-reaction, 2D shallow water), CFO demonstrates superior long-horizon stability and remarkable data efficiency. 
CFO trained on only 25% of irregularly subsampled time points outperforms autoregressive baselines trained on complete data, with relative error reductions up to 87%. 
Despite requiring numerical integration at inference, CFO achieves competitive efficiency, exceeding AR accuracy using only half the AR rollout step budget, while uniquely enabling reverse-time inference and arbitrary temporal querying. Code is available at https://github.com/shannon-hou/CFO_official

Published as a conference paper at ICLR 2026
CFO: LEARNING CONTINUOUS-TIME PDE DYNAM- ICS VIA FLOW-MATCHED NEURAL OPERATORS
Xianglong Hou
Graduate Group in Applied Mathematics and Computational Science University of Pennsylvania Philadelphia. PA 19104
xlhouésas upenn. edu
Xinquan Huang, Paris Perdikaris
Department of Mechanical Engineering and Applied Mechanics
University of Pennsylvania
Philadelphia. PA 19104
[huang26, ,6, pgp) 2 seas.. www.u.ennnnn..m edu
ABSTRACT
Neural operator surrogates for time-dependent partial differential equations (PDEs) conventionally employ autoregressive (AR) prediction schemes, which accumulate error over long rollouts and require uniform temporal discretization. We introduce the Continuous Flow Operator (CFO), framework that learns continuous-time PDE dynamics without the computational burden of standard continuous approaches, c.g., neural ODE. The key insight is repurposing flow matching to directly learn the right-hand side of PDEs without backpropagat- ing through ODE solvers. CFO fits temporal splines to trajectory data, using finite-difference estimates of time derivatives at knots to construct probability paths whose velocities closely approximate the true PDE dynamics. ral operator is then trained via flow matching to predict these analytic velocity fields. This approach inherently time-resolution invariant: training accepts tra- jectories sampled trajectory-specift queries solutions at any temporal resolution through ODE integration. Across four benchmarks (Lorenz, 1D Burgers, 2D diffusion--eaction 2D shallow water), CFO demonstrates superior long-horizon stability and remarkable data efficiency. CFO trained only 25% of irregularly subsampled time points outperforms autore- gressive baselines trained on complete data, with relative error reductions to 87%. Despite requiring numerical integration at inference. CFO achieves compet- itive efficiency, exceeding AR accuracy using only half the AR rollout step budget, while uniquely enabling reverse- everse-time inference and arbitrary temporal querying. Code is available at ? github..com/saanonn Shw.s/sooonooonooooOOO official
neu-
non-uniform time grids while inference
INTRODUCTION
Time-dependent partial differential equations (PDEs) are fundamental to modeling dynamical phe- nomena across the physical, biological, and engineering sciences. Neural PDE solvers have emerged as promising paradigm for learning mappings between function spaces and fast simulation of com- plex PDEs (Azizzadenesheli et al. 2024). However, current approaches face critical limitations: toregressive (AR) methods suffer from error accumulation over long horizons and require uniform time grids Gonzzlez et al. 2020); spatio-temporal methods scale poorly with the space- time volume and struggle with causality (Wang et al.) 2021); continuous approaches like Neural ODEs (Chen et al., (2118) require expensive backpropagation through ODE solvers.
In this study, we propose the Continuous Flow Operator (CFO). continuous-time time neural framework that circumvents these limitations through key insight: repurposing flow match-

---

Published as a conference paper at ICLR 2026
ing (Lipman et al., 2022 Albergo al. 2023) typically used for generative modeling, to learn continuous comnnuuuuuutuuuuueue dynamics without backpropagating through ODE solvers.
We directly learn the right-hand side dynamics of PDE system by matching it with the analytic velocity of spline-based interpolants over trajectories. In doing so, CFO can be trained on arbi- trary, non-uniform time grids to fully exploit available temporal data and queried at any temporal resolution during inference, while retaining the training efficiency of discrete methods.
Our contributions are summarized as:
A continuous- framework for learning PDE dynamics. We propose novel method that learns the inherent continuous-time right-hand side (RHS) by training neural operator to match the analytic velocity (i.e., the time derivative) of spline-based interpolants, thereby avoiding costly backpropagation through ODE solvers via direct derivative supervision.
Time-resolution invariance. The proposed framework inherently handles irregular, trajectory-specific time grids during training and supports arbitrary-resolution inference, which is unattainable with standard autoregressive methods.
State-of-the-art performance with minimal data. Across four benchmarks, CFO achieves superior long-horizon stability; models trained on 25% irregular subsamples out- perform full-data autoregressive baselines, with relative error reductions up to 87%.
2 RELATED WORK
Neural PDE Solvers and Operator Learning. Neural operator learning has emerged as a power-
ful paradigm for solving PDE families by learning mappings between function spaces. Pioneering
frameworks include DeepONet (Lu et al. 2021) with branch-trunk architectures and FNO (Li et al.
2020a) with spectral convolutions. Recent extensions handle irregular geometries (Geo-FNO (Li et al. 2023), multipole graph operators (Li et al. (020b)), multi-scale structures (WNO (Tripura Chakraborty 2022)), improved depth (U-NO (Rahman et al. 2022 and irregular sampling (Galerkin Transformer (Cao) 2021), CViT (Wang et al., 2024)).. For time-dependent PDEs, au-
toregressive prediction suffers from exposure bias and cumulative error, while treating time as a spatio-temporal coordinate scales poorly and complicates causality enforcement (Wang et al. 2021
Neural ODEs for Continuous-time Modeling. Neural ODEs (Chen et al. 2018 and variants (Kidger et al.] 20200 Dupont et al.] (209)) provide continuous-tine modeling for PDE systems (Sholokhov et al. 2023] Wen et al. (023). While theoretically elegant, these methods require back- propagation through ODE solvers via the adjoint method, remaining computationally expensive and memory-intensive even when operating in latent spaces.
Generative Models for PDE Simulation. Recent generative approaches treat solution fields as samples from learned distributions. DiffusionPDE (Huang et al. 2024) uses EDM-style diffusion
(Karras et al. 2022) with physics guidance during sampling; CoCoGen (Jacobsen et al. 2025)
tends score-based models (Song et al. (020)) with physically-consiseent sampling; PBFN M (Baldan)
al. (025) integrates physies directly into flow matching (Lipman al. (022) training via conflict-
free gradient updates. While effective, these methods typically require explicit PDE expressions and
suffer from computational overhead because of multi-step sampling and space-time generation.
Our Approach. CFO bridges continuous-time modeling with computational efficiency through a key innovation: using flow-man matching objectives to learn continuous-tine vector fields directly from trajectory derivatives, avoiding the costly ODE solver backpropagation required by Neural ODEs. Unlike generative PDE solvers that need explicit differential operators, CFO implicitly encodes dynamics through probability paths aligned with PDE evolution, maintaining temporal causality without requiring governing equations.
Three advantages distinguish our framework: (1) Model agnosticism compatible with any neural operator architecture (e.g., FNO, U-Net (Ronneberger et al. 2015)): : (2) Time-resolution invariance
handles irregular, trajectory -speciffi time grids during training while supporting inference; (3) Computational efficiency- achieves autoregressive training speeds while maintaining continuous-time benefits.
2

---

Published as a conference paper at ICLR 2026
While Zhang al. (2024) apply flow matching with piecewise linear interpolation to clinical time series and treat trajectories as generic sequences, CFO instead exploits PDE structure: by fitting high-order splines with inite-difference derivative estimates at knots, construct probability paths whose velocities closely match true PDE dynamics. This physics-aware design keeps the learned flow in high-poobability regions consistent with governing equations, enabling accurate predictions even when trained on only 25% % of irregularly sampled time points.
PRELIMINARIES
Method Lines. The method of lines (MOL) is a numerical method for solving time-dependent PDEs. The core idea is first to discretize the spatial domain, transforming the PDE into a large system ODEs time, which can then be solved using a standard numerical integrator. Formally, after spatial discretization (e.g., grid), the PDE in Eq. (2) becomes system of ODEs:
d (n (un, uy (0) 20, dt
where un(() is vector representing the solution the spatial grid at time t, and N₁. is finite- dimensional approximation of the spatial operator N (c.g., via finite differences or spectral meth- ods). However, when N is unknown or difficult specify, instead of relying on fixed, predefined discretization, Ni. we learn a neural operator Ng that directly approximates the continuous cottuuuuuu-teueeeeee dy- namies from data. This allows us define continuous-time ODE system that can be integrated with arbitrary time steps, inheriting the flexibility of MOL while avoiding the limitations of discrete-time autoregressive models.
While MOL provides conceptual basis for learning a continuous- operator Ng. a naive imple- mentation that matches solution trajectories would require backpropagation through an ODE solver, which is computationally expensive, as seen in Neural ODEs (Chen et al. 2018).
Generative modeling frameworks like flow matching (ipmann et al. 2022) and stochastic inter- polants (Albergo et al.] (2023) have put forth an efficient alternative. These methods learn continuous-time vector field by directly regressing it against the velocity of predefined probability path connecting two distributions. This approach avoids costly ODE integration during training, key feature we leverage in CFO. Flow matching is special case of stochastic interpolants with no noise term. We briefly review these concepts below.
Stochastic Interpolants and Flow Matching. Stochastic interpolants (Albergo et al., 2023) and flow matching (Lipman et al. (202)) provide framework for learning continuous-time vector field that transports one probability distribution another. The core idea is define a probability path Pt that connects source distribution po target p1- A stochastic interpolant is typically defined as 4t = s(t, Xo, X1) + (,)),, where Xo. X1 are drawn from a probability measure q 21), Z ~ AOO, I) is independent noise, and the terms satisfying boundary conditions:
s(0, xo, x1) = Jo, s11, x1) = xı,
70) 711 = 0.
The key insight is that the vector field x)) I I æ] that generates this path can be learned directly by regressing neural network vo (t, x) against the analytic velocity _ Le of the predefined process:
((0) = Et [Iro(t, 1.)
1 dt
(I)
This regression objective avoids costly ODE integration during training. During inference, the learned vector field by defines an ODE, di X, = (t X²), which can be solved with a numerical integrator to map samples from Po P1- CFO adapts this principle by constructing more sophisti- cated probability path based on splines that incorporate derivative information and interpolate entire solution trajectories, rather than only matching the endpoints.
3

---

Published as a conference paper at ICLR 2026
u(G))
Referencce u(t)
u(ss)
Spline interpolant (()
m(ta)
Nolr) dr
15
Beewer fiew vebeace
No
1/
8,7(0) FD ((u(4))L)))
min (oo (
Nold -))
(a) Training
(b) Inference
Figure 1: Overview of the Continuous Flow Operator (CFO) framework. (a) Training: flow match- ing on spline path. For each trajectory with snapshots (u(t:)) ]. we fit a temporal spline s(t) that interpolates the data and matches finite-difference derivative estimates at knots. The spline's s analytic time derivative (()(() provides exact velocity targets for training a neural operator No via the flow matching objective- ODE integration required during training. (b) Inference: continuous-time rollout. The trained operator (t, u) defines continuous vector field dg = (ot uæ). Given initial condition u(t), we compute u(!") by numerical integration; reverse-time inference integrates backward.
4 CONTINUOUS FLOW OPERATOR
4.1 METHODOLOGY
Consider dime-dependent PDE domain [0, Ω2 C R Rd: 0, = (N(u(t,))) x)), u(0,z) = uz)), I TE te [0,7],
(2)
where A is an unknown spatial differential operator and u1(x) is the initial condition.
We suppress spatial coordinates and write u(t) = u(t,)) :) E V. For sampled times T : (t we write the restriction u := (u(to)....,u((( When T is clear from context, abbreviate the snapshot vector to UT.
Notation.
Problem Statement. Let M trajectories (t)) be generated from initial conditions ( For each trajectory j, observe snapshots at (possibly irregular) times T (t)) and write u' == W = (w' (1)).... , w (1x,)) ))- Given the dataset (u {u'} our goal is to learn an operator that maps from any initial condition un to the solution trajectory u(t) for any E [0, For exposition, we drop the superscript j and set the horizon to
Probability Path Formulation. We define a temporal spline s(1; u) that interpolates the snap- shots:
st4; u) = u(t:):
i= 0.....,
The (multi-point) stochastic interpolant is
[(t; u) s((; + (t)
tE [0,1], 1],
(4)
I). independent of u. In the MOL one discretizes space to obtain NN, ~ A with dtu(t) = ~N(u(t)). When N is unavailable, we take a data-driven route: (i) estimate 0,((4) from u using finite-difference stencils along time; (ii) choose the spline s(t; u) so that its derivatives at the knots match these estimates. This allows to control the derivative accuracy via the stencil order and the path smoothness via the spline degree.
with boundary conditions (t.) 0 for O..., and
Finite Difference Approximation. The order of of points used (Strikwerda) 2004). Below are examples for the first and second time derivatives. Let d, Rs (l,( and a; Rs O.e(();) denote derivative estimates at time 4, Define the steps h;; := 4 -1; and AL == max; h4
stencils depends on the number
A two-point stencil is defined as d, v(((+1)-w(t+) which is first-order accurate: 0,u(ts) di + O (At).. A three-point stencil provides second-order accuracy for the first derivative and first-order
4

---

Published as a conference paper at ICLR 2026
accuracy for the second derivative irregular grids. For an interior point t,, h4-1
d(t:)= hi (hi- + hss) u(ti-1) + h4 h4-1 u(ti) 1:(1i=1 + hi)
h4
u(te+1)
14-1h4
(5)
2
u(ti+1) alti+t)-u(ee u(ti) hi
u(ti) u(bs-1))
all4) h1-1 h4
14-1
These yield 0.u(4)) dj + ((Ar2) and On.u(1s) = a; + 0(A1). Endpoint estimates can be obtained with standard one-sided stencils. Higher-order accuracy is achievable with more points (see Appendix A.I.I details).
Segment Interpolation. With the derivative estimates at the knots, we construct the spline s(;; u) segment-wise. On each interval [f;, [i+1], J. we a polynomial that matches the values and estimated derivatives the endpoints. This approach, known as Hermite interpolation, allows us to control the spline's smoothness and the accuracy of its derivatives at the knots. For spline to be C continuous, we must impose conditions on derivatives up to order at each endpoint, which requires polynomial of degree at least 2s + each segment.
Representative Cases.
We consider two cases: linear splines (C") and quintic splines (C²).
LINEAR SPLINE. The simplest choice is a linear spline, which only interpolates the values u(4;) and u((++1). On each segment [1, 11+ the spline is:
s(1 u) = f+ - u(l.) h;
t- u((;+1).
(6)
The resulting spline continuous but not differentiable at the knots. Its one-sided derivative at t; matches the first-order forward-difference estimate, yielding 0,s(t; u) (u((t))) + O(Af).
QUINTIC SPLINE. To incorporate higher-order derivative information, we use a quintic spline. On each segment, this spline matches the values (ut,, Wt, .), first derivatives (de, de,a)) )). and second derivatives (0t;; (t:-1) ) at both endpoints. With the three-point finite difference stencils from Eq. (5)
the spline derivatives the knots satisfy:
(ts(t; u) = (u((t)) + O(At?). Sus(t; u) (,,(t((,)) + O(at).
This construction yields a globally C2 smooth path that more accurately reflects the underlying dynamics. The explicit form is provided in Appendix(A.13. A.1122
Unlike global methods such as natural splines, our segment-wise Hermite interpolation yields closed-form coefficients. It also allows for flexible choices of basis functions, smoothness, and derivative accuracy based on the physics prior. Compared to low-order splines, the higher-order spline here (also with high-order assigned derivatives) benefits from more accurate derivative ap- proximations and higher smoothness: (i) high-order stencils make the probability flow match the exact flow of the physical dynamics as closely as possible, (ii) greater smoothness improves training stability and long-horizon accuracy. We show these effects in Sec. 5.1
In practice, we use quintic splines by default, since they naturally capture acceleration, which ap- pears in many physical laws (e.g., Newton's second law, wave equations). We do not go beyond quintic because of diminishing returns: CFO, the total error comes from both the spline approx- imation and the neural operator training, and beyond quintic we observe only minor improvements in derivative accuracy while the network error dominates (see Sec. A.1.3 for detailed discussion).
For the noise term, can also impose regularity on >(f) with splines. For example, >(t) (1--))"((--11- 70 (fk+1 11641-14)00 4g) 2m l I < (++1 with some constant To > ensures Cm-1 continuity and
vanishing noise at the knots. A detailed study can be found Sec. [.2
4.2 TRAINING AND INFERENCE
Training. We learn time-dependent neural operator Nott, u(t)) that matches the velocity of the basede-based stochastic interpolant T(t; u). Similar to the training objective stochastic inter- polants, the training loss is
(()) = Em.t. [[1 - 0, ul]
(7)
5

---

Published as a conference paper at ICLR 2026
where t ~ Unif(O, 1]. denotes sampled trajectory, Z is the Gaussian noise independent of u. [tt; u) the corresponding stochastic interpolant, and (t u) is its analytic spline derivative. This avoids backpropagating through ODE solver.
Inference. After training, N₁ defines a continuous-tims right-hand side. For a new initial condi- tion No. we predict by integrating
d
dt u = Nott, uo)),, (o( = uo:
(8)
with a standard ODE solver; by default, we use 4th-order Runge-Kutta (RK4). See Sec. 5.2 for details of solver choice and the number of function evaluations (NFE).
Resolutin Invariance. CFO is agnostic to time grids: (i) Training accepts trajectories with trajectory-specife and irregular time stamps. (ii) Inference queries (8) arbitrary time resolutions or schedules (including times unseen training). This allows mixing dense and sparse sequences during training and evaluating any target grid at test time.
Reverse-time Inference. Unlike autoregressive next-step predictors, our learned time-dependent vector field induces locally invertible flow map under standard Lipschitz conditions. This makes CFO suitable for inverse tasks such as inferring earlier states from a state observed at time t.. Given complete state u((+) any time l, E [0, 1J. one can recover the state at earlier times < l, by integrating the learned flow (o u()) backward:
- ( u(t.)+ [." No mo(T)) dr = u((t) - uoll 7)) dr,
We empirically validate this reverse-time inference capability on the dissipative Burgers' equation Sec. 5.2
Probability- Transport. For a fixed (possibly irregular) time grid, if NV minimizes the loss function (7). the ODE (8) transports the distribution across the sequence of data marginals each grid time. It generalizes the two-marginal result of Lipman et al. (2022;; Albergo et al. (2023). (see Proposition A. A.2 for details. Empirical results also show that the induced probability flow closely tracks the true dynamics even with heterogeneous time grids across trajectories, especially for quintic CFO.
EXPERIMENTS
In this section, we evaluate the proposed CFO framework across four canonical benchmarks: Lorenz, ID Burgers', 2D diffusion-reaction. and 2D shallow water equations. Our primary focus is on long-horizon rollouts under irregular subsampling, where compare CFO against autoregres- sive basclines trained with teacher forcing. addition, we perform ablation studies to assess the impact of the noise schedule, inference efficiency, long-horizon temporal extrapolation, and back- bone architectures. All time series are normalized to the range [0, 1]. More results are shared the Appendix A.4
5.1
MAIN RESULTS
We evaluate CFO's time-resolution invariance by training on irregularly subsampled time grids and testing on the full-resolution grid. Models are trained using per-trajectory random time grids with keep ratios of 100%. 50%, and 25%. Autoregressive baselines, which require uniform time steps, are trained only on the full (100%) grid; for autoregressive, we try multiple architectures (including the same backbone as CFO) and report the best. All models learn the mapping from the initial condition to the entire trajectory. Table Treports the rollout accuracy (relative Le error, mean = sd), Notably, Quintic CFO trained on just 25% of irregularly sampled time points still outperforms AR models trained on full-resolution data, achieving relative error reductions of 24.6%, 78.7%, 87.4%, and 82.8% on the four benchmarks, respectively. This highlights the CFO's ability to leverage all available data, regardless of sampling density or regularity, and to generalize effectively to arbitrary time resolutions.

---

Published as a conference paper at ICLR 2026
Table Time-resolution invariant results across four benchmarks: CFO robust to irregular time subsampling. Training uses per-trajectory random time grids at keep ratios 100%, 50%, and 25%; Testing full resolution. Values are relative L2 error, (mean + sd). Bold indicates the best result in each sub-row; "_" means AR not suitable for irregular grids.
Dataset Sampling
Autoregressive
Linear CFO
Quintic CFO
100%
(9.04 # 2.80) 10-2
(6.42 10-2
(4.53 0.68) 10-2
Lorenz
50%
(7.85 x 10
(4.75 0.78) 10
25%
(9.39 x 10
(6.82 + 0.48) 10-2
100%
(3.34 0.14) 10-2
(5.75 0.72) 10-3
(5.89 0.79) 10-3
Burgers 50%
(8.91 # 1.30) x 10-3
(6.97 = 0.73) 10-3
25%
(1.04 x 10
(7.09 1.10) x 10-3
100%
(4.23 # 0.59) 10-1
(4.35
10-2
(4437 x 10-2
DR
50%
(6.88 + 0.49) 10-2
(6.19 0.35) 10-2
25%
(7.25 + x 10
(5.32 10-2
100%
(9.04 + 0.36) 10
+ (..5)) x 10-3
(4.56 0.38) 10-3
SWE
50%
(7.64 + 0.52) x 10-3
(6.57 + 0.66) 10-3
25%
(1.69 + 0.38) x 10-2
(1,55 + 0.22) x 10-2
Lorenz System. We first consider the Lorenz system (Lorenz] 2017). 3D chaotic ODE with spatial-derivative dynamics. In this simple setting, the state is a vector u(t) E R representing system of ODEs rather than PDE with spatial dependence:
t (3 x), à = x(p - 2) y. == 82,
with (a, p) (10, 28, 8/3). We an MLP for both CFOs and baseline autoregressive methods.
Since the Lorenz vector field is derivative-free, we directly evaluate how well the spline-implied velocity stt; u) approximates the true dynamics. We analyze this by regularly downsampling trajec- tories to resolutions of 2 and 4Ar. pigtte the average relative L2 error across normalized trajectory segment. For fixed step size, the quintic spline achieves uniformly lower error. The lin- ear spline is only comparable near the segment's midpoint, where its secant slope approximates a central difference; decreases, both splines improve, but the quintic spline converges faster. Table confirms these con- vergence rates at the endpoints: the linear spline is first-order accurate (O(A)), while the quintic spline is second-order (O(Ar))) matching finite-diffen lifte-difference theory.
error increases near endpoints, where the steneil becomes one-sided. As At
In addition, we study the training stability and long-horizon accuracy under irregular subsam- pling. The spline-implied velocity error (Figure 2) reflects how closely the induced flow aligns with the true dynamics; lower error ensures the learned flow remains in gions. Thus, the high-order spline has two ad- vantages: (i) it improves long-horizon accuracy and is more data-efficient. When only a frac- tion of snapshots available, the Quintic CFO maintains low rollout eтror and degrades more gracefully than the Linear CFO, demonstrating stronger robustness to data scarcity under ir- regular subsampling (Table [). (ii) enhances training stability. As shown in Figure 3 Quin- tic CFO achieves lower training loss and faster convergence compared to Linear CFO, particularly in tasks with irregular time sampling.
Quinti (AI
Luoe: (At)
Relative 12 Errux
Quinte (2.3+)
Linea (2./)
Quintoc (L3))
re-
Liner (4.M)
0.02
0.00
Seymuntt ratio
Figure 2: Spline fitting error on the Lorenz sys- tem. Average velocity- -feld approximation error for quintic and linear splines at different time-step sizes (At, 2At, 4At), with At = 0.005s.
ID Burgers' Equation. We test the CFO on the ID Burgers' equation, a standard benchmark for nonlinear dynamics. The setup follows Wang al. (2021). CFO uses 1D U-Net, while the

---

Published as a conference paper at ICLR 2026
Oparni (ogbte)
0.12
1901 borogaa
10%
1.0
0.133
1.20
4.30
4.2
1800
2002
2040
3000
(a) 100% keep ratio
(b) 50% keep ratio
(c) 25% keep ratio
Figure 3: Training comparison between linear and quintic CFO on the Lorenz system. Training loss
(10k-60) 60k iterations) for linear and quintic CFO under different keep ratios: (a) 100%, (b) 50%, (c)
25%.
Wlep 20
alep
50
alep 75
vlep 1180
66650
lof.
Phod. All Eirel.
6.605-
0.000
-R015
-0050
Figure 4: Worst-case prediction for Burgers' equation. Quintic CFO (trained on 25% irregular data) accurately predicts the full trajectory for the test sample with the highest error. Shown from left to right: solutions at the 25, 50, 75, and 100th time steps.
autoregressive baseline is an MLP, which performed best among several architectures tested (U-Net, MLP, FNO, see Appendix A.4.5.
As shown in Table [ Quintic CFO trained on just 25% of irregularly sampled time points achieves a 78.7% relative еттог reduction compared to the autoregressive baseline trained on full data. Figure visualizes the rollout for the test sample with the highest error (0.0307 relative L2 error); even in this worst-case scenario, the prediction remains highly accurate.
2D Inffusion-Reactinnn Equation. We evaluate CFO on the 2D diffusion-reaction (DR) bench- mark from PDEBench (Takamoto al. 2022). We use 2D U-Net for both CFO and the autore- gressive baseline. As shown in Table CFO significantly outperforms the autoregressive model. Quintic CFO trained on just 25% of irregularly sampled time points achieves an 87.4% relative er- ror reduction compared to the autoregressive baseline trained on full data. Figure 5 visualizes long- horizoorizon rollout, showing that CFO accurately captures the complex pattern evolution, while the autoregressive model quickly diverges.
2D Shallow Water Equation. We evaluate CFO the 2D shallow water equations (SWE) from PDEBench (Takamoto et al.] 2022). CFO uses a 2D U-Net: the strongest autoregressive baseline Diffusion Transtormer (DiT) architecture. As shown in Table I quintic CFO (with 25% data keep ratio) attains 82.8% lower relative error than the autoregressive model trained on the full dataset. Figure 6 compares predictions over time and shows that, even on its worst sample, CFO substantially outperforms the autoregressive model.
5.2 ABLATION STUDY
Impact of ODE solvers and function evaluations. We study how CFO's inference accuracy and computational cost vary with the ODE solver and the number of function evaluations (NFE). We compare Buler, Heun, and RK4. Experiments show that higher-order numerical integrators achieve lower error for a fixed NFE. Table fixes the solver to RK4; additional solver sweeps are provided in Tables g and 9 CFO's continuous-time formulation supports flexible solver choice and NEE budgets, unlike autoregressive (AR) rollouts that use fixed steps. For RK4, across all benchmarks find that: with only 50% of the AR NFE budget, CFO already outperforms AR; increasing the
8

---

Published as a conference paper at ICLR 2026
budget 200% yields further gains and stable high accuracy, while additional NFE up 400% brings diminishing returns with only negligible improvements. To connect NFE to actual runtime, also report wall-clock inference time and accuracy for AR and CFO+RK4 on identical hardware with the same spatial backbone in Table TO and Figure g It confirms that the theoretical reduction in NFE translates directly wall-clock time savings.
Influence of NFE on Quintic CFO (RK4 solver) accuracy (Relative L2 Error) at the final time step, compared the autoregressive baseline.
Table 2:
NFE (% AR)
Lorenz
Burgers
DR
SWE
AR (100%)
0.1481
0.0647
0,3850
0. 0.1048
50%
0.0919
0.0090
0.0876
0.0698
100%
0.0655
0.0089
0.0718
0.0154
200%
0.0674
0.0090
0.0608
0.0060
400%
0.0645
0.0089
0.0593
0.0061
Impact of Noise Schedule. We investigate the impact of the noise schedule model perfor- mance. The noise schedule (t) = 7o 14)-4114))mm
[t (t= for E [tx, fx+1]. where m controls smooth-
ness and "o controls magnitude. Figure shows the impact of noise magnitude "P on the SWE and DR benchmarks. Performance is stable for small noise levels (7o E [0, 1)= but degrades as the noise magnitude increases. For our experiments, we set m = 3 to ensure the noise schedule is C continuous, matching the smoothness of the quintic spline.
Comparison with Other Baselines. Besides AR, we also compare CFO with two other base- lines. The first is continuous-time neural ODE model (Chen et al. 2018) trained with teacher- forcing strategy. CFO attains higher accuracy with substantially less training time (table 3 and the continuous-time supervision of CFO leads more stable optimization (Figure 9.. Second, compare against PDE-Refiner (Lippe et al. 2023), a strong baseline for long-horizon rollouts. CFO outperforms PDE-Refiner on hong-horizon predictions while using fewer computational resources, since PDE-Refiner requires diffusion-stylee denoising simulations during training.
Reverse-time Inference. We study CFO's performance of backward integration the dissipative Burgers' cquation. For Figure 10 a). we fix a terminal time t. and integrate the flow backward varying backward horizons within [0, t.). The panel reports how the error accumulates with increasing backward horizon under different noise levels added to the terminal state. Figure (b)b) visualizes initial-condition recovery by integrating backward from different terminal times to t 0. Despite the ill-posed nature of inverse problems for dissipative PDEs, CFO yields reasonable reverse-time reconstructions over short backward horizons before errors grow.
Temporal Extrapolation Beyond the Training Horizon. To assess temporal extrapolation, we train CFO only on the first half of each trajectory [0, T/2), and evaluate rollouts over the full horizon [0,7] (Table [T]. The error accumulation curves in Figure IT show that, for Lorenz, Burgers, and DR, the extrapolation errors remains nearly identical to full-horizon training results, indicating that CFO has learned the underlying dynamics rather than memorizing the trajectories. For SWE, the error increases 2.28 10-2, but this still remains approximately 4x lower than the full-horizon training autoregressive baseline.
Neural Network Architectures. We assess the effect of the backbone on CFO by swapping the spatial operator while keeping the training recipe fixed. For ID Burgers we use an FNO (Li et al. 2020a); for 2D diffusion-reaction and 2D shallow water (SWE) we use the DiT (Peebles & Xie 2023) Architecture. Results in Table shoo comparable accuracy for backbones, with only modest variations. This indicates that CFO is largely architecture aggnssicc: temporal continuity is handled by CFO, while the spatial inductive bias can be chosen to suit resolution, memory, or hardware constraints.

---

Published as a conference paper at ICLR 2026
Table Comparison of CFO with other baselines: Neural ODE Lippe et al. 2023).
CChnn et al. 2018 and PDE-Refiner
Equation
Method
Rel. L2 Error
Training Time (s/batch)
Lorenz
Neural ODE
0.101
0.133
CFO
0.0453
0.00350
Burgers
Neural ODE
0.0275
3.38
CFO
0.00589
0.00920
DR
PDE-Refiner
0.125
1.38
CFO
0.044
0.40
SWE
PDE-Retiner
0.093
1.40
CFO
0.005
0.40
CONCLUSIONS AND FUTURE WORK
We introduced CFO,, a continuous-time neural operator that learns PDE dynamics by matching the analytic velocity of spline-based interpolants. By repurposing flow matching to avoid ODE solver backpropagation during training. CFO achieves unique combination of capabilities: training on irregular, trajectory-specifi time grids; inference at arbitrary temporal resolutions; and accurate predictions from severely subsampled data. Across four benchmarks, models trained on only 25% of irregularly sampled time points outperform fully-trained autoregressive baselines, demonstrating that physics-aware probability path design can overcome data scarcity.
Several directions warrant future investigation. Sputial resolution: While CFO handles temporal irregularity, coupling with mesh-agnostie operators would enable full spatio-temporal resolution invariance. Inference acceleration: Though CFO matches autoregressive efficiency with modest function evaluations, distilling the learned flow into consistency models (Song et al., (023) could enable single-step predictions. Physics integration: When governing equations partially known, augmenting training with physics-informed constraints (Huang et al.] 2024) could improve back- ward integration stability. Adaptive splines: Current fixed-order splines face trade-off between computational cost and accuracy; learned splines with curvature regularization (Luo et al. 2025) could adapt smoothness local dynamies.
CFO demonstrates that continuous-tin modeling need not sacrifice computational practicality. By bridging the efficiency of discrete methods with the flexibility of continuous formulations, it opens new possibilities for learning from real-world data with irregular temporal measurements, a critical step toward practical neural PDE solvers.
REPRODUCIBILITY STATEMENT
To facilitate reproducibility, we provide extensive implementation details, ablation studies, and the- oretical clarifications throughout the main paper and appendix. Section describes the proposed model architecture and training procedure, while Section outlines the experimental setup across all benchmark tasks. The theoretical results are included Appendices [.J and the full al- gorithms are summarized in Appendix A.3 All hyperparameters and data preprocessing steps are detailed in Appendix A.4 Code is available at httptttss///gtthhee hww.:/uwb...com/aaannnnhhh/e//C- official
ACKNOWLEDGMENTS
This work was supported by the US Department of Energy under the Advanced Scientific Computing Research program (grant DE-SC0024563), and the Penn Fellowship program.
REFERENCES
Michael Albergo, Nicholas M Boffi, and Eric Vanden-Eijjdee Stochastic interpolants: unifying framework for flows and diffusions. arXiv preprint urXiv:2303.08899 2023.
10

---

Published as a conference paper at ICLR 2026
Kamyar Azizzadenesheli, Nikola Kovachki, Zongyi Li, Miguel Liu-Schiaffni,, Jean Kossaifi, and Anima Anandkumar. Neural operators for accelerating scientific simulations and design. Nature Reviews Physics, 6(5):320-328, 2024.
Giacomo Baldan, Qiang Liu, Alberto Guardone, and Nils Thuerey. Flow matching meets pdes: A unified framework for physies generation. arXiv preprint arXiv:2506.0805 2025,
Shuhao Cao. Choose a transformer: Fourier or galerkin. Advances neural information processing systems, 01124924-249492 2021.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
Emilien Dupont, Arnaud Doucet, and Yee Whye Teh. Augmented neural odes. Advances in neural information processing systems, 32, 2019.
Francis Begnaud Hildebrand. Introduction to numerical analysis. Courier Corporation, 1987.
Jiahe Huang, Guandao Yang, Zichen Wang, and Jeong Joon Park. Diffusionpde: Generative pde- solving under partial observation. arXiv preprint arXiv:2406. 17763, 2024.
Christian Jacobsen, Yilin Zhuang, and Karthik Duraisamy. Cocogen: Physically consistent and conditioned score-based generative models for forward and inverse problems. SIAM Journal on Scientific Computing, 47()):C399C4255, 2025.
Tero Karras, Miia Aittala, Timo Aila, and Samuli Laine. Elucidating the design space of diffusion- based generative models. Advances in neural information processing systems, 35:26565-2653 2022.
Patrick Kidger. On Neural Differential Equations. PhD thesis, University of Oxford, 2021.
Patrick Kidger, James Morrill, James Foster, and Tery Lyons. Neural controlled differential equa- tions for irregular time series. Advances in neural information processing systems, 33:6696-6707, 2020.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, An- drew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. arXiv preprint arXi::2010.08894 2020a.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhat- tacharya, and Anima Anandkumar. Multipole graph neural operator for parametric partial differ- ential equations. Advances in Neural Information Processing Systems, 33:6755-6766, 2020b.
Zongyi Li, Daniel Zhengyu Huang, Burigede Liu, and Anima Anandkumar. Fourier neural oper- ator with learned deformations for pdes on general geometries. Journal of Machine Learning Research, 24(388):1-24 2023.
Yaron Lipman, Ricky TQ Chen, Heli Ben Ben-Hamu, Maximilian Nickel, and Matt Le. Flow matching for generative modeling. arXiv preprint arXiw:2210.02747, 2022.
Phillip Lippe, Bas Veeling, Paris Perdikaris, Richard Turner, and Johannes Brandstetter. Pde-refiner: Achieving accurate long rollouts with neural pde solvers. Advances Neural Information Pro- cessing Systems, 36:67398-6748 2023.
Edward Lorenz. Deterministic nonperiodic flow 1. In Universality in Chaos, 2nd edition, pp- 367-378. Routledge, 2017.
Lu Lu, Pengzhan Jin, Guofei Pang. Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature machine intelligence, 3(3):218-229 2021.
Yan Luo, Drake Du, Hao Huang. Yi Fang, and Mengyu Wang. Curveflow: Gurvature-guides flow matching for image generation. arXiv preprint a508v:2508. 15093, 2025.

---

Published as a conference paper at ICLR 2026
Harris Abdul Majid.
suramoto-ssvvshinskt dataset.
httbs://hugghoofcceccoo
attassts//arrsm wwt.sssssssssssssseeeeeae 2024.
William Peebles and Saining Xie. Scalable diffusion models with transformers. In Proceedings of the IEEE/CVF international conference on computer vision, pp. 4195-4205, 2023.
Md Ashiqur Rahman, Zachary Ross, and Kamyar Arizzadenesheli. U-no: U-shaped neural oper-
ators. arXiv preprint arXiv:2204. 1277, 2022.
Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-net: Convolutional networks for biomedi- cal image segmentation. In International Conference on Medical image computing and computer- assisted intervention, pp. 234 241. Springer, 2015.
Alvaro Sanchez-Gonvale,, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovee, and Peter Battaglia. Learning to simulate complex physics with graph networks. In International conference on machine learning. pp. 8459-8468. PMLR, 2020.
Aleksei Sholokhov, Yuying Liu, Hassan Mansour, and Saleh Nabi. Physics-informed neural ode (pinode): embedding physics into models using collocation points. Scientific Reports, 13(1): 10166, 2023.
Yang Song, Jascha Sohl-Dickstein, Diederik Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:201 1. 13456, 2020.
Yang Song, Prafulla Dhariwal, Mark Chen, and Ilya Sutskever. Consistency models. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scar- lett (eds.), Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pp. 32211-32253 PMIR, 23-29 2023.
John Strikwerda. Finite difference schemes and partial differential equations. SIAM, 2004.
Makoto Takamoto, Timothy Praditia, Raphael Leiteritz, Daniel MacKinlay, Francesco Alesiani, Dirk Pflüger, and Mathias Niepert. Pdebench: An extensive benchmark for scientific machine learning. Advances in Neural Information Processing Systems, 15:1566-161,, 2022.
Tapas Tripura and Souvik Chakraborty. Wavelet neural operator: neural operator for parametric partial differential equations. arXiv preprint htt.:22205.0199/ 2022.
Ch Tsitouras. Runge-kutta pairs of order 5 (4) satisfying only the first column simplifying assump- tion. Computers & mathematics with applications, 62(2):770-775, 2011.
Sifan Wang. Hanwen Wang, and Paris Perdikaris. Learning the solution operator of parametric par- tial differential equations with physics-informed deeponets. Science advances, 7(40):eabi8605, 2021.
Sifan Wang, Jacob Seidman, Shyam Sankaran, Hanwen Wang, George Pappas, and Paris Perdikaris. Cvit: Continuous vision transformer for operator learning. arXiv preprint 2:Xiv:2405.13900 2024.
Tianshu Wen, Kookjin Lee, and Youngsoo Choi. Reduced-order modeling for parameterized pdes via implicit neural representations. arXiv preprint arXiv:2311.164100 2023.
Xi Nicole Zhang, Yuan Pu, Yuki Kawamura, Andrew Loza, Yoshua Bengio, Dennis Shung, and Alexander Tong. Trajectory flow matching with applications to clinical time series modelling. Advances in Neural Information Processing Systems, 37:107 198-107224, 2024.
12

---

Published as a conference paper at ICLR 2026
A APPENDIX
A. SPLINE CONSTRUCTION
This section details the two steps of our spline construction: (i) estimating derivatives the knots and (ii) piecewise interpolation between knots. We provide general finitediifferecce framework for derivative estimation on nonuniform grids, and we give the closed-form quintic Hermite interpolant on each interval.
The results follow standard numerical references; see, e.g., (Striwwerda 2004 Hildebrand 1987).
A.1.1 FINITE DIFFERENCE STENCIL ON NONUNIFORM GRID
Let (T) interval T C IR and let z0, #1, .... Fm be pairwise distinct with point of interest zo. Set A == Iz Zo and h := maxosism (A Fix k E {0, m.]
Define the Vandermonde matrix
P:
[Pr Pril rJ
Prj (xj - mo)' = Aj.
Let ej E Rm+ denote the k-th standard basis vector. We choose weights w (Wo...., Wm) T by
solving
P w k! Ck=
Then the 1o-+1)-poinn finite-difference approximation of f(() (ro)
((4) = Σ U1 f(zz),
and the error is controlled by
(zo) (D)) f(zo) = Oh
A.1.2 QUINTIC HERMITE INTERPOLATION
Suppose f E ((T, the value, first-, and second-derivatives are 8, d, a and s d', a' at endpoints respectively. Set == - andr == x that matches values, first-, and second-derivatives at both endpoints can written analytically on [1,2"] as
h E [0, 1]. The unique quintic Hermite interpolant (x)
H(x) Hoo(r) + (dhh) H1)(r) + (a) Hzo(T)
+ Hoo + (d' h) (T) + (a' 12) H2i(T),
where the quintic Hermite basis polynomials are
Hoo = 1 - 1 + 157 67
H1o(T) = 6r + 87 - 37
H2o(T) (T) = (7² - 3 + 37 - -0),
Hoi(r) = 1073 - 15r + 675,
(11(T) = -47 7 375
HT1(T) _ 2 + +
These satisfy the endpoint conditions
(0) = 1, (0) 1, Hy (0) 1,
(11 1, (1 1, (1 1,
and all other values/derivatives up order 2 vanish the opposite endpoint.
A.1.3 SPLINE SELECTION
Spline fitting cost. The cost of spline fitting is negligible. The finite-difference approximation and spline construction involve only simple algebraic operations (no iterative solvers). They are performed once per trajectory tional complexity is Osslline order x time segments # spatial grid points) per trajectory, which negligible compared to neural operator backpropagation.
preprocessing step and reused throughout training. The computa-
13

---

Published as a conference paper at ICLR 2026
Practical guidelines for choosing splines.
1. Given the low cost of spline fitting, we use quintic rather than linear splines as the default, since quintic splines capture acceleration (second derivatives), which central to many physical laws (c.g., Newton's ' second law, wave equations). At the same time, we do not go beyond quintic because higher orders bring diminishing returns accuracy.
2. Linear splines adequate when (i) time sampling is very dense (e.g., 100% keep ratio), (ii) the underlying dynamics are known to be non-smooth (c.g.. shocks or discontinuities), where enforeing higher-order global smoothness may be undesirable.
3. Adaptive spline selection is exciting direction for future work. A dynamic approach could involve learning the optimal spline parameters jointly with the operator, potentially incorporating physics-informed priors to adapt the smoothness.
A.2 PROBABILITY-FLOW TRANSPORT
This section provides the theoretical guarantee of CFO's probability -ffoo transport properties.
Assumption A.1. Let u = (u((o))...., u((x)) ~ 4T with = lo <... < IN 1. Assume:
1. (Knot densities) For each k = ....., N. u(ts) admits a strictly positive Lebesgue density qu R2 (0, ).).
2. (Spline in time) The map + s(l; u) is C' [0, 1}, satisfies s(La; u) u(lx), and O.s is dominated in L.: there exists an integrable random variable G(u) E L,2 such that
[[0,((t; u)ll S G(u) for all t E [0, Eu[G(u)] < <o.
where || denotes the Euclidean norm Ied.
3. (Noise schedule) 7 ([O 1]) with (txx) for all 7t)) for $ and
I 80.
4. (Independent Gaussian) Z NO La) independent of
Proposition A.1. Under Assumption A.1 define
u[t: = s(t; u) + (t) O,I(t:u) = dis(t; + (t)
Then:
1. For every t E. [0, 1}, ut admits strictly positive Lebesgue density plt, ). Moreover, for every knot t.. p(ta, x) = (x(x) for x E IR..
2. Let
" 0t xt)x)) = ([ I(t:u) =].
Then for every test function 1 C ((0, x R
[. JEL (0rb(1,0)+ (1) V(b())) o² (1,
(1, x)) p(l, x)) dr dl 0.
Equivalently, the continuity equation O,p + V (po*) = 0 holds in D'((0, 1) x Rd).
Proof. Fix & From the independence of and Z, the law of
Fix
[(::u) s(; u) +
is the convolution of := (sw(s(t; u)) with the non-degenerate Gaussian 7(00,>(t)1)) Hence I [tt; u) admits the Lebesgue density
pt,,)) (Py(t) Pt.
Because (Pyce) is a strictly positive th kemel,, convolution with PP())) regularizes psi it follows that p(l,.) E t (TR)) and p(t,2r) for all d
14

---

Published as a conference paper at ICLR 2026
a a the 7(tx) 0 and s(tx; u) u(fx), hence
1(tx: u) = u(tx)
By Assumption A.1 1), u((x) has strictly positive Lebesgue density yki therefore p((2,-) = 9c(.)
(2) Take test function / E C. ((0, 1) IR))) and define
F(I) == (1, [(t; u)
LE (0,
Since has compact support in t, F(t) = 0 near = and f 1, and thus
L (t) dt F(1) - F(0) 0,
once F" exists and is integrable.
Since the map / + [1; u) is C² by the chain rule,
d
(t, [t: u)) = (t, [(t; + Vettt, ()) at;(tu
dt
Since E C 1) x R.d, O.bb and Vıù are bounded.
Moreover, by Assumption A. 2), there exists G(u) E L such that [t) S G(u) for all t E [0, Hence, for all
00,
[t; u| S (t0,((t; u)ll + (1) s G(u) + /1/1 G,(u, 2).
The random variable G, is integrable under the joint law of (u, Z) since
Eu,2 [G1] = E., [G(u) + >>> Ez || <0.
Therefore,
(,t Vett, ((; (t; s || + V Gr(u, 2).
and the right-hand side is an L' L'(u, Z) dominating random variable independent of t. Thus, domi- nated convergence justifies differentiating under the expectation, yielding
(t) = (t, [(t: + (t [(t; u)) (t;
Next, by the tower property conditioning on [t: u), we have
E [(Veitt,((tuu) 0,1 u)] E V o utt
Since T(t: u) has density p(,, for E (0, 1), we obtain for LE (0, 1);
F((t) { (are(t. ()) Vv(t, u² x) dz.
Integrating over E (0, 1) and using So F'(t) dt yields
S. J (Brv(),)) Vub(t V+((,, x²' (t. ((t)))((t,, p(t, dx dt 0,
which is exactly the weak form of 0.p + V (pr*) = 0 on (0, x Rd.
Proposition A.2. The conditional mean velocity 1 (t, = E(O, u) I(t: u) 2) is the unique minimizer of
C(v) /' E [[| T(t:u)) a.
15

---

Published as a conference paper at ICLR 2026
Proof.
E ((:))
-B - ?
2R [[O [(;,)) -,,())))), a.((;;))]
The equality holds if and only if u(t, ElO,T((;1) [(t:u) = x]]
Remark A.1. By Proposition A.1 the marginals p(t,-) solve the continuity equation with velocity v². Under standard regularity conditions ensuring existence/unigueness the flow for ü(t) (t, u(t)) (e.g. local Lipschitz in x), this ODE transports qo to p(t, :,, and in particular matches the knot marginals qk t
By Proposition A.2 v* is the minimizer of the regression objective Therefore, if the learned operator No matches e" (e.g., at the minimizer of (9), then the forward ODE 8 induces the same probability-llow transport across the knots.
A.3 ALGORITHM
We summarize the numerical procedures used in the main text. Algorithm describes tion of quintic splines from discrete snapshots. Algorithm update based on these splines; in practice, we sample trajectories and use minibatches. Algorithm| outlines the corresponding inference procedure.
desscrbbes the construc- 2 gives a single-sampl CFO training
Algorithm Quintic Spline Construction
Input: Snapshots (u(tu)..... u(tx)) with time grid T (t=)
2: for = .....,,.. N do
3: Estimate first and second time derivatives dj Rs 0ru((:) and a; s dne((4) using the finite-
difference stencils Appendix A.I.I. for the three-point scheme, use the explicit expression in Eq. (5).
4: end for
5: for O,..., N do
t E [0, 1].
6: Define the local coordinate ti+1 t
7: On [4, 1. construct stt) via the quintic Hermite form in Appendix A.1.2 using the knot data (u(t:), u(te41), d,, d,+1, as a.+1).
8: end for
9: Output: sl) (the collection of per-segment coefficients).
Algorithm CFO Training (One-Sample Example)
1: Input: A trajectory u with time grid T, neural operator (t, u), noise schedule 7().
2: Construct spline st) with (u, T); for the quintic spline, use Algorithm
3: for iter = 1,.... . do
4: Sample
Unif and = 1).
5: x - +t)((t) + (
6: .urget +- 0ts(t) + (t)
7: Update 0 using loss L(0) = [[Nott mangett
8: end for
9: Output: No
16

---

Published as a conference paper at ICLR 2026
Algorithm 3 CFO Inference
1: Input: Trained neural operator Nott, u), initial time to and state un, evaluation times Teval (tx) with to < t1 <... < tk. ODE solver (e.g.. RK4).
2: ((ln) - ug.
3: for k 0...... K do
4: ú(tx+1) <- ú(Ix) flx+1 (t, ((t)) di (approximated by the chosen ODE solver).
5: end for
6: Output: {ū(t)))
A.4 EXPERIMENTAL DETAILS
A.4.1 METRICS
We evaluate model accuracy using the relative L2 norm between the predicted pred and ground-truth Mtrue trajectories, averaged over all test trajectories:
1/2
N
No.
(k) "pred (tj, xii) (k) (t,, (t))
Nisa
i=1
Relative L2
Niest
1/2
Nr
U Merue (ty, 2z))
(K)
A.4.2 NEURAL NETWORK ARCHITECTURES
We evaluate Multi-Layer Perceptrons (MLPs), U-Nets (Ronneberger et al. 2015). Diffusion Trans- formers (DiTs) (Peebles Xie 2023), and Fourier Neural Operators (FNOs) (Li et 2023) CFO models receive sinusoidal time embeddings; autoregressive baselines do not.
MLP A 5-layer MLP with hidden dims [128, 256, 256, 256, 128] and ReLU activations; time- embedding dimension 16 0.202M parameters).
ENO ID Fourier Neural Operator with width 128, 5 spectral blocks, and m retained Fourier modes. Each block applies an rFFT, learned complex-mont mode multiplication with a 1 1x 1 residual conv, and GELU. Inputs are augmented with a normalized grid and lifted to width 128; two-layer head maps output channel 2s 0.591M parameters).
1D U-Net A minimal 1D U-Net with scales (channel 1 mult t [64, 128)) using Conv- GroupNorm- blocks, 4x1 stride-2 downsampling, transposed-convoluue upsampling, and a 1x1 head; time-embedding dimension 64 z= 0.602M parameters).
2D U-Net A 2D U-Net with three scales (channel.mult [64, 128, 256)), Conv- Gron-GrooNNom- Swish blocks, 4x4 stride-2 downsampling, transposed-convolude upsampling, and a 1x1 head; time-embedding dimension 256 8.1M parameters).
DIT A DiT with 8x8 patch embedding, hidden size 384, depth 4, attention heads, and MLP ratio 6. Fixed 2D sinusoidal positional embeddings are added to patch tokens; a linear head reconstructs one output channel. 2s 13.4M parameters).
A.4.3 TEACHER-FORCING TRAINING FOR BASELINE MODELS
We summarize the teacher-forcing schemes used to train the autoregressive and neural ODE base- lines.
Autoregressive baselines. Let (u(4)) denote trajectory sampled uniform time grid. The AR model Fg trained one-step predictor ú(l++1) F& (u(t:)). At each step the input the
17

---

Published as a conference paper at ICLR 2026
ground-truth state u(tz) rather than the model prediction (teacher forcing). The training loss is the average one-step error
N-1
1
(0) = Σ 1=0
Fo (u(t:)) u u(te+1))l2
N
At inference time, we roll out autoregressively by feeding the model's own prediction back input from ground-truth initial condition û(ta) = u(ta):
ú(t++1) Fo(ū(t))), = 0.... ... ,N
Neural ODE baseline. For the continuous-timg baseline, we parameterize the right-hand side as a neural network oe (t, u) and define the ODE
d i(t) fo (t, ū()))
dt
We train this model with a ferchig-forcing strategy on short time intervals. For each pair (14 (++1) in trajectory, initialize the ODE solver the ground-truth state ä(t,) u(t,) and numerically integrate obtain ((+++1) + J, fT( ū(r))dr. The loss the average mismatch between the integrated solution and the next ground-truth snapshot:
1 N-1
(NODE(
Σ
[((t++)) a(t;-1))le
This setup ground-truth states starting points for each integration interval.
A.4.4 LORENZ SYSTEM
We consider the chaotic ODE system with spatial-derivative- dynamics:
dr
dy di x(p y. dl
dz
( x)),
xy Bz,
[0,5].
dl
with (o, p)B (10, 2/38/3).
Dataset. The initial condition (To, 20) is sampled uniformly from the box [-5, 5]S. The dataset comprises 9000 training. 500 validation, and 500 testing trajectories. Each trajectory sampled at 1001 uniformly spaced time points over horizon of 5.0s.
For the irregular sampling experiments, we create subsampled datasets by randomly keeping 50% or 25% of the time points for each trajectory, with the selection performed independently for each trajectory. This procedure consistent across all equations.
Neural Solver Choice. Both the CFO and the autoregressive baseline use the same MLP architec- ture. CFO models, normalized time E (0, 1] is encoded using a 16-band sinusoidal embedding.
Training of CFO and AR. batch size of 2048, a learning rate of 10-4, and betas (51. (2) (0.9, 0.99), and evaluated every 5000 steps. The noise schedule parameter is set to 7o = 10
The model is trained for 200,000 steps using the Adam optimizer with
Training of Neural ODE. We use the Tsit5 adaptive Runge- Kutta solver (Titouras) 201 with absolute tolerance 10-6 and relative tolerance 10-1 Backpropagation through the ODE solu- tion is performed using Diffrax's Recurs veCheckpo wtw.hheeccoiiindjjooooo joint nt method (Kidger 2021 for memory-efficient gradients. The other training setting matches CFO and AR.
A.4.5 BURGERS' EQUATION
Consider the ID Burgers' equation benchmark investigated in (Wang et al. 202
du dl = 1/ dx?
d
du
tr, E (0, x (0. 1), x (0, 1),
dx
u(x, 0) = u0(x),
with the viscosity 0.01.
18

---

Published as a conference paper at ICLR 2026
Table Endpoint error and empirical order p for the spline-implied velocity on the Lorenz system. Errors are average relative L2 at trajectory segment endpoints (r E {0, 1)) over 50000 segments. At 0.005
Spline type
Step size
Error (T=0)
p (T=00
Error (T=1)
p (T==1)
Quintic
At
5.33 x 10
5.32 x 10
2
2.16 10-3
2.02
2.16 10-3
2.02
4
8.63 10-3
2.00
8.63 10
2.00
Linear
At
2.49 10-2
2.48 10
2
4.97 10-2
1.00
4.96 10-2
1.00
4
9.92 x 10-2
1.00
9.88 10
1.00
Dataset. The initial condition (x) is generated from Gaussian Random Field (GRF) ~ N (0,252(-4 + 52 with periodic boundary conditions. The dataset consists 1000 training, 100 validation, and 100 testing trajectories. Each trajectory is sampled at 101 uniformly spaced time points over [0, 1s]. The spatial domain [0, 1] is discretized into 100 points.
the 1D U-Net; results appear in Table 1 To demonstrate flexibility, our ablation study part 5.2 include FNO variant. The autoregressive baseline mirrors these architectures (MLP, U-Net, FNO) with time embeddings disabled; we report only the best- performing variant MLP in Table
Neural Solver Choice. For CFO,
Table Relative Lg Error of autoregressive baselines with different architectures on the 1D Burgers' equation. The MLP architecture, as the best-performing variant, is used for comparison in Table
Architecture
ML.P
U-Net
FNO
Relative La Error
3.34 10-2
2.56
1.16
Training CFO and AR. The model is trained for 60,000 steps using the Adam optimizer with batch size of 256, a learning rate of 10-1, and betas (8),8)) (0.9, 0.99). The model is evaluated every 5000 steps. The noise schedule parameter is set 7o = 10-3
Training of Neural ODE.
We use the Tsit5 adaptive Runge Kutta solver ('sitouras 201 with Runge-Kutta
absolute tolerance 10-6 and relative tolerance 10-4, Backpropagation through the ODE solu- tion is performed using Diffrax'ss 's memory-efficient gradients. The other training setting matches CFO and AR.
method (Kidger 202 for
A.4.6 DIFFUSION- REAITTIO EQUATION
We consider the 2D diffusion-reaction equation from PDEBench (Takamoto al. 2022
0,u Dr u ug + Ro (u,
0, = (D =gyv) + Ro (u,
where 2(t, y) is the activator and v(t, v) denotes the inhibitor the domain (t,x, v) E (0,5) x (-,, The reaction terms are given by
(u( 3 k U,
(u(,,
with &: = 5 x 10 The diffusion coefficients are D., = 1 x - and D. = 5 x 10-3
Dataset. We directly use the dataset provided in PDEBench (Takamoto et al. 2022 The initial condition generated as standard normal random noise u(0, x. y) N (0, 1.0) The dataset con- sists of 900 training, 50 validation, and 50 testing trajectories of activator and inhibitor (2 channels). Each trajectory is sampled at 101 uniformly spaced time points over [0, 5s). but the first time step is truncated ensure stability. The spatial domain is discretized into 128 128 grid.
19

---

Published as a conference paper at ICLR 2026
Neural Solver Choice. For CFO, the 2D U-Net described above; results shown in Table I To demonstrate flexibility, our ablation study part 5.2 include a DiT variant. The autoregressive baseline mirrors these architectures (U-Net, DiT) with time embeddings disabled; we report only the hest-performing variant-UUNNet-iin Table
Table 6: Relative L2 Error of autoregressive baselines with different architectures on the 2D diffusion-reaction Equation. The U-Net architecture, as the best-performing variant, is used for comparison in Table
Architecture
U-Net
DiT
Relative L2 Error
# 10 (9.21 + 0.10) x 10
Training of CFO and AR. The model is trained for 5000 steps using the Adam optimizer with batch size of 256, learning rate of 2 x 10-4, and betas (51, (22) (0.9, 0.99). The model is evaluated every 500 steps. The noise schedule parameter is set to "Yo = 10 10-b.
Training of PDE-Refiner. For the PDE-Refiner baseline, we follow the training setup from the original paper (Lippe et al. 2023) for 2D equations and use the 2D U-Net architecture with pa- rameter count and number of training steps matched to CFO for a fair comparison.
A.4.7 SHALLOW WATER EQUATION
We consider the 2D shallow water equations from PDEBench (Takamoto et al. 2022
0gh + (hw) + +y(hh)) = 0,
1
0,(hu) + az
u2h +
2 9rh
+ (wy(u = h9,
1
(he)) + 0g e +
grh2
+ (uvh) = -9hho,b
2
where h is the water depth, (u, v) are the velocities, is the spatially varying bathymetry, and gr is the gravitational acceleration.
Dataset. We use the dataset from PDEBench TTaamotoo et al. 2022). The initial condition for water depth is generated from GRF, with initial velocities set to zero. The dataset consists of 900 training, 50 validation, and 50 testing trajectories of water depth (one channel). Each trajectory is sampled at 101 uniformly spaced time points [0,1], with the first time step truncated. The spatial domain is 128 X 128 grid.
Neural Solver Choice. For CFO, we use the 2D U-Net; results are shown in Table ! To demon- strate flexibility, ablation study in Section 5.2 includes DiT variant. The autoregressive base- line mirrors these architectures (U-Net, DiT) with time embeddings disabled; we report only the best-performing variant (U-Net) in Table for comparison.
Table Relative L2 Error of autoregressive baselines with different architectures on the Shallow Water Equation. The DiT architecture, the best-performing variant, is used for comparison in Table
Architecture
U-Net
DiT
Relative L2 Error
3.16 x 10-1 9.04 x 10-2
Training of CFO and AR. The model is trained for 5000 steps using the Adam optimizer with batch size 256, a learning rate of 10-5, 5 and betas (81, (2) = (0.9, 0.99). The model is evaluated every 500 steps. The noise schedule parameter is set " 10-6
20

---

Published as a conference paper at ICLR 2026
t- 1,25
1.5.
1.7555
T- 5.0
42
-12
-02
-50
3.000
0.755
000
(a) Activator field u
t- 1.25
t 2.5
3.75
:- 5.0
20
41
-20
200
2.04
3.00
204
110
(b) Inhibitor field
Figure 5: DR trajectory visualization for activator (top) and inhibitor (bottom).
21

---

Published as a conference paper at ICLR 2026
Training of PDE-Refiner. For the PDE-Refine baseline, follow the training setup from the original paper (Lippe et al. 2023) for 2D equations and use the 2D U-Net architecture with pa- rameter count and number of training steps matched to CFO for a fair comparison.
0.25
0.5
t 0.75
t 1.0
1.2
1.1
Ref. h
1.0
0.8
0.0
1.2
CFO Pred.
1.1
1.0
0.9
0.8
1.2
AR Pred.
1.1
0.20
CFO Erro
0.15
0.10
0.05
0.00
0.20
AR Error
0.15
0.10
0.05
0.00
Figure 6: Shallow Water Equation trajectory visualization. Comparison of long-horizon rollout for the test sample with the highest error. Quintic CFO, trained on 25% irregularly sampled data, captures the dynamies, while the autoregressive baseline (trained on full-resolution data) accumu- lates significant error.
A.4.8 ABLATION EXPERIMENTS
To complement the ablation study in the main text, we provide additional results and figures in this subsection.
Impact of noise schedule. We sweep the noise magnitude 7o. See Figure 7 (SWE the left; Diffusion-Reactioon on the right).
Impact of ODE solvers and function evaluations. We study the tradeoff between accuracy and computation by varying the solver and the number of function evaluations (NFE). Errors generally decrease with NEE, and higher-order solvers achieve lower error at fixed budget. See the Euler and Heun sweeps in Tables and for final-time relative L2 errors across tasks. Table IO compares wall-clock inference time and accuracy between the AR bascline and CFO with RK4 solver using 50% of the AR NFE. Figure visualizes the euror-n across different NFE settings.
22

---

Published as a conference paper at ICLR 2026
SWE
Taff. Tiff-Routt.
2.0
6.0
6.0
Rel. 6. Rm.
-36.
6.00
45
0.5
10""
Figure 7: Impact of noise magnitude to model performance for SWE (left) and DR (right) bench- marks.
Table 8: Influence NFE on Quintic CFO (Fuler solver) accuracy (Relative L2 Error) at the final time step, compared the autoregressive (AR) baseline.
NFE (% AR)
Lorenz
Burgers
DR
SWE
AR (100%)
0.1481
0,0647
0.3850
0.1048
50%
0.3032
0,0142
0.0611
0.1345
100%
0.2821
0.0111
0.0603
0.0978
200%
0.2337
0.0099
0.0588
0.0651
400%
0,2076
0.0094
0.0603
0.0389
Table 9: Influence NFE on Quintic CFO (Heun solver) accuracy (Relative L2 Error) the final time step, compared to the autoregressive (AR) baseline.
NFE (% AR)
Lorenz
Burgers
DR
SWE
AR (100%)
0.1481
0.0647
0.3850
0.1048
50%
0.1878
0.0090
0.0694
0.0874
100%
0.1052
0,0090
0.0601
0.0249
200%
0.0705
0.0090
0.0600
0.0076
400%
0.0651
0.0090
0.0594
0.0061
Comparison with other baselines. ODE for Lorenz system: CFO converges faster and more stably, especially in the early stage of training.
Figure9sshon shows the rescaled training loss for CFO and neural
Reverse-time inference. We evaluate this the reverse-time inference capabilities our framework on the dissipative Burgers' equation. Figure 10 visualizes the results.
Temporal extrapolation beyond the training horizon. We evaluate long-horizon temporal ex- trapolation by training the first half of each trajectory, [0, T T/2), and rolling out the full horizon, []. Table 11 reports the resulting errors, and Figure 11 shows error accumulation on the four benchmarks, compared baseline trained the full time horizon.
Neural Network Architectures. CFO is largely architecture-agnostic. swapping the spatial back- bone (FNO for 1D Burgers; DiT for 2D Diffusion-Reaction and SWH) yields comparable accuracy. See Table 12
More Challenging Examples. We further evaluate CFO the one-dimensional Kuramoto- Sivashinsky (KS) equation, a fourth-order chaotic nonlinear PDE,
"r + W0g + Mar + D'lazzaa 0.
We follow the setting and dataset of Majid 2024) (140-step rollout) and use the same neural net- work architecture (1D U-Net) for all methods. As shown in Table T3 CFO achieves roughly 2x
23

---

Published as a conference paper at ICLR 2026
Time (ms/sample)
Rel. L2 Error
Task
Backbone
AR
CFO (50% NFE)
AR
CFO (50% NFE)
Lorenz
MLP
0.16
0.12
1.48 x 10-1
9.19 x 10
Burgers
ID U-Net
0.68
0.31
1.03 x 10a
9.00 x 10-3 3
DR
2D U-Net
117.18
49.36
3.74 10-1
8.76
10
SWE
2D U-Net
111.83
49.63
3.84 x 10-1
6.98 x 10
Table 10: Inference runtime (ms/sample) and relative L2 error at the final time step for the autore- gressive baseline and CFO using 50% of the AR NFE, evaluated single NVIDIA A6000 GPU (JAX, fp32) with the same spatial backbones.
[ Leremu: Fina-timee Error Wall-clock sumple (ms)
Burgess:
Final-time Error Wall-elock per sample (ms)
0.200
10
(RK44 varying steps)
varying steps)
Flou--ime relative La croo
8.173
baselize
Fluni- relative L 4ro0
bareline
0.150
100
0.125
0.100
10-1
0.075
0.2
0.55
3.6
01.8
1.0
Wall-clodk time (rs8/ scmple
Wal-clock sime [ms sample]
(a) Lorenz system
(b) ID Burgers' equation
3.4
DR: Final-time Eroo va Wallellock oee samole (ms)
Final-time Error vs. Wall-cookk per sample (ms)
CFO (RK4, varying stope)
CPO (RK4, varying steps)
baseline
Fina-time relative otsor
AR baceline
Final-Fime velatine La
3.3
9.3
0.2
0.1
100
300
100
200
200
4000
Wall-clock time sample]
Wall-clock time Ims sample]
(c) 2D diffusion reaction system
diffusion-rraattion
(d) 2D shallow water equations
Figure 8: Error versus wall-clock time for CFO with different numbers of function evaluations (NFE). The star marker denotes the AR baseline.
lower error than the neural ODE baseline and 3x lower error than AR, while being vastly more computationally efficient than the neural ODE (comparable training time to AR).
Latent Space Simulation. In addition to direct application, CFO can be extended to operate in latent spaces. For the PDE in (2), consider invertible operator & that maps the function latent representation Ew).
Define the conjugation of A as ENN-! then the dynamics in the latent space can be expressed as:
E(u,) (((.))(((()))..
24

---

Published as a conference paper at ICLR 2026
CFO
15.0
Neural ODE
12.5
Train Loss
10.0
7.5
6.0
2.5
0.0
10000
20003
30000
40000
50000
60000
Step
Figure 9: Training loss versus optimization steps for CFO and the neural ODE baseline the Lorenz system. CFO converges faster and with substantially reduced oscillations, whereas the neural ODE exhibits large fluctuations in the early stages before eventually reaching a comparable loss level.
Burgess: Reseese Time Inference
Predictions the
0.3
0.0001
0.001
0.G
= 0.01
Relative Lg
0.1
0.4
Vidce
0.0
Groand Truth (t=()
-02
Prodicted 1-0.5
Paslicted bum (-1150
0.0
Prodicted t-0.75
100
Paslictiad from tw100
Backward horison (% of obeerved weindow)
0.0
0.2
0.8
0.8
1.0
(a) Relative Le error of reconstructed @(t) for t E 1) under different noise levels added to the terminal state u(1)-
(b) Trajectories reconstructed backward from / 0.75, 0.5, 0.25] to f O.
Figure 10: Reverse-time inference for dissipative Burgers' equation.
Specifically, if & is time-independent linear operator (e.g., Fourier wavelet transform), then commuting with the time derivative yields:
[[()], -(E(w)).
In this case, CFO learns the spatial operator N directly in the latent space. Experiments on the 2D diffusion-reaction and 2D shallow water equations (see Table 14) show that CFO effectively models latent dynamics, achieving performance comparable to direct CFO in the original space.
25

---

Published as a conference paper at ICLR 2026
Equation
Train [0, [']
Train (0,) Test [0, 11
Lorenz
4.53
10-2
3.67 10-2
Burgers
5.89 x 10
7.94 x 10-3
DR
4.37 x 10-2
4.78 x 10-2
SWE
4.56
10-3
2.28 10-2
Table 11: Relative Lg error of temporal extrapolation beyond the training horizon.
Locenz: CEO Triaet on (0. T. vel T/2)
Burgene: CFO Trained ca T vs T/2)
lsaiting sinker 1/2)
0,0170
Eid
CSO Jrain on P. T23))
(PF (toxin on D. T22))
(Loiis D.
0.0150
(taxto on 0 "I)
Iolutte 17 croor
0.00
000022
0.007
0.02
0.000
0.00
0.0000
000
900
1000
Time
Time
(a) Lorenz system
(b) ID Burgers' equation
(0. T/2)
Thiine
T/221
Exst teaining windot T/2)
End. of training window R. 1/2)
CPO Ttuins oo B. 1/2))
C20 [touia on (D. 1/2))
0.00
(triis T
Beladtie 12 61105
01.04
1,
0.03
0.03
0.02
0.00
0.01
0000
0.00
20
20
130
20
00
80
100
Timee
Tiree alog
(c) 2D Diffusion-Reaction
(d) 2D Shallow Water Equation
Figure 11: Long-horizon temporal extrapolation: error accumulation on Lorenz, Burgers, DR, and SWE. Models are trained on E T/2) and evaluated on rollouts over E [0, T7.
Table 12: CFO performance across (Relative L2 Error) neural backbones FNO (ID Burgers) and
DiT (2D Diffusion-RRactionn 2D SWE).
Equation
Burgers
SWE
DR
Rel. Lg Er.
0.0068
0.0083
0.0821
Method
Rel. Lg Err.
Training Time h 40 h 53 min 1 5 min
AR
0.317
Neural ODE
0.197
CFO
0.105
Table 13: Performance on the ID Prramoto-Sivashinsky (KS) equation. CFO attains substantially lower error than AR and Neural ODE while remaining far more computationally efficient than the Neural ODE baseline.
26

---

Published as a conference paper at ICLR 2026
Table 14: CFO in wavelet latent space: Relative L2 Error for 2D Shallow Water and 2D diffusion- reaction equations.
Equation
SWE
DR
Rel. L2 Err.
0.0053
0.0218
27


Reviews:
[{"review_id": "uLgSHXhv3v", "reviewer": "Reviewer_CPTa", "summary": "The paper proposes a method to train spatial neural PDE models without relying on explicit integration-based optimization, as is common in neural ODE frameworks. Instead, temporal derivatives are supervised through stochastically interpolated trajectories derived from discrete training snapshots. The study further evaluates the method’s in-domain predictive performance within the same temporal horizon, under varying temporal sampling densities.", "questions": "See weakness.", "limitations": "", "rating": 2, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "The main strengths of the paper are as follows:\n- It proposes an approach that significantly reduces the computational and memory costs associated with applying directly neural ODE frameworks to PDE prediction tasks.\n- It bridges flow learning and stochastic modeling by introducing a strategy to estimate temporal derivatives at unobserved time points through stochastic interpolation of training snapshots.", "weaknesses": "- Limited baselines: The paper only compares against a discrete-time autoregressive model. While neural ODE-based approaches can indeed be computationally demanding, this alone is not a sufficient reason to omit them from quantitative comparison. Including related continuous-time baselines such as the mentioned continuous-time methods (e.g., Chen et al., 2018 for ODE/PDEs; Yin et al., 2023 for PDEs) or other comparable approaches would strengthen the experimental validation.  \n- Unclear training details for autoregressive models: The paper does not specify how the autoregressive baselines are trained. Please provide more details about their training procedures (e.g., whether the model is rolled out for several time steps and supervised over the full trajectory).  \n- Different training strategies in practice similar to the proposed method: Common strategies such as teacher forcing (supervising one-step-ahead predictions) or temporal curriculum scheduling (progressively extending prediction horizons) can simplify training for both autoregressive and neural ODE models while reducing computational cost. These strategies can be viewed as alternative ways to achieve interpolation. A discussion or comparison of how these techniques relate to the proposed approach would provide valuable insight.  \n- Limited evaluation scope: The experiments focus solely on in-domain prediction. To better assess the model’s ability to learn underlying dynamics, the authors should evaluate extrapolation beyond the training horizon (e.g., training on data within [0, T/2) and testing on both [0, T/2) and [T/2, T)).  \n\nReferences  \n- Chen et al., 2018. Neural Ordinary Differential Equations. *NeurIPS 2018.*  \n- Yin et al., 2023. Continuous PDE Dynamics Forecasting with Implicit Neural Representations. *ICLR 2023.*", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "TcvqL9UIxL", "reviewer": "Reviewer_CE6Q", "summary": "The paper introduces the Continuous Flow Operator (CFO), a framework that learns continuous-time PDE dynamics by flow-matching the analytic time-derivative of a spline interpolant fitted to each training trajectory. Concretely, the authors estimate temporal derivatives on possibly irregular grids, build a piecewise quintic Hermite spline $s(t; u)$ that matches values and derivatives at knots, and then train a neural operator $N_\\theta(t, u)$ to regress $\\partial_t s(t)$ without back propagating through any ODE solver during training. At inference, the learned vector field $du/dt=N_\\theta(t, u)$ is integrated forward or backward for arbitrary time resolutions. This yields time-resolution invariance in training and querying, competitive accuracy-efficiency trade-offs via solver/NFE sweeps, and notable data efficiency: training with only 25 % irregular time points still beats strong autoregressive baselines trained at full resolution across Lorenz, 1D Burgers, 2D diffusion-reaction, and 2D shallow water.", "questions": "1. How far does reverse-time hold before numerical blow-up on dissipative systems like Burgers/SWE, and what step control or regularization is needed in practice? A short curve of backward-horizon vs. $L^2$ error would be enough.\n\n2. For fairness, would the authors include an AR baseline with the same spatial backbone and time embedding that CFO uses, trained on the full grid, plus Latent ODE and Neural CDE implementations adapted to the PDE setting (or at least reported on Lorenz)? This would help isolate the gain due to the training objective rather than the architecture.\n\n3. It would be very valuable to add one harder PDE. Two natural choices are 2D Navier–Stokes (Kolmogorov flow) and Kuramoto–Sivashinsky, both available in community suites such as PDEBench. A single appendix experiment on either would strengthen generality.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "From an originality perspective, the idea of repurposing flow matching, which is originally developed for CNFs/diffusion, to directly learn RHS dynamics is elegant and avoids the heavy simulation-through-solvers that slows Neural ODE/SDE training. The positioning is credible relative to Flow Matching in generative modeling and stochastic interpolants, as CFO leverages fixed probability paths defined by the spline while keeping the operator-learning focus. I appreciate that the authors discuss flow matching lineage and deliver a practical instance specialized to temporal PDEs.\n\nThe time-resolution invariance protocol with per-trajectory random time grids is convincing, and the table shows that quintic-CFO trained on 25% of points outperforms AR trained on 100%, with large relative error reductions. The Lorenz downsampling study that compares linear vs. quintic spline velocities and observes first- vs. second-order endpoint behavior is technically sound and matches finite-difference theory. The NFE sweeps across Euler/Heun provide a clear view of accuracy vs. compute. \n\nCFO serves as a bridge between discrete AR operators and fully continuous Neural ODE training. It inherits irregular-time training, continuous-time querying, and even reverse-time rollouts without the pain of adjoint backprop; in my opinion this is a practical step forward for scientific ML systems that see irregular telemetry or experiment logs. The method also appears architecture-agnostic (FNO / DiT backbones), which increases adoptability in operator-learning pipelines.", "weaknesses": "The reverse-time inference claim is attractive, but dissipative PDEs often make backward integration ill-posed. I would prefer a short quantitative check showing how error grows with backward horizon and perturbation size to avoid over-selling.\n\nBaseline coverage against continuous-time sequence models could be expanded. Since CFO’s training avoids solver backprop similarly to newer flow-matching approaches for time series, a comparison to Trajectory Flow Matching (TFM) would situate the contribution more precisely; likewise Latent Neural CDEs are classic irregular-time references. \n\nFinally, efficiency comparisons mix fixed-step AR with variable NFE CFO. The NFE sweeps are good, but I suggest harmonizing by wall-clock on the same hardware and showing error vs. seconds to eliminate any suspicion of budget mismatch.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "P4qDoLFT9V", "reviewer": "Reviewer_weht", "summary": "This paper proposes the Continuous Flow Operator (CFO), a method that learns continuous-time PDE dynamics from data. Its key innovation is using flow matching and spline interpolation to directly learn the underlying dynamics, avoiding the computational cost of methods like Neural ODEs. CFO is uniquely time-resolution invariant, meaning it can be trained on data with irregular, non-uniform time steps and make predictions at any time resolution.", "questions": "1. When using finite difference approximations for gradients, estimation errors are present. Particularly in cases involving non-uniform time steps and large time intervals, direct flow matching may suffer from insufficient accuracy.\n\n2. During inference, the authors employed the RK4 method, whose computational efficiency remains relatively low.\n\n3. It should be noted that the neural operators referenced by the CFO are distinct from established architectures such as DeepONet or FNO, which may cause potential misunderstandings.\n\n4. See weaknesses.", "limitations": "", "rating": 2, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 1, "strengths": "1. This paper is well-organized and highly readable.\n\n2. The adoption of the flow matching concept can significantly reduce the training complexity of neural networks, and may even enhance the learning accuracy in certain scenarios.\n\n3. The author has conducted experimental analyses across multiple datasets.", "weaknesses": "1. The work lacks originality, as the idea of applying flow matching to train Neural ODEs has been previously explored in several earlier studies.\n\n2. The experimental section lacks validation on real-world datasets, which significantly limits the persuasiveness of the claims.\n\n3. It lacks comparison with existing SOTA baseline methods.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "D3YXO9SQMP", "reviewer": "Reviewer_3Ms9", "summary": "This paper introduces the ​​Continuous Flow Operator (CFO)​​, a novel framework for learning continuous-time PDE dynamics by repurposing flow matching to directly approximate the right-hand side of PDEs. This approach avoids the computational burden of backpropagating through ODE solvers (as in Neural ODEs) while enabling training on irregular time grids and inference at arbitrary temporal resolutions. Evaluated on four benchmarks (Lorenz, 1D Burgers, 2D diffusion-reaction, 2D shallow water), CFO demonstrates exceptional data efficiency and long-horizon stability.", "questions": "See weaknesses.", "limitations": "", "rating": 4, "confidence": 2, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- Novel Continuous-Time Framework​​: CFO bridges flow matching (typically used for generative modeling) with PDE dynamics learning.\n\n- Architecture Agnosticism​​: CFO is compatible with various neural operators (e.g., U-Net, FNO, DiT). Experiments show consistent performance across backbones, allowing users to choose architectures based on spatial inductive biases or hardware constraints.", "weaknesses": "- Spline Order Sensitivity​​: While quintic splines generally outperform linear ones, the optimal spline order may depend on the PDE complexity and data density. The paper does not provide guidelines for selecting spline orders adaptively.\n- High-Dimensional Scaling​​: Experiments are limited to 2D spatial domains. The computational cost of spline fitting and ODE integration in higher dimensions is not discussed.\n- Comparison to Modern Baselines​​: While CFO outperforms autoregressive models and Neural ODEs, comparisons to recent continuous dynamics modeling methods (such as CORAL[1], MARBLE[2]) are limited. These methods integrate Neural ODEs and neural fields, achieves better long-horizon prediction performance.  A broader baseline comparison would better contextualize its advantages.\n\nReference\n\n[1] Serrano, Louis, et al. \"Operator learning with neural fields: Tackling pdes on general geometries.\". NeurIPS 2023. \n\n[2] Wang, Honghui, Shiji Song, and Gao Huang. \"GridMix: Exploring Spatial Modulation for Neural Fields in PDE Modeling.\" ICLR 2025.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)