Improving Long-Range Interactions in Graph Neural Simulators via Hamiltonian Dynamics
Authors: Tai Hoang
Alessandro Trenta
Alessio Gravina
Niklas Freymuth
Philipp Becker
Davide Bacciu
Gerhard Neumann
Venue: ICLR 2026 Poster (2026)
Keywords: Graph Neural Simulators
Long-range interactions
Learning Simulators
AI4Science
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: We propose a novel Graph Neural Simulator that preserves information during propagation, enabling it to model complex physical dynamical systems with long-range dependencies.

Abstract: Learning to simulate complex physical systems from data has emerged as a promising way to overcome the limitations of traditional numerical solvers, which often require prohibitive computational costs for high-fidelity solutions. Recent Graph Neural Simulators (GNSs) accelerate simulations by learning dynamics on graph-structured data, yet often struggle to capture long-range interactions and suffer from error accumulation under autoregressive rollouts. To address these challenges, we propose Information-preserving Graph Neural Simulators (IGNS), a graph-based neural simulator built on the principles of Hamiltonian dynamics. This structure guarantees preservation of information across the graph, while extending to port-Hamiltonian systems allows the model to capture a broader class of dynamics, including non-conservative effects. IGNS further incorporates a warmup phase to initialize global context, geometric encoding to handle irregular meshes, and a multi-step training objective that facilitates PDE matching, where the trajectory produced by integrating the port-Hamiltonian core aligns with the ground-truth trajectory, thereby reducing rollout error. To evaluate these properties systematically, we introduce new benchmarks that target long-range dependencies and challenging external forcing scenarios. Across all tasks, IGNS consistently outperforms state-of-the-art GNSs, achieving higher accuracy and stability under challenging and complex dynamical systems. Our project page: https://thobotics.github.io/neural_pde_matching.

Published as a conference paper at ICLR 2026
IMPROVING LONG-RANGE INTERACTIONS IN GRAPH NEURAL SIMULATORS VIA HAMILTONIAN DYNAMICS
Tai Hoang Alessandro Trenta*: Alessio Gravina":
Niklas Freymuth' Philipp Becker' Davide Bacciu² Gerhard Neumann'
1 Karlsruhe Institute of Technology
University of Pisa
ABSTRACT
Learning simulate complex physical systems from data has emerged promising way to overcome the limitations of traditional numerical solvers, which often require prohibitive computational costs for high-fidelity solutions. Recent Graph Neural Simulators (GNSs) accelerate simulations by learning dynamics on graph-structured data, yet often struggle to capture long-range interactions and suffer from error accumulation under autoregressive rollouts. To address these challenges, we propose Information preseevvig Graph Neural Simulators (IGNS), a graph-based neural simulator built on the principles of Hamiltonian dynamics. This structure guarantees preservation of information across the graph, while ex- tending to Hort-tamiltnnian systems allows the model to capture a broader class of dynamics, including non-conservative effeets. IGNS further incorporates a warmup phase to initialize global context, geometric encoding handle irregular meshes, and a multi-step training objective that facilitates PDE matching, where the trajectory produced by integrating the port-Hamiltonian core aligns with the ground-truth trajectory, thereby reducing rollout error. To evaluate these prop- erties systematically, we introduce new benchmarks that target long-range de- pendencies and challenging external forcing scenarios. Across all tasks, IGNS consistently outperforms state-of-the-arr GNSs, achieving higher accuracy and stability under challenging and complex dynamical systems. Our project page: https: cs. thub.. 1o/ lo/neural _pde_matccinng
INTRODUCTION
Simulating complex physical systems is central to many scientific domains. These systems are typically governed by partial differential equations (PDEs), which are rarely solvable in closed form. To approximate the solution, classical numerical solvers (Dormand) 1996. Ascher 2008 Evans)
2010) discretize the domain, obtaining a mesh in the case the finite element method (Brenner & Scott, 2008). However, most solvers are tied to specific kind of discretization, and sufficiently accurate simulations require tnn-ggrined meshes, which quickly become prohibitively expensive.
In contrast, neural simulators directly learn system dynamics from data (Guo et al. 2016 Bhatna- gar al., [019 Li al., [2099;; Pfaff al.) (202)). enabling simulations that can be several orders of magnitude l'ster than numerical methods. In particular, Graph Neural Simulators (GNSs) (Paaff
let al. 2021 Linkerhägner et al. 2023 Yu al. 2024) leverage Graph Neural Networks (GNNs) (Wu et al. 2020; Gravina & Bacciu 2024) to capture interactions between entities, such as mesh ele- ments, through message passing (Gilmer et al., 2017). By encoding local geometry in graph edges, they achieve accurate and efficient simulations that generalize across meshes in applications ranging from fluid to elastic dynamics (Sanchez-Gonzale Gonzalez et al. 2020 Pfaff et al. 2021 Yu et al. 2025).
GNSs are commonly trained with one-step loss to predict the next state from the current observa- tion, and then rolled out autoregressively to generate long trajectories during inference. While being effective short horizons, even small global errors accumulate rapidly. Implicit (Pfaff et al. 2021 Brandstetter et al. 2022) and explicit (Würth et al.) 2025) denoising objectives alleviate this some extent by reducing local noise, but the loss still fails to capture fow-frequency drift that
1 "Correppondence hoang@ki edu..
'Equal contribution.

---

Published as a conference paper at ICLR 2026
Inpat
Pateps warnuup
Globally Informed initial condittion
Information preserving lseratiive rollout simulation
Port-Ham b= (q.(6)
Port-Ham. (o)
Port-Ham. b= (g)
Eq.(6)
X =
X = X
Conti
1=1
1=2
Figure 1: A high-level overview of the proposed IGNS. The model takes as input the initial node state X and performs an 1-step warmup phase (left), enabling each node to incorporate broader spatial context before the rollout begins. The enriched state is then used to initialize the rollout X x (middle). During the simulation phase (right), the system evolves according to the port-Hamiltonian dynamics of Eq. (6). while the multi-step loss Lmilti-step supervises all intermediate predictions, ensuring stable and accurate simulations.
emerges only after several steps. Multi-step objectives have been proposed to overcome the short- comings of this one-step training (Lienen & Günnemann 2022). In practice, however, they only work over short horizons, since message passing quickly becomes unstable and loses information after a few steps (Arroyo et al. (025), making direct long-horize training infeasible. This limitation fundamental existing GNSs, as their GNN architectures struggle with long-range dependencies (Cai & Wang 2020) Rusch et al. (023), over-squashing (Alon & Yahav) 2021 Di Giovanni al. 2023). and more generally gradient degradation (Arroyo al. (025).
due to over-smoothing
We propose Information-presse preserving Graph Neural Simulators (IGNS), a neural simulator that learns continuous dynamics in a latent space and enables effective information propagation across the graph for accurate predictions. A schematic of IGNS is shown in Fig. | unlike existing approaches with dissipative updates that rapidly degrade information as pressage-passinggg depth increases, IGNS structures information flow through port-Hamiltonian dynamics (Van der Schaft 20177; Galimberti ett al., 2021 Lanthaler al. 2023 Heilig et al. 2025), better preserving long-range dependencies between distant nodes. By training with a multi-step objective, our setup explicitly learns the under- lying PDE by matching the solution of the parameterized Nort-Hamiltonian dynamics to ground-trute trajectories
time window T. enforcing states to remain on-trajectory and thereby reducing ror accumulation (Radler et al. 2025). Each integration step in IGNS thus mirrors a forward update of classical PDE solver, where signals naturally propagate in space time evolves.
In addition, IGNS incorporates a warmup phase to initiate global context before rollout and a geo- metric encoding module to capture local spatial relations irregular meshes. Together, these com- ponents improve stability and generalization, enabling accurate simulation across diverse physical systems while maintaining a physically consistent inductive bias. On benchmarks targeting long- range and complex dynamics, IGNS achieves lower eтror and outperforms strong GNS basclines by staying on- un-trajectory over long horizons with high accuracy up to 7'=200. By contrast, prior work either (i) uses rather short supervision windows (T=)0) at train and test (Lienen & Günnemann Shi et al. 2023b) or trains on short-windows and then rolls out autoregressively (Janny al,,
2022
2023 (ii) applies multi-step losses only in reduced latent space rather than the original graph (Han et al. 2022).
To summarize, we (i) introduce Information-pressevee Graph Neural Simulators (IGNS), graph neural simulator built on the port-Hamiltoiian formalism that supports information-preserving roll- outs with high accuracy; (ii) analyze its theoretical properties, showing that IGNS maintains infor- mation flow across the spatial domain while capturing broad class dynamics; (iii) aim to mimic classical PDE solvers by matching, through a multi-step objective, the trajectory obtained by inte- grating the port-Hane Hamiltonian core of IGNS with the ground-truth trajectory; and (iv) propose new tasks including Plate Deformation, Sphere Cloth, and Wave Balls, specifically designed to evaluate long-range propagation and oscillatory dynamics under external forcing.
2 BACKGROUND
In this section, we provide an overview of PDEs for modeling physical systems, followed by a brief discussion of message-passing. -ased GNSs for learning such systems. An in-depth discussion of
2

---

Published as a conference paper at ICLR 2026
related works GNSs, effective propagation, and Hamiltonian-inspireed neural networks is provided
in Appendix
2.1 PARTIAL DIFFERENTIAL EQUATIONS FOR PHYSICAL MODELING
Many physical phenomena can be described by PDEs, which relate the rates of change of quantity with respect to both space and time. A general form of a time-dependent PDE is given by
deu = F(u, Vu, V',,..;
(æ, E 92 x (0,
(1)
where : f2 x []. + R is the unknown function (without loss of generality, assumed to be scalar field), 52 S I2d is the spatial domain, and (possibly nonlinear) function u and its spatial derivatives.
To solve Eq. ( numerically, one typically discretizes the spatial domain using methods such as finite differences finite elements, resulting in system of ordinary differential equations (ODEs) that can be integrated over time using various time-stepping schemes (Igel 2016). This often referred to the method of lines, where spatial derivatives are approximated at discrete points, leading to a system of ODEs the form
üi = Si(,, u()), i= u ú,
where (t) (t). Un (+)] represents the state of the system at discrete spatial points, and fo encapsulates the spatial discretization of F point i. The initial condition is derived from the initial state of the PDE.
Remark. While first-order systems as in Eq. T are common, many physical phenomena are better captured by second-order dynamical systems (Evans) 2010), which naturally arise in the modeling of oscillatory behavior like propagation of waves and vibrations. A general second-order PDE with damping and external forcing can be written
u,u + B d.uu + L ((,,t),
(z, E Ω (0,
where L and B are linear operators (12). and is prescribed forcing term. A detailed deriva- tion of applying the method of lines to this equation is provided in Appendix C.2 leading to a system of second-order ODEs.
2.2 MESSAGE PASSING FOR MODEL LEARNING
Under the model learning settings, the goal is to leam the underlying dynamics (i.e., the RHS of Eq. (I] from data. Similar to traditional numerical methods, we first need discretize the spatial domain into set of points or elements. This discretization can be naturally represented as graph, G (V, F)) where V is the set of nodes (i.e., the mesh entities) and E S V V is the set of vn edges capturing their interactions. Then, following the method of lines, we obtain the following system ODEs:
1 = X(; 1= ( X.
where (0 E Rdd is the state vector of node at time f, X(t) = ( In UIT E Rud is the node-state matrix, E = [e) E mmmx is the edge-feature matrix, and fü is the unknown function that captures the interactions between nodes based on the graph structure G. Here, also assume the initial condition X is given.
Solving this requires integrating Eq. (3) over time, which can be done using standard ODE solvers (Hairer et al. 1993). However, standard Graph Neural Simulators (GNSs) simplify this process by further discretizing time using fixed step size At and applying the message-passingg paradigm to compute the node states at the next time step x (++1), leading to the following update rule
(++1)
(t) x
At (x) E:G).
(4)
Here, o denotes the function obtained by stacking L message-passinne layers, which aggregate information from neighboring nodes to update the state of each node, defined as
x Onode: Fii Σ e/ts
e/y Oedge (z: By Cit)
(5)
JEN,
3

---

Published as a conference paper at ICLR 2026
Here, Pedge is the message function, Puode is the node update function, and Mt is the set of neighbors of node i. Both Pedge and Pnode are often implemented as MLPs. We can see that under this setting, fa in Eq.
is approximated by a stack of L message-passing layers, i.e., fi(t, X (1); E; 23
98 (x (t), E; G). In practice, one can obtain the whole trajectory by iteratively applying Eq. for time steps, starting from the initial condition X. However, this approach often suffers from error accumulation over long time horizons, limiting long-term accuracy (Brandstetter et al. 2022
In this work, guided by the theory of dynamical systems, we propose a novel GNS that directly models the underlying continuous dynamies in Eq. B built upon the port- Hamiltonian formalism. Under this framework, information is maintained more effectively across the entire space. Combined with the multi-step loss, the learned model demonstrates that it can mitigate the error accumulation issue and achieve greater accuracy on long rollouts. We discuss them in detail in the next section.
3 THE INFORMATION-
GRAPH NEURAL SIMULATORS
Problem Statement. Given a time series of node states {X (4) JT= on an irregular mesh or graph
G (V, E), our goal is to learn and simulate the evolution of each node x; over time. To address this, we propose the Information ceormation-preserving Graph Neural Simulators (IGNS), a GNN-based simu- lator that uses Hort-tamilnnnnan dynamics its latent evolution core. Fig. shows the high-level overview of our method. IGNS is built on four key components:
Port-Hamiltnnian formalism: by parameterizing Hamiltonian H and adding damping and ternal forcing terms, IGNS captures both energy-conserving and non-conservative behavior.
Warmup: the conservation property is further exploited through the use of warmup phase at initialization, enabling the model to obtain more global context from the start of the rollout.
Geometric encoding: by embedding geometric features directly into the edge attributes, IGNS operate irregular meshes and exploit spatial structure for more accurate local interactions. Multi-step objective: it supervises the entire rollout windows, reducing error accumulation, im- proving data efficiency, and better aligning training with long-horizon prediction.
3.1 PORT-HAMILTONIAL FORMALISM
We start by expressing the dynamics on the RHS of the Eq. (3) using the port- nort-Hamiltonian formalism van der Schaft & Jeltsema| 2014 Desai et al. 2021 Heilig et al., 2025), which is generalization of Hamiltonian systems that integrate both conservative and non-conservative dynamics based on energy principles. Specifically, we parameterize the Hamiltonian He (t, X) and evolve the joint state according to
I
à; Jv X)
J
(6)
DoV Ho
ro (1 x)
-1 0
dumping
forcing & residuull
Here, we omit the time index t for clarity, and denote x; pq: with E R and p e Rax $ to represent the generalized coordinates and momenta, respectively. Without forcing and damping, Eq. (6) reduces to the classical Hamiltonian dynamics, which conserves the energy Ha over time. However, most real-world physical systems are non-conservative, where energy can be dissipated (e.g., due to friction) or injected (c.g., due to external forces). To account for this, two additional terms are introduced: Do Vp. p. He(t, and re X), for damping and external forcing, respectively. Next, we parameterize the Hamiltonian Ho as
(o(W(O)) Ho X) = (())): + Σ ) la.
SEV
JEN.
where ō is an anti-derivative of non-linear activation function o, 1, is a row vector of ones of dimension d, and Wtt) and (t) are block diagonal matrices (to ensure the separation into the
and p components), i.e., (t)) diag([W,(s Wp(t)) and V diag([V (t), Vt(t)). The
matrices (t) and (t) can be either shared across time (i.e., W(t)= = w and V(t) = Vt) or
time dependent, as further discussed in Appendix
4

---

Published as a conference paper at ICLR 2026
To justify this choice of dynamies, provide a more theoretical discussion in Section 4 show- ing that Eq. (6) is universal, and via energy conservation, it allows stable long-range information propagation across space. Hence, it is well-suited for modeling complex physical systems.
Symplectic Integrator. To ensure the energy-conseen conggy-cnssrrving property of the Hamiltonian dynamics, we employ symplectic integrator numerically solve Eq. (6). Specifically, we choose the symplectic Euler method (Hairer et al. 2006). which updates the state as follows:
P; (++1) - + At( V Ho tdh) pl")) - DoV,. Ho(t, qtd) p(9) rrott, X ()),
(1+1) 9 + A.V, p. Ho(l, d0) p(tt))
Here, the gradients of the Hamiltonian computed with respect to the generalized coordinates and momenta.
3.2 WARMUP PILASE
Under the message- sastaageepssssigg. passing framework, single update propagates information only to direct neigh- bors. This locality is a fundamental limitation when simulating physical systems, where long-range interactions often play crucial role from the very first timestep. For example, in mesh-based simulations of fluids or elastic materials, the evolution of node may depend not only on its imme- diate surroundings but also on distant regions of the mesh, such boundary conditions of external forces. If the model starts from a completely local representation, it must accumulate long-range information gradually over several rollout steps, which delays the emergence of globally consistent dynamics and may degrade accuracy in the early stages of the simulation.
sources
To alleviate this issue, we introduce warmup phase applied once initialization. Starting from the initial state, we perform additional rounds of message passing without advancing time, as visually summarized on the left of Fig. [ More formally, the initial condition of the simulation is set to ( X (1), where X denotes the state after warmup iterations. This strategy broadcasts information across the graph up to radius from each node, so that for sufficiently large I, each node effectively receives signals from increasingly larger neighborhoods. Thanks to the energy- conserving core of IGNS, this globally informed latent state is preserved throughout the rollout, rather than being dissipated. The effect of is further analyzed in Section 6
3.3 GEOMETRICAL INFORMATION ENCODING
Accurately encoding geometry is essential for modeling physical systems on irregular meshes, where interactions depend on relative positions (i.e., in the mesh space) rather than absolute states (i.e., in the physical world space). To capture this relationship, tem's evolution based on geometrical information, thereby allowing the dynamics to influenced by the spatial structure. We design the external foree following scheme inspired by Praff et al. (2021), where each edge feature Eij encodes both the displacement vectors (Sij pos, )os,) and distances (diy ll pos,ll2). with pos, indicating the position
model the external force to drive the sys-
edge-level pressage-passing
of node the mesh space. More formally, we formulate the external force as re(t, X) = g (.) where each message-passing step is defined as
(node = o W Wnodeeq: + Σ e
= W Wedge (q) - q:) + Eij
(9)
JEN.
with Ei as the concatenation between displacement vectors and distance values, i.c., Eij [Sij, d This formulation explicitly separates the edge in the world space (q) - q.) and the edge in the mesh space (eij) in the edge update. Unlike Pfaff et al. (2021), 1, we do not update edge messages at every time step. Combined with the multi-step loss (discussed below), the edge information serves instead as a static prior that weighs neighbor messages. This makes the model less dependent on explicit geometry, helps prevent overfitting to specific meshes, and still provides sufficiently useful context. We empirically study this hypothesis in Appendix
3.4 MULTI-STEP LOSS
We train our models through multi-step objective. Given window (q ,..., ((+KK) roll out from (()) and optimize both the predicted ((++) and (((+)) to match its corresponding ground truth
5

---

Published as a conference paper at ICLR 2026
fields ((++T) and velocities p(t VT E [1, K],
K
mmulti-stop
- 1 + || -p 11)
(10)
This multi-step loss enhances the consistency between the prediction and ground-trume data on a trajectory level. In contrast, one-step training followed by autoregressive rollout is prone to suffering from the error accumulation problem, especially when the time-series data exhibits a strong temporal correlation, which analyze in Section(] The multi-step loss also fits IGNS well because the dissipative core preserves signals across time; therefore, the gradients from distant steps do not vanish, and the model can make use of the whole trajectory. Standard first-order message passing instead tends to oversmooth, gradients decay over many steps, and the benefit of a long horizon objective is limited. Finally, including p helps improve performance further; when momenta are not explicit in the state, approximate them by finite differences the ground truth q()) and ((++)) A visual exemplification of this multi-step loss strategy is summarized on the right of Fig.
4 THEORETICAL PROPERTIES
In this section, we discuss the theoretical benefits of our IGNS, focusing on universality and effec- tive information propagation during the simulation, underpinning its ability to effectively simulate complex physical systems. A more in-depth discussion is provided in Appendix E
Universality. We start by showing that our IGNS, with state updates in Eq. (6) and Hamiltonian in Eq. 7 yields a universal model. Therefore, IGNS can approximate any given functional + Raxd Raxd that maps the initial condition to the solution at time
Theorem 1. Let Vo be the functional that maps the initial condition To = [qo pol to the solution at time T to the ODE defined Eq. (6) and Eq. (7). In other words, Va represents IGNS. Then, for any F: Raxd + Rnxd with compact support, which represents the evolution of the system after time T, there exist a set of parameters 0 such that (Voo F( S.
The proof is reported in Appendix E.I along with an in-depth discussion and interpretation for phys- ical simulations. As a consequence of the theorem and its proof, IGNS is capable of approximating any functional with compact support that maps the time-dependent signal of external forces and geometrical information into the dynamics of the physical system, making it suitable for accurately modeling complex interactions and dynamics in arbitrary graph-structured systems.
Remark 1. Compactness in the support of the functional a common requirement for Univer- sality theorems. However, this is not restrictive for physical systems, it is sufficient to have limited domain and upper bounds on the physical quantities and energy in the system to satisfy this assumption. All these conditions are perfectly reasonable in contained physical simulations. See Appendix E.1 for a more in-depth discussion on this topic.
Long-Range Information Propagation. Our next focus is on the information propagation capabil- ity of IGNS, a critical aspect for physical simulations, since reliable modeling hinges on the ability to capture long-range interactions in both time and space.
In the absence of damping and external forces, our IGNS reduces to a pure Hamiltonian systems and thus adheres to the corresponding conservation laws (Galimberti et al. 2021 2023 Heilig al.) 2025). Consequently, the vector field induced by the conservative IGNS is divergence- free, ensuring information preservation throughout propagation. In other words, under this setting, the system dynamics can be interpreted as purely rotational, with no information dissipation.
Following Arroyo et al.] (2025). which demonstrates that the inability to perform effective propa- gation is closely tied to the vanishing gradient problem, we investigate the behavior of IGNS with the Hamiltonian proposed in Eq. 7 through a sensitivity analysis. Thus, we study the quantity [|zz(t) (x)))))))))))))))))))))))))))))) [s)]]. as formalized in the following theorem, which follows the sensitivity arguments developed in HHeilig et al. (2025).
Theorem 2. If ax)) is the solution to the ODE h; = JV., Ho that is, the Hamiltonian dynamic defined in Eq. (6 without dampening and external forcing, then dm(t) dz(s)
2 1 for any
0 55a L.

---

Published as a conference paper at ICLR 2026
(a) Plate Deformation
(b) Sphere Cloth
(c) Wave Balls
Three novel long-range simulation tasks. (a) Plate Deformation: a flat plate deformed subject to varying numbers and magnitudes of point forces. (b) Sphere Cloth: a cloth mesh im- pacted by a falling sphere, producing elastic deformation with contact dynamics. (c) Wave Balls: surface waves on water generated by three balls moving linearly from different initial positions.
Figure 2:
The proof is discussed in Appendix E.2 This theorem shows that the sensitivity matrix is lower bounded by a constant, ensuring non-vanisie van-vinishing gradients and thereby enabling effective long-range information propagation; unlike GCNs (Kipf & Welling 2017), which suffer from an exponential decay in propagation rate over time (Gravina et al., 2025). This non-dissipative behavior also en- ables IGNS to propagate and preserve broader context across nodes during the warmup phase (see
Section 3.2)
Although Hamiltonian dynamics provide theoretical guarantees for propagation behavior, additional terms, such as dissipation and external forcing, may counteract this effect. To this end, we recall the discussion in Rusch et al.) (2022); there, they showed that the oscillatory behavior of the system allows the gradients not to decay exponentially even in the presence of these additional terms. More specifically, in the limit of small At, the gradient of the loss function with respect to any learnable weight parameter in the forcing term is independent of the number of updates. We formalize this with statement in Appendix EE.3] inspired from Proposition 3.6 of Rusch et al. (2022). Additionally, as long as the dependence between At and the number of iterations N is polynomial At ~ N gradients do not decay exponentially in N, alleviating the vanishing gradient problem.
5 EXPERIMENTAL SETUP
Datasets. We evaluate our methods on six benchmark datasets in two categories: Lagrangian sys- tems (Plate Deformation, Impact Plate. Sphere Cloth) and Eulerian systems (Wave Balls, Cylinder Flow, Xuramoto-Sivashinsky Sivashinsky). First, we introduce Plate Deformation (Plate Def.) task to evaluate long-range propagation effects on spatial domain, ignoring the evolution over time. Here, given only force positions and magnitudes, the model must learn to predict the final deformed state of an ini- tially flat sheet. Next, the Impact Plate task, proposed in HCMT (Yu et al., 2024), tests the model's ability to capture rapid stress propagation across a large mesh over time. We then introduce two additional new tasks, Sphere Cloth and Wave Balls, both designed to evaluate oscillatory behav- ior under external forcing. Finally, include two complex dynamics tasks, Cylinder Flow (Pfaff ctt al.) 2021) and the fourth-order Kuramoto-Shvvrhinsky (KS) equation, each slightly modified to probe error accumulation under autoregressive training. Fig. 2 highlights the three novel tasks introduced. Full dataset specifications and simulator details are in Appendix (Table
Evaluating Methods. We evaluate broad set of methods, grouped into standard wraph-convolution (Graph Conv.), effective propagation (Eff. Propagation). and our proposed port-Hamiltonian (Port- Нam.) models. Graph Conv.: Our main competitor, MGN (MeshGraphNets) (Pfaff et al.) 2021) is standard message passing GNN trained autoregressively with explicit Buler integration. For the remaining baselines, we use the multi-step loss (Eq. (IO). MP-ODE (Iakovlev al., 2021
Lienen & Günnemann 2022) applies neural ODE integration each step: here, we choose the RK4 integrator (instead of dopri5, 5) ensure a fair training budget compared to other baselines; addi- tionally, we define the message function on nodal points instead of cell elements as in Lienen & Günnemann (2022). Next, GCN+LN is a graph convolutional network with LayerNorm, common technique to mitigate oversmoothing in deep message passing (Luo al. 2024). Eff. Propaga- tion: ADGN (Gravina et al., 2023) adopted an anti-symmetric matrix in the graph convolution to avoid dissipative behavior, CraphCON (Rusch al., 2022) is oscillatory graph network, origi- nally proposed to avoid oversmoothing effects. Port-Hame (Ours): IGNS, (time independent) and IGNS (time varying) are honniio-nninnaee models that follow the state update in Eq. (6). For
7

---

Published as a conference paper at ICLR 2026
Table Per-task test MSE (mean + std, lower is better). Methods are grouped by operator class, (I) Graph Conv., (II) Eff. Propagation and (III) Port-Ham. (Ours). Results are averaged over 4 seeds. Best and second best results are highlighted in orange and teal, respectively.
Plate Def. MSE
Impact Plate MSE
Sphere Cloth MSE (x10-))
Wave Balls MSE (x10
Cylinder Flow MSE (x107
KS
Method
MSE (x107*)
GUN+LN
3.98 + 0.66
3997 98 +: 199 83
26.45 +: 0:28
1.95 +: 0.14
11.71 + 1.03
3:10 +:10+11 0.78
MGN
1.27 + 0.06
3095.75 +:
32.07 2.45
1.78 0.14
12.08 2.60
10.76 9.16
MP-ODE
25.75
1.32
97.04 + 0.70
18.17 + 1.15
59.55
14.82
A-DGN
1.99 + 0.76
3974.56 +: 309.98
30.99 +: 0.80
0.70 +: 0.05
14.72 + 0.94
0.97 0.89
(II)
GraphCON
1.20 4 0.02
2279. 00 +: 720.25
29.00 0.44
1.62 0.03
7.32 4 0.05
0.49 0.02
IGNS&
1.43 - 0.04
75.99 17.17
28.20 + 0.35
0.45 + 0.01
9.12 0.10
0,86 0.12
(III)
IGNS
1.34 = 0.05
67.86 = 14.75
24.16 + 0.67
0.37 + 0.01
8.21 = 0.38
1.20 =: 0.18
fairmess, all baselines, except MGN which updates edge features at each message- paatgggg-pninggggg step, use the same proposed geometric encoding (Section B.3 Appendix Groovdees details on the choice of hyperparameters and training budget.
Our experimental setup highlights the importance of being on-trajectory during rollout via the use of multi-step loss, i.e., to avoid the over-reliance on autoregressive state carry-over (Brandstetter et al.) 2022
Radler et al. 2025). Therefore, a single window T' fixed for both training and evaluation.
RESULTS AND DISCUSSIONS
Table summarizes the main results across all six tasks, reporting mean and standard deviation of test MSE over 4 random seeds. Overall, IGNS, and IGNS consistently outperform the main competitor MGN, as well as other baselines, across all tasks. Notably, GraphCON with the geo- metric encoding also performs strongly on static graph tasks (Plate Def., Cylinder Flow, KS) yet fails severely on the remaining ones. This suggests that while the oscillatory bias in GraphCON im- proves standard graph convolution methods, its formulation limits its expressiveness for more complex dynamics. To gain further insight, we provide an analysis of the relationship between GraphCON and our models in Appendix C Furthermore, introducing time-variation in the weight matrices (IGNS) produces further gains, suggesting that flexible, time-dependens parameterizations better capture complex, nonstationary dynamics. We now analyze per-task behavior below.
Table 2: Comparisons with additional strong baselines on Plate Deformation, Impact Plate and Suramoto-SSvashing tasks. Best results are highlighted in orange.
Method
MSE
Method
MSE
Method
MSE (x 10-3)
MGN-rewiring
3.72 + 0.15
HCMT
646.81 + 27.87
FNO-RNN
6.69 + 0.46
MGN-shared
8.21 + 1.98
MGN
3095.75 + 908.58
MGN
10.76 1 9.16
MGN
1.27 1 0.06
IGNS, IGNS
75.99 17.17
IGNS, IGNS
0.86 + 0.12
IGNS
1.43 0.04
67.86 + 14.75
1.20 0.18
IGNS
1.34 0.05
(a) Plate Deformation
(b) Impact Plate
(c) Seramoto-SSvashinsee
Long-range propagation. In Plate Deformation (final-state prediction only), the Nort-Hamiltniin. based models (IGNS,, IGNS) and GraphCON achieve low error, while A-DGN and GCN+LN per- form substantially worse. Interestingly. MGN also performs well on this task. To better understand this result, we further evaluate two MGN variants: MGN-shared, which uses a shared processor across steps, and MGN-rewiring, which connects each force node to every second node in the mesh. As shown in Table (a), both variants yield higher error than standard MGN, suggesting that MGN's strong performance primarily due to its non-shared processor. Although this design increases the parameter count significantly (x#MP Layers), it allows the model to overfit to geometric details. Additional visualizations supporting this observation provided in Appendix!!
For Impact Plate, standard MGN fails by a large margin, whereas IGNS attains the lowest error, outperforming HCMT by nearly order of magnitude. This result highlights that the energy-
8

---

Published as a conference paper at ICLR 2026
IGNS I MGN
MSE-10
MSE-full
IGNS IINS
1.5
0.1
MSE
0.5
0.05
100
200
400
800
5
10
30
50
100
WB-1 WB-2 CF
Warmup Steps
(a) Data efficiency.
(b) Varying warmup steps.
(c) Longer horizon.
Figure 3: Ablation studies. (a) Data efficiency the final-state-only Plate Deformation task. (b)
Effect of warmup steps on the full-rollout Plate Deformation task. (c) Longer--orrree horizon performance on Sphere Cloth (SC) with T 100. Wave Balls (WB) with two variants: WB-1 (T 100) and
WB-2 (7 = 200), and Cylinder Flow (CF) with T = 100. For better visualization in the bar plots, MSE values for Wave Balls and Cylinder Flow scaled by x100 and x x10, respectively. All results are averaged over 4 seeds.
preserving property of Hamiltonian dynamics enables effective long-range transmission of elastic energy between distant nodes, without requiring specialized hierarchical architectures. Notably, we find that GraphCON's loss consistently diverges after a certain number training iterations across all seeds, indicating that its fixed oscillatory bias may be too restrictive to capture the dynamics present in this task. Note that we omit MP-ODE results on these two tasks, as the method was not proposed for static graphs, and for Impact Plate, it fails converge during training.
Complex dynamics analysis. We now analyze performance across the remaining tasks in Table I that involve complex dynamics and oscillatory behavior. For Sphere Cloth, IGNS achieves the low- est error, while the other models yield errors of similar magnitude, with MGN exhibiting the largest error. In Wave Balls, proposed models (IGNS, IGNS;) significantly outperform all baselines, suggesting that the port-Hamiltoiian structure, viewed as a generalization of wave equations (Ap- pendix C is particularly well-suited for this type of wave dynamics.
and GraphCON models perform comparably, while A-DGN, GCN+LN, and MP-ODE show larger errors under the causal-only setup used here. Notably, in the original paper (Pfaff et al., (02II, MGN performs better due to longer autonomous portions in the dataset; our focusee split reduces that advantage (Appendix Finally, for Kuramoto- Sivashinsky, the chaotic dynamics heavily penalize autoregressive MGN (resulting in large error), while GraphCON, IGNS,/IGNS, and A-DGN remain comparatively strong; GCN+LN struggles to represent the emergence of high- figh-freeuency modes towards the end of the evolution. Although Fourier Neural Operators (FNO-RNN) [i et al. (202I) often perform well on this type of dynamic, they underperform significantly compared to graph-based models, as shown in Table 2 (c). This is likely because, unlike the usual setup with uniformly random initial states, here we initialize with three Gaussian sources randomly placed in the grid to encourage local propagation. Further discussion and visualizations for all tasks are provided in Appendix !
For Cylinder Flow, both port-Hamiltonian-baee
Table 3: Ablations of IGNS SphereCloth-direct. We report test MSE, values scaled by x x103 (a) Core components: geometric encoding, warmup, and multi-step loss. (b) Damping and forcing in IGNS. (c) Damping and forcing in IGNS.
Method
MSE (x10-3)
Method
MSE (x10-3
Method
MSE (x 10 10
IGNS niogoo
36.4 + 3.0
IGNSmood.noof no-d_no-f
36.4 + 1.7
IGNS Saodd,aao
80.6 4.3
IGNS Soowarmmp
27.6 + 0.5
IGNS
25,6 0.3
IGNSu. Suidod
27.8 0.5
IGNS no multi
2388.6 981.0
IG
34.2 + 0.8
IGNS a
60.5 + 16.4
IGNS
25.3 + 1.0
IGNS
28.2 0.8
(a) Core components
(b) IGNS variants
(c) IGNS; variants
Ablation on IGNS components. We ablate IGNS's components on a harder variant, Sphere direct, where the ball connects only the four cloth corners (the original task adds extra connections

---

Published as a conference paper at ICLR 2026
to every-four nodes of the cloth). This stresses long-range propagation: signals must travel from the ball through sparse corners before spreading across the cloth. Table 3 (a) shows that removing geometric encoding or the multi-step objective prevents reliable learning this task, showing these two pieces are necessary for PDE matching; by contrast, dropping the warmup causes only a minor degradation, and we study it further below. We then ablate the nomt-tamitonnian core on Table (b) and (c): removing damping has only minor effect, whereas removing external forcing severely degrades performance. This highlights the central role of forcing in driving information through the graph and supports our universality result in Appendix E.1 We extend this discussion with additional diagnostics and results in Appendix !|
Additional Ablation Studies. Runtime analysis is in Appendix H (a) Data efficiency. We compare
IGNS and MGN on Plate Deformation with varying training set sizes (100, 200, 400 and 800 sam- ples). As shown in Fig. 3a IGNS consistently outperforms MGN across all dataset sizes, with the performance gap widening as the training set decreases. This suggests that the inductive bias from port-Hamillonian dynamics helps IGNS generalize better from limited data, while MGN leads to geometric overfitting when data is scarce. (b) Warmup steps. Next, we study the effect of varying the number of warmup steps used during training for IGNS on Plate Deformation with full-rollout prediction. Fig. shows two metrics: cumulative MSE over the first 10 steps (MSE- and full- rollout validation Toss. Both improve increases. Moving from = 1 to | 5 yields the largest gain (since at / = 1, the task already requires global propagation across the plate); after = 30 the full-rollout loss largely plateaus, while MSE-10 continues to improve with larger indicating larger warmup reduces early error directly, while smaller forces the model to learn dynamics that recover from early mistakes later (visualizations are shown in Fig. 13). (c) Longer horizon. Finally, we compare IGNS and IGNS& on longer rollout horizon setting for Sphere Cloth (T = 100), Wave Balls (T 100, T 200), and Cylinder Flow (T = 100) vs. T 50 (Sphere Cloth, Wave Balls) and T = 30 (Cylinder Flow) in the main experiments. In Sphere Cloth, we additionally increase the number of nodes to 29 x 29 (vs. 19 x 19 in the main experiment) to increase complexity. For Wave Balls, due to the complexity of the task, follow the curriculum scheduling proposed by HHnn et al.) (022):Lieneee & Günnemann (2022). Here, we linearly increase the supervised window by steps from minimal window To to 7 after every epochs. To, and n set as hyperparameters. Mean- while, for Cylinder Flow with T = 100, we adopt the attentional aggregation eLi al. 2019b) in the forcing term to further enhance expressiveness. As shown in Fig. 3c, IGNS outperforms IGNS, on those tasks, with lower std. These results indicate that making weighted matrices time-varying helps increase the model's expressiveness and therefore, enables capture more complex dynamics.
CONCLUSION
In this work, we presented IGNS, graph-based neural simulator that leverages port-Hamiltonian dynamics to improve long-range information propagation and stability in physical modeling prob- lems. By integrating warmup initialization, geometric encoding, and a multi-step training objective, IGNS learns to match the solution trajectories of its underlying port- Hamt-Hamilnnnian dynamics to those of ground-truth PDE solvers, keeping the predicted states on-trajectory and thereby mitigating er- ror accumulation while enhancing long-range dependencies in both spatial and temporal domains. IGNS achieves consistent improvements over strong baselines across six tasks spanning both solid and fluid mechanics, including new benchmarks designed to test complex and long-range dynamics. Our theoretical analysis further shows that the Hamiltonian structure preserves gradient flow, while the nort-Hamiltonia extension enables the model to capture broader class of dynamics through its universality property. Together, these results highlight IGNS as a principled and stable framework for neural simulation of complex physical systems.
Limitations and Future Work. Our method currently operates in an open-loop, causal setting, where sources are specified at the start and the model predicts forward without feedback. However, many real-world applications require closed-loop feedback, in which the model continuously adapts its predictions based on observations obtained from sensor measurements or other exogenous in- puts. Extending IGNS in this direction is therefore natural and meaningful next step. Furthermore, IGNS's ability to propagate information over long ranges naturally opens the door to inverse de- sign (Allen et al. 2022) and control tasks (Hoang et al. (025), where the goal is to identify optimal parameters (e.g., source locations or inputs) that achieve a desired outcome by optimizing inverse objectives. Exploring these directions represents promising avenue for future work.
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
This work is part of the M2 subproject of the DFG AI Research Unit 5339, which focuses on learning dynamical process models based on data and expert knowledge to accelerate manufacturing process maturation. We gratefully acknowledge financial support from the German Research Foundation (DFG) and computational resources provided by bwHPC and the HoreKa supercomputer, funded by the Ministry of Science, Research, and the Arts Baden- Waalen-Württebberg and the German Federal Ministry of Education and Research. We would like to thank Johannes Mitsch and Tobias Würth for developing the data generation platform for the Plate Deformation task in Abaqus. AT, AG, and DB acknowledge funding from EU-EIC EMERGE (Grant No. 101070918).
ETHICS STATEMENT
The research conducted in this paper conforms in every aspect with the ICLR Code of Ethics. Our study does not involve human subjects, sensitive personal data, or applications with foreseeable harmful consequences. No ethical concerns anticipated regarding data usage, methodology, or findings.
REPRODUCIBILITY STATEMENT
We provide all necessary details to implement our IGNS in Section 5 and describe the choices of hyperparameter for each experiment Appendix G and full dataset description in Appendix F therehy ensuring sufficient information to reproduce our results. The full codebase and datasets are released with the camera-ready version to ensure the reproducibility, and are also available at https://thober wtr.bbtttss.g.ht........... 1o/ wtwp//uurllppp. matching
REFERENCES
Kelsey R Allen, Tatiana Lopez-Guavara, Kim Stachenfeld, Alvaro Sanchez-Gonzalez, Peter Battaglia, Jessica Hamrick, and Tobias Pfaff. Inverse design for fluid-structure interac- tions using graph network simulators. In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho (eds.), Advances in Neural Information Processing Systems, 2022. URL https: Tttt::://openreveie net/ www.:neeee/ooeeeeeeeeeee.o Hoorumiid-Hazuqjj0ee
Uri Alon and Eran Yahav. On Bottleneck of Graph Neural Networks and its Practical Im- plications. In International Conference on Learning Representations, 2021. URL https:
//openreview.
Vladimir Igorevich Aroold d. Mathematical methods of classical mechanics, volume 60. Springer Science & Business Media, 2013.
Álvaro Arroyo, Alessio Gravina, Benjamin Gutteridge, Federico Barbero, Claudio Gallicchio, Xi- aowen Dong. Michael Bronstein, and Pierre Vandergheynst. On Vanishing Gradients, Over- Smoothing, and Over-Squashing in GNNs: Bridging Recurrent and Graph Learning. In The Thirty-ninth Annual Conference on Neural Information Processing Systemx, 2025. URL https: www.:et...nmnmm 4=N4cyyMul.y
7 //openreview. wwwprrewwwwww..e.....e..aee
Uri M. Ascher. Numerical Methods for Evolutionary Differential Equations. Society for Industrial and Applied Mathematies, USA, 2008. ISBN 0898716527
Federico Barbero, Ameya Velingker, Amin Saberi, Michael M. Bronstein, and Francesco Di Giovanni. Locality-aware graph rewiring in GNNs. In The Twelfth International Confer- ence on Learning Representations, 2024. URL
net fttt//formmidd
4Ua4hXiAJX
Peter W. Battaglia, Jessica B. Hamrick, Victor Bapst, Alvaro Sanchez-Gonzales Vinicius Zambaldi, Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, Caglar Gulcehre, Francis Song, Andrew Ballard, Justin Gilmer, George Dahl, Ashish Vaswani, Kelsey Allen, Charles Nash, Victoria Langston, Chris Dyer, Nicolas Heess, Daan Wierstra, Pushmeet
11

---

Published as a conference paper at ICLR 2026
Kohli, Matt Botvinick, Oriol Vinyals, Yujia Li, and Razvan Pascanu. Relational inductive biases, deep learning, and graph networks, 2018. URL https:://////... http//awxwvv...oraa.//
012
Saakaar Bhatnagar, Yaser Afshar, Shaowu Pan, Karthik Duraisamy, and Shailendra Kaushik. Pre- diction of aerodynamic flow fields using convolutional neural networks. Computational Mechan- ics, (4(2))525-545, Aug 2019. ISSN 1432-0924. doi: 10. 10.1007760046--01--017776 017-0000
URL
https://dol www./////ooo 10. 1
Johannes Brandstetter, Daniel E. Worrall, and Max Welling. Message passing neural PDE solvers. In International Conference on Learning Representations 2022. UtI.https::// httttttss /otttp:://penrevieer
net/ fo
Susanne C Brenner and Ridgway Scott. The mathematical theory of finite element methods,
volume 3. Springer, 2008.
Chen Cai and Yusu Wang. A note over-smoothing for graph neural networks. arXiv preprint arXiv:2006. 13318, 2020.
Yadi Cao, Menglei Chai, Minchen Li, and Chenfanfu Jiang. Efficient learning of mesh-based phys- ical simulation with bi-stride multi-scale graph neural network. In International Conference on Machine Learning. 2023. | http:: h/tptpes///ppeneviee net Eorummirrmi 2idermmmid-MMbooo IEt
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary dif- ferential equations. In S. Bengio, H. Wallach, H. Larochelle. K. Grauman, N. Cesa-Bianchi, and R. Garnett (eds.), Advances in Neurul Information Processing Systems, volume 31. Cur- Associates, Inc., 2018 http::ss///roeeeeggg...eeeee neurips wwtp:eerer.......ppppree..e files/ paper/2018 http://001//11///////////ee/eeee http:////......teeeeeeeeeaaaaaa.
Philipp Dahlinger, Tai Hoang, Denis Blessing, Niklas Freymuth, and Gerhard Neumann. MaNGO adaptable graph network simulators via meta-learning. In The Thirty-ninth Annual Confer- ence on Neural Information Processing Systems, 2025. URL /ttp://eeeeeeeeeeeeeeeeeeee net./
Huayu Deng, Xiangming Zhu, Yunbo Wang, and Xiaokang Yang. EvoMesh: Adaptive physical simulation with hierarchical graph evolutions. In Aarti Singh, Maryam Fazel, Daniel Hsu, Si- mon Lacoste-Julien, Felix Berkenkamp. Tegan Maharaj. Kiri Wagstaff, and Jerry Zhu (eds.), Proceedings of the 42nd International Conference on Machine Learning, volume 267 of Pro- ceedings of Machine Learning Research, pp. 13327-13346. PMLR, 13-19 Jul 2025. URL httpp::/prroeeedings. mlr prw.s/pesssee..................e Seress/v26777 w2/ /etp///den2...... .html
Shaan A. Desai, Marios Mattheakis, David Sondak, Pavlos Protopapas, and Stephen J. Roberts. hort-tamilnnian neural networks for learning explicit time-dependent dynamical systems. Phys. Rev. E. 104:034312, Sep 2021. doi: 10. 10.1103/PhysRev.. w.w.ReEE......44333. 104.034312. URL hhttps://iim 104. 104.00.33......33122
aps. wes.oorgddo///
Francesco Di Giovanni, Lorenzo Giusti, Federico Barbero, Giulia Luise, Pietro Liò, and Michael Bronstein. On over-squashing in message passing neural networks: the impact of width, depth, and topology. In Proceedings of the 40th International Conference on Machine Learning, ICML 23. JMLR.org, 2023.
J.R. Dormand. Numerical Methods for Differential Equations: A Computational Approach. Engi- neering Mathematics. Taylor& & Francis, 1996. ISBN 9780849394331.
Moshe Eliasof, Eldad Haber, and Eran Treister. PDE-GCN: Novel Architectures for Graph Neu- ral Networks Motivated by Partial Differential Equations. In Advances in Neural Information Processing Systems, volume 34, pp. 3836-3849. Curran Associates, Inc., 2021.
Moshe Eliasof, Alessio Gravina, Andrea Ceni, Claudio Gallicchio, Davide Bacciu, and Carola- Bibiane Schönlieb. Graph Adaptive Autoregressive Moving Average Models. In Proceed- ings of the 42nd International Conference Machine Learning, volume 267 of Proceedings of Machine Learning Research, pp. 15232 15232-15265. PMLR, 13-19 Jul 2025. URL https: //proceed ngs. ml http:ssmw..pppas////emmm wwwppsawee..7
asof25a .ttml
12

---

Published as a conference paper at ICLR 2026
Lawrence C Evans. Partial Differential Equations. Graduate studies in mathematics, American Mathematical Society, Providence, RI, 2 edition, March 2010.
Clara Lucia Galimberti, Liang Xu, and Giancarlo Ferrari Trecate. unified framework for hamilto- nian deep neural networks. In Ali Jadbabaie, John Lygeros, George J. Pappas, Pablo A. P.rrario,, Benjamin Recht, Claire J. Tomlin, and Melanie N. Zeilinger (eds.), Proceedings of the 3rd Con- ference Learning for Dynamics and Control, volume 144 of Proceedings of Machine Learning .mlr.
Research. pp. 275-286. 286. PMLR, 07 08 June 2021. URL press/v144 44 ga imbert 21a html
Clara Lucia Galimberti, Luca Furieri, Liang Xu, and Giancarlo Ferrari-Trecate. Hamiltonian deep neural networks guaranteeing nonvanishing gradients by design. JEEE Transactions on Automatic Control, 68(5):3155-3112) 2023. doi: 10. 10.11/9A..2222333399999
Johannes Gasteiger, Stefan Weiß enberger, and Stephan Günnemann. Diffusion Improves Graph Learning. In Advances in Neural Information Processing Systems, volume 32. Curran Associates,
Inc., 2019.
Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl. Neural message passing for Quantum chemistry. In Proceedings of the 34th International Conference on Machine Learning, volume 70 of ICML'17, pp. 1263-1272. JoLR.org, 2017.
Alessio Gravina and Davide Bacciu. Deep Learning for Dynamic Graphs: Models and Benchmarks. IEEE Transactions on Neural Networks and Learning Systems, pp. 1-14, 2024. doi: 10.1
NNNNS..202.3379755
Alessio Gravina, Davide Bacciu, and Claudio Gallicchio. Anti-Symmetric DGN: a stable architec- ture for Deep Graph Networks. In The Eleventh International Conference on Learning Represen- lutions, 2023, URL htt.htttt///peeeee/iiee
net (tte///oommmndd33/o2o0000
Alessio Gravina, Claudio Gallicchio. and Davide Bacciu. Non-Dissipative Propagation by Random- ized Anti-Symmetric Deep Graph Networks. In International Workshops of ECML PKDD 2023, Turin, Italy, September 18-22, 2023, Revised Selected Papers, Part V, 2024a.
Alessio Gravina, Giulio Lovisotto, Claudio Gallicchio, Davide Bacciu, and Claas Grohnfeldt. Long Range Propagation on Continuous- Time Dynamic Graphs. In Proceedings of the 41st Interna- tional Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Re- search, pp. 16206-16225. PMLR, 21-27 Jul 2024b. URL http::s://pocceednngg mmlr
press/ 00533/v235/42 2ra.i///raiia222
html
Alessio Gravina, Daniele Zambon, Davide Bacciu, and Cesare Alippi. Temporal graph odes for irreularty-sampled time series. In Kate Larson (ed.), Proceedings of the Thirty-Third Inter- national Joint Conference on Artificial Intelligence, IJCA/-24, pp. 4034-4034. International Joint
Conferences on Artificial Intelligence Organization, 2024c. doi: 22249963/jcai.22224444 URL https://dol ww..orr//// / 10. 24963/291413 icai. 2024/445)
Alessio Gravina, Moshe Eliasof, Claudio Gallicchio, Davide Bacciu, and Carola-Bibiana Schönlicb. On oversquashing in graph neural networks through the lens of dynamical systems. In The 39th Annual AAAI Conference on Artificial Intelligence, 2025.
Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks. Advances in neural information processing systems, 32, 2019.
Xiaoxiao Guo, Wei Li, and Francesco Torio. Convolutional neural networks for steady flow ap- proximation. In Proceedings of the 22nd ACM SIGKDD International Conference on Knowl- edge Discovery and Data Mining, KDD '16, pp. 481 4-4-4 New York, NY, USA. 2016. Associ- ation for Computing Machinery. ISBN 9781450342322.. doi: 10.11452299967222999788. URL
Ernst Hairer, Gerhard Wanner, and Syvert P Nørsett. Solving ordinary differential equations I:
Nonstiff problems. Springer, 1993.
Emst Hairer, Marlis Hochbruck, Arieh Iserles, and Christian Lubich. Geometric numerical integra- tion. Oberwolfach Reports, 3(1):055-822, 2006.
13

---

Published as a conference paper at ICLR 2026
Xu Han, Han Gao, Tobias Pfaff, Jian-Xum Xun Wang, and Liping Liu. Predicting physies in mesh- reduced space with temporal attention. In Interaational Conference on Learning Representations, 2022.
Ali Hariri, Alvaro Arroyo, Alessio Gravina, Moshe Eliasof, Carola-Bibiane Schönlieb, Davide Bacciu, Xiaowen Dong, Kamyar Azizzadenesheli, and Pierre Vandergheynst. Return of Cheb- Net: Understanding and Improving an Overlooked GNN on Long Range Tasks. In The Thirty- ninth Annual Conference on Neural Information Processing Systems, 2025. URL https:
//openreview.
Simon Heilig, Alessio Gravina, Alessandro Trenta, Claudio Gallicchio, and Davide Bacciu. Port- Hamiltonian Architectural Bias for Long-Range Propagation in Deep Graph Networks. In The Thirteenth International Conference on Learning Representations, 2025. URL https:
//opnrevveww http:///oooommmoo000ooeee33
Tai Hoang, Huy Le, Philipp Becker, Vien Anh Ngo, and Gerhard Neumann. Geometry-aware RL for manipulation of varying shapes and deformable objects. In The Thirteenth International Confer- ence on Learning Representations, 2025. URL httph:://///........ http///ppeeeeeeeeeeeeeeeee www.:nnnttfoororooe i /
7BLXhmWvwFF)
Valerii Iakovlev, Markus Heinonen, and Harri Lähdesmäki. Learning continuous-time (pde)s from sparse data with graph neural networks. International Conference on Learning Representa- tions, 2021. URL /ttttps://ppereeeeet
nnt/ http:/ooummmmaaaaam add
Heiner Igel. Computational Seismology: Practical Introduction. Oxford University Press, edition, 2016. ISBN 9780198717409. URL httptttstt/:glll.. www.:uwpp.ccccacaaaeeaa c/
prrouctt/computatte
see otw..ooyyyy777111177777777771111 97 www.7//1188117777ccdddd.........
Steeven Janny, Aurélien Benetteau, Madiha Nadri, Julie Digne, Nicolas Thome, and Christian Wolf. Eagle: Large-scale learning of turbulent fluid dynamics with mesh transformers. In International Conference Learning Representations (ICLR), 2023.
Thomas N. Kipf and Max Welling. Supi-Supe Classification with Graph Convolutional
Networks. In International Conference on Learning Representations, 2017. URL https:
7/openreview htt.://toooommmdaaS/tttaaaseee
Samuel Lanthaler, Konstantin Rusch, and Siddhartha Mishra. Neural oscillators are universal. In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine (eds.), Advances in Neural Information Processing Systems, volume 36, pp. 46886. Curran Associates, Inc., 2023. URL h pttpps://prooceeedngg.... neurips http://pprrfflllepppaee///////// fi e/ 911/9/22222d80eeee 4eat
pddf
Guohao Li, Mathias Müller, Ali Thabet, and Bernard Ghanem. Deepgens: Can gens go deep as cnns? In 2019 IEEE/CVF International Conference on Computer Vision (ICCV), . pp. 0266-9275, 2019a. doi: 10.1109/CCC22201900009
Yujia Li, Chenjie Gu, Thomas Dullien, Oriol Vinyals, and Pushmeet Kohli. Graph matching net- works for learning the similarity of graph structured objects. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 3835-3845. PMLR, 09-15 Jun 2019b. URL
www.tpress/aee wttpss/////11191/0/ html
Yunzhu Li, Jiajun Wu, Russ Tedrake, Joshua B. Tenenbaum, and Antonio Torralba. Learning particle dynamics for manipulating rigid bodies, deformable objects, and fluids. In International Confer- ence on Learning Representations, 2019c. URL :ttttt:s///ppereviiie net/ / forum? i=
rJgbSn09Ym
Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede liu, Kaushik Bhat- tacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, 2021. URL
wttprrewwwwwwwweeeooooeeee https://openreviie net wttp:/nettooooooooooooooooooooo.
14

---

Published as a conference paper at ICLR 2026
Marten Lienen and Stephan Günnemann. Learning the dynamics of physical systems from sparse observations with finite element networks. In International Conference on Learning Representa- tions, 2022. httL:ttss///peeeeeeaait wttp://wwwweeeeeeeeeeeeeeaa
net/ forum? Lor.k
Jonas Linkerhägner, Niklas Freymuth, Paul Maria Scheikl, Franziska Mathis-Ullrichh and Ger- hard Neumann. Grounding graph network simulators using physical sensor observations. In The Eleventh International Conference on Learning Representations, 2023. URL https:
//openrevien m d=tnnndddjsssssvs87
Yuankai Luo, Lei Shi, and Xiao-Ming Wu. Classic GNNs strong baselines: Reassessing GNNs for node classification. In The Thirly-eighi Conference on Neural Information Processing Systems Datasets and Benchmarks Track, 2024. URL https: httt:p::///peereviiieee net wwwt///ooo/rrte forum? di
xk1 [x1jKdGe4B 4E
Miles Macklin. Warp: A high-performance python framework for gpu simulation and graphics. March 2022. NVIDIA GPU Technology Confer-
https: //gitt/b//ttbbb./ ence (GTC).
Mayank Mittal, Calvin Yu, Qinxi Yu, Jingzhou Liu, Nikita Rudin, David Hoeller, Jia Lin Yuan, Ritvik Singh, Yunrong Guo, Hammad Mazhar, Ajay Mandlekar, Buck Babich, Gavriel State, Marco Hutter, and Animesh Garg. Orbit: A unified simulation framework for interactive robot learning environments. IEEE Robotics and Automation Letters, 86)):3740-3747, 2023. doi:
RQ.1009/LRA..223.22700
Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzalee and Peter Battaglia. Learning mesh-based simulation with graph networks. In International Conference on Learning Representations, 2021. URL
net www.:/et/ooooreee id=roNgY1. XP
Michael Poli, Stefano Massaroli. Junyoung Park, Atsushi Yamashita, Hajime Asama, and Jinkyoo
Park. Graph neural ordinary differential equations. arXiv preprint arXiv:1911.07553 2019.
Andreas Radler, Vincent Seyfried, Stefan Pirker, Johannes Brandstetter, and Thomas Lichtenegger. Paint: Parallel-in-timg time neural twins for dynamical system reconstruction, 2025. URL https:
Andreas Roth and Thomas Liebig. Rank collapse causes svorthhogting and over-correlation in graph neural networks. In Proceedings of the Second Learning on Graphs Conference, volume 231 of Proceedings of Machine Leurning Reseurch, pp. 35:1-35:23. PMLR, 27-30 Nov 2024.
T. Konstantin Rusch and Daniela Rus. Oscillatory state-space models. In The Thirteenth Intern- tional Conference on Learning Representationss 2025. URL
forum? d=GRMFXcAFFh
Konstantin Rusch, Ben Chamberlain, James Rowbottom, Siddhartha Mishra, and Michael Bron- stein. Graph-coupled oscillator networks. In International Conference on Machine Learning, pp- 18888-18909. PMLR, 2022.
T. Konstantin Rusch, Michael M. Bronstein, and Siddhartha Mishra. A Survey on Oversmoothing in Graph Neural Networks. arXiv preprint arXiv:2303. 2023.
Alvaro Gonzalz, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovee, and Peter Battaglia. Learning to simulate complex physics with graph networks. In Proceedings of the 37th International Conference on Machine Learning, ICML'20. JMLR.org, 2020.
Dai Shi, Andi Han, Lequan Lin, Yi Guo, and Junbin Gao. Exposition on over-squashing problem on GNNs: Current Methods, Benchmarks and Challenges. arXiv preprint arXiv:231 2023a.
Haochen Shi, Huazhe Xu, Samuel Clarke, Yunzhu Li, and Jiajun Wu. Robocook: Long-horizor elasto-plastic object manipulation with diverse tools. In Conference on Robot Learning, pp. 642- 660. PMLR, 2023b.
Michael Smith. ABAQUS/Standard User's 's Manual, Version 6.9. Dassault Systèmes Simulia Corp, United States, 2009.
15

---

Published as a conference paper at ICLR 2026
Jake Topping, Francesco Di Giovanni, Benjamin Paul Chamberlain, Xiaowen Dong, and Michael M. Bronstein. Understanding over-squashing and bottlenecks on graphs via curvature. Interna- tional Conference Learning Representations, 2022. URL hattps://openrevieey net/
forum? dd:mmmidd-mmmm J=7Um
Alessandro Trenta, Alessio Gravina, and Davide Bacciu. SONAR: Long-Range Graph Propaga- tion Through Information Waves. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, 2025. Uttp:t:s: wwtp:///e.ppeeeeeee.oeeeoe.
jmc95rl
Arjan Van der Schaft. L2-gain and passivity techniques in nonlinear control. Springer, third edition, 2017.
Arjan van der Schaft and Dimitri Jeltsema. Port-Hamiltonian Systems Theory: An Introductory Overview. 2014. doi: 10.1561/2600000002
Petar Veličković. Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pictro Liò, and Yoshua Bengio. Graph Attention Networks. In International Conference on Learning Representations, 2018. URL BttL:tttst///peeveeleet
net/ o
Amaury Wei and Olga Fink. Integrating physics and topology in neural networks for learning rigid body dynamics. Nature Communications, 16(1), July 2025. ISSN 2041-1723. doi: 10. 1038/ 051466-025-62250- URL/tttp://xx.co// org/ http://1.1111///11 025 62250-7
Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S Yu Philip. A Comprehensive Survey on Graph Neural Networks. IEEE Transactions Neural Networks and Learning Systems, 32(1):4-24, 2020.
Tobias Würth, Niklas Freymuth, Gerhard Neumann, and Luise Kärger. Diffusion- based hi- erarchical graph neural networks for simulating nonlinear solid mechanics. arXiv preprint 2:Xiv:25006.06055, 2025.
Youn-Yea Yu, Jeongwhan Choi, Woojin Cho, Kookjin Lee, Nayong Kim, Kiseok Chang, ChangSe- ung Woo, ILHO KIM, SeokWoo Lee, Joon Young Yang, SOOYOUNG YOON, and Noseong Park. Learning flexible body collision dynamics with hierarchical contact mesh transformer. In The Twelfth International Conference on Learning Representations, 2024. URL https:
Youn- Yu, Jeongwhan Choi, Jaehyeon Park, Kookjin Lee, and Noseong Park. PIORF: Physics-
informed ollivier-ricci flow for long-range interactions in mesh graph neural networks. In The Thirteenth International Conference on Learning Representations, 2025. URL https:
/oopeneview net. www.:/neee/tiimiieeoeee q
Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty. Symplectic ode-net: Learning
hamiltonian dynamics with control. International Conference on Learning Representations, 2020. URL htt.:httt///oeeeet
wttp:/eet/.oouttietee.mmm wtt.:/ooommmmmmdmmmmmma rKDS
David Zwicker. py-pde: python package for solving partial differential equations. Journal of Open Source Software, 5(48), 2020. doi: 10.21 105/j0ss.022158.. URL https: /ttts:::/oii.. //doi.
02158
16

---

Published as a conference paper at ICLR 2026
A LLMS USAGE
Large Language Models (LLMs) were used as general-purpose assistive tools to improve the writing quality of this paper. Specifically, we used LLMs help with grammar correction, rephrasing for clarity, and suggesting some improvements to the overall structure of the text. All LLM- -ggeerrtt text was carefully reviewed and edited by the authors to ensure that it accurately reflects the authors' intentions and scientific content. No LLMs were used to generate scientific content, including but not limited to research direction, hypothesis formulation, experimental design, data analysis, or interpretation of results.
B RELATED WORK
Graph Neural Simulators. Graph Neural Networks have been widely adopted in machine learn-
ing for physical simulation tasks (Battaglia et al. [0118] Plaff et al. 2021). Notable examples include modeling particle interactions (Sanchez-Gonzate,
al. 2020 Li al. 2019c; Sanzaez-Gonzalez 2020). fluid dynamics (Pran et al. 2021 Yu et al. 2025), and deformable materials (Yu
et al.
al.
2024 Linkerhägner
aa.3) Here, given the observations at the current time step, one
can predict the next states by applying few iterations of message- passsegaasing steps to incorporate neighborhood information, and output the next state. This underlying message-passinn
step of-
ten implemented with graph convolutional layers (Kipf & Welling (017)) or attention-based meth- ods (Veličković et al.) 2018), or more generally, message-passing framework (Battaglia et al. 2018)) During rollout, the model can generate the whole trajectories in an autoregressive manner (Plafi et al. 2021 Wei & Fink 2025). However, these models often struggle with long-term accuracy due error accumulation over long time horizons. In contrast, another line of work focuses on directly learning the underlying continuous dynamics using neural ODEs (Chen et al. (018). which helps mitigate the error accumulation problem by enabling the model to output whole trajectories in a single forward pass and better capture the continuous nature of physical systems. Recent works extend this idea for graph-based architectures (Iakovlev al., 202 Lienen & Günnemann (022).
These models are often trained with multi-step objectives, showing accuracy improvement and with the additional ability to handle irregular time series and improve generalization. However, these models are often limited to short-term predictions (Lienen & Günnemann, (202) Janny et al.) (2023)
or supervision in a reduced latent space (Han et al. (0022) due to the challenges in propagating information over both long spatial and temporal horizons.
Effective propagation in Graph Neural Networks. Effectively propagating information across the graph, and thus effectively modeling long-range dependencies, remain critical challenge for GNNs (Shi et al. 2023a; Rusch et al. 2023 Roth & Liebig 2024). Stacking many GNN layers to capture long- long-range dependencies often leads to issues such as over-smoothing (Cai & Wang. 2020 Rusch et al.] 2023), when node states become indistinguishable as the number of layers in- creases. To address this, methods like residual connections (Li et al. (019a),, normalization tech- niques (Luo et al. (2024), and advanced architectures (Rusch et al.) 2022, Eliasof et al.) [222] have been proposed. Another issue that prevent effective long-range propagation is the over-squashing phenomenon (Alon & Yahav 102 Topping al. 2022: Di Giovanni et al. 2023), which arises from topological bottlenecks that squash exponentially increasing amounts of information into node states, causing loss of information between distant nodes. To overcome this, variety of strategies have been proposed. For example, graph rewiring methods (Topping et al. 2022 Gasteiger et al.) 2019 Barbero et al. (2024) modify the graph topology facilitate communication. In physical sys- tem Simulation, hierarchical (Deng et al. 2025 Cao et al.) 2023 Janny et al., 2023) and rewiring methods, including HCMT (Yu et al. 2024) and PIORF Yu et al. 2025). improve the global in- formation flow by introducing shortcuts or multi-scale pooling. However, these approaches can make models overly tied to the mesh structure, which can result in geometric overfitting. More recently, AArroyo et al. (2025) demonstrated that both over-smoothing and over-squashing problems are closely linked to the problem of vanishing gradients in message passing, and works like Eliasof et al. (202I): Gravina et al. (2023 2024a; 2025); Heilig et al. (2025): ; Arroyo al. (2025); Trenta
et al. (2025): Hariri et (2025 Blasof et (2025) have explored this perspective to design GNNs
that preserve the gradient flow over many layers. However, most of these methods focus on static
graphs and node classification tasks, and it is still unclear how to extend them to physical simulation
tasks.
17

---

Published as a conference paper at ICLR 2026
Neural ODEs and Hamiltonian Dynamics. Neural Ordinary Differential Equations (ODEs) have been widely used to model continuous- time dynamics (Chen et al., 2018). offering flexibility in time integration and improved generalization, and have also been extended to graph-structured data Poli et al. (211)); Gravina et al. c). There has been significant interest extending neu- ral ODEs capture physical systems by drawing inspiration from Hamiltonian dynamics (Amol'' 2013). Most closely related to our work is HOGN (Greydanus et al.) 2019). which parameter- izes Hamiltonian dynamics with graph-based message passing, but does explicitly study long- range spatio-temporal dependencies. Other notable examples include Hamiltonian Neural Networks (HNNs) (Greydanus al.) 2 and Symplectic ODE-Net (Zhong et al.) 2020p, which achieve en- ergy conservation and stable long-term predictions by learning parameterized Hamiltonian func- tion. This conservation property has also been explored as way to mitigate gradient vanishing in deep neural networks (Galimberti et al., 2023 Lanthaler et al. 2023 Rusch & Rus 2025). and it has inspired new architectures for enhancing long-range propagation in GNNs (Heilig et al., 2025) as well as long-horizon forecasting in state-space-model (Rusch & Rus, 2025). However, these works remain different from ours, as they focus either on static graphs or purely temporal sequences.
In this work, propose novel GNS that leverages dynamical systems theory to preserve the gradient flow across distant nodes in spatial domains. Trained with multi-step loss, this model can mitigate error accumulation and hence, enable accurate long-term physical simulation.
C DERIVATIONS
In this section, establish the connection between the port-Hamiltonian formalism in Eq. (6) and the PDE wave equation defined on a space-time domain. In particular, we show that both can be expressed in the form of mass-spring damneer system. These two links make clear that port-Hamiltonian formalism induces a second-order inductive bias, which can be interpreted as a weak form of underlying
equation.
C.1 IGNS MASS-SPRING- WORRING-DAMPER DAMPER ODE
To begin with, we first consider the following Hamiltonian function:
1 ーー M-'p-tt 2 -g"Kg Kg
(t, q.
where define the generalized coordinates q and momenta p Mq. and let =1g.p] !q,p/=
In if we set ō(x) x and
fact, this definition is consistent with the Hamiltonian defined Eq. 7
W diag ([K1/2 M 1 and V = 0. The update of and is then given by:
q(t) = V vp (t, q.p) = -mt)((t)
p(t) (t. -t).
We now eliminate p(t) by substituting p(t) M q(f). and add the damping term D q(f) and the external force rt,, q). yielding:
M @(t) + D q(t) + K q(t) = rt)
(1I)
This recovers the mass-spring-damper system.
For more general case, consider Hamiltonian H(),q.) E C² in all arguments. Taking the time derivative of q (t, q, p) gives
2 11 2 pq a.VHH. pp
Assuming H exists and is invertible (which holds for certain activation functions, e.g. a(:)
tanh(-) and regularity conditions for invertibility), define M(,, q.p) H. -1 Solv- ing for p.
P ( q. p) q 2 H 0, V, We now use the Hamiltonian equation for p with damping and external forcing from Equation (6)
(12)
p v D4 - 9)
(13)
18

---

Published as a conference paper at ICLR 2026
Substituting Equation (2) into Equation 13 and rearranging, arrive at the general second-order
form:
Mt..q.p)
D ((t....)) g+ VgH = r(t. + ( q.P..
(14)
mass
damping
stfflness
external force & residual
This analysis demonstrates that, for general Hamiltonians, the resulting second-order dynamics exhibit state- and time-dependent mass, damping, and stiffness terms. In particular, for separa- ble Hamiltonians, have V H 0, and for time-invariant Hamiltonians, O, V. 0. In these cases, the dynamics reduce the form of Eq. with Mtt,q, q. pp H) and V qp)..))
Comparison with GraphCON. We can further observe that the baseline GraphCON (Rusch et al. 2022 is in fact simpler instantiation Eq. (TT). Its update takes the form
q(t) + (t( + (t) = o (R&(t, (t)
where Rg implemented as a standard graph convolution or attention operator. This corresponds directly mass-spring-dant dariner system with M I, D al, K = and external forcing r-: G) = a (Ro (t, q)). Notably, the mass matrix M is the identity, meaning that coupling between nodes occurs only through the external forcing term r, whereas in our formulation coupling also appears on the left-hand side through M, D, and K.
C.2 MAAS-SPRINGGDAMAAE ODE AS A SECOND PDE
In this subsection, we show how a second-order PDE with damping and forcing reduces to a finite- dimensional mass-spring-damie damper ODE via Galerkin projection. Interested readers are referred to Igel 2016) for further details.
Let Z C
Rd be bounded domain and consider the following damped, forced equation
u. + 8. Lu g(#,(), (æ: E Ω (0,T),
(15)
where the solution is u x (0, 7] - R. i.e., scalar field. We set - with constant wave speed e > 0, where A is the spatial Laplacian, and define B as damping operator (e.g., a7) acting on the time-derivative term d.u.
Let V C (S) be the test space. To obtain the weak form of Eq. v E V, multiply the strong equation by v. and integrate over S2:
15 take test function
[ Ow u dx + S( BBquddx - 2 AAuudx 1/ g(: t) v dx.
The first two terms already appear in weak form. For the Laplacian term, we apply integration by parts and assume that u(t) vanishes at the boundary. This gives the following bilinear form:
alu,v) v) = -c2 Auvda c² / Vu Vudæ.
Altogether, the weak form reads: for all VV
(dnu, v) + (Bā,u, u, v) + a(u, v) = (g(. t), v)
(16)
where (S.g) = Ja (xx) g(x) dæ denotes the (22 inner product.
Choose finite-dimensional subspace V,, = span(on......... Φn) V and rewrite the solution as
un un( 1)
Σ 4 (w))
(17)
where (t) time-dependent coefficients to be determined. Since Eq. (16) holds for all v di.
i 1,... k. rewriting using the Eq. (17) yields
Σ (e,, (t) + BB oi) (t(t) + (()) (l(.
(18)
19

---

Published as a conference paper at ICLR 2026
To show the equivalent mass-spring-danme damper form, now define the following matrices and vectors:
My = (Øj, di),
(mass matrix)
Dij = (B di),
(damping matrix)
Ky a( Di);
(stiffnesss matrix)
ri(:t) ()= (o(: 1), Qi);
(external forcing)
with vector q = 1qyly thus Eq. (18) becomes
M q()) + D (()) + q(1) = r(;
Interpretation. Because our latent dynamics (Eqs. O and (14) have the same form as the sys- tem above, IGNS can be seen as learning effective mass, damping. and stiffness operators, consis- tent with weak-form second-order PDE (in a PI nodal basis, coefficients equal nodal values, i.e.. ul ((x,t) = (t))). This method-of-lime view is related to Finite Element Networks (FEN) Lienen & Günnemann 2022), which were the first introduce the connection between message passing and the finite element method. However, while FEN focuses on diffusion/heat dynamics (parabolic, first order in time), here model more general equation with forcing and damping, yielding a more general second-order hyperbolic PDE.
D TIME-VARYING WEIGHT MATRICES
To further enhance the model's expressiveness, we propose the following parameterization to make (t) and (t) in Eq. (7) time-varying, allowing the Hamiltonian to adapt over time, which is particularly useful for modeling more complex, non-stationary dynamics.
of (t),, as (tt) can be treated similarly. Assuming W(t) is a square matrix of size d x d, we decompose it into symmetric and skew-symmetric parts, i.c.
Here, consider only the
((I) S(1 (t), where
S11
812
8id
Q12
Gid
812
822
82d
012
0
02d
S
A
1idd
82d
sdd_
L-aid
-
02d
with 8ij 8ji. Oij Ojjs and as time-varying components
0.
We then factor these parts through time- invariant and
5(1) = U, diag(Ta(())...
A(1) - diag((To()) - (i))
where U,, cient vectors output by MLPs with parameters 0. This two-step construction helps keeping all the learned parameters unconstrained while reducing the parameter counts to O(d), compared to the naive approach of letting MLP output all the matrix entries, which leads to d2 parameters. In practice, we implement (t), (t) as small MLPs that take time-embedding t as input and output the corresponding coefficient vectors.
Ua, Pa are learned (time-invariant) bases, and (t),,, (t( are two time-varying coeffi-
E PROOFS OF THE THEORETICAL STATEMENTS AND FURTHER DISCUSSION
In this section, prove the theoretical statements in this paper together with more in-depth discussion.
E. NEURAL OSCILLATORS AND UNIVERSALITY
We start this section by recalling the definition of neural oscillators and their universality property (Lanthaler et al. )023). The general setting involves learning an operator Φ u + (uu), that maps time-dependent input signal u(l)
output function (())(() Ra. Here, & defined
continuous operator (w.r.t. the 00 norm) that maps the space
Co( 7] RP) {u [0, T]
TRP vu(t) is continuous and (0))= 0)
(19)
20

---

Published as a conference paper at ICLR 2026
into itself. The assumption u(0) 0 is not restrictive, as the main paper showed that the proposed neural oscillators can operate even in the case u(0) is non-zero. Here, by adding an arbitrarily small time interval to "warmup", the oscillator starting from the rest state synchronize with this non-zero input signal. The only remaining assumption is that O needs to be causal, in the sense that the value of (t))(t) does not depend future values (u(r)]7 > t), which matches our setup considered in this paper well. Formally speaking, & is causal if giving two input signals and assuming that u/lo.1 v/fa.j] then (t) (t))
Then, we recall the general form of the neural oscillator, given in the following
i(t) (Ay(t) + Bu(t) + cc,,
( (0) =0, (t) = (t(() +e,
(20)
where y E Rd and d is C (TR activation function with (0) and o'()) 1, e.g., tanh sin. Here, Eq. 220 defines an input-output mapping u(t) + (t) E R as a solution of a second-order
dynamics. Furthermore, Lanthaler et al. (2223) (Section 2) also considers another form of neural
oscillator, which is given by
((()) (Ay(t) (c) y(0) 5(0)=
(21)
(1(1) = Dy(()- +e.
In this case, the input is provided via r(t), an external forcing term which in their setting is linear transformation the inputs. For the universality result presented shortly below, we will consider the neural oscillator in Eq. (2I). which is the closest setting to ours. Here, the universality result remains same, as the original proof does not change if the forcing is inside or outside the activation (see the discussion at the end of Section 2 in Lanthaler et al. (2023)).
The main result of Lanthaler et al. (2023), for Neural Oscillators the form of Equation 21 can be summarized in the following statement:
Theorem (Theorem 3.1 of Lanthaler et al. 2023 ). Let 4 Co(OOTT);P) Co(lo, T): RP) be causal and continuous operator and let K C Co( T); RP) be compact. Then, for any > O, there
exists latent dimension d, weights A, D, and biases such that the output of the multi-layer neural oscillator with input satisfies:
sup (t)) z(t) S e, Yu K. 1E(0,1]
(22)
While this theorem is stated for maps between time-dependent signals utt) and <t)) Lanthaler et al. (2023) provides an additional result for continuous functions F IR + IR9, i.e. neural oscillators can learn to approximate arbitrary precision any such map F between vector spaces. The full
discussion is provided in Appendix A of Lanthaler et al. (202)).
To show that model is universal in either sense, need cast model, defined by Eq. (6) and Eq. 7.. into the definition in Eq. (2T). First, we provide a proof of the universality of our model, cast into the neural oscillator form, and then we provide interpretation of our results.
Theorem (Universality of IGNS, Theorem Of Section D. Let V be the functional that maps the initial condition xo = [go. Pol to the solution at time T to the ODE defined Eq. (6) and Eq. (7). In other words, Ve represents IGNS. Then, for any F i - Ruxd d with compact support, which represents the evolution of the system after time T, there exist set of parameters 9 such that I o F(xu)]l c.
Proof. Our goal for the proof is to reconcile our definition of the port- nort-Hamiltonian system in Eq. (6) with the definition of the oscillatory model Eq. (2I). We start from the derivation of IGNS as
mass-spring-damper system in Eq. (14). Here, we consider the separable Hamiltonian as in Eq. (7) hence V2.4 0. We can additionally set the dissipation term to be null D = 0: as we will see, without this term, our model can be cast more simply into the form of the Neural Oscillator in Equation (2T). Adding this term allows for further flexibility to learn a wider class of functionals, but does not break the proof (D can be inserted into the definition of the auxiliary function
21

---

Published as a conference paper at ICLR 2026
later). Conversely, IGNS can simply learn to have null dissipation term, if needed by the particular task. With this in mind, Eq. (14 can be written as
== - 1=-M- (Volt + q)) + V.
(23)
For this proof, we consider the time-invariant case, while the general case is very similar, as O p.V.H plays the same role as the external forcing r(t, q). By using the definition of the Hamiltonian in Eq. (
we rewrite Eq. (23) for each node as
q. =
Wg
Woq. +
Vgq +b +b) + tr((, + +r(t,9)+ JEA(()
Σ
(24)
where g(q); = EjENt Wo + EKENC) Vq + b) accounts for the terms in which node d is a neighbor of node j (there exists a k (() for which h = i). Stacking all nodes, we arrive at the form
q= (Wao (We9+ + + f(t. 9)):
(25) where 6 = (b, ......... b] contains one copy of for each node and W is a block matrix, where each diagonal block has copy of W and block j has a copy of V iff j E N()). In fact, W absorbs the effects of node aggregation into one bigger block matrix.
Since the definition from Lanthaler et al. 2023) has a null initial condition, we adopt a trick with the forcing terms to reconcile the two forms finally. Let q.P be auxiliary functions, solutions to & p and p = (q, (t, q). which lead to q M-1 (g. q). and initial conditions qo and po-
Then, the function y q q evolves through the following equations
ÿ M 1 y + W) (Way+ + WG
M - 4
1) @) -
(26)
00) 0,
(0) = 0.
Here, we can see that Eq.
26 is indeed in a similar form as Eq. 21 The only difference is in
the additional matrix M
2 H. which depends only on W. (since H is separable in p, Ø)
and makes Eq. (66) more general than Eq. 21. Finally, we define z(t) as either x(t) yt)
or as obtained via a decoder z (t) = (ec(y(t)) Hence, IGNS is a universal neural oscillator by Theorem 3 To conclude our proof, we need two extra observations. First, since p(t) Mt))).
IGNS is universal also in terms of œ = [q. p|__ Second, this result shows that IGNS is an operator input signal rtt, q) into any possible evolution, and it can be further extended to any map between the space of our interest Rand to itself. For this, it is sufficient to consider the
that can convert
proof given in Appendix A Lanthaler al. (2023) and recall the definition of q, which adds the information on the initial condition through r ( q)
O
Remark 2. This proof also shows that IGNS can learn any possible transformation with compact support between the forcing terms and the dynamics of the physical system. Specifically, the Hamil- tonian in Eq. 7 can learn any transformation between continuous functions in C(O,, T), Baxd).
Remark 3. Like most Universality theorems for Neural networks or operators, our result requires the support of the operator to compact subspuce of R We argue that this assumption is not restrictive for our purposes, as we work with physical simulations. In particular, subset K C Raxd is compact if it is closed and bounded. In our case, this means that there exists A >> 0 such that ||gol| < A and ||poll < A. It is reasonable assume that in contained physical simulations, the scale of physical quantities is limited and that there is a limit to the energy inside or added to the system. Hence, these are fair assumptions for our purposes.
There is an interesting interpretation of this result, which links the definition of IGNS in Eq. (6) and the one of neural oscillators. In particular, looking at the formulation in Eq. (66 we see that all the terms that are outside the activation function play the same role as the input signal r(t), defined in Eq. (2I). In our case, these terms contain the geometrical information and the external or driving forces acting the system, well as the initial condition. Hence, IGNS learns to transform these
22

---

Published as a conference paper at ICLR 2026
inputs into the desired dynamics. In this sense, the Hamiltonian learns how propagate through the mesh the information from the initial conditions, geometry, and external forcing.
Another interpretation comes directly from Eq. (24) by looking at the state of cach node qı- In particular, we can interpret the state of neighboring nodes qz. with j E (t),. as the input function (t)). Hence, IGNS learns how each node should behave in the simulation as consequence of the influence of the neighboring nodes, dissipation, and forcing.
E.2 VANISHING GRADIENTS AND LONG-RANGE INTERACTIONS
The radius of propagation exploited by the GNN has a large influence on its performance (Li al, Cai & Wang. 2020 Alon & Yahav 2021 Rusch et al.) 2023 Giovanni al. 2023 Roth & Liebig 2024), as it is often related to over-smoothing, which causes node representations collapse, and over-squashing, which causes information dissipation. Recently, it has been observed that these
2019a Cai
behaviors are strongly related to the sensitivity of the Jacobians 0x.(++),, /dx!!) (or dz. /0x,) ()) of the GNN updates and their norms (Arroyo et al., 2025), as low sensitivity values can cause contractive behavior that induces information loss. Hence, it is necessary to employ models that provide non-contractive updates to avoid these issues. Additionally, real-world simulations often require modeling long-range interactions, as physical information can propagate over long distances. In this discussion, aim to show that all these problems be solved or limited by adopting an oscillatory dynamic at the core of the model. Particularly, the Hamiltonian formalism provides further guarantees vanishing gradients (Galimberti et al. 2027 (2023) and long-range interactions (Heilig et al. (0255, strengthening performance. From a dynamical system perspective, a key aspect of oscillatory systems lies in their purely imaginary eigenvalues:
Theorem 4 (Eigenvalues of Oscillatory System). Consider second-order system of the form Mã + Kx = 0, where M and K are positive definite. Then, the eigenvalues of the system are pure imaginary.
Proof. Since M is positive definite, it invertible, and we can consider the following system in- stead:
# + Kr. = 0.
(27)
We decompose this second-order system into a system of first-order ODEs:
i = v, ú K s.
(28)
Calling y = (ar, v] the system can be represented as
1 0
I 0
K
ÿ K 0 1
0
y.
(29)
That is, the evolution is determined by a matrix which is the product of an antisymmetric matrix and a positive definite one. Hence, as per Lemma 2 of Galimberti et al. (202T). the eigenvalues of the system are all imaginary.
In terms of sensitivity matrices, similar property has been derived for general model dynamics defined by Æ = f(((t) (Heilig et al., 2025 Galimberti et al., 2023 2021). The starting point is understanding how the sensitivity matrix dæ effect of past states on the current one:
(T)/ax t) evolves through time, measuring the
Lemma 1 (Lemma I of Galimberti et al. 2021 ) Let æ be the solution to the first-order ODE = (t)())) Then, calling Φ(T, T t)= (x(T) it holds that
d dt
öf
0 (T T T - t) dx (TT-t)
(T,TT-
(30)
Hence, the sensitivity matrix and, most importantly, the vanishing gradient effect are strictly related to the eigenvalues of O /0x which, for oscillatory models, are purely imaginary. In reality, Heilig al. (2025) shows much stronger property for IGNSg. In particular:
23

---

Published as a conference paper at ICLR 2026
Theorem (Theorem 2 of Section and Theorem 2.3 of Heilig et al. (2025) ), If is the solution
to the ODE i; = JV æ. (t, X) that is, the Hamiltonian dynamic with Eq. (6) without dampening
and external forcing, then (m(())
3 2 for any 885L
This result represents lower bound on the sensitivity matrix norm for any time f. This ensures that gradients never vanish, even over extended periods. We emphasize that not all oscillatory dynamics exhibit the same desirable properties as IGNS. For instance, while GraphCON (Rusch et al. 2022) models a mass- seass-spring--ee -aaiinggddaerraaaee system (see Section C.I) and alleviates some of the propagation issues presented. A Hamiltonian-lie like dynamic (Heilig et al. 2025), as adopted in our IGNS, provides stronger theoretical guarantees on long-range interactions and sensitivity matrices. These advantages also translate into improved empirical performance, as demonstrated in Section [ GraphCON adds node coupling only in its foreing term, while employing diagonal matrices in the left-hand side of Eq. (II.. See Section for additional details the interpretation of these models mass-spring- damper systems.
First-order models and the Heat Equation. Without inductive biases introduced in oscillatory and second-order models, there is no theoretical guarantee of preventing vanishing gradients. On the contrary, there are both practical (Li et al., 2019c) and theoretical (Arroyo et al. (025) insights that show how GNNs are prone to these exponentially decaying gradients. Similarly, first-order differentalleeuationn-nnspp
GNNs are based on the heat equation 0pu-Au Au 0, which leads again to exponentially decaying gradient and information propagation (Gravina et al. 2025), provable also through the lenses of over-smoothing (Eliasof al. 2021
E.3 FURTHER DETAILS ON THE VANISHING GRADIENTS WITH DISSIPATION AND EXTERNAL. FORCING
In this Section, we provide further details how IGNS limits the vanishing gradient effect. This discussion inspired by the results in Rusch et al. (2022), in particular Propositions 3.5 and 3.6. We formalize this in the following statement:
Theorem 5. For =
T, let (5) the prediction of IGNS at time t. Then, for sufficiently small At << 1, the leading order of At the gradient of the multi-step loss term time T, with respect ag
to a parameter at time t. ((() (that is Powies ), is independent of the number of iterations. Thus, although can be small, the gradient will not vanish even with many iterations of the model.
Proof. We follow the steps for the proof of Proposition 3.6 in Rusch et al. 2022). ). First, we use the chain rule to calculate the gradient of the loss with respect to the model parameters. We will consider only one term of the multistep loss, the one for the T-th step, which we call L for simplicity. As will see at the end of the proof, this will not be issue. We interested in the influence of model parameter at step s, which we indicate (s),, on the loss at step T.
OC OC (x) (m) ac (x(++1) (t)(() Ow dx(() (z(() Owl() (T)
II
(31)
Ox(()) (w(
Thus, need to calculate the one-step Jacobian On (++1) (z)( for this analysis, which requires calculating
the Jacobians of the combinations of the p and q terms. The calculations for the port- Nort-Hamiltonian
dynamics of IGNS are very similar to those Heilg et al. (025)) The main differences are that weight matrices are time-dependent (which does not influence the calculations partial derivatives
with respect to (+)) and that our dissipative term does not depend on q. Since we are interested in the leading order terms (as in Rusch et al.] (2022)), ), report the four quantities of interest up
to order At. which we adapt from Heilig et al. (2025) (equations (48) and (49), here we exchange
between numerator and denominator notation). We omit the explicit time-dependence of the matrix
24

---

Published as a conference paper at ICLR 2026
and recall it when needed.
Öp. (++1)
1ij - Ac (1)
θp (1)
o2(K)V? )
A,
+ b,
E
RENiNW,
Op. (++1)
= A1 1,5 (1 ((()) E (H())
3qg (1)
KEN
AA (,(,),,,, W,,,, +
E ob(())3) +
aF(q.1)
(32)
dq,
kEAGNUN
0q. (1+1)
-At IN (2OW3 + E C
0p) (1)
KEN:
+ Ai (,,,VV + ,,,,
E
kENTTAG
0q: (1+1) 0q, (1)
1, + O
where we used the same notation in Heilig et al. 2025 for d'(()) d' (W,p, EKEN, Vpps) and the other similar terms.
Now, repeat the process of Propositions 3.5 and 3.6 from Rusch al. (2022) and calculate
(x(++2) (z(+))
BatFFy Bxll" up to the first-order, to see how the terms in the chain rule multiply. To do so, we need to calculate the four submatrices of this product. Considering that only the first and fourth
terms in the above calculations have leading term of order O(1), we can simplify as
@p(t+2)
@p((+2) @p(++)
@p(t+2) (++1)
0p(2+2)
@p(++)
1 + 0(213)
0p(()
Sp((+T) 0p(t)
tq(t+T) 0p(()
0p(++3)
0p(()
@p(++2
0p((+2)
(+1++
0p(++2 0p
@p(++2)
@p(++1)
+ O(AP)
0q(()
0q(++1)
0q(()
Sp((+1)
0p(0)
dq(t+1)
0q(1)
(33)
(tt2)
(112+
(t
0q(()2)
(p-1-1)
(1
(1(++
+OAFF)
0p(() (1+1) Opl()
Ep((+1) 0pl1)
Sp(c+1)
0p(()
(1+2) (1+2) (q
lq(+22 (p((-1) 0g
0q(1) 0q(++-1) 0q(() + 0p(r+1) dq(²) (q(++1)
0
0q(0)
This is very important, as it shows that up to the first order,
@x((+2) (x+(++1 (x(++2)
(1+1)
I + O(A?)
(34)
(x+1+) (x(()
(x(++
ax(()
which is very similar to the results in Rusch et al. 2022). With similar steps, we arrive at
(T)
(z(++1)
+ O(dt)²
(x(s)
dz((t)
(0x(++1))
8- +
(35)
dar(t))
T-1 I
Σ E
where, in our case, 1,
Bm (t
and has leading order At. Finally, we note from Equa-
tion (32) that (m((t)
Pwiee + that is, the gradient with respect to parameters, has leading order O(1) or
ac O(Ar), and, naturally BalTy is of order (1)) Hence, in Rusch et al. 2022), the leading order of Bw is O(A(2) independent the number of iterations between and T. n
the gradients octr)
25

---

Published as a conference paper at ICLR 2026
(a) Plate Deformation
(h) Sphere Cloth
(c) Imppct Plate
(d) Wave Balls
(e) Kuramoto-SSvashing
(I) Cylinder Flow
Figure 4: Dataset overview. (a) Plate Deformation: flat plate deformed subject to varying num- bers and magnitudes of point forces, simulated with linear elasticity. (b) Sphere Cloth: a cloth mesh impacted by falling sphere, producing elastic deformation with contact dynamics. (c) Impact Plate: an elastic plate under impact with rigid ball. (d) Wave Balls: surface waves on water gen- erated by three balls moving linearly from different initial positions. (e) Kuramoto-Sivashinsk chaotic spatio-temporal evolution of scalar field governed by the KS equation. (f) Cylinder Flow: incompressible fluid flow around cylindrical obstacles with vortex shedding.
F DATASETS DESCRIPTIONS
Table 4: Dataset summary: average number of nodes, average number of edges, total horizon H,
window length T (prediction sub-horizon), and train/val/test split (number of trajectories).
Dataset
Avg. Nodes
Avg. Edges
Horizon H
Window T
Train/Val/Test
Plate Deformation
851
3279
48
48
800/100/1
Sphere Cloth
401
3020
49
49
400/100/100
Impact Plate
2201
13031
51
51
2000/200/200
Wave Balls
1596
6048
50
50
250/100/50
Cylinder Flow
1885
10843
180
30
200/100/10
Kuramoto-Sivashinsky
1600
6240
300
100
250/100/50
Table reports the dataset details, including the average number of nodes/edges, as well horizon length and window size. For the full-rollout Plate Deformation task variant used in the ablation in
Section 6 we predict the whole trajectory with length = 48. In addition, for Cylinder Flow and Kuramoto-Sivassiiskk, we follow the suggestion from Lienen & Günnemann (002) to only split the sub-trajectories into equal-size, non-overlapping windows. Compared to the Cylinder Flow used in the Lienen & Günnemann (002 paper,
the window size of 30 instead of 10 to further
increase the complexity.
Next, report the node-level inputs and outputs used for each task in Table 5 For Lagrangian systems (Plate Deformation, Sphere Cloth, Impact Plate), the state consists of displacement q: and velocity di together with statie properties n, such material parameters, forces magnitude (for the
26

---

Published as a conference paper at ICLR 2026
Table 5: Task Inputs/Outputs used in our experiments. q: displacement, q:: velocity, nı; static properties (e.g., node-type), P:: density/scalar field, oz: velocity field (for Cylinder Flow). For methods training with multi-step loss ([MS)), they also receive initial absolute positions æ. (notated inline).
Tasks
Type
Node Inputs
Node Outputs
Plate Deformation
Lagrangian
Th Th: q:, (q.
q
Sphere Cloth
Lagrangian
[MS] ) [MS]]
qi q.
Impact Plate
Lagrangian
The q: (q. (0)
qt q.
Wave Balls
Eulerian
Thi P:
Pt. P:
Cylinder Flow
Eulerian
ni, Ui, 0s
vi, ti
Surramoto-Sivashine
Eulerian
Yeg, ps
Pi: Ps
Plate Deformation task) and node types. For these tasks, methods trained with the multi-step objec- tive additionally receive the initial absolute position qf. which represents the rest configuration. For Eulerian systems (Wave Balls, Suraamto-Sivashinskyy field Pt, and the outputs include both P: and its temporal derivative يي. In Cylinder Flow, the state instead contains velocity vs with ù; provided to capture the temporal evolution of the velocity field.
the state is represented by a scalar density
In the following, we give the full details for each dataset, with examples given in Fig.
Solid mechanics. This category includes Plate Deformation, Sphere Cloth, and Impact Plate (c.f. Yu et al. (2024)).
Plate Deformation: The task is to predict the final state (or full horizon) of an initially flat plate subjected to external forces of varying positions and magnitudes. These forces are encoded as special nodes connected locally to the sheet within a small radius, and their number ranges from 1 to 19. This setup primarily evaluates the model's ability to capture long-range spatial interactions. Ground-truth simulations are obtained using lincar-elastic dynamics with shell elements in the commercial Abaqus software (Smith) 2009).
Sphere Cloth: ball dropped from random positions and heights onto a cloth mesh of size 19 x 19, and the system is simulated for T = 50 steps. To ensure fair training and compatibility with dissipative models, we further enhance the cloth mesh connectivity by introducing edges between the ball and every fourth cloth node. The dynamics are governed by elasticity, and the initial ball position strongly influences the outcome. Simulations are performed using multi-body physics in NVIDIA Isaac Sim (Mittal et al. (023). We also consider a more challenging variant of this task, used to ablate the importance of the time- varying component, with larger cloth size of 29 x 29 and rollout length T 100.
Note that this dataset has been redesigned from a similar scene introduced in Dahlinger et al.] (2025), which was originally developed for a meta-learning setup, to instead probe mangrraneean dependencies under direct supervised training.
Impact Plate: This is standard benchmark from the HCMT paper (Yu et al. 2024 to test long-range interaction, where the model must predict both stress and displacement propagation from distant nodes. The ground-truth simulation is Ansys.
Fluid dynamics. This category includes Wave Balls, Cylinder Flow, and the Kuramoto- -Svassiiisyy (KS) equation.
Cylinder Flow: This task is based on the standard MeshGraphNets task (Pfaff et al. 2021), but we use only H 180 time steps (vs. 600) and 200 trajectories (vs. 1000), windowed into non-overlapping segments of length 30. After 180 steps, the system converges to a pe- riodic steady state. This setup highlights error accumulation under autoregressive training. explaining the lower MGN performance in our results.
Wave Balls: We let three balls, starting with random positions, move from left to right linearly along a water surface in various grid shapes (cross, L, U, T). The dynamics are governed by a second-order hyperbolic PDE with extermal forcing generated by those balls:
Oxu 2 2
Sextemal (2r, 1),
27

---

Published as a conference paper at ICLR 2026
The simulation is implemented in NVIDIA Warp (Macklin 2022). An extended version of this task, used in the ablation study, runs for T 100 steps. For the first 50 steps the balls move acToss the surface, and for the next 50 steps they stop. The waves, however, keep oscillating since there is no damping. This creates an interesting case where the model must learn not only the forced response but also continue the free oscillation, and figure which phase the system is in.
Kuramoto Siiassiiii (KS): A 4th-order PDE commonly used neural operator bench- marks (Li et al. (021). generating chaotic behavior. Instead of random initial states, we use three Gaussian sources on a 40 x 40 grid and apply Neumann boundary conditions (not periodic) to focus on local rather than global behaviors. Simulated with the py-pde library (wwicker 2020). The KS equation in 2D is given by
Ocu = 2
Vul Vuu V4u.
G HYPERPARAMETERS AND TRAINING TIME
Table 6: Hyperparameters used in all experiments. "AR" is autoregressive methods (e.g., MGN).
"MS" is multi-step methods (others).
Hyperparameter
Common
AR
MS
Node feature dimension
Latent representation dimension
128
Decoder hidden dimension
128
Optimizer
Adam
Learning rate
x 10
possage-passing blocks
15"
Window size none 8
Training noise (std)
x 10-2
Task-specific overrides (MS unless noted):
Plate Deformation
message blocks =
Sphere Cloth
Warmup 20.
Impact Plate
Warmup 15.
Wave Balls
Warmup
Prramoto-SSvvshinne
Warmup = 10; MS training noise = 10 -0
Cylinder Flow
Warmup 10: MS training noise 1 x 10
MGN uses pressage-passingg blocks by default; Plate Deformation uses 50. For MS methods we disable training noise except on Cylinder Flow and Buramoto-Sivashinee where we use 1 x 10-2 due to varying initial conditions.
(ummmriiu the shared and task casc-specific hyperparameters. Most settings are common across models, with feature sizes fixed 128 and Adam the optimizer. For MGN, we adopt mean aggre- gation with Leaky ReLU, use 15 mnessage-passing blocks by default (50 for Plate Deformation), and apply Gaussian training noise with
x 10 In contrast, MS methods use ReLU activations in the encoder and decoder, while the graph-convolution blocks rely on tanh in A-DGN, IGNS and IGNS, Training noise is generally disabled for MS methods, but we enable it with in Cylinder Flow and Kuramoto- Siviottoo--vvvvsk...... where the initial conditions vary. In all cases, we nor- malize inputs and unnormalize outputs for every training batch, which stabilizes optimization and allows us to use the same learning rate and training noise across tasks.
1 10
Table (a) shows the training time allocated to each task. ensure fairness, all methods were trained for the same wall-clock time on single NVIDIA A100 GPU. For the default settings, the budgets were chosen based on validation curves, where extending training further did not yield no- ticeable improvements. For the extended version of Sphere Cloth and Wave Balls tasks, we instead fixed a hard budget to ensure sufficient training for fair comparison in the ablation studies, indepen- dent of validation plateaus. In Table 7 (b), we report the number of parameters for each method. It clear that the standard MGN with a non-shared processor requires by far the largest number of parameters.
28

---

Published as a conference paper at ICLR 2026
Table (a) Training time budget for each lask, where all methods were trained with the same wall-
clock time. (b) Parameter counts (measured on the Kuramoto-Sivashinsky task).
Task
Default (hours)
Long (hours)
Method
#Parameters
Plate Deformation
12
IGNS
216,000
Sphere Cloth
12
24
IGNS;
80,100
Impact Plate
48
GCN+LN
79,700
Wave Balls
16
24 / 48
MP-ODE
101,000
Cylinder Flow
24
GraphCON
67,800
Kuramoto Sivashinsky
24
A-DGN
100,000
MGN (15 MP)
1,800,000
(a) Training time
(b) Parameter counts
IGNS
MGN
Method
Runtime (s)
MGN
0.6690 0.0845
Forward time (s)
GCN_LN
0.0381 + 0.0007
MP-ODE RK4
0.1414
0.0068
MP-ODE dopri5
0.3062
0.0064
4-
ADGN
0.0721 0.0008
GraphCON
0.0433 0.0005
IGNS; IGNS
0.1253 0.0011
0.1703 0.0205
24x24
48x48
96x96
192x192
Resolution
Runtime per forward pass on WaveBalls (7-200, batch size 1). Mean + std over samples (seconds) on different methods (left) and with varying resolutions (right).
Figure 5:
H COMPLEXITY AND RUNTIMES
We discuss the theoretical complexity our method, followed by a comparison of runtimes with other methods.
IGNS preserves the complexity of its underlying backbone GNN. For instance, when instantiated with a GCN backbone, each iteration is linear in the number of nodes VV| and edges [E], achieve a time complexity of [(VV] + [E]) Considering warmup steps and IGNS iterations, the total time computational cost is therefore O(LL + T)(V) El)). For comparison, MGN performs message- sasssagg-pssing updates per simulation step. yielding O(ST(VV)+ ([V]+([+))))) in our setup $=15, which explains why it is the slowest. GraphCON is simpler: it applies single GCN only in the external-forcing branch, so its per-step cost remains ([V[+1E)) with small constant. By contrast, IGNS/IGNS, uses three GCNs for (q.p) and forcing, raising the constant factor but increasing expressiveness, while still staying linear-time and much faster than MGN.
Runtimes. Fig. 5 reports the inference runtimes for one full-trajectory forward pass (in seconds) on Wave Balls with 7-200 (batch size 1), measured on a single NVIDIA RTX 4080 (16GB). On the left, IGNS is competitive with other GNN backbones and substantially faster than adaptive ODE solvers and MGN: the latter is slowest because it performs inner spatial updates at every time step. On the right, we compare IGNS and MGN as we increase the Wave Balls resolution from 24 x24 up to 192 x 192. This corresponds to graphs with roughly 450 nodes at 24x24, about 6k nodes at 96x96, and around 23k nodes and 93k edges at 192x 192. While the runtime of IGNS remains al- most constant up to 96 96x96 and grows only moderately at 192x 192, the runtime of MGN increases sharply and becomes an order of magnitude slower at the highest resolution. This behavior consis- tent with the architectural design: MGN's additional spatial update, which is applied sequentially at every forward step, makes its cost scale with both the time horizon and the graph size, and prevents efficient GPU utilization compared IGNS, which avoids these extra per-step spatial updates.
29

---

Published as a conference paper at ICLR 2026
Ground Truth
MGN-Rewiring
MGN
IGNS
Figure 6: Plate Deformation qualitative comparison on three test samples with predictions from MGN-rewiring, MGN and IGNS. The arrows indicate the forces (with magnitude proportional to its length) applied to the plate at specific positions.
I ADDITIONAL RESULTS
Table 8: Comparison of train and test losses, shown as mean std, on Plate Deformation. Best results highlighted with orange.
Method
Test loss
Train loss
MGN (m=50)
1.27 + 0.06
0.09 + 0.01
MGN-rewiring (m=5)
3.72 + 0.15
0.20 + 0.01
IGNS
1.34 + 0.05
0.71 + 0.05
Plate Deformation. Fig. shows qualitative results three test samples, comparing IGNS, MGN, and a rewired variant of MGN. In MGN-rewiring, we improve connectivity by directly link- ing the force nodes the entire mesh (through every second node), a common technique to mitigate over-squashing. We also reduce the number of message passing steps from 50 (standard MGN) to rewG-newiring). The results highlight clear differences: only IGNS produces smooth and phys- ically consistent surfaces, while MGN predictions are less smooth despite lower eπтоr, and MGN- rewiring generates the roughest outputs. Quantitative results Table & further confirm this trend, showing that both MGN variants severely overfit by relying too heavily on geometric features, while IGNS achieves a better trade-off between train and test performance.
Impact Plate. Table 9 reports quantitative results the Impact Plate task. comparing IGNS and IGNS against the HCMT model from Yu et al. (2024) and MGN. Here, we take the results from Yu et al. (224) for both HCMT and MGN, and report with the metrics RMSE for position
30

---

Published as a conference paper at ICLR 2026
Table Impact Plate quantitative results. Reported are root mean squared errors (RMSE) for posi- tion and stress, shown as mean + std. Best results are highlighted with orange.
Method
Position RMSE
Stress RMSE (x 109)
HCMT
20.71 + 0.57
14.74 + 5.02
MGN
40.73 2.94
35.87 11.89
IGNSa
7.99 + 1.03
3.25 + 0.69
IGNS
7.65 = 0.94
2.85 + 0.55
IGNS/GGNS
geometric encoding no damping
no warmap
no forcing
no forcing, no damping
MSE
0.12
Validation MSE
0.12
0.12
0.10
0.10
0.10
0.08
0.08
0.08
0.00
0.00
0.00
0.04
0.04
0.04
0.02
0.02
0.02
200
800
200
400
800
Epochs
Epochs
Epochs
Figure 7: Validation MSE SphereCloth-diiee (a) IGNS core components: geometric encoding, warmup (multi-step loss omitted as curve explodes beyond the limit). (b) Damping and forcing in IGNS. (c) Damping and forcing in IGNSS
and stress. Both IGNS and IGNSa significantly outperform the baselines in both position and stress prediction, with IGNS achieving the best performance overall. This demonstrates the effectiveness of our approach in capturing long-range interactions and complex dynamics in this challenging task.
Sphere Cloth quantitative results. Reported Table 10:
squared errors (MSE) for overall loss, position, and velocity, shown as mean + std (x10-3 3). Best results are highlighted with orange.
Model
Loss
Position
Velocity
GCN+LN
26.45 + 0.23
4.14 + 0.04
22.31 + 0.20
MGN
32.07 2.45
6.53 + 0.76
25.54 1.73
MP-ODE
25.75 + 1.32
4.26 + 0.26
21.49 + 1.06
ADGN
30.99 + 0.80
4.57 + 0.25
26.42 + 0.66
GraphCON
29.00 0.44
4.19 + 0.07
24.82 0.40
IGNS
28.20 0.35
3.74 + 0.08
24.46 + 0.30
IGNS 24.16 0.67
3.55 + 0.15
20.60 + 0.55
Sphere Cloth. In this task, the model must jointly predict the motion of rigid sphere (represented only by its center position) and deformable cloth, creating strong coupling dependencies. Figure 9 illustrates two representative test samples. In the top sample. MGN produces wrong predictions due to error accumulation, while in the bottom sample it manages produce more reasonable result. This highlights the instability of first order methods, in contrast IGNS, which consistently yields results closest to the ground truth. The quantitative results in Table I0 confirm this trend: although IGNS, achieves slightly lower position error, its surfaces remain rougher, while MGN prone to enror accumulation and overfitting. Other baselines (GraphCON, ADGN, and GCN) fail generate smooth surfaces.
Sphere Cloth direct. In Fig. 7 the validation curves show that even without forcing, the time- varying model (IGNS) can reach a reasonable solution after a longer training phase, which explains the larger variance across seeds. This suggests that time variation partly substitute for external
31

---

Published as a conference paper at ICLR 2026
IGNSa
IGNS mod Ae 10
IGNS f At
IGNS S Ar 10
Normalized Energy
2.5
Normalized Energy
0.05
Ammma
0.00
1.5
-0.05
1.0
0.10
20
30
40
50
60
10
30
Steps
Steps
Figure 8: Energy behavior on SphereClott-direct Left: IGNSa (with damping and forcing) shows an increasing energy profile, as expected for non-conservative model that learns external actua- tion/dissipation. Right: removing damping and forcing yields approximately conservative lution; smaller time steps (At E (10-3, 10 10-2)) keep energy nearly constant, while a larger step (At 10-1) produces mild oscillations due discretization error.
foreing by mediating energy exchange and therefore can help increasing the model's expressiveness. By contrast, the time-invariant model (IGNSa) depends heavily on forcing: when it is removed, optimization plateaus at a poor local optimum for a long training time (until 400 epochs), and then starts decaying with slower consistent with the Test MSE reported in Table 3 Energy conservation analysis. Next, to assess energy behavior, we load the best validation checkpoint of IGNS, on this task. Fig. the full port- Hamiltonian model with all the terms (left) increases normalized energy over time, reflecting learned non- non-conservative effects. When we disable forcing and damping (right) and vary the integrator step, the dynamics become nearly conservative: the energy stays flat for small At and oscillates slightly for larger At. This aligns with our ablations, where forcing acts as the main driver of energy injection while the time-invariant core without con-conservative terms maintains near-constant energy.
compared the time-varying versions. These behaviors
Wave Balls. Fig. I0 compares IGNS with the baselines MGN, GraphCON, and A-DGN. The absolute error maps show that IGNS is the only method producing consistently reasonable results. MGN is able to track the ball positions but overestimates wave amplitude, leading to large errors. near the sources.
GraphCON captures the wave pattern but mislocalizes the balls, causing A-DGN fails to learn this task.
Wave Balls (long) (T=100 and T=200). We generate trajectories with high-resolution simulator using a ground-truth horizon 400 and form two settings by subsampling: T=100 (4x) and T-200 (2x). This probes whether models learn under larger step sizes while remaining accurate over long horizons. As shown in Table II MGN fails severely on this task due to over-reliance on autoregressive state. contrast, both IGNS and IGNS, perform well; IGNS attains the lowest error and captures finer wave-density details than its time-invariant counterpart.
Comparison of test losses on WaveBall at Te(100, 200} (mean + std; lower is better). Values are reported in x 10-3 Best results highlighted with orange. MGN is clearly worse, and IGNS surpasses the time-invariant variant IGNS,, confirming the benefit of the time-varying nort-Hamiltonian core.
Table 11:
Method
T-100
T-200
MGN
1.644 0.096
2.469 = 0.221 2.011 + 0.063
GraphCON
2.476 + 0.082
IGNS#
0,675 + 0.014
0.498 = 0.022
IGNS
0.549 + 0.014
0.402 + 0.037
32

---

Published as a conference paper at ICLR 2026
Ground Truth
MGN
MP-ODE
ADGN
Ground Truth
GraphCON
IGNSt
IGNS
Ground Truth
MGN
MP-ODE
ADGN
Ground Truth
GraphCON
IGNS
IGNS
Figure 9: Sphere Cloth qualitative comparison on two test samples (separated by a horizontal line) with predictions at the final step /. T from first-order models (MGN, MP-ODE) and second-order models (ADGN, GraphCON, IGNS, and IGNS).
Kuramoto-Sivashinsky. To illustrate performance, we show predictions and absolute errors for FNO-RNN, MGN, GraphCON, and IGNS on two samples (Figs. [] and 12 The first sam- ple reflects an early stage where local behavior dominates; the second shows later, more chaotic regime. Two patterns emerge: FNO-RNN produces more globally distributed responses, while the graph-based models remain more localized. In both samples, MGN shows larger absolute errors than the second-order models (GraphCOO and IGNS). Fig. [T] MGN also misses one mode, con- sistent with error accumulation under autoregressive rollout. Overall, GraphCON performs best on this task, with IGNS and IGNS, close competitors.
Warmup Phase. Figure 13 illustrates the effect of different warmup steps / on the Plate Deforma- tion task over a full rollout. With almost no warmup (I 1), the model produces sharp local patterns that gradually recover towards the correct global shape, matching the ground truth only after about
33

---

Published as a conference paper at ICLR 2026
20
20
10
10
10
20
10
30
Initial
Final ground-tree -truth
30
30
30
20
20
20-
20-
10
MGN
GraphCON
A-DGN
IGNS
40
40
40-
40
30
30
30
10
19
0
20
30
30
20
30
MGN сгrоr
GraphCON сттоr
A-DGN eгror
IGNS error
0.00
0.05
0.10
0.15
Absolute Eroor
Figure 10: Wave Balls qualitative comparison on one test sample. Top: initial condition and ground- truth final state at time (-T. Middle: predictions at 1=T from MGN, GraphCON, A-DGN and IGNS. Bottom: absolute error maps, /ppeed - Ugrl- The horizontal colorbars indicate the error mag- nitude in the bottom row.
t = 10. Increasing accelerates this recovery: with = 5, the deformation aligns well already by t = 5, and with / 10 by / 3. For larger values, / 30 and / 50, the global shape is already close to the ground truth at t = 1. These results highlight the role of the warmup phase in seeding long-rang range information before rollout, enabling the model to capture the dynamics more effectively from the very first step.
34

---

Published as a conference paper at ICLR 2026
10
20
Initial
Final ground-truth
30
30
30
30
10
20
30
10
20
30
10
FNO-RNN
MGN
GraphCON
IGNS
30
30
35
20
20
20
20
10
10
10
10
20
10
20
10
20
30
10
FNO-RNN error
MGN error
GraphCON error
IGNS error
0.00
0.05
0.10
0.15
Absolute
Saramoto-SSvvshinsee (KS) qualitative comparison one test sample. Top: initial condition and ground-truth final state at time t=T'. Middle: predictions at t=T from FNO-RNN, MGN, GraphCON, and IGNS. Bottom: absolute error maps, |Mpred Mgal The horizontal colorbars indicate the error magnitude in the bottom row.
Figure 11:
35

---

Published as a conference paper at ICLR 2026
10
20
Initial
Final ground-truth
30
30
30
30
10
20
30
10
20
30
10
FNO-RNN
MGN
GraphCON
IGNS
30
30
30
20
20
20
20
10
10
10
10
20
10
10
20
30
10
FNO-RNN error
MGN error
GraphCON error
IGNS error
0.00
0.20
0.40
0.60
Absolute
Figure 12: Huramoto-Sivashinsk (KS) qualitative comparison one test sample. Top: initial
condition and ground-truth final state at time t=7.. Middle: predictions t=7 from FNO-RNN MGN, GraphCON, and IGNS. Bottom: absolute error maps, |Mpred Ugal The horizontal colorbars indicate the error magnitude in the bottom row.
36

---

Published as a conference paper at ICLR 2026
Figure 13:
Comparison between different warmup steps 1 from top to bottom {1,5, 10, 30, 50) with ground truth in the last row, on Plate Deformation task.
37


Reviews:
[{"review_id": "hPPz1S8Gb7", "reviewer": "Reviewer_Rdux", "summary": "This paper introduces Information-preserving Graph Neural Simulators (IGNS) for learning physical dynamics of  both conservative and non-conservative port-Hamiltonian systems. Key features of the work are: 1) Hamiltonian formalism to eliminate part of the model systematic drifts; 2) a warmup phase for global context initialization; 3) a multi-step training loss for long-horizon stability.\nThe authors also present new benchmarks (Plate Deformation, Sphere Cloth, Wave Balls) to test long-range dependencies.\nIGNS achieves good results across all datasets, showing higher or comparable accuracy and stability than MeshGraphNets and GraphCONs together with some other methods.", "questions": "* What are the computational expensies of your model training and inference compared to competitors and numerical simulator?\n* What are the limits of your model generalization? (out of distribution, new geometries, etc.)\n* Why for Kuramoto-Sivashinsky equation the results for IGNS_ti are so much better than for IGNS?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 1, "contribution": 2, "strengths": "* novelty in introducing a port-Hamiltonian formalism to graph simulators\n* theoretical proofs of universality and non-vanishing gradients\n* data efficiency, importance of warm-up steps and length of confident prediction horizon were investigated\n* code attached", "weaknesses": "* the advantages of Hamiltonian dynamics simulation were not properly studied (like the energy conservation for the conservative systems)\n* the generalizability is under question\n* the file dataset.zip in supplimentary link is corrupt\n* the paper lacks qualitative discussion of the distinction between two versions of the algorithm is not clear (IGNS, IGNS_ti (time-independent))", "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": "qG7Y8ly2gj", "reviewer": "Reviewer_KM4o", "summary": "This paper presents Information-preserving graph neural simulators (IGNS). It is a novel approach aimed to improve long-range propagation of information and to reduce error accumulation in task related with physical modelling. By using port-Hamiltonian dynamics, IGNS preserves information across graphs and captures conservative and non-conservative dynamics. Authors use some key innovations, for example: a warmup phase for initializing global context, geometric encoding for irregular meshes, and a multi-step training objective for stable long-horizon predictions.\n\nAdditionally, the authors provide strong theoretical guarantees for IGNS's universality and gradient preservation. Authors provide comprehensive experiments on six benchmarks, including new tasks. IGNS consistently outperforms many baselines.", "questions": "1. See W2 and W3. It is interesting to see results of inference time and scale analysis.\n\n2. In the message-passing update you use only q_{i} and q_{j} (see eq. 9). Does this exclude the momentum variables p_{i} and p_{j} from exchange between neighbors, or is q there meant to denote the full node state [q,p]? Please clarify.\n\n3. Did you conduct an ablation study to directly compare static geometric encoding of IGNS against a dynamic geometric encoding strategy (e.g., updating edge features per time step) within the IGNS framework (not changing the model)? It is needed to empirically validate its claimed benefit in preventing overfitting. If so, where are these results presented? If not, please conduct this ablation study.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. The idea is novel and interesting. \n2. The paper has solid theoretical foundations (universality, non-vanishing gradients) and robust experimental design against many strong baselines across different physical systems.\n3. The paper is generally well-written, logically structured, with clear problem statements and architecture descriptions. Good use of figures, tables, and appendices.\n4. It offers a substantial advancement in GNSs by solving critical long-standing limitations.", "weaknesses": "1. A direct ablation for authors' static geometric encoding feels like missed. The authors make an argument that model's specific geometric encoding helps avoid \"overfitting,\" which is important for generalizing well. However, authors don't really show us an experiment where they directly compare their own model with and without developed static encoding. Instead, they rely an indirect comparison to other architectures.\n2. No inference time analysis. It it is hard to understand, beyond just training time, how fast do the models predict one or multiple steps during inference? This is crucial for \"accelerate simulations\". \n3. No scale analysis. How does the computational cost scale with increasing graph size (nodes, edges) or rollout horizon?", "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": "reD1qXiKA4", "reviewer": "Reviewer_eSDm", "summary": "The paper introduces a graph neural simulator built on a port-Hamiltonian system. Unlike the existing relevant literature, the core idea of the proposed method is to introduce the symplectic integrator while proposing a port Hamiltonian involving non-conservative energy terms in order to closely align the proposed emulator with the ODE dynamics. The warmup iteration is also proposed in an attempt to enhance the long-range message propagation. The paper also conducts theoretical analysis about the universality and the sensitivity of the model. The proposed framework is evaluated on a range experiment including a couple of new scenarios designed to assess the long-range propagation capability under external forcing.\n\n\nOverall, I think the paper is not ready for the publication in the current form, because 1) the main arguments about the theoretical analysis and proposed architecture are over-claimed and loose, and 2) the experiments miss some relevant baselines and essential ablation study. In particular, Theorem 2 is almost completely identical to a theoretical result in an existing literature. A detailed clarification is necessary to highlight the difference. The detail of the warmup iteration is also missing, which is another factor that causes the difficulty to assess the significance and soundness of the contribution of the paper. The details are given in \"Weaknesses\" and \"Questions\" columns.", "questions": "**Regarding the warmup iteration:**\n- How exactly does the warmup iteration work? I cannot find its detail including the update formula for each iteration.\n- Is the port-Hamiltonian used in the warmup phase shared with the time-evolving forward model? \n\n**Proof of Theorem 1:**\n- Why is it fine to set $D$ to be $0$ without the loss of generality? It is not obvious at first glance.\n- What does B represent at line 989. Is it a parameter involved in $r(t)$?", "limitations": "", "rating": 2, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "- The universality result (although I do not fully understand the proof yet) is novel. \n- Proposal of new tasks, specifically designed to test long-range propagation and oscillatory dynamics under external forcing.", "weaknesses": "**Over-claimed assertions:**\n- A sentence starting at line 290 regarding Theorem 1 is over-claiming. Being with/without compact supports makes a huge difference in the significance.\n- The multi-step objective is a pretty common loss function for training auto-regressive models. The authors need to cite relevant papers or argue that this is a pretty common approach. The use of symplectic integrator in this context is also not novel. \n- Theorem 2 is almost identical to a sensitivity result in [1], and it is unclear if this theorem is different enough from the result of [1] to claim it as novel and/or original.\n\n**Misdirected experiments:**\n- The core idea of the paper is closely related to the idea in [1], but apparently the baselines in the main experiments do not include the model, which makes the experiments setting look unfair and the result unfairly unconvincing.\n- The ablation version of damping and forcing/residual terms is missing. This ablated model is essentially equivalent to the idea in [1] adopted to symplectic integrator with well-adopted multi-step loss, which I believe is valuable to be compared.\n- The other fundamental difference is use of symplectic integrator and warmup iteration, but the ablation experiment also misses these aspects.\n- The paper addresses the long-range propagation problem, but the experiment misses the impact of increasing the resolution of the space, which controls the difficulty of the long-range message propagation. \n\n**Inaccurate description on the assertion and proof of Theorem 1:**\n- Line 287: $\\dot{x}_{0}$ should not be included.\n- Theorem 1 should need a compact support on which $F$ is approximated by $\\Psi_{\\theta}$.\n- Line 855: The second Hamiltonian equation. This should be introduced formally.\n- Line 1028 misses the definition of $\\bar{p}$ and the relation between $p$ and $q$ is unclear, so the derivation of (24) is non-trivial.\n\n**Minor:**\n- Typo: the map $\\Phi$ is duplicated in line 987\n- Diameter is important metric in this context, but it is missing from Table 3\n\n\n[1] Heilig, et.al., \"Port Hamiltonian Architectural Bias for Long-Range Propagation in Deep Graph Networks.\", 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": {}}, {"review_id": "HcyQAmUy8g", "reviewer": "Reviewer_S8R7", "summary": "This paper introduces Information-preserving Graph Neural Simulators (IGNS), a graph-based neural simulator that improves modeling of complex physical systems. IGNS enforces Hamiltonian dynamics to preserve long-range interactions and extends to non-conservative systems. It includes warmup initialization, geometric encoding, and multi-step training to enhance stability. Evaluated on new benchmarks with long-range dependencies and external forces, IGNS outperforms state-of-the-art methods in accuracy and robustness for dynamic systems.", "questions": "1. In L201-203 and in Appendix D, the paper states that $\\gamma_\\theta(t)$ and $\\tau_\\theta(t)$ are time-varying coefficient vectors produced by MLPs with parameters θ. Could the authors clarify what exactly is used as the input t to these MLPs? Is t normalized to a fixed range (e.g., [0, 1]) or directly represented as the raw timestep index? Additionally, if the model is trained on trajectories with 400 steps, can the learned time-dependent MLPs generalize to longer rollouts (e.g., 1000 steps) without retraining, or does the model rely on an absolute temporal scale?\n2. L231-233 states: \"Thanks to the energy conserving core of IGNS, this globally informed latent state is preserved throughout the rollout, rather than being dissipated.\" Could the authors clarify why this happens? Can you provide qualitative or quantitative analysis showing how the latent state is preserved over time？\n3. Regarding the multi-step loss: how many time steps were included in the loss computation? Are the reported results based on single-step MSE or on the rollout of the entire sequence? Please include these experiment details.\n4. Why does WaveBall not require a warmup phase?\n5.  In the supplementary code (`igns.py`), L651: `x = self.one_step(x, edge_index, edge_weight, batch, t=i)` passes `i` as `t`. Here, `i` corresponds to the layer index rather than the time step. Can the authors explain why this is done?", "limitations": "", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The proposed Information-preserving Graph Neural Simulator introduces a principled integration of port-Hamiltonian dynamics into graph-based simulators, marking a significant step beyond existing message-passing and oscillatory GNN frameworks.\n- Theoretical analyses are thorough and provide a clear justification for the model’s ability to capture complex and long-range physical interactions. \n- Experimental evaluation is comprehensive, spanning six datasets and consistently demonstrating the superior accuracy and stability of IGNS compared to strong baselines.\n- The paper is clearly written and well structured: the motivation for each component (port-Hamiltonian core, warmup phase, geometric encoding, and multi-step loss) is clearly articulated, and the accompanying figures effectively convey both the methodology and the empirical findings.", "weaknesses": "- Although the theoretical analysis establishes information preservation and universality, it remains largely qualitative in linking these properties to the observed empirical improvements. A more quantitative or ablation-based verification (e.g., measuring gradient norms or energy conservation over rollouts) would provide stronger evidence for the theoretical claims.\n- The training/testing computational overhead of the port-Hamiltonian formulation and the warmup phase is not explicitly analyzed; reporting runtime or memory costs relative to standard GNSs would clarify the practical trade-offs. \n- Although the benchmarks are diverse, most tasks are synthetic or controlled simulations. It would strengthen the paper’s significance to include or discuss applications in more realistic or large-scale physical systems.\n- The geometric encoding used to map edges to features follows the same formulation as previous works (e.g., MGN) and therefore cannot be considered a novel contribution.\n- Several related approaches are not cited or compared [1–4], which limits the contextual positioning of this work within recent advances in graph-based physical simulation. \n- Introducing \"warmup phase\" in GNNs is not novel. Eagle [2] employs a warmup-like phase in its encoder, using multiple message-passing blocks to aggregate local and global context before rollout, which parallels the proposed initialization strategy.\n- The separation of state variables into coordinates and momenta, as well as the coordinate–momentum supervision in Eq. (10), are established techniques already used in [2–3].\n\n[1] EvoMesh: Adaptive Physical Simulation with Hierarchical Graph Evolutions. ICML 2025\n\n[2] Eagle: Large-Scale Learning of Turbulent Fluid Dynamics with Mesh Transformers. ICLR 2023\n\n[3] Efficient Learning of Mesh-Based Physical Simulation with BSMS-GNN. ICML 2023\n\n[4] Physics meets Topology: Physics-informed topological neural networks for learning rigid body dynamics", "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)