Learning coarse-step dynamics and internal mechanical response with graph networks

# Learning coarse-step dynamics and internal mechanical response with graph networks
Source: [https://arxiv.org/html/2609.30344](https://arxiv.org/html/2609.30344)
Vinay SharmaOlga FinkAffiliation:Intelligent Maintenance and Operations Systems, EPFL, Lausanne, SwitzerlandAffiliation:Corresponding Author

## Abstract

Modern sensing records the motion of physical systems, but often leaves the forces and mechanical response governing that motion unobserved\. Inferring these quantities from discretely sampled trajectories is especially difficult at coarse time scales, when mechanical response evolves between observations and interactions propagate across the system\. Here we introduceNewmark\-β\\beta\-DGN, a graph neural network\-based framework that combines two structures inspired by computational mechanics\. First, a semi\-implicit update inspired by the Newmark\-β\\betamethod uses learned momentum fluxes and matrix\-valued response operators to advance the state over each observed interval\. Second, an operator\-weighted virtual hub provides system\-wide coupling through a sparse set of connections\. The learned quantities thus determine the predicted motion and remain accessible for mechanical analysis\. Across a deformable beam, human motion and protein dynamics,Newmark\-β\\beta\-DGNsupports long\-horizon prediction at time steps for which explicit learned simulators deteriorate\. Without force, moment or constitutive relation supervision, forces inferred from walking kinematics track independently derived hip and knee joint moments, while response operators learned on the beam recover the relative spatial and directional structure of its finite\-element stiffness tangent\.Newmark\-β\\beta\-DGNtherefore links coarse\-step prediction to the inference of mechanical quantities that were never observed during training\.

## 1Introduction

Motion provides a direct window into the evolution of a physical system, but not necessarily into the mechanics that produce it\. Structural behavior is governed by internal stresses and load transfer; biomechanical function by forces and joint moments; and interacting systems by forces transmitted across contacts and interfaces\. Yet these quantities are often difficult or impossible to measure continuously, and different internal mechanical states can produce similar observable motion\. Accurate trajectory prediction alone, therefore, does not establish that a model has recovered the mechanics underlying the observed dynamics\. Bridging this gap requires models that infer unobserved mechanical quantities from kinematic observations while using those same quantities to generate the predicted evolution\. In such a model, forces, stresses, and mechanical operators are not auxiliary explanations attached to a prediction, but explicit components of the learned dynamics itself\.

Inferring these mechanical quantities becomes increasingly challenging as the observations become more widely spaced in time\. At coarse temporal resolution, the observations no longer resolve wave propagation, contact, or vibration over their native timescales\. They instead reveal only the net state change produced by these processes over the observation interval\. Over the same interval, a disturbance can travel through a larger part of the physical system\. Predicting a single coarse transition, therefore, requires both a stable representation of the response accumulated over the interval and a mechanism for propagating that response across the corresponding spatial range\.

Graph\-based models offer a natural starting point for addressing these requirements in discretized mechanical systems\. Material points, particles, or components become nodes, while their interactions become edges, allowing information to be propagated across the system through message passing\. Graph Network\-based Simulators \(GNS\)\[[1](https://arxiv.org/html/2609.30344#bib.bib2)\]and MeshGraphNets \(MGN\)\[[2](https://arxiv.org/html/2609.30344#bib.bib3)\]use this organization to predict trajectories across changing geometries and connectivities\.E⁡\(n\)E\(n\)\-Equivariant Graph Neural Networks \(EGNNs\) make these predictions consistent under changes of reference frame\[[3](https://arxiv.org/html/2609.30344#bib.bib15)\], Equivariant Graph Hierarchy\-based Neural Networks \(EGHNs\) shorten communication paths between distant parts of a system\[[4](https://arxiv.org/html/2609.30344#bib.bib4)\], and the Equivariant Graph Neural Operator \(EGNO\) represents evolution over wider temporal intervals\[[5](https://arxiv.org/html/2609.30344#bib.bib5)\]\. These developments extend the spatial and temporal reach of graph\-based prediction, but the quantities exchanged within these networks generally lack an explicit mechanical interpretation\. They are optimized through training to predict the future state without requiring the learned interactions to correspond to the forces, moments, or stiffness and damping like response operators responsible for that evolution\.

A complementary line of work introduces mechanical meaning through either direct supervision or explicit physical structure specified a priori\. Neural Equivariant Interatomic Potentials \(NequIP\)\[[6](https://arxiv.org/html/2609.30344#bib.bib16)\]and the Message\-Passing Atomic Cluster Expansion \(MACE\)\[[7](https://arxiv.org/html/2609.30344#bib.bib17)\]are directly supervised using interatomic energies and forces from quantum\-mechanical calculations\. This provides a direct mechanical interpretation of the learned internal quantities, but requires those quantities to be available as training labels\.

When direct interaction labels are unavailable, different types of interactions can instead be inferred from observed trajectories\. Neural Relational Inference \(NRI\)\[[8](https://arxiv.org/html/2609.30344#bib.bib18)\]represents interactions as latent categorical edge types, while Collective Relational Inference \(CRI\)\[[9](https://arxiv.org/html/2609.30344#bib.bib20)\]extends this construction to heterogeneous interactions\. Both infer these latent relations through trajectory prediction rather than direct supervision of the interactions themselves\. These relations characterize how components interact but do not themselves correspond to mechanically defined quantities such as forces or moments\. A stronger mechanical structure is introduced by the Physics\-Induced Graph Network for Particle Interaction \(PIG’N’PI\)\[[10](https://arxiv.org/html/2609.30344#bib.bib19)\]by representing the trajectory\-inferred interactions as pairwise Newtonian forces\. In these approaches, however, the inferred interactions either remain latent relational quantities or acquire mechanical meaning through an interaction form specified in advance\.

The mechanical structure can also be provided through the physical equations governing trajectory evolution\. Port\-Hamiltonian neural networks represent conservative and dissipative effects through distinct Hamiltonian and dissipation terms, with the port\-Hamiltonian equations specifying how these terms contribute to the state evolution\[[11](https://arxiv.org/html/2609.30344#bib.bib21)\]\. The Information\-preserving Graph Neural Simulator \(IGNS\) extends this port\-Hamiltonian formulation to graph\-based dynamics, using GNNs to parameterize the Hamiltonian and dissipative terms while retaining their roles within the port\-Hamiltonian state evolution\[[12](https://arxiv.org/html/2609.30344#bib.bib6)\]\. Mechanical information can also be supplied directly from physics\-based simulations or models\. Equi\-Euler GraphNet\[[13](https://arxiv.org/html/2609.30344#bib.bib26)\]learns internal interaction forces using force labels generated by high\-fidelity multiphysics simulations, whereas the Physics\-encoded Time Integrator Graph Network \(PeTIGN\)\[[14](https://arxiv.org/html/2609.30344#bib.bib22)\]uses mass, damping, and stiffness operators obtained from a finite\-element model of the system as known inputs to a Newmark\-based nodal update\. These approaches therefore require mechanical quantities to be available independently of the observed trajectories, either as supervision or as prescribed model inputs\.

Across these approaches, learned internal quantities acquire mechanical meaning from physical information or model structure beyond the observed trajectories themselves\. Direct mechanical labels specify what these quantities represent, whereas assumed interaction forms or governing equations define the mechanical form they take and the role they play in the dynamics\. Each of these requirements can be restrictive in real systems\. Internal forces, moments, and load\-transfer quantities are often difficult to measure directly because they arise at contacts, interfaces, or along structural load paths, while reconstructing them with physics\-based models requires constitutive relations, mechanical parameters, boundary conditions, and system\-specific calibration\. Predefining the mechanical formulation not only constrains the learned dynamics to an assumed interaction or response structure, but also requires uncertain and difficult\-to\-estimate system\-specific quantities, such as mass, damping, and stiffness, may themselves be uncertain or unavailable\. These limitations motivate learning mechanically interpretable internal quantities from observed kinematics without direct supervision or assuming any predefined form for those quantities, using general physical principles to define their mechanical role while using the same inferred quantities to generate the predicted state evolution\.

Momentum conservation provides such a general physical principle\. Dynami\-CAL GraphNet \(DGN\) uses this universal law to assign a mechanical interpretation to interactions inferred from trajectories\[[15](https://arxiv.org/html/2609.30344#bib.bib1)\]\.Each edge represents a pairwise mechanical interaction between two components\. DGN represents this interaction as a momentum exchange, aggregates the resulting momentum fluxes at each node, and uses them directly to advance the state\. At coarse temporal resolution, however, the interaction represented over one observed state\-transition is an effective exchange accumulated over unresolved bending, torsion, contact, or articulated motion, and its resultant force need not act along the line joining the interacting components\. Such a non\-central interaction cannot be described by linear\-momentum exchange alone\. Non\-central forces generate a moment and thereby contribute to angular\-momentum exchange, even when equal and opposite forces balance linear momentum\. DGN therefore adopts a Cosserat\-type representation in which each node carries a rotational degree of freedom and each edge carries an antisymmetric angular\-momentum flux\. This allows the orbital moment generated by the force and the intrinsic rotational contribution to be accounted for jointly in the angular\-momentum balance\. When rotation is not observed, the rotational state provides a latent representation of the angular\-momentum exchange underlying the observed translational dynamics\.

DGN’s interaction representation still leaves three limitations in the finite\-time update, all of which become more pronounced as the observation interval increases\. First, the interaction fluxes are evaluated from the current state and then used to advance it explicitly\. Over a larger interval, the state can change substantially, and therefore the mechanical response can evolve during the transition rather than remain fixed at its initial value\. Second, explicit integration of stiff interactions or increasingly fine spatial discretizations requires progressively smaller time steps, thereby increasing the number of integration steps needed to span a single observation interval\. Third, mechanical influence is communicated only along edges connecting directly interacting components\. As the observation interval grows, a disturbance can propagate across a larger part of the system within a single observed state\-transition, requiring information to traverse an increasing number of physical interaction edges\. This third limitation is illustrated by the clamped rod in Fig\.[1](https://arxiv.org/html/2609.30344#S1.F1)a\. A disturbance introduced at the loaded tip propagates along the rod through successive local interactions\. Once the observation interval is sufficiently long for the disturbance to reach the distant clamp, the state reached within that single transition must reflect the clamp’s influence on the system\-wide mechanical response\. Representing this finite\-time coupling through local interactions alone would require information to traverse the full chain of physical edges within one observed transition\. The update therefore must capture the system\-wide mechanical coupling within single state\-transition\.

In this paper, we introduceNewmark\-β\\beta\-DGNto address these limitations through two coupled components: a semi\-implicit nodal update for stable finite\-time integration and an operator\-weighted virtual hub that approximates system\-wide mechanical coupling within a single observed state transition\. Both components are motivated by the structure of the conventional implicit average\-acceleration Newmark scheme\[[16](https://arxiv.org/html/2609.30344#bib.bib7)\]\.

In an explicit update, mechanical forces are evaluated from the state at the beginning of the step and are then used to advance the system\. An implicit Newmark scheme instead evaluates the mechanical response at the updated state, such that the forces depend on displacement and velocity increments that are themselves unknown\[[16](https://arxiv.org/html/2609.30344#bib.bib7)\]\. In a discretized mechanical system, this makes the update inherently coupled: the response of each component depends not only on its own state increments but also on those of mechanically connected components\. The increments therefore cannot be determined independently\. Their equations must be assembled into a global matrix system and solved jointly\.

This global formulation is what allows the implicit formulation to address both the temporal and spatial limitations of the explicit update\. Evaluating the response at the updated state allows the mechanics to evolve over the step and permits stable integration at step sizes that would be too large for the corresponding explicit scheme\. The coupled equations are assembled into a global system matrix whose sparsity reflects the local connectivity of the underlying mechanical interactions\. Although this matrix is sparse, its inverse is generally dense\. Consequently, a force or constraint acting on one component can affect distant components within the same update\. The global solve thereby converts local mechanical interactions into a system\-wide finite\-time response\. The global solve, however, is computationally expensive\. Replacing it with independent nodal solves preserves efficiency but removes precisely the nonlocal coupling induced by the global solve\.Newmark\-β\\beta\-DGNtherefore separates these two roles: it retains the local finite\-time response through independent semi\-implicit nodal updates and introduces a separate mechanism to approximate the system\-wide coupling omitted by those local solves\.

To recover the system\-wide coupling without solving the full global system, we consider a rank\-one approximation of the dense coupling implied by the implicit solve\. Specifically, we approximate the dense all\-to\-all coupling induced by the implicit solve with a rank\-one operator\. We show that, under this approximation, the long\-range pairwise coupling can be represented through a single shared state\. Each component contributes to this shared state through its learned response operator and responds relative to it, resulting in operator\-weighted equilibrium\. This induces a star topology in which long\-range interactions are mediated through the shared state, while local interactions remain represented by direct physical connections between components\. The hub, therefore, approximates system\-wide mechanical coupling withO⁡\(N\)O\(N\)virtual edges, without constructing or inverting a global matrix\. The resulting augmented graph has an effective diameter of two\. The operator\-weighted hub is defined in Section[4\.5](https://arxiv.org/html/2609.30344#S4.SS5)\(Eq\. \([21](https://arxiv.org/html/2609.30344#S4.E21)\)\); its derivation, which eliminates the internal degrees of freedom of the implicit system and approximates the resulting dense coupling by a rank\-one form, is given in Supplementary Information, Sections 9\.2\.1–9\.2\.2\.

Figure 1:Observation\-learned mechanics within a coarse graph update\.a, Over a coarse observation interval, the loaded tip of a clamped rod must respond to the distant constraint within the same transition\.b, Condensation of an implicit mechanical system produces dense coupling\. Its rank\-one approximation is represented by an operator\-weighted virtual hub, giving a graph of diameter two with𝒪⁡\(N\)\\mathcal\{O\}\(N\)edges\.c, Each message\-passing round is one semi\-implicit Newmark substep: edge quantities are decoded and aggregated, the state is advanced and the resulting state enters the next round\.d, Physical edges exchange non\-central linear\- and angular\-momentum fluxes antisymmetrically, whereas virtual hub edges satisfy collective zero\-sum balance\.e, Positive\-definite learned response operators enter the local translational and rotational Newmark systems\.f, The quantities generating the transition remain accessible as internal\-force, torque andresponse\-operatorreadouts\.We implement this rank\-one\-motivated star construction inNewmark\-β\\beta\-DGNby augmenting the physical graph with a single virtual hub\. The physical nodes represent the system components, and bidirectional physical edges represent their direct interactions, as illustrated for the clamped rod in Fig\.[1](https://arxiv.org/html/2609.30344#S1.F1)a\. The hub is connected to allNNphysical nodes through𝒪⁡\(N\)\\mathcal\{O\}\(N\)bi\-directional virtual edges, giving the augmented graph a diameter of two \(Fig\.[1](https://arxiv.org/html/2609.30344#S1.F1)b\)\. Unlike hierarchical or long\-range graph architectures introduced primarily to shorten communication paths, such as EGHN\[[4](https://arxiv.org/html/2609.30344#bib.bib4)\], this virtual connectivity is motivated by the rank\-one approximation of the dense system\-wide coupling implied by the implicit formulation\. The resulting topology thus has a mechanical interpretation beyond merely facilitating message propagation\.

Newmark\-β\\beta\-DGNdoes not assume that the true underlying system\-wide coupling is exactly rank one\. Instead, it retains the star topology revealed by the rank\-one derivation while replacing the scalar couplings with learned matrix\-valued response operators\. The hub in the star topology is neither a physical body nor an independently integrated mechanical degree of freedom\. Rather, it is a shared state whose position, velocity, and angular velocity are recomputed at every substep as operator\-weighted equilibria of the physical\-node states\. The hub, therefore, provides a tractable approximation to the system\-wide coupling that would otherwise arise through the global implicit solve\.

Newmark\-β\\beta\-DGNadvances each observed state transition throughSSsemi\-implicit sub time\-steps, whereSSis a chosen hyperparameter and each substep corresponds to one complete round of message passing \(Fig\.[1](https://arxiv.org/html/2609.30344#S1.F1)c–e\)\. This recurrent update allows the mechanical response and the corresponding system\-wide coupling to evolve within the observation interval\. At the beginning of each sub\-step, the hub state is recomputed from the current physical\-node states using the aggregated nodal response operators decoded during the preceding round\. For the first substep, these operators are initialized as identity matrices, so the hub position, velocity, and angular velocity reduce to the corresponding arithmetic averages over the physical nodes\.

The resulting physical\-node and hub states form the augmented graph state for the current substep\. Following the DGN construction, an edge\-local reference frame is first constructed for each physical and virtual edge\. Vector quantities associated with the edge and its incident nodes are projected onto this frame, yielding invariant scalars that are passed to the edge encoders\. From the resulting edge representations, the decoders produce linear and angular\-momentum fluxes\. On each bidirectional physical edge, the antisymmetric DGN construction ensures that these momentum fluxes are equal and opposite in the two directions \(Fig\.[1](https://arxiv.org/html/2609.30344#S1.F1)d\)\. For the virtual\-edges, the linear and angular\-momentum fluxes are projected onto the subspace of zero collective sum\. Because the hub carries no mass, inertia, or external load, this projection ensures that it redistributes internal forces and torques without introducing a net external momentum contribution to the system\.Newmark\-β\\beta\-DGNadditionally decodes symmetric positive\-definite matrix\-valued response operators:Ki​jK\_\{ij\}andDi​jD\_\{ij\}for the translational channel, andKi​jrotK^\{\\mathrm\{rot\}\}\_\{ij\}andDi​jrotD^\{\\mathrm\{rot\}\}\_\{ij\}for the rotational channel\.

The decoded edge quantities are then aggregated at each physical node\. The aggregated linear and angular\-momentum fluxes provide the internal force and torque contributions, whereas the edge response operators are summed to obtain the nodal response operators\. Together with the learned inverse mass and inertia from the node decoder, these quantities define the translational and rotational semi\-implicit Newmark updates in Section[4\.6](https://arxiv.org/html/2609.30344#S4.SS6)Eqs\. \([27](https://arxiv.org/html/2609.30344#S4.E27)\) and[28](https://arxiv.org/html/2609.30344#S4.E28), respectively\.

The external force is either provided as an input when observed or decoded from the current state when unobserved\. In Eq\. \([27](https://arxiv.org/html/2609.30344#S4.E27)\), the aggregated internal force is combined with the external force to form the net force drive on the right\-hand side, while in Eq\. \([28](https://arxiv.org/html/2609.30344#S4.E28)\), the aggregated torque provides the corresponding rotational drive\. The nodal response operators enter the3×33\\times 3coefficient matrices on the left\-hand sides of the respective updates, allowing each node to account semi\-implicitly for the local evolution of the mechanical response\. Solving these small nodal systems avoids the global matrix solve of conventional implicit integration\. The operators are then carried forward for the next time substep, where they determine the operator\-weighted hub state\. The nodal update and hub construction are therefore coupled recurrently: the hub mediates the system\-wide response used in the current update, while the newly decoded response operators determine the hub state used in the next\.

Each unclamped physical node then solves separate3×33\\times 3systems for its translational and rotational increments\. Fig\.[1](https://arxiv.org/html/2609.30344#S1.F1)e illustrates this replacement for the translational update\. These local3×33\\times 3systems require neither assembly nor factorization of a global3​N×3​N3N\\times 3Nmatrix and can be solved in parallel\. The learned response operators on each edge are constrained to be symmetric positive definite, ensuring that each nodal solve is well posed for every substep\. Because the independent nodal solves approximate rather than reproduce the fully assembled system, exact conservation after integration is not assumed; the resulting finite\-substep conservation residuals are derived and quantified in Supplementary Information, Section 11\.

The proposed architecture is trained solely from observed kinematics, together with boundary and load inputs when available\. When boundary conditions are unobserved, their effects are learned from the observed motion; when external forces are unobserved, they are decoded from the current graph state\. The model receives no supervision from internal forces, moments, or mechanical operators such as stiffness and damping\. The decoded linear and angular momentum fluxes provide internal force and torque contributions for the translational and rotational updates; the translational update also incorporates the observed or decoded external force\. The learned response operators determine how the translational and rotational states respond within each time substep\. Because these quantities participate directly in the state update, they remain accessible after training as internal forces, moments, and response operators \(Fig\.[1](https://arxiv.org/html/2609.30344#S1.F1)f\)\. We evaluate mechanical quantities derived from the decoded interactions against independent references that were never provided during training\. The position\-response operator indicates where the structure is locally stiff, while the velocity\-response operator indicates where it is locally dissipative\. Both are learned coefficients of the coarse\-step update, interpretable as relative mechanical indicators rather than uniquely identified material or constitutive properties\.

We evaluate the proposed framework on four systems\. Three are trajectory prediction benchmarks: a clamped finite\-element beam, human motion capture and protein dynamics\. The fourth is a benchmark for human walking biomechanics that tests the recovery of internal mechanics\. These four settings test coarse\-step stability, prediction on sparse irregular graphs, generalization across loads, geometries, and spatial resolutions, and recovery of internal forces and moments from observed kinematics\. Across the three trajectory\-prediction benchmarks,Newmark\-β\\beta\-DGNmaintains bounded autoregressive rollouts at coarse temporal resolution\. On the beam, it remains accurate under extrapolation in load, geometry, and mesh resolution\. On human motion and protein dynamics, it maintains bounded autoregressive predictions using the sparse hub\-augmented graph, while requiring fewer edges than the dense distance\-based or multi\-hop interaction graphs used by the compared methods\.

A staged ablation examines the contributions of the hub and semi\-implicit Newmark update\. The mean whole\-body error decreases from3\.00%3\.00\\%for the twelve\-substep hub\-free explicitDGNbaseline to1\.26%1\.26\\%for the four\-substep hub\-augmented explicit update, and further to0\.70%0\.70\\%with the semi\-implicit response and0\.54%0\.54\\%with the full Newmark update\. These results support complementary roles for the two components: the hub communicates mechanical response across the system within an observed transition, while the semi\-implicit update supports integration of stiff dynamics at a coarse temporal resolution\.

Beyond predicting trajectories, we ask whether the quantities that participate in the learned update recover internal mechanics never observed during training\. In the human biomechanics experiment, joint\-moment estimates derived from the inferred forces track the reference hip and knee moments, reaching correlations ofr=0\.94r=0\.94andr=0\.87r=0\.87, respectively\. In the beam experiment, the learned response operators recover the spatial and directional organization of the finite\-element tangent stiffness, while the predicted dynamics reproduce the fundamental frequency and dominant vibration mode of the structure\. These results indicate that the learned interactions and response operators used for prediction also provide accessible representations of the mechanical processes underlying the observed motion\.

These results indicate that the learned interactions and response operators do more than support accurate trajectory prediction: they retain mechanical structure that can be evaluated against independent references, despite receiving no supervision on internal forces, moments, or response operators\. By coupling an operator\-weighted virtual hub to semi\-implicit nodal integration,Newmark\-β\\beta\-DGNadvances dynamics at coarse observation intervals without a global implicit solve\. The forces and response operators are part of that update itself\. This makes it possible to learn from observed motion both how a mechanical system evolves and which internal mechanical responses govern the learned evolution\.

## 2Results

### 2\.1Overview of experiments

We evaluateNewmark\-β\\beta\-DGNon four systems: a clamped finite\-element beam, human motion capture, a solvated protein, and human walking biomechanics\. The four settings test coarse\-step stability, spatial generalization, prediction on sparse irregular graphs, and recovery of internal mechanical quantities from kinematics alone\. For the biomechanics case, we introduce a new graph\-based benchmark curated from the instrumented\-treadmill recordings of van der Zee et al\.\[[17](https://arxiv.org/html/2609.30344#bib.bib10)\]\. Each walking trial is converted into a sequence of graphs in which the motion\-capture markers form the nodes and the anatomical connections form the edges\. The source recordings contain synchronized kinematics, ground\-reaction forces, and inverse\-dynamics joint moments\. Only the graph kinematics are provided toNewmark\-β\\beta\-DGN; the joint moments are excluded from training and used solely to evaluate whether the decoded internal forces recover independently derived mechanical structure\. Graph construction, observed inputs, and data splits for all four systems are specified in Supplementary Information, Sections 4\.1–7\.1\.

Across all four systems,Newmark\-β\\beta\-DGNrepresents the observed components as physical nodes connected by their native interaction graph and augments this graph with a single virtual hub connected bidirectionally with virtual edges to every physical node\. The physical nodes carry position and velocity together with case\-specific attributes, while an edge\-type attribute distinguishes physical from virtual edges\.

We compareNewmark\-β\\beta\-DGNwith six graph\-based simulation approaches spanning explicit mechanics, learned particle and mesh simulators, hierarchical equivariant models, and temporal neural operators\.DGN\[[15](https://arxiv.org/html/2609.30344#bib.bib1)\]is the closest comparison for the proposed integration scheme: it shares the edge\-local frames, antisymmetric momentum fluxes, and rotational channel ofNewmark\-β\\beta\-DGN, but integrates the decoded fluxes explicitly\. GNS\[[1](https://arxiv.org/html/2609.30344#bib.bib2)\]uses an encode, process, and decode graph network in which particles are nodes and Cartesian relative displacements describe their interactions\. MGN\[[2](https://arxiv.org/html/2609.30344#bib.bib3)\]applies the same general architecture to simulation meshes, with messages passed along mesh edges and, where required, between points that are close in physical space\. EGHN\[[4](https://arxiv.org/html/2609.30344#bib.bib4)\]uses equivariant hierarchical pooling, while EGNO\[[5](https://arxiv.org/html/2609.30344#bib.bib5)\]applies an equivariant temporal convolution in the Fourier domain\. EGHNO, introduced in the same work, combines the temporal operator of EGNO with the hierarchical backbone of EGHN\. On the clamped beam, we additionally compare against IGNS\[[12](https://arxiv.org/html/2609.30344#bib.bib6)\], an explicit port\-Hamiltonian graph integrator, to determine whether increasing the number of internal integration steps is sufficient for an explicit method spanning the same coarse observation interval\. Supplementary Information, Section 2 documents the baseline objectives and shared adaptations, while Sections 4\.2–6\.2 give the benchmark\-specific implementations and adaptations required for autoregressive evaluation\.

To make the comparison controlled, all models are trained on identical trajectories, data splits, and use the same optimizer, early\-stopping criterion, and rollout protocol\. We retain their published training objectives except where an adaptation is required for the common kinematic task; Supplementary Information, Section 2 records which objectives are retained or adapted, and Section 3\.3 defines the common autoregressive rollout\. Performance is evaluated through full autoregressive rollouts, with each predicted state fed back as input to the next step\. Errors in both the learned dynamics and the state update, therefore, accumulate over the complete prediction horizon\. We measure error over all body nodes rather than at a single point, preventing locally accurate predictions from masking deformation elsewhere in the system\. EGNO predicts a window of future states, whereas the remaining models predict one state at a time\. Because the published EGHN, EGHNO and EGNO architectures do not return velocity, we add an equivariant velocity head following the construction of each model’s position decoder; Supplementary Table 2 specifies each added head and the representation it reads\.

The clamped beam provides a controlled test of the proposed coarse\-step integration\. Its material parameters and reference time step are known; a temporally converged reference removes any dependence on the reference discretization\. Mesh refinement then systematically increasesωmax\\omega\_\{\\max\}, pushing the system beyond the frequency range encountered during training while progressively tightening the stability limit of explicit integration\. Human motion capture andthe protein transitiontest the same architecture on sparse, irregular graphs \(Sections[2\.3](https://arxiv.org/html/2609.30344#S2.SS3)and[2\.4](https://arxiv.org/html/2609.30344#S2.SS4)\)\. Finally, we test whether the quantities learned for prediction also recover mechanical quantities that are never provided during training\. In human walking, inverse\-dynamics joint moments provide an independent reference for joint moments assembled from the internal forces inferred from kinematics alone\. In the beam, the learned response operators are compared with the finite\-elementtangent and with the modal structure recovered from the rollout\(Section[2\.5](https://arxiv.org/html/2609.30344#S2.SS5)\)\. Conservation and computational cost are examined last \(Section[2\.6](https://arxiv.org/html/2609.30344#S2.SS6)\)\.

### 2\.2Stable beam elastodynamics under extrapolation

Figure 2:Newmark\-β\\beta\-DGNremains bounded over the clamped\-beam rollout\.Rollout position error against step, as a percentage of the beam lengthLL, averaged over the 12 held\-out test configurations\.a, Tip\-position error over the full rollout range \(left\) and zoomed to≤10%\\leq 10\\%\(right\)\.b, Whole\-body error, the root\-mean\-square error over all body nodes, over the full range \(left\) and zoomed to≤3%\\leq 3\\%\(right\)\.Newmark\-β\\beta\-DGN\(red\) stays bounded and non\-monotone, tracking the beam’s oscillation, whereas every baseline grows or diverges; a cross marks the first non\-finite value\. Line colour denotes architecture and line style the number of message\-passing steps or the cluster count\. Curves are means over the 12 test configurations at a single training seed \(seed 42\); no across\-seed dispersion is shown\.![Refer to caption](https://arxiv.org/html/2609.30344v1/figures/main/beam_fea/beamfea_main1_rollout_tipdisp_Tc0.4.png)Figure 3:Beam rollout under load\-duration extrapolation\.a, Reference and predicted configurations at four instants of free vibration for a held\-out case whose load is removed at44\\,unit\-time, against cut\-off times of22to33\\,units in training\.b, Components of the tip displacement over the same trajectory: the transverse response follows the phase and decaying envelope of the reference while a small axial drift accumulates\. Single seed \(n=1n=1, seed 42\); no dispersion is shown\.Figure 4:Whole\-body rollout error on the nine extrapolation cases\.Each case extrapolates outside the training range in load, cut\-off time, cross\-section, length or mesh resolution\. Bars show the whole\-body error averaged over the 95\-step rollout \(00–9\.5​time\-units9\.5\\,\\text\{time\-units\}\), expressed as a percentage of the beam lengthLL; the y\-axis is logarithmic\.Newmark\-β\\beta\-DGNuses four Newmark sub\-steps,DGNand MGN twelve message\-passing steps, and IGNS twenty\-four internal steps, six timesNewmark\-β\\beta\-DGN’s budget\. Single seed \(n=1n=1, seed 42\);Figure 5:Modal structure recovered from the rollout\.a, Natural\-frequency spectrum of the tip transverse displacement for the finite\-element reference \(black\),Newmark\-β\\beta\-DGN\(red\) andDGN\(blue dashed\), with the fundamental frequencies annotated and the true fundamental marked\.b–d, The first three proper\-orthogonal mode shapes, finite element againstNewmark\-β\\beta\-DGN, each labelled with its energy fraction and Modal Assurance Criterion\.Newmark\-β\\beta\-DGNmatches the finite\-element fundamental and dominant mode;DGNpeaks at a spurious lower frequency\. No frequency or mode is supervised\. Single seed \(n=1n=1, seed 42\)\.The beam case study considers a three\-dimensional cantilever clamped at one end and driven by a transverse load at the other\. The load is removed at a cut\-off timeTcT\_\{c\}, after which the beam vibrates freely under Rayleigh damping\.All quantities in this beam benchmark are nondimensional; geometric, load, material and temporal values are therefore reported without physical units\.Reference trajectories are generated with a finite\-element solver\[[18](https://arxiv.org/html/2609.30344#bib.bib8),[19](https://arxiv.org/html/2609.30344#bib.bib9)\]on a tetrahedral mesh atΔ​t=0\.1\\Delta t=0\.1\\,for 100 steps\. The in\-distribution dataset comprises 108 trajectories spanningL∈\{1\.0,1\.5,2\.0\}L\\in\\\{1\.0,1\.5,2\.0\\\},W,D∈\{0\.5,1\.0\}W,D\\in\\\{0\.5,1\.0\\\},F∈\{1\.5,2\.0,2\.5\}F\\in\\\{1\.5,2\.0,2\.5\\\}andTc∈\{2\.0,2\.5,3\.0\}T\_\{c\}\\in\\\{2\.0,2\.5,3\.0\\\}\\,at four elements per unit length, divided into 86 training, 10 validation and 12 test configurations\. Nine additional configurations are held out entirely from these sets to test extrapolation beyond the training domain in operating conditions, geometry, and spatial discretizations\. They include a 20% increase beyond the maximum training load \(F=3\.0F=3\.0versus2\.52\.5\); load removal outside the trained temporal range, occurring 25% earlier and 33% later than the respective training bounds \(Tc=1\.5T\_\{c\}=1\.5and4\.04\.0\\,versus−3\.02\.0\\\!\-\\\!3\.0\\,\); cross\-sections for which one dimension is 20% below the trained minimum while the other is 50% above the trained maximum,\(W,D\)=\(0\.4,1\.5\)\(W,D\)=\(0\.4,1\.5\)and\(1\.5,0\.4\)\(1\.5,0\.4\); beams that are 50% and 75% longer than the longest training beam \(L=3\.0L=3\.0and3\.53\.5versus2\.02\.0\); and a twofold increase in mesh resolution \(res=8\\mathrm\{res\}=8versus44\) atL=1\.5L=1\.5and2\.02\.0\. Supplementary Information, Section 4\.1 specifies the finite\-element model, material parameters, observed inputs and data splits\.

The prescribed step is large relative to the mesh time scale\. With non\-dimensionalE=1000E=1000andρ=1\\rho=1the elastic wave speed isc=E/ρ=31\.6c=\\sqrt\{E/\\rho\}=31\.6, so a wave crosses a unit\-length beam in0\.0320\.032\\,and one outer step spans roughly three such transits\.Newmark\-β\\beta\-DGNdoes not advance this interval in a single step, but divides it into four equal semi\-implicit Newmark substeps ofΔ​ts=Δ​t/4=0\.025\\Delta t\_\{s\}=\\Delta t/4=0\.025\\,\. A complementary spectral estimate is obtained from the largest natural angular frequency of the discretized finite\-element system,ωmax\\omega\_\{\\max\}, computed from the generalized eigenvalue problemK​ϕ=ω2​M​ϕK\\phi=\\omega^\{2\}M\\phi\. For the corresponding explicit update, the critical step size isδ​tcrit=2/ωmax\\delta t\_\{\\mathrm\{crit\}\}=2/\\omega\_\{\\max\}, giving approximately0\.00450\.0045\\,on the training mesh and0\.00230\.0023\\,on the refined mesh\. The outer stepΔ​t=0\.1\\Delta t=0\.1\\,exceeds the explicit stability limitΔ​t​ω≤2\\Delta t\\,\\omega\\leq 2by a factor of2222on the training mesh and4444on the refined mesh; even with fourNewmark\-β\\beta\-DGNsubsteps, the corresponding factors remain5\.55\.5and1111\. Supplementary Information, Section 4\.7 derives this spectral estimate, and Supplementary Table7reports the frequencies, stability multiples and amplification factors\. Both the finite\-element solver andNewmark\-β\\beta\-DGNuse the average\-acceleration Newmark parametersγ=12\\gamma=\\tfrac\{1\}\{2\}andβ=14\\beta=\\tfrac\{1\}\{4\}\.

The finite\-element mesh defines the graph representation: each mesh vertex is represented as a node carrying its current position and velocity, together with a free/clamped boundary indicator and the applied load at that node, while bidirectional edges follow the mesh connectivity\.Newmark\-β\\beta\-DGNaugments this graph with one virtual hub connected to every mesh node through directed virtual edges, identified by an edge attribute of−1\-1to distinguish them from the mesh edges\. Each node additionally carries an angular\-velocity state, initialized to zero and updated latently across the Newmark substeps\. The hub carries no external load, and no material properties or finite\-element matrices are supplied to the model \(Supplementary Table 3 summarizes the beam graph representation\)\.

Each outer time step comprises four Newmark substeps, each implemented as a single message\-passing round with a linearized nodal solve\. We roll the model out autoregressively for 95 time steps \(9\.5​time\-units9\.5\\,\\text\{time\-units\}\), covering the free\-vibration response until it is substantially attenuated by Rayleigh damping\. Over this horizon, we evaluate both the loaded\-tip trajectory and the deformation of the full mesh\. Because a metric evaluated at a single point can remain small while the interior of the mesh over\-deforms, we report the whole\-body error, defined as the root\-mean\-square position error over all body nodes as a percentage of the beam length\. Supplementary Information, Section 4\.3 lists the model hyperparameters\. All beam results use one trained seed \(n=1n=1, seed 42\), so no dispersion is reported\.

On the test cases, the error ofNewmark\-β\\beta\-DGNrelative to the finite element trajectory remains below approximately1%1\\%of the beam length throughout the 95\-step autoregressive rollout\. Figure[2](https://arxiv.org/html/2609.30344#S2.F2)shows the evolution of the tip and whole\-body errors, averaged over the test cases: forNewmark\-β\\beta\-DGN, shown in red, both errors are non\-monotone, rising in the first half and falling afterwards, indicating that it tracks the beam’s oscillation rather than accumulating the error in position\.

The same rollout behavior extends beyond the training domain\. We define the whole\-body error as the root\-mean\-square position error over all body nodes, expressed as a percentage of the beam length\. Across the nine held\-out configurations, its mean value over the 95\-step rollout ranges from0\.21%0\.21\\%to3\.45%3\.45\\%\(red bars in Fig\.[4](https://arxiv.org/html/2609.30344#S2.F4)\)\. A representative load\-duration extrapolation case is shown in detail in Fig\.[3](https://arxiv.org/html/2609.30344#S2.F3)\. Panelashows the finite\-element and predicted beam configurations at four instants of the rollout, colored by displacement magnitude, and panelbshows the correspondingxx,yyandzzcomponents of the tip displacement\. When the load is removed at44\\,time\-units, one unit beyond the latest training cut\-off time, the predicted tip displacement continues to follow the finite\-element free\-vibration response, with only a small axial drift developing\. Supplementary Figure 2 shows the same rollout visualization for four additional held\-out cases: a refined mesh atL=1\.5L=1\.5andres=8\\mathrm\{res\}=8, a longer beam withL=3\.0L=3\.0, and the two cross\-section extrapolations\(W,D\)=\(1\.5,0\.4\)\(W,D\)=\(1\.5,0\.4\)and\(0\.4,1\.5\)\(0\.4,1\.5\)\. Across all four cases, the predicted response remains close to the finite\-element trajectory, although larger but still bounded deviations appear for the longer beam\.

We next compare theNewmark\-β\\beta\-DGNwith the baseline models introduced in Section[2\.1](https://arxiv.org/html/2609.30344#S2.SS1)\. All receive the same node features and mesh graph but without the virtual hub\. Supplementary Information, Section 4\.2 details their inputs and adaptations\. Figure[2](https://arxiv.org/html/2609.30344#S2.F2)shows that the baselines exhibit distinct forms of autoregressive rollout error growth across the same test cases\. EGHN and EGHNO reach non\-finite values at steps 11 and 19, respectively\. Their behavior shows that shortening spatial communication paths through hierarchical processing, even when combined with a temporal operator in EGHNO, is not sufficient by itself to prevent rollout error amplification\. The error of EGNO remains finite, but its tip and whole\-body errors grow steadily, showing that a temporal operator over the prediction interval likewise does not prevent error accumulation during autoregressive rollout\. Increasing the message\-passing depth of MGN from 4 to 12 rounds slows this growth, yet the 12\-round model still reaches approximately70%70\\%tip error and61%61\\%whole\-body error by step 95, indicating that deeper local propagation improves the rollout without stabilizing it\. ForDGN, increasing the number of explicit substeps from 8 to 12 reduces the final tip error from about5\.6%5\.6\\%to2\.2%2\.2\\%, and the whole\-body error from about3\.6%3\.6\\%to1\.8%1\.8\\%\. IGNS shows the same general dependence on temporal resolution: increasing its internal\-step budget from 4 to 24 symplectic\-Euler steps substantially reduces both errors, yet the 24\-step variant still reaches about4\.3%4\.3\\%tip error and2\.6%2\.6\\%whole\-body error by step 95\. Thus, finer temporal subdivision improves bothDGNand IGNS, but their errors continue to accumulate over the rollout rather than exhibiting the low, non\-monotone behavior of the proposedNewmark\-β\\beta\-DGN\.

The extrapolation comparison in Fig\.[4](https://arxiv.org/html/2609.30344#S2.F4)shows the same performance gap betweenNewmark\-β\\beta\-DGNand the baselines\.DGNis the closest competitor, with a lower mean whole\-body error thanNewmark\-β\\beta\-DGNonly for the two longest beams,L=3\.0L=3\.0andL=3\.5L=3\.5; IGNS is less accurate throughout and MGN has the largest errors\. Mesh refinement is the sharpest test of stiffness: it lowers the critical explicit time stepΔ​tcrit\\Delta t\_\{\\mathrm\{crit\}\}the most and is the only extrapolation that makes the mesh stiffer than any training case\. On the refined mesh case,Newmark\-β\\beta\-DGNstays at0\.5%0\.5\\%–1\.2%1\.2\\%whole\-body error, against15%15\\%–16%16\\%forDGN,54%54\\%–69%69\\%for MGN and126%126\\%–149%149\\%for IGNS, even at 24 internal steps, six timesNewmark\-β\\beta\-DGN’s budget\. Adding internal steps for IGNS does not close this gap: a larger IGNS budget lowers its single\-step validation loss without stabilizing the rollout \(Supplementary Information, Section 4\.7 and Supplementary Table 8\)\. The reason is spectral: at the four\-substep budget, each explicit symplectic\-Euler substep multiplies any error in the stiffest vibration mode by109109on the training mesh and by455455on the refined mesh, so the error grows within a few substeps\. In contrast, for the corresponding classical undamped linear mode, the average\-acceleration Newmark update has unit amplification, so its amplitude is neither amplified nor numerically damped, although phase error remains; Supplementary Information, Section 10\.4 derives this unit spectral radius and the associated phase error\.

On the two longest beams in the extrapolation regime,L=3\.0L=3\.0andL=3\.5L=3\.5,DGNwith twelve message\-passing rounds attains a slight advantage compared toNewmark\-β\\beta\-DGN, for example, the whole body error of1\.9%1\.9\\%versus2\.1%2\.1\\%atL=3\.0L=3\.0\.Newmark\-β\\beta\-DGNcouples distant nodes through a single operator\-weighted hub, a rank\-one approximation of the global coupling; for longer beams, whose bending response spans the full length, this coupling may require a higher\-rank representation, whereasDGNreaches the same range through progressive local message passing\. Increasing beam length is therefore the regime in which the advantage of the hub\-based approximation gets diminished\.

To identify which components ofNewmark\-β\\beta\-DGNsustain its rollout accuracy across both the test and extrapolation cases, we perform a staged ablation, reporting the whole\-body error averaged over all 21 test and extrapolation configurations\. Starting fromDGN, which is hub\-free and uses twelve explicit substeps, the mean whole\-body error is3\.00%3\.00\\%\. Adding the operator\-weighted virtual hub with an explicit per\-node update \(𝐀i=𝐈\\mathbf\{A\}\_\{i\}=\\mathbf\{I\}\) and only four substeps reduces the error to1\.26%1\.26\\%\. Replacing this update with the semi\-implicit Newmark solve lowers it further to0\.70%0\.70\\%withβ=0\\beta=0and to0\.54%0\.54\\%withβ=14\\beta=\\tfrac\{1\}\{4\}for the fullNewmark\-β\\beta\-DGN\. Supplementary Information, Section 4\.5 and Supplementary Table 5 report this four\-stage ablation\. The operator\-weighted hub provides the largest accuracy gain, with further improvements from the semi\-implicit solve and stiffness term\.

The ablation establishes how the hub and semi\-implicit update reduce whole\-body rollout error\. We next ask whether the resulting trajectory also preserves the beam’s modal structure, despite receiving no frequency or mode supervision \(Fig\.[5](https://arxiv.org/html/2609.30344#S2.F5)\)\. From the tip\-displacement spectrum of the9595\-step rolloutNewmark\-β\\beta\-DGNreproduces the finite\-element fundamental frequency, and a proper\-orthogonal decomposition of the transverse motion recovers the dominant mode shape at a Modal Assurance Criterion of1\.001\.00\. The explicitDGNrollout instead peaks at a spurious lower frequency, so a bounded trajectory needs not carry the correct spectrum\. This is the empirical counterpart of the unit amplification of the average\-acceleration update: the period elongation of a direct integrator is a property of the scheme rather than of its one\-step accuracy\[[20](https://arxiv.org/html/2609.30344#bib.bib24),[21](https://arxiv.org/html/2609.30344#bib.bib25),[22](https://arxiv.org/html/2609.30344#bib.bib23)\]\. Supplementary Information, Section 4\.8 reports the frequency and mode recovery for every configuration \(Supplementary Table 9 and Supplementary Figure 3\)\.

### 2\.3Human motion prediction benchmark

Figure 6:Motion\-capture rollout accuracy and rotation robustness\.a, Position rollout MSE on the unrotated test trajectories\. Each group represents one model, and its five bars correspond to successive rollout steps of 30 frames each\.DGN\(global reference frame\) usesNewmark\-β\\beta\-DGN’s global reference frame inDGN’s external\-force channel; both are rotation\-equivariant\.b, Position rollout MSE after fixed rotations of the test trajectories about the vertical axis\. The singleNewmark\-β\\beta\-DGNgroup applies at every angle because its prediction is rotation\-equivariant\. Both panels are capped at55;trianglesmark values beyond this limit, with the true value printed above\. A step is marked as non\-finite and omitted if any seed diverges\. Bars show means over seeds\{0,42,100\}\\\{0,42,100\\\}and error bars show one standard deviation \(n=3n=3\)\.The beam tests coarse\-step prediction on a discretized structure with a defined physical interaction graph\. Articulated human motion presents a different setting: observations are sparse markers linked by a skeletal graph, and the model must predict coordinated movement across the body over an interval containing many recorded frames\. We use this setting to test whether the virtual hub supports long\-range coordination without a dense interaction graph\. We evaluateNewmark\-β\\beta\-DGNon the walking sequences of subject 35 from the CMU Motion Capture Database\[[23](https://arxiv.org/html/2609.30344#bib.bib11)\]and follow the trial\-level split and prediction interval of EGHN\[[4](https://arxiv.org/html/2609.30344#bib.bib4)\]\. One rollout step predicts the marker configuration 30 frames ahead, and every model is trained from seeds\{0,42,100\}\\\{0,42,100\\\}\. The body is represented by 31 marker nodes and 60 directed skeletal edges\. Each node carries its current position and a velocity computed as the backward difference of positions\.Newmark\-β\\beta\-DGNadds a single virtual hub connected to all markers through 62 directed edges; these virtual edges are separated from skeletal edges by an edge attribute\. While, the baselines add 70 two\-hop edges in addition to the skeletal edges and separate them by an edge attribute\. Supplementary Information, Section 5\.1 defines the observed graph and data split, while Section 5\.3 lists the training configurations\.

Because the external forces acting on the walking body are unobserved,Newmark\-β\\beta\-DGNinfers them from the current state\. To express their directions without introducing a fixed Cartesian reference, it constructs a parameter\-free spatial frame from the current marker positions and the two leading non\-trivial eigenvectors of the body\-graph Laplacian\. The resulting frame rotates with the body, is invariant to translation, while the eigenvectors defining its construction depend only on the graph topology\. We compare our model with GNS\[[1](https://arxiv.org/html/2609.30344#bib.bib2)\], EGHN\[[4](https://arxiv.org/html/2609.30344#bib.bib4)\], EGHNO\[[5](https://arxiv.org/html/2609.30344#bib.bib5)\], EGNO\[[5](https://arxiv.org/html/2609.30344#bib.bib5)\],DGN\[[15](https://arxiv.org/html/2609.30344#bib.bib1)\], andDGN\(global reference frame\), aDGNvariant that uses the same graph\-level equivariant reference frame asNewmark\-β\\beta\-DGNto represent the external force channel\. Because the underlyingDGNinteraction model remains unchanged, theDGN\(global reference frame\) variant isolates the effect of expressing the external force in the graph\-level equivariant frame used byNewmark\-β\\beta\-DGNinstead ofDGN’s original reference frame\. The original published configurations of EGHN, EGHNO and EGNO receive absolute height, whileDGNconstructs it internally; neither GNS norNewmark\-β\\beta\-DGNreceives an absolute coordinate\. The rotations used in the following experiment preserve height, so this difference in available information remains constant\.

On unrotated trajectories,Newmark\-β\\beta\-DGNachieves the lowest step\-five MSE among the equivariant models, at1\.3351\.335\(Fig\.[6](https://arxiv.org/html/2609.30344#S2.F6)a\)\. GNS reaches a lower error of0\.50060\.5006on these trajectories, but is not rotation\-equivariant: although translation\-invariant through its relative inputs, its encoder operates on Cartesian vector components and its decoder predicts Cartesian vectors directly\. We therefore examine whether this advantage persists when the observation frame is rotated\. A rotation of just5∘5^\{\\circ\}, which is small enough to arise from a slight misalignment of the motion\-capture setup, raises the step\-five error of GNS to2\.50972\.5097, above the unrotated error ofNewmark\-β\\beta\-DGN’s1\.3351\.335\. At15∘15^\{\\circ\}and30∘30^\{\\circ\}, the GNS error is35\.635\.6and128128times its unrotated value, respectively\. Figure[6](https://arxiv.org/html/2609.30344#S2.F6)b shows the rotation sweep against the angle\-invariant error ofNewmark\-β\\beta\-DGN\. Supplementary Table 13 reports the full rotation sweep for GNS,Newmark\-β\\beta\-DGNis rotation\-equivariant so its error is unchanged at every angle, while Supplementary Information, Section 5\.4 provides the unclipped rollout results and rotation analysis\.

### 2\.4Protein dynamics prediction benchmark

We next evaluateNewmark\-β\\beta\-DGNon a large collective conformational change of the protein adenylate kinase, from its closed to its open conformation\. We use the trajectory dataset of Seyler and Beckstein\[[24](https://arxiv.org/html/2609.30344#bib.bib12)\], accessed through MDAnalysis\[[25](https://arxiv.org/html/2609.30344#bib.bib13),[26](https://arxiv.org/html/2609.30344#bib.bib14)\], comprising 200 trajectories of the protein moving from the closed to the open state\. Each trajectory contains the 855 backbone atoms over 90–106 frames, during which the two mobile domains move away from the rigid core\. At each prediction step, the model forecasts the protein configuration 15 frames into the future\. Repeating this autoregressively for five steps yields a total prediction horizon of 75 frames, covering most of the closed\-to\-open transition\. We split the data by trajectory into 140 training, 30 validation and 30 test trajectories\. Each node represents a backbone atom and carries its position, backward\-difference velocity, and atom type\.Newmark\-β\\beta\-DGNrepresents local interactions using the covalent\-bond graph and adds a single virtual hub connected to all atoms for system\-wide communication\.DGNuses the same covalent\-bond graph but without the hub\. In contrast, EGHN and EGNO augment the covalent bonds with their published10​Å10\\,\\text\{\\AA\}distance\-based contact graph, producing a substantially denser interaction graph\. Supplementary Information, Section 6\.1 specifies the dataset and data split, while Section 6\.3 lists the training configurations\.

All four models use the same trajectory\-level data split and are trained with a single seed \(n=1n=1, seed4242\), retaining the checkpoint with the lowest validation loss\. Their training objectives follow their respective formulations:Newmark\-β\\beta\-DGNpredicts position and velocity increments,DGNuses the same increment objective on the hub\-free graph, EGHN predicts the end\-state position and velocity together with its auxiliary link\-prediction objective, and EGNO retains its multi\-frame prediction objective\. Despite these differences in training objective, all models are evaluated identically using unscaled position MSE over the same autoregressive rollout\. Because each result comes from a single trained seed, no across\-seed dispersion is reported\.

Figure 7:Autoregressive rollout on the adenylate\-kinase closed–open transition\.Position mean squared error over five autoregressive steps of1515frames each, forNewmark\-β\\beta\-DGNand the three baselines from their minimum\-validation checkpoints\.Newmark\-β\\beta\-DGNstays bounded across the horizon, whileDGNdiverges on the hub\-free covalent graph and EGHN and EGNO operating on dense contact\-distance based graphs diverge: the ordinate is capped at1\.21\.2\\,Å2, triangles mark values above the cap with the true value printed, and crosses mark a non\-finite \(diverged\) step\. One seed \(n=1n=1, seed4242\)\. The full unclipped table and the mean\-trajectory floor \(the average\-displacement baseline\) are given in the Supplementary Information \(Section 6\.4\)\.Newmark\-β\\beta\-DGNremains bounded across the full horizon, rising from0\.2210\.221to a maximum of0\.3460\.346and settling near0\.280\.28\\,Å2\(Fig\.[7](https://arxiv.org/html/2609.30344#S2.F7)\)\. TheDGNremains finite at one step \(0\.3860\.386\), but its rollout climbs monotonically to6\.126\.12; EGHN and EGNO trackNewmark\-β\\beta\-DGNto the third step and then diverge, reaching361361and11\.011\.0at step four before becoming non\-finite\. The hub supplies the global coupling through17101710virtual edges, against the roughly55 61055\\,610directed edges the contact\-distance based baselines carry: providing better stability at a fraction of the graph density\.

Figure[8](https://arxiv.org/html/2609.30344#S2.F8)renders the predicted structural transition for a held\-out trajectory\. For visualization, we extend the rollout to six autoregressive steps, reaching frame\+90\+90\.Newmark\-β\\beta\-DGNreproduces the collective opening in which the two mobile domains, the NMP arm and the LID, swing away from the core; the predicted backbone tracks the ground\-truth conformation with a root\-mean\-square Cα\\alphaposition error of1\.071\.07\\,Å at the midpoint of the transition and0\.780\.78\\,Å at the open state\.

![Refer to caption](https://arxiv.org/html/2609.30344v1/dims_cartoon_rollout.png)Figure 8:Structural rollout of the closed–open transition\.Newmark\-β\\beta\-DGN’s predicted backbone, coloured by domain \(rigid core grey, NMP arm orange, LID blue\), overlaid on the ground\-truth conformation \(light grey\) at the closed start \(frame\+0\+0\), the transition midpoint \(\+45\+45\) and the open end \(\+90\+90\) of a held\-out trajectory, superposed on the core\. The two mobile domains swing away from the core through a collective motion of about66\\,Å, tracked with a root\-mean\-square Cα\\alphaposition error of1\.071\.07\\,Å at the midpoint and0\.780\.78\\,Å at the open state\.To put this prediction accuracy in perspective, we find that the margin over a trivial predictor is small\. A state\-independent baseline that applies the training\-averaged displacement at each step reaches0\.230\.23to0\.370\.37\\,Å2across the horizon, andNewmark\-β\\beta\-DGNstays within44to10%10\\%of it, so the benchmark separates the models by bounded rollout and graph sparsity rather than by a large accuracy margin \(Supplementary Table 16\)\. We use the closed\-to\-open transition trajectories rather than the established equilibrium benchmark because, at its 15\-frame horizon, the equilibrium motion is largely decorrelated from its preceding trajectory and is dominated by fluctuations around the mean conformation\. A simple linear restoring rule therefore outperforms all learned models, indicating that the benchmark mainly measures equilibrium fluctuation statistics rather than predictive collective dynamics \(Supplementary Information, Section 6\.5\)\.

Together, these two benchmarks, in addition to the finite\-element beam case, demonstrate the applicability ofNewmark\-β\\beta\-DGNacross systems with markedly different structures and dynamics\.

### 2\.5Internal mechanics inferred from observations

A model can reproduce observed motion while assigning mechanically implausible forces to the interactions that generate it\. We therefore test the decoded quantities against mechanical references withheld from training\. Human walking provides a stringent setting: marker trajectories describe the motion, while synchronized ground\-reaction measurements and inverse\-dynamics estimates of joint moments provide separate evidence about the forces and moments involved\. To make this comparison, we introduce a graph\-based benchmark using the instrumented\-treadmill recordings of van der Zee et al\.\[[17](https://arxiv.org/html/2609.30344#bib.bib10)\]\. For subject p2, we convert 33 walking trials, containing 18,631 marker frames, into graph sequences sampled at120120\\,Hz\. Each graph contains 37 marker nodes and 63 anatomical edges, augmented by a single virtual hub connected to all markers through bi\-directional virtual edges\. The measured joint moments and ground\-reaction forces are retained as independent references and are not supplied toNewmark\-β\\beta\-DGNduring training\. Supplementary Information, Section 7\.1 describes the graph construction, gait\-condition split, deterministic evaluation windows and moment readout; Supplementary Table 17 lists the observed inputs\.

The training and validation sets contain preferred walking over a range of speeds and cadence variants at1\.25​m​s−11\.25\\,\\mathrm\{m\\,s^\{\-1\}\}\. We withhold every constant\-step\-length and constant\-step\-frequency trial, together with preferred walking at2\.0​m​s−12\.0\\,\\mathrm\{m\\,s^\{\-1\}\}\. The first two groups change the gait strategy, whereas the last changes the speed within a familiar strategy\. Each trial is evaluated using deterministic, non\-overlapping windows that cover the recorded gait cycle rather than randomly selected starting frames\.Newmark\-β\\beta\-DGNreceives only marker positions and velocities obtained by finite differencing\. Its loss contains the normalised errors in position and velocity increments and no force, moment or ground\-reaction term\.

The joint\-moment readout uses the internal segment forces decoded at each step of a five\-step autoregressive rollout\. These are the same forces thatNewmark\-β\\beta\-DGNuses to advance the markers\. For a jointJJ, we sum the orbital moment of the forces acting on the markers distal to that joint \(the markers on the limb segments below the joint\),

𝐌J\(t\)=−∑i∈𝒟⁡\(J\)\(𝐱i\(t\)−𝐱J\(t\)\)×𝐅i\(t\)\.\\mathbf\{M\}\_\{J\}\(t\)=\-\\sum\_\{i\\in\\mathcal\{D\}\(J\)\}\\left\(\\mathbf\{x\}\_\{i\}\(t\)\-\\mathbf\{x\}\_\{J\}\(t\)\\right\)\\times\\mathbf\{F\}\_\{i\}\(t\)\.\(1\)This produces a joint\-moment vector throughout the rollout\. We retain its flexion\-extension component, which corresponds to sagittal\-plane motion, and compute its root\-mean\-square value over the gait cycle\. Applying the same reduction to the inverse\-dynamics measurement gives one joint\-demand value for each joint, leg and walking condition\. We exclude the decoded spin contribution because the nodal spin is an unsupervised angular\-momentum ledger rather than an identifiable joint torque\. The readout is therefore obtained directly from the forces used by the update, without a separately trained inference head\.

Figure 9:Joint\-moment demand inferred from marker kinematics alone\.a, Representative autoregressive rollout of the marker graph for subject p2 for 8 steps, showing the recorded and predicted configurations\. The joint\-moment readout is evaluated over the five\-step rollout\.b, Normalised root\-mean\-square sagittal knee moment inferred from the decoded internal forces using Eq\. \([1](https://arxiv.org/html/2609.30344#S2.E1)\), compared with the inverse\-dynamics reference\. Constant step length and constant step frequency are withheld gait strategies\. Preferred walking at2\.0​m​s−12\.0\\,\\mathrm\{m\\,s^\{\-1\}\}is a withheld speed within a familiar strategy\. The shaded step\-frequency sweep at1\.25​m​s−11\.25\\,\\mathrm\{m\\,s^\{\-1\}\}contains training and validation conditions\. Demand is normalised to the preferred1\.25​m​s−11\.25\\,\\mathrm\{m\\,s^\{\-1\}\}trial at unit subject mass and is therefore a ratio rather than a calibrated torque\. Error bars show variation across deterministic evaluation windows for one trained model and one subject \(n=1n=1\), not variation across training seeds\. Neither joint moments nor ground\-reaction forces enter the training objective\.Across five evaluation experiments and for both legs of the subject, the inferred demand profiles correlate with the inverse\-dynamics measurements atr=0\.943r=0\.943for the hip andr=0\.869r=0\.869for the knee\. Figure[9](https://arxiv.org/html/2609.30344#S2.F9)shows the knee result across familiar and withheld gait conditions\. Under the withheld constant\-step\-length strategy, the inferred demand rises with walking speed\. Under the withheld constant\-step\-frequency strategy, it follows the measured rise to a maximum at1\.8​m​s−11\.8\\,\\mathrm\{m\\,s^\{\-1\}\}\. The shaded step\-frequency sweep at1\.25​m​s−11\.25\\,\\mathrm\{m\\,s^\{\-1\}\}contains training and validation conditions and provides an interpolation control\. There, the readout reproduces the measured minimum near the preferred cadence\. Because no moment enters the objective, this non\-monotonic response supports the mechanical interpretation of the readout, although it is not an extrapolation result\.

At1\.8​m​s−11\.8\\,\\mathrm\{m\\,s^\{\-1\}\}under constant step frequency, the inferred right\-knee demand is1\.2171\.217, compared with1\.5321\.532in the reference, an under\-prediction of21%21\\%\. At the withheld2\.0​m​s−12\.0\\,\\mathrm\{m\\,s^\{\-1\}\}preferred\-speed condition, the corresponding difference is5%5\\%\. The model therefore extrapolates speed within the familiar preferred\-walking strategy more accurately than the unseen constrained strategy\.

The pooled hip, knee, and ankle results across the different experiments are reported in Supplementary Information, Section 7\.4, and Supplementary Figure6\. The ankle profile is not recovered: its pooled correlation isr=0\.492r=0\.492, its amplitude is approximately one order of magnitude below the reference, and one profile is anticorrelated\. This can be attributed to the fact that the ankle readout contains only two distal markers\.

A separate in\-distribution control using walking conditions seen during training is shown in Supplementary Figure7\. The inferred hip and knee profiles again follow the inverse\-dynamics reference across speed and cadence\. At the knee, the correlation isr=0\.81r=0\.81\. The control is therefore no more accurate than the evaluation containing withheld conditions\.We also assess the contribution of the angular\-momentum channel to the physical fidelity of the inferred moments\. Removing this channel leaves rollout accuracy essentially unchanged but degrades the joint\-moment estimates \(see ablation experiments detailed in Supplementary Information, Section 8\.1 and Supplementary Figure 8\.\)

The learned response operators provide a second mechanical readout\. The per\-node position\- and velocity\-response operators used by the Newmark solve align with the corresponding finite\-element matrices in their dominant direction, with cosine similarities of approximately0\.850\.85and0\.920\.92, and recover their spatial organisation, with correlations of approximately0\.630\.63and0\.940\.94between nodal traces\. They do not recover the directional anisotropy or absolute scale\. Moreover, the position\-response operator is nearly orthogonal to the pointwise Jacobian of the decoded force, but predicts the change in the model’s momentum response to a velocity perturbation with cosine similarity0\.990\.99to1\.001\.00\. Supplementary Information, Section 4\.9 defines and reports these probes; Supplementary Figure 4 summarizes the comparison with the finite\-element tangent\. These quantities are therefore response operators of the learned update, not identified material tangents\.They give a relative, interpretive readout of where the structure is stiff or soft, not a calibrated material tangent\.

Thus,Newmark\-β\\beta\-DGNreturns more than a trajectory\. On the walking biomechanics benchmark, the joint\-moment readouts assembled from the internal forces used to advance the marker graph follow independently measured hip and knee demand across familiar and withheld gait conditions\.On the beam, the same response operators recover the structure of the finite\-element tangent, and the rollout reproduces the beam’s fundamental frequency and dominant vibration mode \(Section[2\.2](https://arxiv.org/html/2609.30344#S2.SS2)\)\.None of these mechanical targets is supplied during training\. The model therefore exposes aspects of the mechanics underlying its predictions, rather than providing trajectories alone\.

### 2\.6Conservation, scaling and computational cost

Figure 10:Edge scaling and computational cost\.a, Directed edges per graph for the beam mesh \(left\) and the protein backbone \(right\)\.Newmark\-β\\beta\-DGNreplaces a distance\-based contact graph with2​N2Nhub edges, which reduces the protein edge count by a factor of16\.316\.3\. On the fixed beam mesh, no model uses a distance cutoff and no corresponding reduction arises\.b, Peak memory against wall\-clock time per training iteration, using a batch size of four and latent widths of128128for the protein and6464for the beam\. Two protein configurations are exceptions: EGNO andDGNwith global edges require9622​MB9622\\,\\mathrm\{MB\}and10 048​MB10\\,048\\,\\mathrm\{MB\}, respectively, at width128128and therefore do not fit on the11​GB11\\,\\mathrm\{GB\}device; their plotted points are measured at width6464\. Five forward\-and\-backward iterations are timed after one discarded warm\-up iteration, with each model profiled in a separate process pinned to a single device and retaining its published edge topology\. Single measurement per configuration\.We examine two consequences of the proposed formulation separately: its momentum\-conservation properties and the computational cost of the virtual hub\.

##### Momentum conservation\.

Before time integration, the physical edges conserve linear and angular momentum pairwise, as inherited fromDGN\[[15](https://arxiv.org/html/2609.30344#bib.bib1)\]\. For each interacting pair, the two directed edges carry equal\-and\-opposite linear\- and angular\-momentum fluxes and use the same force\-application point\. Their contributions therefore cancel exactly when summed over the pair \(Methods, Section[4\.4\.1](https://arxiv.org/html/2609.30344#S4.SS4.SSS1); Supplementary Information, Section 11\.1\.1\)\.

The virtual hub satisfies a different, collective balance\. Its decoded linear\- and angular\-momentum fluxes are projected so that the corresponding sums over all hub edges are zero \(Methods, Section[4\.4\.1](https://arxiv.org/html/2609.30344#S4.SS4.SSS1), Equation \([13](https://arxiv.org/html/2609.30344#S4.E13)\); Supplementary Information, Section 9\.3\.2\)\. The force projection therefore guarantees that the hub contributes no net linear momentum\. For angular momentum, however, zero summed angular\-momentum flux is not sufficient to guarantee zero total angular contribution: the hub forces act at different application points and can therefore generate a non\-zero net moment about a common origin\. The hub thus enforces collective linear balance and a zero summed angular\-momentum flux, but not exact total angular\-momentum balance \(Supplementary Information, Section 11\.1\.2\)\.

A separate conservation error is introduced by the independent semi\-implicit nodal solves\. At nodeii, the total force drive𝐛i\\mathbf\{b\}\_\{i\}is transformed by the local3×33\\times 3coefficient matrix𝐀i−1\\mathbf\{A\}\_\{i\}^\{\-1\}, and the velocity update contains an additional position\-response term proportional to𝐊i​𝐯i,t\\mathbf\{K\}\_\{i\}\\mathbf\{v\}\_\{i,t\}\(Methods, Section[4\.6](https://arxiv.org/html/2609.30344#S4.SS6)\)\. Before this solve, the physical\-edge forces cancel pairwise and the projected hub\-edge forces sum to zero, so

∑i𝐛i=𝐅ext,\\sum\_\{i\}\\mathbf\{b\}\_\{i\}=\\mathbf\{F\}^\{\\mathrm\{ext\}\},where reactions at prescribed degrees of freedom are included in𝐅ext\\mathbf\{F\}^\{\\mathrm\{ext\}\}\.

Exact preservation of this balance would give∑iΔ​𝐩i=δ​t​𝐅ext\\sum\_\{i\}\\Delta\\mathbf\{p\}\_\{i\}=\\delta t\\,\\mathbf\{F\}^\{\\mathrm\{ext\}\}\. Instead, the semi\-implicit nodal update gives the residual

𝐑P=δ​t​∑i\(𝐀i−1−𝐈3\)​𝐛i−12​δ​t2​∑i𝐀i−1​𝐊i​𝐯i,t\.\\mathbf\{R\}\_\{P\}=\\delta t\\sum\_\{i\}\\left\(\\mathbf\{A\}\_\{i\}^\{\-1\}\-\\mathbf\{I\}\_\{3\}\\right\)\\mathbf\{b\}\_\{i\}\-\\frac\{1\}\{2\}\\delta t^\{2\}\\sum\_\{i\}\\mathbf\{A\}\_\{i\}^\{\-1\}\\mathbf\{K\}\_\{i\}\\mathbf\{v\}\_\{i,t\}\.\(2\)The two terms have distinct origins\. The first appears because the balanced nodal drives𝐛i\\mathbf\{b\}\_\{i\}are transformed by node\-dependent matrices𝐀i−1\\mathbf\{A\}\_\{i\}^\{\-1\}, so their sum is not generally preserved\. The second is the contribution of the position\-response term𝐊i​𝐯i,t\\mathbf\{K\}\_\{i\}\\mathbf\{v\}\_\{i,t\}\. In a fully assembled system, such internal responses are accompanied by cross\-node reactions that cancel when the coupled equations are summed\. The independent nodal solves omit these cross\-node reactions, leaving a finite residual \(Supplementary Information, Section 11\.2\)\.

For bounded masses, states and response operators, both terms are𝒪⁡\(δ​t2\)\\mathcal\{O\}\(\\delta t^\{2\}\)\. The absolute one\-substep momentum residual is therefore second order in the substep size, and its value relative to a non\-zero external impulse is first order \(Supplementary Information, Section 11\.3\)\. The twelve\-node synthetic test confirms this behaviour\. The complete update gives an observed relative\-residual order of0\.840\.84; retaining the two contributions separately gives orders of0\.970\.97and0\.810\.81, respectively; and removing both reduces the residual to approximately1\.7×10−81\.7\\times 10^\{\-8\}, at numerical round\-off \(Supplementary Information, Section 11\.4 and Supplementary Table 20\)\. The residual also vanishes in the explicit\-response limit𝐊i,𝐃i→0\\mathbf\{K\}\_\{i\},\\mathbf\{D\}\_\{i\}\\rightarrow 0, where the linear\-momentum update reduces to the conserving explicit force update ofDGN, apart from the trapezoidal position advance \(Supplementary Information, Section 11\.3\)\.

##### Computational cost\.

The virtual hub adds2​N2Ndirected edges forNNphysical nodes, so its additional graph cost scales linearly with system size\. This is advantageous when the hub replaces a much denser long\-range interaction graph\. On the protein benchmark,Newmark\-β\\beta\-DGNuses 1708 covalent\-bond edges and 1710 hub edges, giving 3418 directed edges in total, whereas the distance\-based baselines use approximately55 61055\\,610directed edges, or16\.316\.3times as many \(Fig\.[10](https://arxiv.org/html/2609.30344#S2.F10)a\)\. This reduction is reflected in memory use: at the common profiling setting,Newmark\-β\\beta\-DGNrequires1067​MB1067\\,\\mathrm\{MB\}of peak memory, compared with3843​MB3843\\,\\mathrm\{MB\}for EGHN, while EGNO andDGNwith global edges exceed the available memory at the common latent width \(Fig\.[10](https://arxiv.org/html/2609.30344#S2.F10)b\)\.

This advantage is not expected when the underlying graph is already sparse\. On the human skeleton, the edge budgets are similar \(122 versus 130 directed edges\), while on the beam no distance\-based contact graph is used and the hub therefore provides no meaningful reduction in edge count\. Consistently, MGN and EGHN are cheaper thanNewmark\-β\\beta\-DGNper training iteration on the beam benchmark\. The computational benefit of the hub therefore arises specifically when system\-wide communication would otherwise require substantially denser connectivity\.

## 3Discussion

This study introducesNewmark\-β\\beta\-DGN, a mechanically structured framework for learning coarse\-step dynamics and examining the mechanical quantities that generate them from discretely sampled trajectories\. Coarse observations create two coupled challenges: the mechanical response can evolve substantially within one observed transition, while mechanical influence can propagate across much of the system over the same interval\.Newmark\-β\\beta\-DGNaddresses the first through a semi\-implicit update inspired by the average\-acceleration Newmark method, in which conserved linear and angular\-momentum fluxes, observed or decoded external forces, learned inverse mass and inertia, and matrix\-valued response operators determine the finite\-time state increment\. The second is addressed through an operator\-weighted virtual hub motivated by the star topology obtained from a rank\-one approximation of implicit non\-local coupling\. Because these learned mechanical quantities are used to generate the predicted state update through a mechanically structured formulation, they remain accessible after training as mechanical readouts\.

We test these two architectural components and their mechanical readouts across four settings\. The finite\-element beam \(Section[2\.2](https://arxiv.org/html/2609.30344#S2.SS2)\) provides a controlled test of coarse\-step rollout and extrapolation across loading, geometry, stiffness and mesh resolution\. Human motion capture \(Section[2\.3](https://arxiv.org/html/2609.30344#S2.SS3)\) and protein dynamics \(Section[2\.4](https://arxiv.org/html/2609.30344#S2.SS4)\) test the formulation on systems with different graph structures and scales, while walking biomechanics \(Section[2\.5](https://arxiv.org/html/2609.30344#S2.SS5)\) tests whether quantities learned only through trajectory supervision retain independently verifiable mechanical meaning\.

The beam results show why finite\-time response and system\-wide coupling must be considered together\.Newmark\-β\\beta\-DGNmaintains low, bounded error across the test and extrapolation cases and reproduces the fundamental vibration frequency and dominant mode without direct supervision\. Increasing explicit temporal resolution does not reproduce this behavior:DGNand IGNS improve with finer sub\-time stepping, but their errors continue to accumulate\. with the mean whole\-body error decreasing from3\.00%3\.00\\%for twelve\-substepDGNto0\.54%0\.54\\%forNewmark\-β\\beta\-DGN\. Hierarchical EGHN and EGHNO shorten communication paths across the graph, yet both become non\-finite during coarse\-step rollout\. These comparisons show that neither finer explicit integration nor broader communication alone is sufficient\. The controlled ablation makes this separation explicit: adding the operator\-weighted hub to the explicit update reduces the mean whole\-body error from3\.00%3\.00\\%to1\.26%1\.26\\%, introducing the semi\-implicit solve reduces it further to0\.70%0\.70\\%, and the full formulation reaches0\.54%0\.54\\%\. These results support a formulation in which finite\-time response and system\-wide coupling are incorporated jointly through structure inspired by computational mechanics\.

The same formulation also applies beyond the continuum setting\. On human motion,Newmark\-β\\beta\-DGNretains rotation\-equivariant prediction, while on the 855\-node protein system the hub provides system\-wide connectivity with 16\.3 times fewer edges than the distance\-based interaction graphs used by the baselines\. These results establish the reach of the formulation across different topologies and scales\. The hub’s edge saving, however, is specific to settings in which long\-range interactions would otherwise require a dense contact graph; it does not arise on a fixed mesh\.

The stronger test is whether a model trained only on trajectories can infer forces and response operators that were never observed during training, and whether those inferred quantities agree with independent mechanical references\. Joint\-moment estimates assembled from internal forces inferred solely from walking kinematics track independently derived hip and knee demand trends, with correlations ofr=0\.94r=0\.94andr=0\.87r=0\.87, despite receiving no force or moment supervision\. Removing the angular\-momentum channel leaves rollout accuracy essentially unchanged while degrading the inferred joint moments, further separating predictive and mechanical fidelity\. On the beam, the learned response operators recover relative directional and spatial organization of the finite\-element response, but not its anisotropy or absolute scale\. These results show that mechanically informative quantities can emerge from trajectory supervision when they participate directly in generating the predicted dynamics\.

The beam results also define the present limitations and directions for extension\. A single rank\-one hub becomes less effective for the longest beams, motivating higher\-rank or multiple\-hub representations that retain more of the non\-local interaction\. Independent nodal solves avoid a global implicit system but do not guarantee exact finite\-step momentum conservation; future work could therefore consider more strongly coupled updates that improve update\-level conservation while preserving coarse\-step efficiency\. Observed external interactions or limited calibration data could additionally anchor the scale of the learned mechanical quantities, extending the present relative readouts towards absolute forces and moments\.

Overall,Newmark\-β\\beta\-DGNestablishes a mechanically structured approach to coarse\-step learning in which conserved interactions, finite\-time response, and system\-wide coupling are learned jointly within the state transition, linking accurate dynamical prediction with access to internal mechanics underlying the observed motion\. By making inferred forces and response operators accessible from trajectory data,Newmark\-β\\beta\-DGNopens a route to studying mechanical behavior when direct measurements are unavailable\. These inferred mechanical quantities can reveal how mechanical demand is distributed across the system and how it changes with operating conditions, paving the way for unsupervised virtual sensing for condition monitoring applications\. With its demonstrated generalization across unseen loads, operating conditions, and configurations,Newmark\-β\\beta\-DGNcan also support what\-if scenario evaluation by showing how mechanical demand redistributes under changed loads, configurations, or operating conditions\.

## 4Methods

### 4\.1Graph representation and model inputs

Newmark\-β\\beta\-DGNtakes as input a graph representation of the observed system state\. We denote the physical graph by𝒢=\(𝒱,ℰ\)\\mathcal\{G\}=\(\\mathcal\{V\},\\mathcal\{E\}\), where theNNnodes in𝒱\\mathcal\{V\}represent the discretised system components and the edges inℰ\\mathcal\{E\}represent their direct physical connections\. Each physical connection is represented by two directed edges,i→ji\\\!\\to\\\!jandj→ij\\\!\\to\\\!i\. The physical nodes correspond to finite\-element vertices for the beam, anatomical markers for the human systems, and backbone atoms for the protein; the corresponding physical connections follow the mesh connectivity, anatomical skeleton, and covalent bonds, respectively\.

Each physical nodeiicarries its position𝐫i\\mathbf\{r\}\_\{i\}and velocity𝐯i\\mathbf\{v\}\_\{i\}\. Velocities are obtained by causal finite differencing when they are not directly observed\. The model additionally carries a spin𝝎i\\bm\{\\omega\}\_\{i\}as a latent angular\-momentum state, initialized to zero at the beginning of each prediction interval\. Additional node and edge inputs are case\-specific and are given in Supplementary Information, Sections 4\.1–7\.1\.

The physical graph is augmented with one virtual hub connected to every physical node through bi\-directional virtual edges, adding2​N2Nvirtual edges \(Fig\.[1](https://arxiv.org/html/2609.30344#S1.F1)b\)\. Physical and virtual edges are distinguished by an edge\-type attribute\. The hub is not a physical body, carries no external load, and is not advanced as an independent degree of freedom\. Its position, velocity, and spin are recomputed during the model update as described in Section[4\.5](https://arxiv.org/html/2609.30344#S4.SS5)\.

### 4\.2Model update overview

Given the augmented graph at the current observation time,Newmark\-β\\beta\-DGNpredicts the physical\-node state after an intervalΔ​t\\Delta t\. The interval is divided intoSSsubsteps of size

δ​t=Δ​tS,\\delta t=\\frac\{\\Delta t\}\{S\},\(3\)with one complete message\-passing round performed at each substep \(Fig\.[1](https://arxiv.org/html/2609.30344#S1.F1)c\)\. We useS=4S=4throughout this work\.

Each substep follows the same sequence\. The hub state is first recomputed from the current physical\-node states using the nodal response operators retained from the preceding substep; for the first substep, these operators are initialized as the identity\. The resulting augmented graph state is then encoded using the equivariant representation of Section[4\.3](https://arxiv.org/html/2609.30344#S4.SS3)\. From this representation, the edge decoders produce the linear\- and angular\-momentum fluxes and the translational and rotational response operators, while the node decoder produces the inverse mass and inverse inertia and, when the external force is unobserved, its components in the global frame, as described in Section[4\.4](https://arxiv.org/html/2609.30344#S4.SS4)\. The virtual\-edge fluxes are collectively projected and the edge response operators are aggregated at the physical nodes\. Together with the observed or decoded external force, these quantities enter the translational and rotational semi\-implicit Newmark updates of Section[4\.6](https://arxiv.org/html/2609.30344#S4.SS6)\. The updated physical state and the newly aggregated response operators are then carried to the following substep\.

AfterSSsubsteps, the updated physical\-node positions and velocities define the prediction at the next observation time\. Section[4\.7](https://arxiv.org/html/2609.30344#S4.SS7)describes how this transition is learned from kinematic supervision and how the mechanical quantities used to generate it remain accessible after training\.

### 4\.3Equivariant interaction encoding and processing

We retain the edge\-local scalarisation–vectorisation construction ofDGN\[[15](https://arxiv.org/html/2609.30344#bib.bib1)\]\. Each directed edgei​jijis equipped with an orthonormal frame𝐑i​j=\[𝐚i​j​\|𝐛i​j\|​𝐜i​j\]\\mathbf\{R\}\_\{ij\}=\[\\mathbf\{a\}\_\{ij\}\|\\mathbf\{b\}\_\{ij\}\|\\mathbf\{c\}\_\{ij\}\]whose first axis lies along the edge,

𝐚i​j=𝐫j−𝐫i∥𝐫j−𝐫i∥\.\\mathbf\{a\}\_\{ij\}=\\frac\{\\mathbf\{r\}\_\{j\}\-\\mathbf\{r\}\_\{i\}\}\{\\lVert\\mathbf\{r\}\_\{j\}\-\\mathbf\{r\}\_\{i\}\\rVert\}\.\(4\)The remaining axes follow the construction ofDGNand make the complete frame rotation\-equivariant, translation\-invariant, and antisymmetric under interchange of the two connected nodes,𝐑j​i=−𝐑i​j\\mathbf\{R\}\_\{ji\}=\-\\mathbf\{R\}\_\{ij\}\.

The construction of the node embedding depends on whether the external force is observed\. When the external force is observed, as for the beam, no global reference frame is required for force inference\. The case\-specific scalar node features are passed directly through the node encoder to obtain the invariant node embedding𝐡i\\mathbf\{h\}\_\{i\}, while the observed force is provided separately to the mechanical update\.

When the external force is unobserved, as for the human and protein systems, we construct a second, graph\-level global frame𝐑glob\\mathbf\{R\}^\{\\mathrm\{glob\}\}\. This frame is used both to express the vector\-valued node features as invariant scalars and to reconstruct the decoded external force\. From the physical\-node graph, excluding the virtual hub, we form the combinatorial Laplacian𝐋=𝐃−𝐀\\mathbf\{L\}=\\mathbf\{D\}\-\\mathbf\{A\}, where𝐀\\mathbf\{A\}is the adjacency matrix of the physical graph and𝐃\\mathbf\{D\}is its diagonal degree matrix, withDi​i=∑jAi​jD\_\{ii\}=\\sum\_\{j\}A\_\{ij\}\. Its eigenvectors satisfy

𝐋​ϕk=λk​ϕk,0=λ0≤λ1≤⋯,\\mathbf\{L\}\\bm\{\\phi\}\_\{k\}=\\lambda\_\{k\}\\bm\{\\phi\}\_\{k\},\\qquad 0=\\lambda\_\{0\}\\leq\\lambda\_\{1\}\\leq\\cdots,\(5\)whereϕk\\bm\{\\phi\}\_\{k\}is thekkth eigenvector andλk\\lambda\_\{k\}its corresponding eigenvalue\. The eigenvalueλ0=0\\lambda\_\{0\}=0corresponds to the constant mode of the connected physical graph\.

We retain the first two non\-trivial eigenvectors,ϕ1\\bm\{\\phi\}\_\{1\}andϕ2\\bm\{\\phi\}\_\{2\}\. Because−ϕk\-\\bm\{\\phi\}\_\{k\}is an equally valid eigenvector, its sign is fixed by locating the component with the largest magnitude and multiplying the entire eigenvector by−1\-1when that component is negative\.

At each substep, the eigenvectors are combined with the current physical\-node positions to form two three\-dimensional anchor vectors,

𝐪k=∑i∈𝒱ϕk​\(i\)​𝐫i,k∈\{1,2\},\\mathbf\{q\}\_\{k\}=\\sum\_\{i\\in\\mathcal\{V\}\}\\phi\_\{k\}\(i\)\\,\\mathbf\{r\}\_\{i\},\\qquad k\\in\\\{1,2\\\},\(6\)where𝐫i∈ℝ3\\mathbf\{r\}\_\{i\}\\in\\mathbb\{R\}^\{3\}is the current position of nodeiiandϕk​\(i\)\\phi\_\{k\}\(i\)is the component ofϕk\\bm\{\\phi\}\_\{k\}associated with that node\. Since the non\-trivial eigenvectors are orthogonal to the constant mode,

∑i∈𝒱ϕk​\(i\)=0,k∈\{1,2\},\\sum\_\{i\\in\\mathcal\{V\}\}\\phi\_\{k\}\(i\)=0,\\qquad k\\in\\\{1,2\\\},\(7\)translating every node by a common vector𝐜\\mathbf\{c\}leaves the anchors unchanged,

∑i∈𝒱ϕk​\(i\)​\(𝐫i\+𝐜\)=𝐪k\+𝐜​∑i∈𝒱ϕk​\(i\)=𝐪k\.\\sum\_\{i\\in\\mathcal\{V\}\}\\phi\_\{k\}\(i\)\(\\mathbf\{r\}\_\{i\}\+\\mathbf\{c\}\)=\\mathbf\{q\}\_\{k\}\+\\mathbf\{c\}\\sum\_\{i\\in\\mathcal\{V\}\}\\phi\_\{k\}\(i\)=\\mathbf\{q\}\_\{k\}\.\(8\)
The anchors are orthonormalised to define the frame axes,

𝐞1=𝐪1∥𝐪1∥,𝐞2=𝐪2−\(𝐪2⋅𝐞1\)​𝐞1‖𝐪2−\(𝐪2⋅𝐞1\)​𝐞1‖,𝐞3=𝐞1×𝐞2,\\mathbf\{e\}\_\{1\}=\\frac\{\\mathbf\{q\}\_\{1\}\}\{\\lVert\\mathbf\{q\}\_\{1\}\\rVert\},\\qquad\\mathbf\{e\}\_\{2\}=\\frac\{\\mathbf\{q\}\_\{2\}\-\(\\mathbf\{q\}\_\{2\}\\\!\\cdot\\\!\\mathbf\{e\}\_\{1\}\)\\mathbf\{e\}\_\{1\}\}\{\\left\\\|\\mathbf\{q\}\_\{2\}\-\(\\mathbf\{q\}\_\{2\}\\\!\\cdot\\\!\\mathbf\{e\}\_\{1\}\)\\mathbf\{e\}\_\{1\}\\right\\\|\},\\qquad\\mathbf\{e\}\_\{3\}=\\mathbf\{e\}\_\{1\}\\times\\mathbf\{e\}\_\{2\},\(9\)where⋅\\cdotdenotes the Euclidean inner product and×\\timesthe three\-dimensional cross product\. The resulting global frame is

𝐑glob=\[𝐞1​\|𝐞2\|​𝐞3\]\.\\mathbf\{R\}^\{\\mathrm\{glob\}\}=\[\\,\\mathbf\{e\}\_\{1\}\\,\|\\,\\mathbf\{e\}\_\{2\}\\,\|\\,\\mathbf\{e\}\_\{3\}\\,\]\.\(10\)The frame is translation\-invariant by Equation \([8](https://arxiv.org/html/2609.30344#S4.E8)\); under a global rotation, the node positions, anchors, and frame axes rotate together, making the construction rotation\-equivariant\.

For systems using this frame, the node vector features, including velocity and hub\-relative position𝐫irel=𝐫i−𝐫H\\mathbf\{r\}^\{\\mathrm\{rel\}\}\_\{i\}=\\mathbf\{r\}\_\{i\}\-\\mathbf\{r\}\_\{H\}, are projected onto𝐑glob\\mathbf\{R\}^\{\\mathrm\{glob\}\}to obtain invariant scalar components\. These are combined with the scalar node features and passed through the node encoder to obtain𝐡i\\mathbf\{h\}\_\{i\}\.

Given the node embeddings, each directed edge forms its interaction embedding\. The vector features of the sender and receiver are projected onto𝐑i​j\\mathbf\{R\}\_\{ij\}and−𝐑i​j\-\\mathbf\{R\}\_\{ij\}, respectively, to obtain invariant scalar edge features\. FollowingDGN\[[15](https://arxiv.org/html/2609.30344#bib.bib1)\], these scalars are combined with the symmetric node representation𝐡i\+𝐡j\\mathbf\{h\}\_\{i\}\+\\mathbf\{h\}\_\{j\}, the edge length, and the edge type to obtain the invariant interaction embeddingϵi​j\\bm\{\\epsilon\}\_\{ij\}\.

The recurrent edge processor ofDGN, including its edge\-level latent memory, is retained across theSSsubsteps\. The first substep initializes the edge latent representation, while subsequent substeps combine the current interaction with the latent state retained from the preceding round\. We refer to\[[15](https://arxiv.org/html/2609.30344#bib.bib1)\]for the complete edge\-frame construction, scalarization and recurrent processing\.

### 4\.4Decoded interactions and response operators

#### 4\.4\.1Edge\-wise Interaction Fluxes

From the edge interaction embeddings and the node embeddings,Newmark\-β\\beta\-DGNdecodes the quantities required by the semi\-implicit update\. On each edge, in addition to the linear and angular\-momentum fluxes ofDGN\[[15](https://arxiv.org/html/2609.30344#bib.bib1)\], we decode translational and rotational response operators\. From each node embedding, a node decoder produces the node’s inverse mass and inverse inertia and, when the external force is unobserved, its external force components in the global frame of Section[4\.3](https://arxiv.org/html/2609.30344#S4.SS3)\.

For each edge, invariant scalar coefficients are decoded from the invariant interaction embeddingϵi​j\\bm\{\\epsilon\}\_\{ij\}and combined with the edge\-local basis to reconstruct the linear\-momentum flux𝐟i​j\\mathbf\{f\}\_\{ij\}and the angular\-momentum flux𝐀i​j\\mathbf\{A\}\_\{ij\}\. By the antisymmetric edge\-frame construction and the coefficients shared between the two traversals,

𝐟i​j=−𝐟j​i,𝐀i​j=−𝐀j​i\.\\mathbf\{f\}\_\{ij\}=\-\\mathbf\{f\}\_\{ji\},\\qquad\\mathbf\{A\}\_\{ij\}=\-\\mathbf\{A\}\_\{ji\}\.\(11\)The two directions of each physical connection also share a common reference point,𝐫i​j0=𝐫j​i0\\mathbf\{r\}^\{0\}\_\{ij\}=\\mathbf\{r\}^\{0\}\_\{ji\}, about which the angular\-momentum exchange is balanced\. The spin torque delivered to nodeiiis obtained by removing the orbital contribution of the force,

𝝉i​j=𝐀i​j−\(𝐫i−𝐫i​j0\)×𝐟i​j\.\\bm\{\\tau\}\_\{ij\}=\\mathbf\{A\}\_\{ij\}\-\(\\mathbf\{r\}\_\{i\}\-\\mathbf\{r\}^\{0\}\_\{ij\}\)\\times\\mathbf\{f\}\_\{ij\}\.\(12\)Together with the antisymmetric fluxes, these relations give pairwise linear and angular\-momentum balance on physical interactions before time integration, as inherited fromDGN\[[15](https://arxiv.org/html/2609.30344#bib.bib1)\]\(Supplementary Information, Section 11\.1\.1\)\.

The virtual edges connecting the hub to the physical nodes decode the same fluxes as the physical edges\. For these edges, the decoded force and angular\-momentum fluxes are projected onto the subspace of zero collective sum, so that the hub only redistributes internal force and torque contributions among the physical nodes,

𝐟~H​j=𝐟H​j−1N​∑k𝐟H​k,𝐀~H​j=𝐀H​j−1N​∑k𝐀H​k,\\widetilde\{\\mathbf\{f\}\}\_\{Hj\}=\\mathbf\{f\}\_\{Hj\}\-\\frac\{1\}\{N\}\\sum\_\{k\}\\mathbf\{f\}\_\{Hk\},\\qquad\\widetilde\{\\mathbf\{A\}\}\_\{Hj\}=\\mathbf\{A\}\_\{Hj\}\-\\frac\{1\}\{N\}\\sum\_\{k\}\\mathbf\{A\}\_\{Hk\},\(13\)with∑j𝐟~H​j=𝟎\\sum\_\{j\}\\widetilde\{\\mathbf\{f\}\}\_\{Hj\}=\\mathbf\{0\}and∑j𝐀~H​j=𝟎\\sum\_\{j\}\\widetilde\{\\mathbf\{A\}\}\_\{Hj\}=\\mathbf\{0\}\. Unlike the pairwise balance on physical edges, this balance is collective over all virtual edges\. Because the virtual edges have distinct reference points, the zero\-sum projection does not by itself guarantee zero total moment about a common origin; the corresponding residual couple is derived in Supplementary Information, Sections 9\.3\.2 and 11\.1\.2\.

#### 4\.4\.2Edge\-wise Response Operators

Alongside the momentum fluxes, each edge interaction embedding is also decoded into four state\-dependent matrix\-valued response operators that represent the finite\-time response to position and velocity increments within the corresponding semi\-implicit substep: the translational position\- and velocity\-response operators𝐊i​j\\mathbf\{K\}\_\{ij\}and𝐃i​j\\mathbf\{D\}\_\{ij\}, and their rotational counterparts𝐊i​jrot\\mathbf\{K\}^\{\\mathrm\{rot\}\}\_\{ij\}and𝐃i​jrot\\mathbf\{D\}^\{\\mathrm\{rot\}\}\_\{ij\}\. These operators are learned coefficients of the update rather than derivatives of the decoded force or identified material stiffness and damping matrices\.

Each response operator is parameterized independently by six invariant scalars\. For a generic operator

𝐖i​j∈\{𝐊i​j,𝐃i​j,𝐊i​jrot,𝐃i​jrot\},\\mathbf\{W\}\_\{ij\}\\in\\left\\\{\\mathbf\{K\}\_\{ij\},\\mathbf\{D\}\_\{ij\},\\mathbf\{K\}^\{\\mathrm\{rot\}\}\_\{ij\},\\mathbf\{D\}^\{\\mathrm\{rot\}\}\_\{ij\}\\right\\\},\(14\)the six scalars define a lower\-triangular factor

𝐋i​j\(𝐖\)=\[d000s1d10s3s4d2\],\(d0,d1,d2\)=softplus⁡\(s0,s2,s5\)\+10−4,\\mathbf\{L\}^\{\(\\mathbf\{W\}\)\}\_\{ij\}=\\begin\{bmatrix\}d\_\{0\}&0&0\\\\ s\_\{1\}&d\_\{1\}&0\\\\ s\_\{3\}&s\_\{4\}&d\_\{2\}\\end\{bmatrix\},\\qquad\(d\_\{0\},d\_\{1\},d\_\{2\}\)=\\operatorname\{softplus\}\(s\_\{0\},s\_\{2\},s\_\{5\}\)\+10^\{\-4\},\(15\)and the operator is reconstructed in global coordinates as

𝐖i​j=𝐑i​j​𝐋i​j\(𝐖\)​𝐋i​j\(𝐖\)⊤​𝐑i​j⊤\.\\mathbf\{W\}\_\{ij\}=\\mathbf\{R\}\_\{ij\}\\mathbf\{L\}^\{\(\\mathbf\{W\}\)\}\_\{ij\}\\mathbf\{L\}^\{\(\\mathbf\{W\}\)\\top\}\_\{ij\}\\mathbf\{R\}\_\{ij\}^\{\\top\}\.\(16\)The positive diagonal of𝐋i​j\(𝐖\)\\mathbf\{L\}^\{\(\\mathbf\{W\}\)\}\_\{ij\}makes each operator symmetric positive definite, while its construction in the edge\-local frame makes it rotation\-equivariant\. Because the decoder coefficients are shared between the two traversals and𝐑j​i=−𝐑i​j\\mathbf\{R\}\_\{ji\}=\-\\mathbf\{R\}\_\{ij\}, the quadratic construction is unchanged under edge reversal,

𝐖j​i=𝐖i​j\.\\mathbf\{W\}\_\{ji\}=\\mathbf\{W\}\_\{ij\}\.\(17\)
The same construction is used on physical and virtual edges; unlike the momentum fluxes, the virtual\-edge response operators are not subjected to the zero\-sum projection of Equation \([13](https://arxiv.org/html/2609.30344#S4.E13)\)\. The edge operators are summed over all edges incident on each physical node,

𝐊i\\displaystyle\\mathbf\{K\}\_\{i\}=∑j∈𝒩⁡\(i\)𝐊i​j,\\displaystyle=\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\mathbf\{K\}\_\{ij\},𝐃i\\displaystyle\\qquad\\mathbf\{D\}\_\{i\}=∑j∈𝒩⁡\(i\)𝐃i​j,\\displaystyle=\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\mathbf\{D\}\_\{ij\},\(18\)𝐊irot\\displaystyle\\mathbf\{K\}^\{\\mathrm\{rot\}\}\_\{i\}=∑j∈𝒩⁡\(i\)𝐊i​jrot,\\displaystyle=\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\mathbf\{K\}^\{\\mathrm\{rot\}\}\_\{ij\},𝐃irot\\displaystyle\\mathbf\{D\}^\{\\mathrm\{rot\}\}\_\{i\}=∑j∈𝒩⁡\(i\)𝐃i​jrot,\\displaystyle=\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\mathbf\{D\}^\{\\mathrm\{rot\}\}\_\{ij\},where𝒩⁡\(i\)\\mathcal\{N\}\(i\)contains the physical neighbours of nodeiiand the virtual hub\. These nodal operators enter both the operator\-weighted hub of Section[4\.5](https://arxiv.org/html/2609.30344#S4.SS5)and the semi\-implicit update of Section[4\.6](https://arxiv.org/html/2609.30344#S4.SS6)\.

#### 4\.4\.3Node\-wise inverse mass, inverse inertia, and external force

In addition to the edge\-wise quantities, the node embedding𝐡i\\mathbf\{h\}\_\{i\}is decoded into the node\-wise quantities required by the semi\-implicit update\. The node decoder produces positive scalar inverse mass and inverse inertia,

𝐌i−1=mi−1​𝐈3,𝐈i−1=Ii−1​𝐈3,mi−1\>0,Ii−1\>0,\\mathbf\{M\}^\{\-1\}\_\{i\}=m\_\{i\}^\{\-1\}\\mathbf\{I\}\_\{3\},\\qquad\\mathbf\{I\}^\{\-1\}\_\{i\}=I\_\{i\}^\{\-1\}\\mathbf\{I\}\_\{3\},\\qquad m\_\{i\}^\{\-1\}\>0,\\quad I\_\{i\}^\{\-1\}\>0,\(19\)where𝐈3\\mathbf\{I\}\_\{3\}is the3×33\\times 3identity matrix\. Their isotropic form preserves rotation equivariance\. The inverse mass enters the translational update, while the inverse inertia enters the rotational update in Section[4\.6](https://arxiv.org/html/2609.30344#S4.SS6)\.

The treatment of the external force depends on whether it is observed\. When it is unobserved, the node decoder maps the invariant embedding𝐡i\\mathbf\{h\}\_\{i\}to three invariant components𝐬i∈ℝ3\\mathbf\{s\}\_\{i\}\\in\\mathbb\{R\}^\{3\}, which are reconstructed in global coordinates using the frame of Section[4\.3](https://arxiv.org/html/2609.30344#S4.SS3),

𝐟iext=𝐑glob​𝐬i\.\\mathbf\{f\}^\{\\mathrm\{ext\}\}\_\{i\}=\\mathbf\{R\}^\{\\mathrm\{glob\}\}\\mathbf\{s\}\_\{i\}\.\(20\)Because𝐬i\\mathbf\{s\}\_\{i\}is invariant and𝐑glob\\mathbf\{R\}^\{\\mathrm\{glob\}\}rotates with the system, the decoded force is rotation\-equivariant and translation\-invariant\. When the external force is observed, no force decoder is used and the observed force is provided directly to the semi\-implicit update\.

### 4\.5Operator\-weighted virtual hub

The response operators also set the virtual\-hub state at the start of each substep\. At the first substep, the nodal operators are initialized as the identity, so the hub position, velocity, and spin reduce to arithmetic averages of the corresponding physical\-node states\. At each subsequent substep, the hub is recomputed from the updated physical\-node states using the operators aggregated in the preceding substep; it therefore carries no independently propagated state\.

The hub state is

𝐫H\\displaystyle\\mathbf\{r\}\_\{H\}=\(∑i𝐊i\)−1​∑i𝐊i​𝐫i,\\displaystyle=\\left\(\\sum\_\{i\}\\mathbf\{K\}\_\{i\}\\right\)^\{\-1\}\\sum\_\{i\}\\mathbf\{K\}\_\{i\}\\mathbf\{r\}\_\{i\},\(21\)𝐯H\\displaystyle\\mathbf\{v\}\_\{H\}=\(∑i𝐃i\)−1​∑i𝐃i​𝐯i,\\displaystyle=\\left\(\\sum\_\{i\}\\mathbf\{D\}\_\{i\}\\right\)^\{\-1\}\\sum\_\{i\}\\mathbf\{D\}\_\{i\}\\mathbf\{v\}\_\{i\},𝝎H\\displaystyle\\bm\{\\omega\}\_\{H\}=\(∑i𝐃irot\)−1​∑i𝐃irot​𝝎i\.\\displaystyle=\\left\(\\sum\_\{i\}\\mathbf\{D\}^\{\\mathrm\{rot\}\}\_\{i\}\\right\)^\{\-1\}\\sum\_\{i\}\\mathbf\{D\}^\{\\mathrm\{rot\}\}\_\{i\}\\bm\{\\omega\}\_\{i\}\.Thus, the hub position, velocity, and spin are weighted by the corresponding nodal response operators\. Because the operators are symmetric positive definite, their sums are positive definite and the inverses in Equation \([21](https://arxiv.org/html/2609.30344#S4.E21)\) are well defined\. With identity operators at the first substep, these weighted states reduce to arithmetic averages\.

The rank\-one motivation for this operator\-weighted star construction is derived in Supplementary Information, Sections 9\.2\.1–9\.2\.2\.

### 4\.6Semi\-implicit nodal Newmark update

At each substep, the quantities defined in Sections[4\.4\.1](https://arxiv.org/html/2609.30344#S4.SS4.SSS1)–[4\.4\.3](https://arxiv.org/html/2609.30344#S4.SS4.SSS3)are combined at each physical node\. Using the same incident\-edge set𝒩⁡\(i\)\\mathcal\{N\}\(i\)as in Equation \([18](https://arxiv.org/html/2609.30344#S4.E18)\), the translational force drive and rotational torque drive are

𝐛i=𝐟iext\+∑j∈𝒩⁡\(i\)𝐟i​j,𝝉i=∑j∈𝒩⁡\(i\)𝝉i​j\.\\mathbf\{b\}\_\{i\}=\\mathbf\{f\}^\{\\mathrm\{ext\}\}\_\{i\}\+\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\mathbf\{f\}\_\{ij\},\\qquad\\bm\{\\tau\}\_\{i\}=\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\bm\{\\tau\}\_\{ij\}\.\(22\)Here𝐟iext\\mathbf\{f\}^\{\\mathrm\{ext\}\}\_\{i\}is the observed or decoded external force of Section[4\.4\.3](https://arxiv.org/html/2609.30344#S4.SS4.SSS3), while𝐟i​j\\mathbf\{f\}\_\{ij\}and𝝉i​j\\bm\{\\tau\}\_\{ij\}are the interaction force and spin\-torque contributions of Section[4\.4\.1](https://arxiv.org/html/2609.30344#S4.SS4.SSS1)\. For a virtual edge, the projected force and angular\-momentum flux of Equation \([13](https://arxiv.org/html/2609.30344#S4.E13)\) are used to form these contributions\.

Letttandt\+1t\+1denote the beginning and end of the current substep, and define

Δ​𝐯i=𝐯i,t\+1−𝐯i,t,Δ​𝐱i=𝐫i,t\+1−𝐫i,t\.\\Delta\\mathbf\{v\}\_\{i\}=\\mathbf\{v\}\_\{i,t\+1\}\-\\mathbf\{v\}\_\{i,t\},\\qquad\\Delta\\mathbf\{x\}\_\{i\}=\\mathbf\{r\}\_\{i,t\+1\}\-\\mathbf\{r\}\_\{i,t\}\.\(23\)The aggregated position\- and velocity\-response operators of Equation \([18](https://arxiv.org/html/2609.30344#S4.E18)\) define the change in translational mechanical response over the substep as

Δ​𝐟iresp=−𝐃i​Δ​𝐯i−𝐊i​Δ​𝐱i,\\Delta\\mathbf\{f\}^\{\\mathrm\{resp\}\}\_\{i\}=\-\\mathbf\{D\}\_\{i\}\\Delta\\mathbf\{v\}\_\{i\}\-\\mathbf\{K\}\_\{i\}\\Delta\\mathbf\{x\}\_\{i\},\(24\)whereΔ​𝐟iresp\\Delta\\mathbf\{f\}^\{\\mathrm\{resp\}\}\_\{i\}denotes the response associated with the state increment\. Equation \([24](https://arxiv.org/html/2609.30344#S4.E24)\) specifies how the learned operators enter the update; it does not identify them as derivatives of the decoded force\.

Using the learned inverse mass𝐌i−1\\mathbf\{M\}^\{\-1\}\_\{i\}of Equation \([19](https://arxiv.org/html/2609.30344#S4.E19)\), the trapezoidal response balance is

Δ​𝐯i=𝐌i−1​𝐛i​δ​t\+12​δ​t​𝐌i−1​Δ​𝐟iresp\.\\Delta\\mathbf\{v\}\_\{i\}=\\mathbf\{M\}^\{\-1\}\_\{i\}\\mathbf\{b\}\_\{i\}\\,\\delta t\+\\frac\{1\}\{2\}\\delta t\\,\\mathbf\{M\}^\{\-1\}\_\{i\}\\Delta\\mathbf\{f\}^\{\\mathrm\{resp\}\}\_\{i\}\.\(25\)For the average\-acceleration Newmark parametersγ=12\\gamma=\\tfrac\{1\}\{2\}andβ=14\\beta=\\tfrac\{1\}\{4\}, the corresponding position increment is

Δ​𝐱i=δ​t2​\(𝐯i,t\+𝐯i,t\+1\)=δ​t​\(𝐯i,t\+12​Δ​𝐯i\)\.\\Delta\\mathbf\{x\}\_\{i\}=\\frac\{\\delta t\}\{2\}\\left\(\\mathbf\{v\}\_\{i,t\}\+\\mathbf\{v\}\_\{i,t\+1\}\\right\)=\\delta t\\left\(\\mathbf\{v\}\_\{i,t\}\+\\frac\{1\}\{2\}\\Delta\\mathbf\{v\}\_\{i\}\\right\)\.\(26\)Substituting Equations \([24](https://arxiv.org/html/2609.30344#S4.E24)\) and \([26](https://arxiv.org/html/2609.30344#S4.E26)\) into Equation \([25](https://arxiv.org/html/2609.30344#S4.E25)\) and collecting the unknown velocity increment gives the translational system

\[𝐈3\+12​δ​t​𝐌i−1​𝐃i\+14​δ​t2​𝐌i−1​𝐊i\]​Δ​𝐯i=𝐌i−1​𝐛i​δ​t−12​δ​t2​𝐌i−1​𝐊i​𝐯i,t\.\\boxed\{\\left\[\\mathbf\{I\}\_\{3\}\+\\frac\{1\}\{2\}\\delta t\\,\\mathbf\{M\}^\{\-1\}\_\{i\}\\mathbf\{D\}\_\{i\}\+\\frac\{1\}\{4\}\\delta t^\{2\}\\,\\mathbf\{M\}^\{\-1\}\_\{i\}\\mathbf\{K\}\_\{i\}\\right\]\\Delta\\mathbf\{v\}\_\{i\}=\\mathbf\{M\}^\{\-1\}\_\{i\}\\mathbf\{b\}\_\{i\}\\,\\delta t\-\\frac\{1\}\{2\}\\delta t^\{2\}\\mathbf\{M\}^\{\-1\}\_\{i\}\\mathbf\{K\}\_\{i\}\\mathbf\{v\}\_\{i,t\}\.\}\(27\)Here𝐈3\\mathbf\{I\}\_\{3\}is the3×33\\times 3identity matrix\. The translational state is then advanced as𝐯i,t\+1=𝐯i,t\+Δ​𝐯i\\mathbf\{v\}\_\{i,t\+1\}=\\mathbf\{v\}\_\{i,t\}\+\\Delta\\mathbf\{v\}\_\{i\}and𝐫i,t\+1=𝐫i,t\+Δ​𝐱i\\mathbf\{r\}\_\{i,t\+1\}=\\mathbf\{r\}\_\{i,t\}\+\\Delta\\mathbf\{x\}\_\{i\}\.

The rotational channel follows the same construction\. DefiningΔ​𝝎i=𝝎i,t\+1−𝝎i,t\\Delta\\bm\{\\omega\}\_\{i\}=\\bm\{\\omega\}\_\{i,t\+1\}\-\\bm\{\\omega\}\_\{i,t\}, and using the inverse inertia𝐈i−1\\mathbf\{I\}^\{\-1\}\_\{i\}from Equation \([19](https://arxiv.org/html/2609.30344#S4.E19)\), the aggregated torque𝝉i\\bm\{\\tau\}\_\{i\}from Equation \([22](https://arxiv.org/html/2609.30344#S4.E22)\), and the rotational response operators from Equation \([18](https://arxiv.org/html/2609.30344#S4.E18)\), gives

\[𝐈3\+12​δ​t​𝐈i−1​𝐃irot\+14​δ​t2​𝐈i−1​𝐊irot\]​Δ​𝝎i=𝐈i−1​𝝉i​δ​t−12​δ​t2​𝐈i−1​𝐊irot​𝝎i,t\.\\boxed\{\\left\[\\mathbf\{I\}\_\{3\}\+\\frac\{1\}\{2\}\\delta t\\,\\mathbf\{I\}^\{\-1\}\_\{i\}\\mathbf\{D\}^\{\\mathrm\{rot\}\}\_\{i\}\+\\frac\{1\}\{4\}\\delta t^\{2\}\\,\\mathbf\{I\}^\{\-1\}\_\{i\}\\mathbf\{K\}^\{\\mathrm\{rot\}\}\_\{i\}\\right\]\\Delta\\bm\{\\omega\}\_\{i\}=\\mathbf\{I\}^\{\-1\}\_\{i\}\\bm\{\\tau\}\_\{i\}\\,\\delta t\-\\frac\{1\}\{2\}\\delta t^\{2\}\\mathbf\{I\}^\{\-1\}\_\{i\}\\mathbf\{K\}^\{\\mathrm\{rot\}\}\_\{i\}\\bm\{\\omega\}\_\{i,t\}\.\}\(28\)The spin is then advanced as𝝎i,t\+1=𝝎i,t\+Δ​𝝎i\\bm\{\\omega\}\_\{i,t\+1\}=\\bm\{\\omega\}\_\{i,t\}\+\\Delta\\bm\{\\omega\}\_\{i\}\. The spin remains a latent angular\-momentum state and receives no direct supervision\.

Equations \([27](https://arxiv.org/html/2609.30344#S4.E27)\) and \([28](https://arxiv.org/html/2609.30344#S4.E28)\) are solved independently at every free physical node as two3×33\\times 3systems; the translational solve is illustrated in Fig\.[1](https://arxiv.org/html/2609.30344#S1.F1)e\. Because the aggregated response operators are symmetric positive definite and the learned inverse mass and inverse inertia are positive isotropic operators, the coefficient matrices are symmetric positive definite with minimum eigenvalue at least one\. Each local solve is therefore nonsingular for everyδ​t\\delta t\(Supplementary Information, Section 10\.3\)\.

We refer to the update as*semi\-implicit*because the unknown state increments enter the mechanical response being solved at each node, while the decoded drives and response operators are evaluated from the current graph state and held fixed during that solve\. Unlike a classical global implicit Newmark step,Newmark\-β\\beta\-DGNperforms one local solve per substep, without an inner Newton iteration or an assembled3​N×3​N3N\\times 3Nsystem\. The connection to the average\-acceleration Newmark method and the scope of its classical stability properties are given in Supplementary Information, Sections 10\.2–10\.4\.

The physical\-edge interactions satisfy the pairwise balance of Section[4\.4\.1](https://arxiv.org/html/2609.30344#S4.SS4.SSS1), while the virtual\-edge interactions satisfy the collective balance of Equation \([13](https://arxiv.org/html/2609.30344#S4.E13)\), before integration\. Because these balanced drives are subsequently processed by different nodal response matrices, the independent local solves do not preserve total linear momentum exactly at finiteδ​t\\delta t\. Supplementary Information, Section 11 derives the resulting one\-substep residual and shows that its absolute magnitude is𝒪⁡\(δ​t2\)\\mathcal\{O\}\(\\delta t^\{2\}\)under the stated boundedness assumptions\.

### 4\.7Training and mechanical readouts

Newmark\-β\\beta\-DGNis trained solely from observed kinematic transitions, with boundary and load information included only when available\. Before encoding, vector inputs are normalized using statistics computed from the training partition\. The target position and velocity increments are likewise standardized using training\-set statistics\.

For one observed transition, the training objective is

ℒ=1\|𝒱\|​∑i∈𝒱\(‖Δ​𝐱^i−Δ​𝐱i‖2\+‖Δ​𝐯^i−Δ​𝐯i‖2\),\\mathcal\{L\}=\\frac\{1\}\{\|\\mathcal\{V\}\|\}\\sum\_\{i\\in\\mathcal\{V\}\}\\left\(\\left\\\|\\widehat\{\\Delta\\mathbf\{x\}\}\_\{i\}\-\\Delta\\mathbf\{x\}\_\{i\}\\right\\\|^\{2\}\+\\left\\\|\\widehat\{\\Delta\\mathbf\{v\}\}\_\{i\}\-\\Delta\\mathbf\{v\}\_\{i\}\\right\\\|^\{2\}\\right\),\(29\)where the increments in Equation \([29](https://arxiv.org/html/2609.30344#S4.E29)\) denote their standardised values and𝒱\\mathcal\{V\}contains only physical nodes\. The hub is excluded from the objective\.

No internal force, moment, stress, ground reaction force, stiffness or damping operator, finite\-element matrix, or constitutive residual enters the loss\. The latent spin, inverse mass, inverse inertia, momentum fluxes, and response operators are therefore learned through their role in generating the observed state transition\. Training uses Adam with early stopping on the validation objective, and the checkpoint with the lowest validation loss is retained\. The learning rate, weight decay, batch size, training budget, and other case\-specific settings are reported in Supplementary Information, Sections 4\.3–7\.3; the settings shared across all four implementations are listed in Supplementary Table 1\.

The quantities used to generate the transition remain accessible after training as mechanical readouts \(Fig\.[1](https://arxiv.org/html/2609.30344#S1.F1)f\)\. The decoded physical\-edge interactions provide internal\-force and spin\-torque contributions\. For the walking\-biomechanics experiment, the moment about a joint centre is assembled from the decoded internal forces acting on the distal markers,

𝐌=−∑i∈distal𝐫i×𝐅i,\\mathbf\{M\}=\-\\sum\_\{i\\in\\mathrm\{distal\}\}\\mathbf\{r\}\_\{i\}\\times\\mathbf\{F\}\_\{i\},\(30\)where𝐫i\\mathbf\{r\}\_\{i\}is measured from the joint center\. These inferred moments are compared only at evaluation time with the inverse\-dynamics references withheld from training\.

The learned position\- and velocity\-response operators are also retained directly\. They are interpreted as relative indicators of the position\- and velocity\-dependent response encoded by the coarse\-step update, rather than as identified material or constitutive tangents\. Their relation to the finite\-element tangent structure is evaluated independently in Supplementary Information, Section 4\.9\.

### 4\.8Datasets and evaluation protocol

The clamped\-beam trajectories are generated with FEniCS\[[18](https://arxiv.org/html/2609.30344#bib.bib8)\], adapting the elastodynamics formulation of Bleyer\[[19](https://arxiv.org/html/2609.30344#bib.bib9)\]\. The reference trajectories use the average\-acceleration Newmark parametersγ=12\\gamma=\\tfrac\{1\}\{2\}andβ=14\\beta=\\tfrac\{1\}\{4\}; Supplementary Information, Section 4\.6 additionally compares the trained model with a temporally converged reference\.

The human\-motion benchmark uses the walking sequences of subject 35 from the CMU motion\-capture database\[[23](https://arxiv.org/html/2609.30344#bib.bib11)\], following the prediction span and data split used by EGHN\[[4](https://arxiv.org/html/2609.30344#bib.bib4)\]\. The protein benchmark uses the driven closed–open transition ensemble of adenylate kinase protein\[[24](https://arxiv.org/html/2609.30344#bib.bib12)\]\. The walking\-biomechanics benchmark is constructed from the instrumented\-treadmill recordings of van der Zee et al\.\[[17](https://arxiv.org/html/2609.30344#bib.bib10)\]\. Dataset construction, graph inputs, prediction spans, train–validation–test splits, and case\-specific metrics are given in Supplementary Information, Sections 4–7\.

Evaluation is autoregressive: each predicted state becomes the input to the next prediction interval\. Errors are evaluated without clipping\. The beam uses whole\-body position error relative to beam length together with the deformation measure defined in Supplementary Information, Section 4\.5; human motion and protein dynamics use position mean\-squared error; and the walking\-biomechanics experiment evaluates the inferred joint\-moment demand against the withheld inverse\-dynamics reference\. The common preprocessing and rollout protocol is given in Supplementary Information, Section 3\.

### 4\.9Baseline implementations

We compareNewmark\-β\\beta\-DGNwithDGN\[[15](https://arxiv.org/html/2609.30344#bib.bib1)\], GNS\[[1](https://arxiv.org/html/2609.30344#bib.bib2)\], MGN\[[2](https://arxiv.org/html/2609.30344#bib.bib3)\], EGHN\[[4](https://arxiv.org/html/2609.30344#bib.bib4)\], EGNO\[[5](https://arxiv.org/html/2609.30344#bib.bib5)\], and EGHNO, and additionally with IGNS\[[12](https://arxiv.org/html/2609.30344#bib.bib6)\]on the clamped beam\. Comparisons are controlled within each benchmark: the models use the same trajectories, data partitions, and observed physical state, while retaining their stated architectures and training objectives unless an adaptation is required for autoregressive evaluation\.

DGNis the closest controlled comparison for the proposed update because it uses the same edge\-local conserved interaction representation and rotational channel but integrates the decoded interactions explicitly\. For EGHN, EGHNO, and EGNO, equivariant velocity heads are added because their published architectures do not return the velocity required as input to the next rollout step\. IGNS is evaluated only on the beam using its released port\-Hamiltonian architecture and symplectic\-Euler integrator; stress supervision from its published solid\-mechanics experiment is removed so that it receives the same kinematic supervision as the other beam models\. The baseline objectives and shared comparison protocol are documented in Supplementary Information, Sections 2–3, while case\-specific adaptations, hyperparameters, and evaluation settings are reported in Sections 4\.2–7\.3\.

## Acknowledgments

This research was funded by the Swiss National Science Foundation \(SNSF\) Grant Number200021−200461200021\-200461\.

## References

- \[1\]A\. Sanchez\-Gonzalez, J\. Godwin, T\. Pfaff, R\. Ying, J\. Leskovec, and P\. W\. Battaglia\(2020\)Learning to simulate complex physics with graph networks\.InProceedings of the 37th International Conference on Machine Learning \(ICML\),Proceedings of Machine Learning Research, Vol\.119,pp\. 8459–8468\.Cited by:[§1](https://arxiv.org/html/2609.30344#S1.p3.1),[§2\.1](https://arxiv.org/html/2609.30344#S2.SS1.p3.1),[§2\.3](https://arxiv.org/html/2609.30344#S2.SS3.p2.1),[§2](https://arxiv.org/html/2609.30344#S2a.p1.1),[§4\.9](https://arxiv.org/html/2609.30344#S4.SS9.p1.1)\.
- \[2\]T\. Pfaff, M\. Fortunato, A\. Sanchez\-Gonzalez, and P\. W\. Battaglia\(2021\)Learning mesh\-based simulation with graph networks\.InInternational Conference on Learning Representations \(ICLR\),Cited by:[§1](https://arxiv.org/html/2609.30344#S1.p3.1),[§2\.1](https://arxiv.org/html/2609.30344#S2.SS1.p3.1),[§2](https://arxiv.org/html/2609.30344#S2a.p1.1),[§4\.9](https://arxiv.org/html/2609.30344#S4.SS9.p1.1)\.
- \[3\]V\. G\. Satorras, E\. Hoogeboom, and M\. Welling\(2021\)E\(n\) equivariant graph neural networks\.InProceedings of the 38th International Conference on Machine Learning \(ICML\),Proceedings of Machine Learning Research, Vol\.139,pp\. 9323–9332\.Cited by:[§1](https://arxiv.org/html/2609.30344#S1.p3.1)\.
- \[4\]J\. Han, W\. Huang, T\. Xu, and Y\. Rong\(2022\)Equivariant graph hierarchy\-based neural networks\.InAdvances in Neural Information Processing Systems 35 \(NeurIPS\),Cited by:[§1](https://arxiv.org/html/2609.30344#S1.p14.1),[§1](https://arxiv.org/html/2609.30344#S1.p3.1),[§2\.1](https://arxiv.org/html/2609.30344#S2.SS1.p3.1),[§2\.3](https://arxiv.org/html/2609.30344#S2.SS3.p1.1),[§2\.3](https://arxiv.org/html/2609.30344#S2.SS3.p2.1),[§2](https://arxiv.org/html/2609.30344#S2a.p1.1),[§4\.8](https://arxiv.org/html/2609.30344#S4.SS8.p2.1),[§4\.9](https://arxiv.org/html/2609.30344#S4.SS9.p1.1),[§5\.1](https://arxiv.org/html/2609.30344#S5.SS1.p1.1),[§6\.5](https://arxiv.org/html/2609.30344#S6.SS5.p1.1.1)\.
- \[5\]M\. Xu, J\. Han, A\. Lou, J\. Kossaifi, A\. Ramanathan, K\. Azizzadenesheli, J\. Leskovec, S\. Ermon, and A\. Anandkumar\(2024\)Equivariant graph neural operator for modeling 3D dynamics\.InProceedings of the 41st International Conference on Machine Learning \(ICML\),Proceedings of Machine Learning Research, Vol\.235,pp\. 55015–55032\.Cited by:[§1](https://arxiv.org/html/2609.30344#S1.p3.1),[§2\.1](https://arxiv.org/html/2609.30344#S2.SS1.p3.1),[§2\.3](https://arxiv.org/html/2609.30344#S2.SS3.p2.1),[§2](https://arxiv.org/html/2609.30344#S2a.p1.1),[§4\.9](https://arxiv.org/html/2609.30344#S4.SS9.p1.1),[§6\.5](https://arxiv.org/html/2609.30344#S6.SS5.p1.1.1),[§6\.5](https://arxiv.org/html/2609.30344#S6.SS5.p4.1.1.1)\.
- \[6\]S\. Batzner, A\. Musaelian, L\. Sun, M\. Geiger, J\. P\. Mailoa, M\. Kornbluth, N\. Molinari, T\. E\. Smidt, and B\. Kozinsky\(2022\)E\(3\)\-equivariant graph neural networks for data\-efficient and accurate interatomic potentials\.Nature Communications13,pp\. 2453\.External Links:[Document](https://dx.doi.org/10.1038/s41467-022-29939-5)Cited by:[§1](https://arxiv.org/html/2609.30344#S1.p4.1)\.
- \[7\]I\. Batatia, D\. P\. Kovács, G\. N\. C\. Simm, C\. Ortner, and G\. Csányi\(2022\)MACE: higher order equivariant message passing neural networks for fast and accurate force fields\.InAdvances in Neural Information Processing Systems 35 \(NeurIPS\),Cited by:[§1](https://arxiv.org/html/2609.30344#S1.p4.1)\.
- \[8\]T\. Kipf, E\. Fetaya, K\. Wang, M\. Welling, and R\. Zemel\(2018\)Neural relational inference for interacting systems\.InProceedings of the 35th International Conference on Machine Learning \(ICML\),Proceedings of Machine Learning Research, Vol\.80,pp\. 2688–2697\.Cited by:[§1](https://arxiv.org/html/2609.30344#S1.p5.1)\.
- \[9\]Z\. Han, O\. Fink, and D\. S\. Kammer\(2024\)Collective relational inference for learning heterogeneous interactions\.Nature Communications15,pp\. 3191\.External Links:[Document](https://dx.doi.org/10.1038/s41467-024-47098-7)Cited by:[§1](https://arxiv.org/html/2609.30344#S1.p5.1)\.
- \[10\]Z\. Han, D\. S\. Kammer, and O\. Fink\(2022\)Learning physics\-consistent particle interactions\.PNAS Nexus1\(5\),pp\. pgac264\.External Links:[Document](https://dx.doi.org/10.1093/pnasnexus/pgac264)Cited by:[§1](https://arxiv.org/html/2609.30344#S1.p5.1)\.
- \[11\]S\. A\. Desai, M\. Mattheakis, D\. Sondak, P\. Protopapas, and S\. J\. Roberts\(2021\)Port\-hamiltonian neural networks for learning explicit time\-dependent dynamical systems\.Physical Review E104,pp\. 034312\.External Links:[Document](https://dx.doi.org/10.1103/PhysRevE.104.034312)Cited by:[§1](https://arxiv.org/html/2609.30344#S1.p6.1)\.
- \[12\]T\. Hoang, A\. Trenta, A\. Gravina, N\. Freymuth, P\. Becker, D\. Bacciu, and G\. Neumann\(2025\)Improving long\-range interactions in graph neural simulators via hamiltonian dynamics\.Note:Accepted at ICLR 2026External Links:2511\.08185Cited by:[§1](https://arxiv.org/html/2609.30344#S1.p6.1),[§2\.1](https://arxiv.org/html/2609.30344#S2.SS1.p3.1),[§2](https://arxiv.org/html/2609.30344#S2a.p1.1),[§4\.9](https://arxiv.org/html/2609.30344#S4.SS9.p1.1)\.
- \[13\]V\. Sharma, R\. T\. Oddon, P\. Tesini, J\. Ravesloot, C\. Taal, and O\. Fink\(2025\)Equi\-euler graphnet: an equivariant, temporal\-dynamics informed graph neural network for dual force and trajectory prediction in multi\-body systems\.Mechanical Systems and Signal Processing241,pp\. 113533\.External Links:[Document](https://dx.doi.org/10.1016/j.ymssp.2025.113533)Cited by:[§1](https://arxiv.org/html/2609.30344#S1.p6.1)\.
- \[14\]Z\. Liang, J\. Li, L\. Deng, and X\. Kong\(2026\)PeTIGN: a physics\-encoded graph network for dynamic response computation of large\-scale structures\.Engineering Structures364,pp\. 123094\.External Links:[Document](https://dx.doi.org/10.1016/j.engstruct.2026.123094)Cited by:[§1](https://arxiv.org/html/2609.30344#S1.p6.1)\.
- \[15\]V\. Sharma and O\. Fink\(2026\)Dynami\-CAL GraphNet: a physics\-informed graph neural network conserving linear and angular momentum for dynamical systems\.Nature Communications17,pp\. 1045\.External Links:[Document](https://dx.doi.org/10.1038/s41467-025-67802-5)Cited by:[§1](https://arxiv.org/html/2609.30344#S1.p8.1),[§2\.1](https://arxiv.org/html/2609.30344#S2.SS1.p3.1),[§2\.3](https://arxiv.org/html/2609.30344#S2.SS3.p2.1),[§2\.6](https://arxiv.org/html/2609.30344#S2.SS6.SSS0.Px1.p1.1),[§2](https://arxiv.org/html/2609.30344#S2a.p1.1),[§4\.3](https://arxiv.org/html/2609.30344#S4.SS3.p1.1),[§4\.3](https://arxiv.org/html/2609.30344#S4.SS3.p8.1),[§4\.3](https://arxiv.org/html/2609.30344#S4.SS3.p9.1),[§4\.4\.1](https://arxiv.org/html/2609.30344#S4.SS4.SSS1.p1.1),[§4\.4\.1](https://arxiv.org/html/2609.30344#S4.SS4.SSS1.p2.3),[§4\.9](https://arxiv.org/html/2609.30344#S4.SS9.p1.1),[§8\.1](https://arxiv.org/html/2609.30344#S8.SS1.p1.1.1.2)\.
- \[16\]N\. M\. Newmark\(1959\)A method of computation for structural dynamics\.Journal of the Engineering Mechanics Division85\(3\),pp\. 67–94\.Note:Issue EM3Cited by:[§1](https://arxiv.org/html/2609.30344#S1.p10.1),[§1](https://arxiv.org/html/2609.30344#S1.p11.1)\.
- \[17\]T\. J\. van der Zee, E\. M\. Mundinger, and A\. D\. Kuo\(2022\)A biomechanics dataset of healthy human walking at various speeds, step lengths and step widths\.Scientific Data9,pp\. 704\.External Links:[Document](https://dx.doi.org/10.1038/s41597-022-01817-1)Cited by:[§2\.1](https://arxiv.org/html/2609.30344#S2.SS1.p1.1),[§2\.5](https://arxiv.org/html/2609.30344#S2.SS5.p1.1),[§4\.8](https://arxiv.org/html/2609.30344#S4.SS8.p2.1),[§7\.1](https://arxiv.org/html/2609.30344#S7.SS1.p1.1)\.
- \[18\]M\. S\. Alnæs, J\. Blechta, J\. Hake, A\. Johansson, B\. Kehlet, A\. Logg, C\. Richardson, J\. Ring, M\. E\. Rognes, and G\. N\. Wells\(2015\)The FEniCS project version 1\.5\.Archive of Numerical Software3\(100\),pp\. 9–23\.External Links:[Document](https://dx.doi.org/10.11588/ans.2015.100.20553)Cited by:[§2\.2](https://arxiv.org/html/2609.30344#S2.SS2.p1.1),[§4\.1](https://arxiv.org/html/2609.30344#S4.SS1a.p1.1),[§4\.8](https://arxiv.org/html/2609.30344#S4.SS8.p1.1)\.
- \[19\]J\. Bleyer\(2018\)Numerical tours of computational mechanics with FEniCS\.Zenodo\.External Links:[Document](https://dx.doi.org/10.5281/zenodo.1287832),[Link](https://doi.org/10.5281/zenodo.1287832)Cited by:[§2\.2](https://arxiv.org/html/2609.30344#S2.SS2.p1.1),[§4\.1](https://arxiv.org/html/2609.30344#S4.SS1a.p1.1),[§4\.8](https://arxiv.org/html/2609.30344#S4.SS8.p1.1)\.
- \[20\]T\. J\. R\. Hughes\(2000\)The finite element method: linear static and dynamic finite element analysis\.Dover Publications,Mineola, NY\.Note:Reprint of the 1987 Prentice\-Hall editionCited by:[§2\.2](https://arxiv.org/html/2609.30344#S2.SS2.p11.1.1)\.
- \[21\]K\. Bathe\(1996\)Finite element procedures\.Prentice Hall,Englewood Cliffs, NJ\.Cited by:[§2\.2](https://arxiv.org/html/2609.30344#S2.SS2.p11.1.1)\.
- \[22\]H\. M\. Hilber, T\. J\. R\. Hughes, and R\. L\. Taylor\(1977\)Improved numerical dissipation for time integration algorithms in structural dynamics\.Earthquake Engineering & Structural Dynamics5\(3\),pp\. 283–292\.External Links:[Document](https://dx.doi.org/10.1002/eqe.4290050306)Cited by:[§2\.2](https://arxiv.org/html/2609.30344#S2.SS2.p11.1.1)\.
- \[23\]CMU Graphics Lab\(2003\)CMU graphics lab motion capture database\.Note:[http://mocap\.cs\.cmu\.edu](http://mocap.cs.cmu.edu/)External Links:[Link](http://mocap.cs.cmu.edu/)Cited by:[§2\.3](https://arxiv.org/html/2609.30344#S2.SS3.p1.1),[§4\.8](https://arxiv.org/html/2609.30344#S4.SS8.p2.1),[§5\.1](https://arxiv.org/html/2609.30344#S5.SS1.p1.1)\.
- \[24\]S\. L\. Seyler and O\. Beckstein\(2014\)Sampling large conformational transitions: adenylate kinase as a testing ground\.Molecular Simulation40\(10\-11\),pp\. 855–877\.External Links:[Document](https://dx.doi.org/10.1080/08927022.2014.919497)Cited by:[§2\.4](https://arxiv.org/html/2609.30344#S2.SS4.p1.1.1),[§4\.8](https://arxiv.org/html/2609.30344#S4.SS8.p2.1),[§6\.1](https://arxiv.org/html/2609.30344#S6.SS1.p1.1.1)\.
- \[25\]N\. Michaud\-Agrawal, E\. J\. Denning, T\. B\. Woolf, and O\. Beckstein\(2011\)MDAnalysis: a toolkit for the analysis of molecular dynamics simulations\.Journal of Computational Chemistry32\(10\),pp\. 2319–2327\.External Links:[Document](https://dx.doi.org/10.1002/jcc.21787)Cited by:[§2\.4](https://arxiv.org/html/2609.30344#S2.SS4.p1.1.1)\.
- \[26\]R\. J\. Gowers, M\. Linke, J\. Barnoud, T\. J\. E\. Reddy, M\. N\. Melo, S\. L\. Seyler, J\. Domański, D\. L\. Dotson, S\. Buchoux, I\. M\. Kenney, and O\. Beckstein\(2016\)MDAnalysis: a python package for the rapid analysis of molecular dynamics simulations\.InProceedings of the 15th Python in Science Conference \(SciPy\),pp\. 98–105\.External Links:[Document](https://dx.doi.org/10.25080/Majora-629e541a-00e)Cited by:[§2\.4](https://arxiv.org/html/2609.30344#S2.SS4.p1.1.1),[§6\.1](https://arxiv.org/html/2609.30344#S6.SS1.p1.1.1)\.

Supplementary Information Learning coarse\-step dynamics and internal mechanical response with graph networks Vinay Sharma1and Olga Fink1 1Intelligent Maintenance and Operations Systems, EPFL, Lausanne, Switzerland

## Contents

## 1Proposed Model: Training and Shared Configuration

This section specifies the objective and implementation settings used byNewmark\-β\\beta\-DGNin all four benchmarks\. Dataset construction, graph inputs, prediction spans and case\-specific hyperparameters are given in the corresponding case sections\.

### 1\.1Model state and update

Newmark\-β\\beta\-DGNrepresents each system as a graph of physical nodes and edges augmented by one virtual hub\. Over one observed transition, it performsSSsemi\-implicit Newmark substeps\. At each substep, the network decodes edge\-resolved exchanges of linear and angular momentum together with symmetric positive\-definite translational and rotationalresponseoperators\. The edge quantities are aggregated at the physical nodes and used in the Newmark update\. The operator\-weighted hub supplies the shared graph state required by this update without constructing a dense inverse\.

The model receives no prescribed mass, damping or stiffness matrix\. Its inverse mass, inverse inertia, fluxes andresponseoperators are learned through the state transition\. Whether an external load is observed or inferred is specified for each benchmark\.

### 1\.2Learning objective

For every system,Newmark\-β\\beta\-DGNis trained to predict the normalised increments of position and velocity over one prediction span,

ℒ=1\|𝒱\|​∑i∈𝒱\(‖Δ​𝐱^i−Δ​𝐱i‖2\+‖Δ​𝐯^i−Δ​𝐯i‖2\),\\mathcal\{L\}=\\frac\{1\}\{\|\\mathcal\{V\}\|\}\\sum\_\{i\\in\\mathcal\{V\}\}\\left\(\\left\\lVert\\hat\{\\Delta\\mathbf\{x\}\}\_\{i\}\-\\Delta\\mathbf\{x\}\_\{i\}\\right\\rVert^\{2\}\+\\left\\lVert\\hat\{\\Delta\\mathbf\{v\}\}\_\{i\}\-\\Delta\\mathbf\{v\}\_\{i\}\\right\\rVert^\{2\}\\right\),\(1\)where each increment is standardised using statistics computed from the training split\. The set𝒱\\mathcal\{V\}contains only physical nodes; the hub is excluded from the loss\.

Equation \([1](https://arxiv.org/html/2609.30344#S1.E1)\) is the complete training objective\. No force, moment, stress or ground\-reaction\-force labels are used\. Theresponseoperators𝐊i​j\\mathbf\{K\}\_\{ij\}and𝐃i​j\\mathbf\{D\}\_\{ij\}, their rotational counterparts, and the angular\-velocity channel are also not supervised directly\.

### 1\.3Optimisation and shared model settings

In every case,Newmark\-β\\beta\-DGNis trained with the Adam optimiser and early stopping on the validation value of Equation \([1](https://arxiv.org/html/2609.30344#S1.E1)\)\. We retain the checkpoint with the lowest validation loss\. Table[1](https://arxiv.org/html/2609.30344#S1.T1)lists the settings shared by all four implementations; each case table reports only the quantities that differ\.

Supplementary Table 1:Settings shared byNewmark\-β\\beta\-DGNacross the four benchmarks\.

## 2Baseline Implementations and Adaptations

The baselines retain the training objectives of their original publications unless an adaptation is stated explicitly\.Dynami\-CAL GraphNet\[[15](https://arxiv.org/html/2609.30344#bib.bib1)\]uses a position–velocity increment loss\. The published GNS\[[1](https://arxiv.org/html/2609.30344#bib.bib2)\]formulation predicts normalised acceleration and adds Gaussian noise to the input velocities; our implementation instead uses Equation \([1](https://arxiv.org/html/2609.30344#S1.E1)\) while retaining this noise injection\. MGN\[[2](https://arxiv.org/html/2609.30344#bib.bib3)\]is trained on normalised acceleration\. EGHN\[[4](https://arxiv.org/html/2609.30344#bib.bib4)\]augments a position–velocity mean\-squared error with a linkage\-prediction term weighted byλlink\\lambda\_\{\\mathrm\{link\}\}\. EGNO\[[5](https://arxiv.org/html/2609.30344#bib.bib5)\]uses a position–velocity loss over the predicted trajectory, and EGHNO\[[5](https://arxiv.org/html/2609.30344#bib.bib5)\]adds the same linkage term\. IGNS\[[12](https://arxiv.org/html/2609.30344#bib.bib6)\]is trained in its published form by matching trajectories over a multi\-step window; the adaptation required for the beam comparison is specified in Section[4\.2](https://arxiv.org/html/2609.30344#S4.SS2a)\.

Section[3](https://arxiv.org/html/2609.30344#S3a)defines the data splits, checkpoint selection and rollout rules shared within each comparison\. The graph, inputs, capacity settings and case\-specific adaptations are reported in the case sections rather than here\.

### 2\.1Velocity heads for autoregressive rollout

Autoregressive rollout requires each model to provide every state variable used as input at the next step\.Newmark\-β\\beta\-DGNandDynami\-CAL GraphNetprovide updated positions and velocities, while GNS and MGN predict the acceleration from which both are advanced using their respective published updates\. The published EGHN, EGNO and EGHNO architectures predict positions but not the velocities required by our rollout state\. We therefore add one velocity head to each architecture\. The added head follows the construction of the corresponding position decoder and reads the representation already produced by the backbone\. The message\-passing backbone, pooling, temporal convolution and training objective remain unchanged\. Table[2](https://arxiv.org/html/2609.30344#S2.T2)summarises the three heads\.

Supplementary Table 2:Velocity heads added to the position\-only baselines\.Each head uses the same network class, hidden width and activation as the corresponding position decoder\.𝐥n​f\\mathbf\{l\}\_\{nf\}and𝐧𝐟\\mathbf\{nf\}are the pooled and unpooled equivariant vectors of the hierarchical block;𝐥X\\mathbf\{l\}\_\{X\}and𝐥V\\mathbf\{l\}\_\{V\}are the pooled coordinate and velocity;hhand𝐥H\\mathbf\{l\}\_\{H\}are the low\- and high\-level invariant features; andhouth\_\{\\mathrm\{out\}\}is the hidden state after temporal convolution\.For EGHN and EGHNO, the velocity and position heads read the same equivariant vectors and invariant scalars but have independent parameters\. For EGNO, the velocity head reads the equivariant backbone velocity, displacement from the initial position and invariant hidden state at each predicted time\. Each head follows the equivariant scalarisation–vectorisation construction of the corresponding position decoder and is therefore rotation\-equivariant and translation\-invariant\.

The position and velocity heads are trained jointly but do not share weights\. In theEGMNimplementation, the velocity decoder receives a shallow copy of its vector\-input list because the decoder appends intermediate outputs to that list at every layer\.

Velocity is decoded independently because it is defined at one\-frame spacing, whereas a position increment spans the complete prediction interval\. Their ratio would therefore be a chord velocity over the prediction span, not the one\-frame velocity consumed by the next rollout step\.

### 2\.2Scope of the adaptations

Only adaptations required to place a published architecture under the stated comparison protocol are introduced\. Baseline processor, pooling and temporal\-operator components are otherwise retained\. A capacity count has different meanings across architectures: it denotes Newmark substeps forNewmark\-β\\beta\-DGN, processor depth forDynami\-CAL GraphNetand MGN, temporal or graph layers for EGNO and IGNS, and clusters for EGHN and EGHNO\. These quantities are therefore reported without treating them as equivalent\. Published results that were not reproduced are identified explicitly in the relevant case\.

## 3Data Processing and Evaluation Protocol

This section defines the conventions shared by all models within a benchmark\. It does not impose a common architecture or training objective\. The objective ofNewmark\-β\\beta\-DGNis given in Section[1\.2](https://arxiv.org/html/2609.30344#S1.SS2), and the objectives and adaptations of the baselines are given in Section[2](https://arxiv.org/html/2609.30344#S2a)\.

### 3\.1Data splits and preprocessing

Within each benchmark, all models use the same trajectories and training, validation and test partitions\. Normalisation statistics are computed from the training partition only\. The construction of each partition and the inputs supplied to each architecture are stated in the corresponding case section\.

Across all four systems, velocities are computed from recorded positions using the causal one\-frame backward difference

𝐯⁡\[t\]=𝐱⁡\[t\]−𝐱⁡\[t−1\]\.\\mathbf\{v\}\[t\]=\\mathbf\{x\}\[t\]\-\\mathbf\{x\}\[t\-1\]\.\(2\)The same position and velocity records are used by every model evaluated within a case\.

### 3\.2Optimisation and checkpoint selection

All trained models use the Adam optimiser and early stopping on their respective validation objectives\. We evaluate the checkpoint attaining the lowest validation loss\. Learning rates, weight decay, batch size, stopping patience and training budget are reported for each case\.

### 3\.3Autoregressive rollout

Rollouts are autoregressive: each call advances the state by one prediction span, and the predicted state becomes the input to the next call\. Errors are evaluated separately at each rollout step\. An aggregate over a complete rollout is reported only where it is defined explicitly in the corresponding case\.

Errors are not clipped\. If any seed produces a non\-finite value at a given step, that step is reported as non\-finite rather than averaged over the remaining seeds and is marked rather than plotted\.

### 3\.4Comparison protocol

Comparisons are controlled within each benchmark, not across benchmarks\. Models use the same trajectories, partitions and observed physical state, while retaining their stated objectives and update rules\. Any additional input required by a published architecture, any withheld mechanical label and every implementation adaptation are identified in the corresponding case section\. Seed sets, error definitions and whether a result is reproduced or quoted from the original publication are also stated there\.

## 4Elastodynamics of a Clamped Beam

### 4\.1Case description

The system is a three\-dimensional cantilever beam of lengthLLand rectangular cross\-sectionW×DW\\times D, clamped at one end\. A transverse load of magnitudeFFis applied from rest and removed at a cut\-off timeTcT\_\{c\}, after which the beam vibrates freely and its oscillation decays through Rayleigh damping\. Reference trajectories are produced with the finite\-element solver FEniCS\[[18](https://arxiv.org/html/2609.30344#bib.bib8)\], adapting the elastodynamics demonstration of Bleyer\[[19](https://arxiv.org/html/2609.30344#bib.bib9)\], on an unstructured tetrahedral mesh whose density is set by a resolution parameterres, the number of elements per unit length\.

The material is linear elastic and nondimensional, with Young’s modulusE=1000E=1000, Poisson ratioν=0\.3\\nu=0\.3and densityρ=1\\rho=1\. Damping is of Rayleigh type,𝐂=ηm​𝐌\+ηk​𝐊\\mathbf\{C\}=\\eta\_\{m\}\\mathbf\{M\}\+\\eta\_\{k\}\\mathbf\{K\}withηm=ηk=0\.01\\eta\_\{m\}=\\eta\_\{k\}=0\.01\. Each trajectory spansT=10T=10\\,in100100steps, so one observed transition isΔ​t=0\.1\\Delta t=0\.1\\,\. The finite\-element trajectories are advanced with the average\-acceleration Newmark scheme,γ=12\\gamma=\\tfrac\{1\}\{2\}andβ=14\\beta=\\tfrac\{1\}\{4\}\. These parameters are fixed, not fitted\. Section[4\.6](https://arxiv.org/html/2609.30344#S4.SS6a)separately evaluates the trained model against a temporally converged reference to ensure that its accuracy is not an artefact of sharing the coarse reference discretisation\.

The elastic wave speed isc=E/ρ=31\.6c=\\sqrt\{E/\\rho\}=31\.6, so a wave crosses the base beam inL/c=0\.032L/c=0\.032\\,and one outer step spans roughly three wave transits\. Measured against the largest natural frequency of the model graph, the outer step exceeds the explicit stability limitδ​t​ω≤2\\delta t\\,\\omega\\leq 2by more than an order of magnitude at the outer step, and by 5\.5–11× even when resolved into the four Newmark sub\-steps; the spectral analysis, with the exact factors, is given in Section[4\.7](https://arxiv.org/html/2609.30344#S4.SS7a)\. A coarser element\-size estimate,δ​tcrit≈h/c\\delta t\_\{\\mathrm\{crit\}\}\\approx h/c, gives the same order of magnitude\.

Observed inputs and graph\.Every model receives the current nodal positions and velocities, the mesh graph, the free/clamped boundary mask and the current applied load\. No material property, constitutive parameter, finite\-element matrix, stress, internal force or moment is supplied\. Nodes are the mesh vertices and physical edges follow the mesh connectivity\.Newmark\-β\\beta\-DGNaugments this graph with one hub joined to allNNphysical nodes by2​N2Ndirected virtual edges, identified by an edge attribute of−1\-1; the baselines retain their native hub\-free topology\. Graphs are preprocessed once and are not rebuilt during training\.

Supplementary Table 3:Beam representation supplied toNewmark\-β\\beta\-DGN\.NNis the number of mesh vertices and depends on the configuration\. Baselines receive the same observed physical state on the hub\-free graph\.Splits\.The training distribution is the Cartesian productL∈\{1\.0,1\.5,2\.0\}L\\in\\\{1\.0,1\.5,2\.0\\\},W∈\{0\.5,1\.0\}W\\in\\\{0\.5,1\.0\\\},D∈\{0\.5,1\.0\}D\\in\\\{0\.5,1\.0\\\},F∈\{1\.5,2\.0,2\.5\}F\\in\\\{1\.5,2\.0,2\.5\\\},Tc∈\{2\.0,2\.5,3\.0\}T\_\{c\}\\in\\\{2\.0,2\.5,3\.0\\\}\\,atres=4\\texttt\{res\}=4, giving108108configurations and one trajectory each, divided into training, validation and test partitions\. Nine further configurations are held out entirely and never seen in any form; they vary the load amplitude, the cut\-off time, the cross\-section, the beam length and the mesh resolution, one axis at a time\.

Metric\.At each rollout step, the whole\-body error is the root\-mean\-square position error over the physical nodes, expressed as a percentage of the beam lengthLL\. Where one value summarises a trajectory, we report its mean over the9595autoregressive steps\. The deformation ratio, predicted over reference total displacement, measures over\- or under\-deformation\. Section[4\.5](https://arxiv.org/html/2609.30344#S4.SS5a)explains why this whole\-body measure is preferred to a tip\-only error\. All results use one seed,4242, so no dispersion is reported\.

### 4\.2Evaluated baselines

We compareNewmark\-β\\beta\-DGNwithDynami\-CAL GraphNet, MGN, EGHN, EGHNO, EGNO and IGNS\. All receive the observed inputs listed above, and none receives the finite\-element mass, damping or stiffness matrices\. The baselines use the physical mesh graph without the hub\. EGHN, EGHNO and EGNO use the velocity heads of Section[2\.1](https://arxiv.org/html/2609.30344#S2.SS1a)\.Dynami\-CAL GraphNetis evaluated with eight and twelve message\-passing rounds, and MGN with four, eight and twelve processor rounds, to test whether additional local propagation closes the gap to the hub\-mediated update\. MGN hasno rotational channeland is included only in the rollout comparison\. No global frame is used because the applied load is observed\.

IGNS is evaluated only in this case\. Its published solid\-mechanics experiment uses multi\-step trajectory matching and stress supervision\. Here it receives no stress target and is trained on Equation \([1](https://arxiv.org/html/2609.30344#S1.E1)\), so that the comparison tests the integration mechanism under the same one\-step kinematic supervision asNewmark\-β\\beta\-DGN\. We use the released port\-Hamiltonian architecture and the symplectic\-Euler integrator specified in the paper\. Because the input edge features are standardised and may be signed, the released positive\-distance adjacency weighting is replaced by symmetric degree normalisation; distance remains in the edge features\. All other inputs, normalisation, optimisation and training budgets match the hub\-free beam protocol\. The number of internal integration steps is swept over\{1,4,12,24\}\\\{1,4,12,24\\\}; four matches the substep budget ofNewmark\-β\\beta\-DGN\.

For MGN, we corrected the output normalisation before evaluation; without this correction, its initial loss was of order10510^\{5\}and rollout drift began immediately\. All reported MGN results use the corrected implementation\.

### 4\.3Hyperparameters

Settings forNewmark\-β\\beta\-DGNnot listed in Table[4](https://arxiv.org/html/2609.30344#S4.T4)follow Table[1](https://arxiv.org/html/2609.30344#S1.T1)\.

Supplementary Table 4:Beam elastodynamics: hyperparameters\.Capacity parameters have architecture\-specific meanings, as explained in Section[2\.2](https://arxiv.org/html/2609.30344#S2.SS2a)\.All models train for at most20002000epochs with early stopping, evaluated every second epoch, and the checkpoint minimising the validation loss is retained\. Seed4242throughout\.

### 4\.4Additional results

![Refer to caption](https://arxiv.org/html/2609.30344v1/figures/si/beam_fea/training/beamfea_training_curves_seed42.png)Supplementary Figure 1:Training and validation loss forNewmark\-β\\beta\-DGNon the beam\.Both terms of Equation \([1](https://arxiv.org/html/2609.30344#S1.E1)\) are shown\. Early stopping selects the checkpoint at the minimum of the validation curve\. Seed 42\.![Refer to caption](https://arxiv.org/html/2609.30344v1/figures/si/beam_fea/combined_rollout/L1.5_W0.75_D0.75_F2.0_Tc0.25_res8.png)

![Refer to caption](https://arxiv.org/html/2609.30344v1/figures/si/beam_fea/combined_rollout/L3.0_W0.75_D0.75_F2.0_Tc0.25_res4.png)

![Refer to caption](https://arxiv.org/html/2609.30344v1/figures/si/beam_fea/combined_rollout/L1.75_W1.5_D0.4_F2.0_Tc0.25_res4.png)

![Refer to caption](https://arxiv.org/html/2609.30344v1/figures/si/beam_fea/combined_rollout/L1.75_W0.4_D1.5_F2.0_Tc0.25_res4.png)

Supplementary Figure 2:Rollout ofNewmark\-β\\beta\-DGNon four held\-out configurations\.Top left, the refined mesh \(res=8\\texttt\{res\}=8\), the only extrapolation that leaves the training envelope of the stiffness spectral radius\. Top right, the beam extrapolated to three times the training length, where the rank\-one hub begins to be an inadequate surrogate for the long\-range coupling\. Bottom, the two cross\-sections outside the training product,W=1\.5W=1\.5,D=0\.4D=0\.4on the left andW=0\.4W=0\.4,D=1\.5D=1\.5on the right, which exchange the bending stiffness of the two transverse axes\. The remaining five configurations behave comparably and are not reproduced here\.
### 4\.5Whole\-body error and the role of the semi\-implicit solve

The metric of a single point like tip can under\-report the errors since a bounded tip motion can coexist with an interior that over\-deforms: a mesh that swells or buckles between the clamp and the tip can still return a small tip error\. We therefore also report the*whole\-body*error, the root\-mean\-square position error over all body nodes as a percentage of the beam length, together with a deformation ratio, the predicted total displacement divided by the reference total displacement, whose departure from unity measures over\- or under\-deformation\. Both are accumulated over the same9595steps and evaluated on the same seed\.

Table[5](https://arxiv.org/html/2609.30344#S4.T5)decomposes the model into the ingredients that distinguish it fromDynami\-CAL GraphNet: the operator\-weighted hub, the semi\-implicit solve, and the stiffness termβ\\betaof the Newmark left\-hand side\. Adding the hub toDynami\-CAL GraphNet—with a purely explicit per\-node update \(𝐀i=𝐈\\mathbf\{A\}\_\{i\}=\\mathbf\{I\}\) and only four sub\-steps in place of twelve—reduces the mean whole\-body error from3\.00%3\.00\\%to1\.26%1\.26\\%, a factor of2\.42\.4: the hub is the dominant ingredient, improving accuracy even as the sub\-step count is reduced\. Restoring the semi\-implicit solve then refines the result, to0\.70%0\.70\\%atβ=0\\beta=0and0\.54%0\.54\\%with the stiffness term, further factors of1\.81\.8and1\.31\.3\. The gap toDynami\-CAL GraphNetis widest on the refined mesh, whereDynami\-CAL GraphNetreaches1515–16%16\\%whole\-body error \(for each of the two refined\-mesh configurations, all nodes, mean over 95 steps\) and over\-deforms the mesh by factors of2929–8585, against0\.50\.5–1\.2%1\.2\\%and a deformation ratio near two forNewmark\-β\\beta\-DGN\.

Supplementary Table 5:Whole\-body error decomposition \(mean over2121configurations, % ofLL\)\.The hub is the dominant ingredient; the semi\-implicit solve and the stiffness termβ\\betaare refinements\. Cells are cumulative: each row adds one ingredient to the row above\. The first hub row also reduces the sub\-step count from1212to44, so its factor bundles the two changes\. Seed 42\.
### 4\.6Accuracy against a temporally converged reference

Because the finite\-element reference uses the same Newmark parameters and the same step asNewmark\-β\\beta\-DGN, its accuracy could in principle reflect a shared discretisation rather than the physics\. To rule this out, we regenerate the reference with the same solver at a100100times finer step,δ​t=10−3\\delta t=10^\{\-3\}\\,, and subsample it back to the0\.10\.1\\,grid; halving this step again changes the whole\-body trajectory by0\.005%0\.005\\%, soδ​t=10−3\\delta t=10^\{\-3\}\\,is the converged solution\. The mesh is deterministic, so the fine and coarse body nodes coincide exactly\. Table[6](https://arxiv.org/html/2609.30344#S4.T6)evaluates the unchanged seed\-42 checkpoint against this converged reference on eight representative configurations spanning the training distribution, both refined meshes, a longer beam and a cut\-off\-time extrapolation\.

Newmark\-β\\beta\-DGNreproduces the converged solution to between0\.19%0\.19\\%and2\.06%2\.06\\%whole\-body error, essentially identical to its error against the coarse reference \(means0\.72%0\.72\\%versus0\.67%0\.67\\%\)\. The reference’s own coarse\-step discretisation error is only0\.090\.09–0\.58%0\.58\\%, so the comparison is not integrator\-matching; the model is accurate to the physics, not to a co\-discretised reference\. On the stiffest and longest configurations the model is marginally closer to the converged solution than the single\-step coarse reference is, consistent with its four sub\-steps ofδ​t=0\.025\\delta t=0\.025\\,resolving the interval more finely than a single0\.10\.1\\,step\.

Supplementary Table 6:Whole\-body error against a temporally converged reference\(δ​t=10−3\\delta t=10^\{\-3\}\\,, mean over 95 steps, % ofLL\)\. The first two error columns compareNewmark\-β\\beta\-DGNwith the indicated reference; the final column reports the coarse\-reference discretisation error\. Seed 42\.
### 4\.7Spectral radius and the IGNS baseline

The outer step ofΔ​t=0\.1\\Delta t=0\.1\\,lies far beyond the reach of any explicit stability limit\. Linearizing the reference dynamics about its trajectory gives a largest natural frequencyωmax≈442\\omega\_\{\\max\}\\approx 442on the training mesh and≈876\\approx 876on the refined mesh; at the outer stepΔ​t=0\.1\\Delta t=0\.1\\,the productΔ​t​ωmax\\Delta t\\,\\omega\_\{\\max\}is4444and8888, that is2222and4444times the explicit stability limitΔ​t​ω≤2\\Delta t\\,\\omega\\leq 2, and even resolved into the four sub\-steps ofNewmark\-β\\beta\-DGNit remains5\.55\.5and1111times that limit\. The corresponding symplectic\-Euler update therefore amplifies the stiffest mode by a factor of109109per sub\-step on the training mesh and455455on the refined mesh, whereas the corresponding classical undamped average\-acceleration Newmark update has unit amplification \(Table[7](https://arxiv.org/html/2609.30344#S4.T7), and Section[10](https://arxiv.org/html/2609.30344#S10),Equation \([66](https://arxiv.org/html/2609.30344#S10.E66)\)\)\. Mesh refinement is the decisive axis because it is the only held\-out change that raisesωmax\\omega\_\{\\max\}beyond the training range\.

Supplementary Table 7:Spectral stability of the beam operator\.ωmax\\omega\_\{\\max\}is the largest natural frequency of the linearised system; the explicit limit isΔ​t​ω≤2\\Delta t\\,\\omega\\leq 2\. The Newmark amplification is unity at every step size; the symplectic\-Euler amplification is per sub\-step \(S=4S=4\)\.This prediction is borne out by IGNS, evaluated under the protocol of Section[4\.2](https://arxiv.org/html/2609.30344#S4.SS2a)as the method\-specific analogue of the identity\-matrix ablation\. At its native budget of roughly one internal step per frame the rollout diverges to10610^\{6\}–107%10^\{7\}\\%on every configuration\. At four internal steps, the substep budget ofNewmark\-β\\beta\-DGN, it fails on all configurations\. Even at2424internal steps, six times that budget, it reaches231231–338%338\\%whole\-body error on the two refined meshes, against0\.50\.5–1\.5%1\.5\\%forNewmark\-β\\beta\-DGN, and its per\-case error on the refined meshes does not decrease as the budget grows \(Table[8](https://arxiv.org/html/2609.30344#S4.T8)\)\. Its single\-step validation loss meanwhile improves monotonically with the budget, from0\.210\.21at one step to0\.0080\.008at2424, so accurate fitting of one transition does not imply a stable repeated application\. Increasing the number of explicit internal steps therefore does not reproduce the coarse\-step behaviour of the semi\-implicit response\.

Supplementary Table 8:IGNS rollout error under an increasing internal\-step budget\(whole\-body, final step, % ofLL\)\. Four internal steps matches the sub\-step budget ofNewmark\-β\\beta\-DGN; the refined\-mesh error does not fall as the budget grows\.Newmark\-β\\beta\-DGN\(four sub\-steps\) is shown for reference\.The trained one\-step maps confirm this mechanism empirically\. Off the reference trajectory the amplification factor ofNewmark\-β\\beta\-DGNis near unity, between0\.970\.97and1\.081\.08across configurations, a few percent short of the physical Rayleigh contraction rather than expansive, so its rollout stays bounded and non\-monotone rather than growing\.Dynami\-CAL GraphNet’s map is likewise near\-neutral \(0\.9990\.999–1\.0061\.006, including the refined mesh\), which is why it drifts rather than diverges; EGHN’s map instead amplifies, by1\.11\.1–1\.6×1\.6\\timeson the training mesh and by roughly11×11\\timeson the refined mesh, which is why it reaches a non\-finite value within tens of steps\.Newmark\-β\\beta\-DGNis accurate on the trajectory and near\-neutral off it, and that is what keeps its rollout bounded\.

Consistent with the absence of numerical dissipation in the linear analysis, the model’s free\-vibration ring\-down matches the reference’s physical Rayleigh decay to within about ten percent in\-distribution \(decay ratios0\.920\.92–1\.091\.09\) and shifts the vibration frequency by under one percent, rather than damping more strongly as an over\-dissipative scheme would\. On the refined mesh this breaks down: the model over\-damps by roughly2\.5×2\.5\\timesand mis\-predicts the frequency by−27%\-27\\%, a genuine extrapolation limitation of the learned operators rather than of the integrator\.

### 4\.8Modal structure recovered from the rollout

The spectral analysis of Section[4\.7](https://arxiv.org/html/2609.30344#S4.SS7a)concerns the integrator alone: the average\-acceleration Newmark update neither amplifies nor damps a linear mode, whereas the explicit update amplifies the stiff modes it cannot resolve\. A model can fit each single transition accurately and still reproduce the wrong frequency, because frequency error accumulates only over the rollout\. We therefore read the free\-vibration content directly from the autoregressive trajectory, with no frequency or mode supervision\.

For each held\-out configuration, the tip transverse displacement over the9595rollout steps is detrended to remove the slow loading envelope and its power spectrum is taken; the peak is the fundamental frequency\. Independently, a proper\-orthogonal decomposition of the transverse displacement along the beam returns the dominant spatial mode shape, compared with the finite\-element mode by the Modal Assurance Criterion \(MAC\), which equals one for identical shapes\. Figure[3](https://arxiv.org/html/2609.30344#S4.F3)shows both for three representative configurations and Table[9](https://arxiv.org/html/2609.30344#S4.T9)reports them for every configuration\.

Supplementary Figure 3:Fundamental frequency and dominant mode recovered from the rollout,on three held\-out configurations \(top to bottom: a load\-amplitude, a load\-duration and a cross\-section extrapolation, identified in each title\)\. Left, the natural\-frequency spectrum of the tip transverse displacement, finite element \(black\),Newmark\-β\\beta\-DGN\(red\),Dynami\-CAL GraphNet\(blue dashed\), with the fundamental frequencies annotated; right, the dominant proper\-orthogonal mode shape, finite element againstNewmark\-β\\beta\-DGN, with its MAC\.Newmark\-β\\beta\-DGNmatches the finite\-element fundamental and mode;Dynami\-CAL GraphNetpeaks at a spurious lower frequency\. No frequency or mode is supervised\.Supplementary Table 9:Natural frequency and dominant mode recovered from the rollout, every configuration\.FE is the finite\-element fundamental; theNewmark\-β\\beta\-DGNvalue is read from its rollout tip spectrum; MAC compares the dominantNewmark\-β\\beta\-DGNmode shape with the finite\-element mode\. A frequency counts as recovered when it falls in the same spectral bin \(≈0\.1\\approx 0\.1\\,Hz atres=4\\texttt\{res\}=4\)\. Seed 42\.The dominant mode shape is recovered on every configuration, atMAC≥0\.97\\mathrm\{MAC\}\\geq 0\.97\. The fundamental frequency ofNewmark\-β\\beta\-DGNmatches the finite\-element value to the spectral resolution of the record on the loading and thin\-cross\-section extrapolations, and departs from it only on the length and mesh\-refinement extrapolations, the axes that raiseωmax\\omega\_\{\\max\}beyond the training range \(Section[4\.7](https://arxiv.org/html/2609.30344#S4.SS7a)\)\. The explicitDynami\-CAL GraphNetrollout never recovers the fundamental, peaking at a spurious0\.30\.3–0\.40\.4\\,Hz on every configuration, so a bounded rollout does not imply a correct spectrum\. This recovery is the empirical counterpart of the integrator’s unit amplification; it is inherited from the linear analysis and is not proved for the learned nonlinear rollout\.

### 4\.9Learned operators and the finite\-element tangent

We ask how much of the true finite\-element tangent the decoded operators recover, since a readout is trustworthy only to the extent that its content is identifiable from the observed motion\. At each free body node, the position\-response operator𝐊i=∑j𝐊i​j\\mathbf\{K\}\_\{i\}=\\sum\_\{j\}\\mathbf\{K\}\_\{ij\}and velocity\-response operator𝐃i\\mathbf\{D\}\_\{i\}, the3×33\\times 3symmetric positive\-definite blocks used by the learned update, are compared with the corresponding finite\-element matrices: the stiffness𝐊i​i\\mathbf\{K\}\_\{ii\}, which is state\-independent for the linear beam, and the Rayleigh damping𝐂i​i=ηk​𝐊i​i\+ηm​𝐌i​i\\mathbf\{C\}\_\{ii\}=\\eta\_\{k\}\\mathbf\{K\}\_\{ii\}\+\\eta\_\{m\}\\mathbf\{M\}\_\{ii\}\. The comparison is made at three stages of the trajectory \(loaded, free vibration and rest\) for three held\-out configurations\.

A symmetric operator splits orthogonally, in the Frobenius inner product⟨𝐀,𝐁⟩F=∑a​bAa​b​Ba​b\\langle\\mathbf\{A\},\\mathbf\{B\}\\rangle\_\{F\}=\\sum\_\{ab\}A\_\{ab\}B\_\{ab\}, into an isotropic part13​\(tr⁡𝐌\)​𝐈\\tfrac\{1\}\{3\}\(\\operatorname\{tr\}\\mathbf\{M\}\)\\mathbf\{I\}and a trace\-free deviatoric part, so that∥𝐌∥F2=∥𝐌iso∥F2\+∥𝐌dev∥F2\\lVert\\mathbf\{M\}\\rVert\_\{F\}^\{2\}=\\lVert\\mathbf\{M\}\_\{\\mathrm\{iso\}\}\\rVert\_\{F\}^\{2\}\+\\lVert\\mathbf\{M\}\_\{\\mathrm\{dev\}\}\\rVert\_\{F\}^\{2\}; physically the isotropic part is the overall stiffness magnitude, which is direction\-independent, and the deviatoric part is the anisotropy, that is, which directions are stiffer than others\. We separate the two with three measures \(Fig\.[4](https://arxiv.org/html/2609.30344#S4.F4)\)\. The*direction*is the Frobenius cosine between the decoded and exact operators,cosi=⟨𝐊i,𝐊i​i⟩F/\(∥𝐊i∥F∥𝐊i​i∥F\)\\cos\_\{i\}=\\langle\\mathbf\{K\}\_\{i\},\\mathbf\{K\}\_\{ii\}\\rangle\_\{F\}/\(\\lVert\\mathbf\{K\}\_\{i\}\\rVert\_\{F\}\\,\\lVert\\mathbf\{K\}\_\{ii\}\\rVert\_\{F\}\): it is scale\-invariant, equals one when the two are proportional, and because it is dominated by the large isotropic part it reports the mean direction and not the anisotropy\. The*anisotropy*is the same cosine taken on the deviatoric partsdev⁡\(𝐌\)=𝐌−13​\(tr⁡𝐌\)​𝐈\\operatorname\{dev\}\(\\mathbf\{M\}\)=\\mathbf\{M\}\-\\tfrac\{1\}\{3\}\(\\operatorname\{tr\}\\mathbf\{M\}\)\\mathbf\{I\}, which isolates the direction\-dependent content\. The*spatial pattern*is the Pearson correlation across nodes between the tracestr⁡𝐊i\\operatorname\{tr\}\\mathbf\{K\}\_\{i\}andtr⁡𝐊i​i\\operatorname\{tr\}\\mathbf\{K\}\_\{ii\}; the trace is the sum of the eigenvalues, a per\-node stiffness magnitude, so this correlation tests whether the model places stiffness where the finite\-element tangent does, independently of any global scale\. A fourth quantity, the ratio of the median traces, measures the absolute scale\.

Supplementary Figure 4:The learned response operators recover part of the finite\-element tangent structure, but not its anisotropy or magnitude\.Agreement of the decoded per\-node position\-response operator𝐊i\\mathbf\{K\}\_\{i\}\(blue\) and velocity\-response operator𝐃i\\mathbf\{D\}\_\{i\}\(red\) with the corresponding finite\-element matrices under three scale\-separated measures:*direction*is the Frobenius cosine between the operators;*anisotropy*is the cosine between their trace\-free \(deviatoric\) parts;*spatial pattern*is the correlation of their traces across nodes\. Each dot is one of three held\-out configurations at one of three trajectory stages \(loaded, free vibration, rest\); bars are the mean and error bars one standard deviation\. The operators align in direction \(cosine≈0\.85\\approx 0\.85for stiffness,0\.920\.92for damping\) and reproduce where the material is stiff or soft \(trace correlation≈0\.63\\approx 0\.63for stiffness,0\.940\.94for damping\), but do not recover the directional anisotropy \(deviatoric cosine≈0\\approx 0\) and are off in absolute scale by22–9×9\\timesfor stiffness and roughly103×10^\{3\}\\timesfor damping\. They are therefore operational operators, usable as relative readouts, rather than identified material tangents\.The decoded operators align with the corresponding finite\-element matrices in direction \(cosine≈0\.85\\approx 0\.85for stiffness and0\.920\.92for damping\) and reproduce their spatial pattern \(trace correlation≈0\.63\\approx 0\.63and0\.940\.94, respectively\), but they do not recover the directional anisotropy \(the deviatoric cosine is approximately zero\) or the absolute scale, which differs by factors of two to nine from the stiffness and by roughly10310^\{3\}from the damping\. The position\-response operator also changes across the three stages although the true stiffness is constant, so it adapts to the state rather than reproducing a fixed material tangent\. Two further probes locate the operators mechanically\. The position\-response direction is nearly orthogonal to the pointwise Jacobian of the decoded internal force \(cosine−0\.05\-0\.05to\+0\.09\+0\.09\), confirming that it is not the derivative of the force channel\. It nevertheless predicts the change in the model’s momentum response to a velocity perturbation with a cosine of0\.990\.99to1\.001\.00\. The operators therefore describe the response produced by the learned update, not an identified constitutive tangent\.

The recovered content therefore determines what the operators can be used for\. They recover the direction and the spatial organisation of the finite\-element tangent, which a*relative*readout requires, but not its anisotropy or its absolute scale, which a calibrated material tangent would need\.They are therefore usable as a relative, interpretive readout of where the structure is stiff or soft; they cannot be read as a calibrated material or constitutive tangent, as the derivative of the decoded force, or as a predictor of an absolute stress or of the response to a change of material\.

## 5Human Motion Capture

### 5\.1Case description

We use the walking sequences of subject 35 of the CMU motion capture database\[[23](https://arxiv.org/html/2609.30344#bib.bib11)\]and follow the prediction span and data\-split protocol introduced by EGHN\[[4](https://arxiv.org/html/2609.30344#bib.bib4)\], so that our numbers are directly comparable with the values published for that benchmark\. The task is to predict the marker configuration one prediction span,delta\_frame=30\\texttt\{delta\\\_frame\}=30raw frames, ahead of the observed state\.

Observed inputs and graph\.All models receive marker positions, the causal velocities of Equation \([2](https://arxiv.org/html/2609.30344#S3.E2)\) and the anatomical graph\. The3131markers are connected by6060directed skeletal edges\.Newmark\-β\\beta\-DGNaugments this graph with6262directed virtual edges to one hub, giving122122edges in total\. The baseline graph contains two\-hop skeletal shortcuts and has130130directed edges\. The edge budgets are therefore similar, while the topologies differ\.

The external forcing is not observed\.Newmark\-β\\beta\-DGNrepresents it using a parameter\-free spectral frame formed from the two leading non\-trivial eigenvectors of the body\-subgraph Laplacian\. The frame is rotation\-equivariant, translation\-invariant and determined by graph topology rather than pose\.

Supplementary Table 10:Human\-motion representation supplied toNewmark\-β\\beta\-DGN\.Baseline\-specific features are described in Section[5\.2](https://arxiv.org/html/2609.30344#S5.SS2)\.Splits and seeds\.Training, validation and test partitions follow the reference protocol\. Every model is trained from three random seeds,\{0,42,100\}\\\{0,42,100\\\}, and results are reported as the mean and standard deviation over those seeds,n=3n=3\. This is the only case in this work for which a dispersion measure is available\.

Rotation evaluation\.To assess whether a model’s predictions are covariant with a rotation of the observation frame, we rotate the entire test trajectory about the vertical walking axis by a fixed angle of5∘5^\{\\circ\},15∘15^\{\\circ\}or30∘30^\{\\circ\}and re\-run the identical rollout evaluation\. The0∘0^\{\\circ\}condition reproduces the unrotated row exactly and serves as a control\.

### 5\.2Evaluated baselines

We compareNewmark\-β\\beta\-DGNwith GNS, EGHN, EGHNO, EGNO,Dynami\-CAL GraphNetandDynami\-CAL GraphNet\(global reference frame\),which givesDynami\-CAL GraphNetthe global reference frame used byNewmark\-β\\beta\-DGN\. EGHN, EGHNO and EGNO use the velocity heads of Section[2\.1](https://arxiv.org/html/2609.30344#S2.SS1a)\. The corrected rollout variant of EGHNO is used: its velocity is decoded from the hierarchical representation rather than from the low\-level block alone\.

The architectures retain their published input representations\. GNS receives the current and previous velocities at each node; its edges carry relative position, relative velocity and their norms\. These inputs are translation\-invariant, but its Cartesian multilayer perceptrons are not rotation\-equivariant\. EGHN, EGHNO and EGNO receive absolute height as an invariant scalar\.Dynami\-CAL GraphNetconstructs a normalised height internally at every message\-passing round\.Newmark\-β\\beta\-DGNand GNS receive no absolute height\.

Dynami\-CAL GraphNet\(global reference frame\) givesDynami\-CAL GraphNet’s external\-force channel the same parameter\-free spectral frame thatNewmark\-β\\beta\-DGNuses;Dynami\-CAL GraphNetis itself rotation\-equivariant, and the originalDynami\-CAL GraphNetis retained to separate the effect of this frame\.Because the rotation test is performed about the vertical axis, height is unchanged by the applied transformation\.

### 5\.3Hyperparameters

Settings forNewmark\-β\\beta\-DGNnot listed in Table[11](https://arxiv.org/html/2609.30344#S5.T11)follow Table[1](https://arxiv.org/html/2609.30344#S1.T1)\.

Supplementary Table 11:Human motion\-capture hyperparameters\.Capacity parameters have architecture\-specific meanings, as explained in Section[2\.2](https://arxiv.org/html/2609.30344#S2.SS2a)\.All models train for at most10 00010\\,000epochs with early stopping\. Prediction spandelta\_frame=30\\texttt\{delta\\\_frame\}=30frames for every model\.

### 5\.4Additional results

Table[12](https://arxiv.org/html/2609.30344#S5.T12)reports the per\-step rollout error for the seven models, unclipped\.Two behaviours visible in the table are omitted from the main text\. First, the equivariant baselines EGHN, EGHNO and EGNO stay finite on the canonical split but grow large and unstable by the five\-step horizon, reaching mean squared errors of roughly four to six, with EGHNO carried to a large finite value by a single seed; their equivariant construction keeps them bounded but does not control the coarse\-step growth\. Second,Newmark\-β\\beta\-DGN, whose learned response operators enter a semi\-implicit solve, is the most accurate of the equivariant models at every step, whereas the non\-equivariant GNS is more accurate still on the unrotated set but loses that advantage under the rotation of Table[13](https://arxiv.org/html/2609.30344#S5.T13)\.

Supplementary Table 12:Per\-step rollout error on the human\-walk benchmark\.Unscaled position mean squared error, mean±\\pmstandard deviation over seeds\{0,42,100\}\\\{0,42,100\\\},n=3n=3, on the canonical split\. Values are not clipped\. Every model stays finite across the five\-step rollout on this split; EGHNO reaches a large finite value at step five, driven by seed 0 \(1\.7×1071\.7\\times 10^\{7\}\) while its other seeds remain near1010\.Supplementary Table 13:Degradation of GNS under a rotation of the test set\.Per\-step rollout error, mean±\\pmstandard deviation over seeds\{0,42,100\}\\\{0,42,100\\\},n=3n=3\. The0∘0^\{\\circ\}row reproduces the GNS row of Table[12](https://arxiv.org/html/2609.30344#S5.T12)\.Newmark\-β\\beta\-DGNis rotation\-equivariant by construction and its error is unchanged at every angle, so it is not tabulated per angle\.The error of GNS increases sharply with rotation angle and compounds over the rollout horizon, which is the behavior expected of a model that has no built\-in rotational symmetry and was trained on trajectories captured in a single orientation\. At5∘5^\{\\circ\}, a misalignment small enough to arise from the calibration of a capture volume, GNS already predicts less accurately thanNewmark\-β\\beta\-DGNat the five\-step horizon\.

## 6Protein Dynamics in Solvent

### 6\.1Case description

We evaluateNewmark\-β\\beta\-DGNon a large collective conformational change of adenylate kinase: the driven transition ensemble of Seyler and Beckstein\[[24](https://arxiv.org/html/2609.30344#bib.bib12)\], accessed through the MDAnalysis toolkit\[[26](https://arxiv.org/html/2609.30344#bib.bib14)\]\. The ensemble contains200200dynamic\-importance\-sampling paths, each following the855855backbone atoms through the conformational change over9090to106106frames\. Along every path the radius of gyration and the inter\-domain distances increase monotonically, so the trajectories are closed–open transitions in which the two mobile domains swing away from the rigid core\. These are enhanced\-sampling paths rather than unbiased kinetics: one step is a fixed fraction of the transition, not a physical time, and we read the benchmark as reproduction of the collective pathway\.

Unlike the equilibrium trajectory of the same protein \(Section[6\.5](https://arxiv.org/html/2609.30344#S6.SS5)\), which samples confined fluctuation within a single basin and is bounded by a linear floor, this ensemble carries a large directed motion: a root\-mean\-square displacement of about7​Å7\\,\\text\{\\AA\}between the closed and open states\. Its predictable content is nonetheless dominated by a single canonical route\. A state\-independent mean\-trajectory baseline, which moves every test path along the same average closed\-open transition of the training paths regardless of its configuration, already reaches a position error of0\.230\.23to0\.37​Å20\.37\\,\\text\{\\AA\}^\{2\}across the five\-step horizon;Newmark\-β\\beta\-DGNstays only44to10%10\\%below it while remaining bounded\. The benchmark therefore tests bounded, stable rollout of the collective transition and graph sparsity, not a large accuracy margin over this trivial predictor\.

Task and split\.We predict the configurationΔ=15\\Delta=15frames ahead over five autoregressive steps, a7575\-frame horizon spanning most of the transition\. Because each path is a distinct trajectory, the split is by path:140140training,3030validation and3030test paths \(241241five\-step test windows\)\. All models receive atomic positions and the causal backward\-difference velocities of Equation \([2](https://arxiv.org/html/2609.30344#S3.E2)\)\.Newmark\-β\\beta\-DGNuses a one\-hot atom\-type feature, the17081708directed covalent bonds and one hub with17101710virtual edges \(34183418directed edges\), identical to its equilibrium configuration\.Dynami\-CAL GraphNetis the controlled ablation with the same bonds and no hub\. EGHN and EGNO augment the bonds with their published10​Å10\\,\\text\{\\AA\}radius graph, adding approximately55 61055\\,610directed contact edges\.

Atom type and partial charge are distinct physical descriptors and are reported explicitly rather than treated as equivalent inputs\. Adding raw partial charge toNewmark\-β\\beta\-DGNprevented convergence in our implementation, soNewmark\-β\\beta\-DGNis evaluated with atom type alone\.EGHN and EGNO retain their published charge input\.

Supplementary Table 14:Protein representation supplied toNewmark\-β\\beta\-DGN\.The input and contact graph used by EGHN and EGNO aredescribed in the text\.Metric\.Unscaled position mean squared error, per rollout step, over five autoregressive steps of1515frames each\.

### 6\.2Evaluated baselines

All four models are trained on this ensemble under the by\-path split, one seed each\.Newmark\-β\\beta\-DGNuses Equation \([1](https://arxiv.org/html/2609.30344#S1.E1)\);Dynami\-CAL GraphNetapplies the same normalised\-increment objective on the hub\-free bond graph; EGHN uses the velocity head of Section[2\.1](https://arxiv.org/html/2609.30344#S2.SS1a)with itsλlink\\lambda\_\{\\mathrm\{link\}\}\-weighted linkage term; EGNO uses its multi\-frame objective\. All are evaluated by the same unscaled position error under autoregressive rollout and retain the checkpoint with minimum validation loss\. The contact\-graph models carry the cost of their10​Å10\\,\\text\{\\AA\}graph: approximately55 61055\\,610directed contact edges against the17101710hub edges ofNewmark\-β\\beta\-DGN, so the hub supplies the global coupling required for the collective transition at a fraction of the graph density\.

### 6\.3Hyperparameters

Settings forNewmark\-β\\beta\-DGNnot listed in Table[15](https://arxiv.org/html/2609.30344#S6.T15)follow Table[1](https://arxiv.org/html/2609.30344#S1.T1)\. This is the one case in which the latent width ofEGHN and EGNOdeparts from ours: the reference protocol trains EGHN and EGNO at128128, and we retain their settings rather than reduce them\.

Supplementary Table 15:Protein dynamics: hyperparameters\.Newmark\-β\\beta\-DGNandDynami\-CAL GraphNetare trained atn​f=64nf=64, as on every other system; EGHN and EGNO retain their published width of128128\.All four models train with early stopping at a patience of5050validation checks\. Prediction spanΔ=15\\Delta=15frames throughout\. Seed4242\.

### 6\.4Additional results

Table[16](https://arxiv.org/html/2609.30344#S6.T16)gives the full rollout for all evaluated models\.Newmark\-β\\beta\-DGNremains bounded across the full horizon, rising from0\.2210\.221to a maximum of0\.3460\.346and settling near0\.28​Å20\.28\\,\\text\{\\AA\}^\{2\}\.Dynami\-CAL GraphNet, the hub\-free ablation, is accurate at one step \(0\.3860\.386\) but its rollout diverges to6\.126\.12; EGHN and EGNO track to the third step and then diverge, reaching361361and non\-finite values\. The state\-independent mean\-trajectory baseline applies the training\-averaged displacement at each step, regardless of the test path \(last row of Table[16](https://arxiv.org/html/2609.30344#S6.T16)\)\.

Supplementary Table 16:Per\-step rollout error on the transition ensemble\.Unscaled position mean squared error \(Å2\)\. One step is1515frames, evaluated on the3030held\-out test paths \(241241five\-step windows,855855atoms\) from each model’s minimum\-validation checkpoint\. NaN entries mark a diverged rollout\.
### 6\.5Why we do not use the equilibrium benchmark

Adenylate kinase also has an established prediction benchmark on its*equilibrium*trajectory \(apo, explicit solvent,300300\\,K and11\\,bar, recorded every240240\\,ps over approximately1​μ1\\,\\mus, the same855855backbone atoms, a1515\-frame horizon\), introduced by EGHN\[[4](https://arxiv.org/html/2609.30344#bib.bib4)\]and adopted by EGNO\[[5](https://arxiv.org/html/2609.30344#bib.bib5)\]\. We do not use it to rank dynamical models: its single\-step error is bounded below by a trivial predictor that no learned model beats, so it scores thermal\-fluctuation statistics rather than learnable, collective dynamics\. All numbers below are measured on the trajectory, using the reference benchmark split and using the unscaled single\-step position mean squared error\.

The motion is a confined fluctuation\.Over the trajectory each backbone atom fluctuates about a fixed mean position with a root\-mean\-square amplitude of1\.9​Å1\.9\\,\\text\{\\AA\}, from0\.8​Å0\.8\\,\\text\{\\AA\}in the rigid core to about6​Å6\\,\\text\{\\AA\}in the flexible loops \(Fig\.[5](https://arxiv.org/html/2609.30344#S6.F5)a\), with no net translation or change of fold\. The position autocorrelation at the3\.63\.6\\,ns \(1515\-frame\) horizon isC=0\.40C=0\.40, a relaxation time ofτ≈2\.7\\tau\\approx 2\.7\\,ns\.

A linear restoring map captures the predictable component\.The per\-atom position increments are approximately Gaussian, with an excess kurtosis of0\.710\.71\. The larger kurtosis obtained by pooling all atoms results from combining rigid\-core and flexible\-loop atoms with substantially different fluctuation amplitudes\. Under a stationary Gaussian approximation, the minimum\-mean\-square predictor based on the present position is

𝐱^t\+1=𝐱¯\+C⁡\(𝐱t−𝐱¯\),C=0\.40\.\\widehat\{\\mathbf\{x\}\}\_\{t\+1\}=\\bar\{\\mathbf\{x\}\}\+C\\left\(\\mathbf\{x\}\_\{t\}\-\\bar\{\\mathbf\{x\}\}\\right\),\\qquad C=0\.40\.Equivalently, the predicted displacement is

Δ​𝐱^=\(C−1\)​\(𝐱t−𝐱¯\)=−0\.60​\(𝐱t−𝐱¯\)\.\\Delta\\widehat\{\\mathbf\{x\}\}=\(C\-1\)\\left\(\\mathbf\{x\}\_\{t\}\-\\bar\{\\mathbf\{x\}\}\\right\)=\-0\.60\\left\(\\mathbf\{x\}\_\{t\}\-\\bar\{\\mathbf\{x\}\}\\right\)\.This one\-coefficient rule obtains a single\-step mean squared error of1\.7252​Å21\.7252\\,\\text\{\\AA\}^\{2\}, while fitting one coefficient per atom reduces the error to1\.6915​Å21\.6915\\,\\text\{\\AA\}^\{2\}, the lowest value obtained by the predictors examined here\. Relative to the assume\-no\-motion error of2\.6127​Å22\.6127\\,\\text\{\\AA\}^\{2\}, the per\-atom map explains0\.9212​Å20\.9212\\,\\text\{\\AA\}^\{2\}of the displacement error\. The remaining1\.6915​Å21\.6915\\,\\text\{\\AA\}^\{2\}shows no detectable correlation with the present position, velocity, or neighbouring\-atom coordinates\.

![Refer to caption](https://arxiv.org/html/2609.30344v1/protein_floor.png)Supplementary Figure 5:The single\-step benchmark is bounded by a data floor\.\(a\) One backbone atom’s position at every recorded frame, relative to its mean; its root\-mean\-square fluctuation is1\.3​Å1\.3\\,\\text\{\\AA\}\. \(b\) Single\-step position error for assume\-no\-motion, the published and reproduced learned models, and the one\-coefficient restoring rule \(Δ​𝐱=−0\.60​\(𝐱−𝐱¯\)\\Delta\\mathbf\{x\}=\-0\.60\(\\mathbf\{x\}\-\\bar\{\\mathbf\{x\}\}\), labelled the spring rule\), against the lowest error any predictor of the current state can reach \(1\.6921\.692, dashed\)\. Every learned model lies above the restoring rule and the floor\. Lower is better\.No model beats the floor\.Newmark\-β\\beta\-DGN\(1\.88681\.8868\), EGHN \(2\.08012\.0801\) and the published\[[5](https://arxiv.org/html/2609.30344#bib.bib5)\]EGNO \(2\.2312\.231\),Dynami\-CAL GraphNet\(2\.3012\.301\) and EGHNO \(1\.8011\.801\) all lie above both the1\.72521\.7252of the one\-coefficient restoring rule and the1\.69151\.6915floor \(Fig\.[5](https://arxiv.org/html/2609.30344#S6.F5)b\)\. A linear map using all atom coordinates jointly reaches29\.9%29\.9\\%skill, no better than the per\-atom rule, and after the restoring prediction is removed a regularised nonlinear model explains none of the remaining error out of sample\. The ranking among the learned models therefore measures how closely each reproduces one linear coefficient, not the accuracy of a dynamical model\.

Why the horizon leaves only statistics\.The atomic trajectory decorrelates on a Lyapunov timescale of about0\.20\.2\\,ps, so the3\.63\.6\\,ns horizon is roughly18,00018\{,\}000times longer and only the2\.72\.7\\,ns statistical relaxation survives\. Shortening the horizon lowers the attainable skill, to21%21\\%at a single frame, so no choice of horizon converts this benchmark into a test of dynamics\. We therefore evaluate the collective dynamics on the driven transition ensemble of Section[6\.1](https://arxiv.org/html/2609.30344#S6.SS1), where the models are separated by their rollout stability rather than by a shared linear floor\.

## 7Joint\-Moment Inference in Human Walking

### 7\.1Case description

We use a public dataset of healthy human walking recorded on an instrumented treadmill\[[17](https://arxiv.org/html/2609.30344#bib.bib10)\], which provides marker trajectories together with joint moments and ground reaction forces obtained from an inverse\-dynamics pipeline\. We use subject p2, comprising3333trials and18 63118\\,631frames recorded at120120\\,Hz\. The measured moments are read only for evaluation and never enter Equation \([1](https://arxiv.org/html/2609.30344#S1.E1)\)\.

Observed inputs and graph\.The3737markers are the nodes and carry positions and the causal velocities of Equation \([2](https://arxiv.org/html/2609.30344#S3.E2)\)\. The6363anatomical connections between markers are the physical edges, andNewmark\-β\\beta\-DGNadds7474directed virtual edges to one hub\. The model receives no joint moment or ground reaction force\. One prediction span isdelta\_frame=30\\texttt\{delta\\\_frame\}=30frames, matching the motion\-capture case\.

Supplementary Table 17:Graph representation of the biomechanics system\.Gait conditions and splits\.The split is by gait condition rather than by time\. Training covers preferred walking at0\.70\.7,1\.11\.1,1\.251\.25,1\.41\.4,1\.61\.6and1\.8​m​s−11\.8\\,\\mathrm\{m\\,s^\{\-1\}\}, together with two cadence variants at1\.25​m​s−11\.25\\,\\mathrm\{m\\,s^\{\-1\}\}in which step frequency is raised while speed is held fixed\. Validation covers preferred walking at0\.90\.9and1\.25​m​s−11\.25\\,\\mathrm\{m\\,s^\{\-1\}\}and the remainder of the1\.25​m​s−11\.25\\,\\mathrm\{m\\,s^\{\-1\}\}cadence sweep\. Held out entirely are the two constrained gait strategies, constant step length and constant step frequency, at every recorded speed, together with preferred walking at2\.0​m​s−12\.0\\,\\mathrm\{m\\,s^\{\-1\}\}\.

Cadence and gait strategy are distinct axes\. Varying step frequency at a fixed preferred speed is a condition on which the model is trained and validated\. Holding step length or step frequency fixed while speed varies is not\. The non\-monotonic cadence response reported in the main text therefore concerns supervision rather than generalisation, whereas the constant\-step\-length, constant\-step\-frequency and2\.0​m​s−12\.0\\,\\mathrm\{m\\,s^\{\-1\}\}conditions are the extrapolation result\.

Moment readout\.From the decoded internal forces the moment about each joint centre is formed by summing over the markers distal to that joint,𝐌=−∑i∈distal𝐫i×𝐅i\\mathbf\{M\}=\-\\sum\_\{i\\in\\text\{distal\}\}\\mathbf\{r\}\_\{i\}\\times\\mathbf\{F\}\_\{i\}, with𝐫i\\mathbf\{r\}\_\{i\}measured from the joint centre\. Only this contribution, the moment about the joint centre, is used\. No term of the loss constrains the spin channel, so it is not an identifiable physical quantity and is not read out\. The summation runs over thirteen markers at the hip, seven at the knee and two at the ankle\.

Evaluation sampling\.Trend metrics tile each trial with non\-overlapping windows rather than sampling start frames at random\. The two legs of a walking human are close to half a gait cycle out of phase, so a short randomly placed window can capture the stance\-phase peak of one leg while clipping the other, and averaging over many such windows does not cancel the bias because the sampling is not phase\-aware\. Deterministic tiling covers every phase of every trial exactly once\. The trial supplying the normalisation denominator is tiled at unit stride so that the divisor is low\-noise\. The measured moments are read at the native120120\\,Hz rate rather than at the sparse prediction\-span rate: sampling the ground truth sparsely aliases the≈0\.8\\approx\\\!0\.8\\,s stride period and returns a left–right ankle RMS ratio of0\.470\.47where the densely sampled signal gives0\.990\.99\.

Normalisation\.Demand is a ratio rather than a calibrated moment\. For each joint the root\-mean\-square of the moment profile is divided by the corresponding value in the preferred1\.25​m​s−11\.25\\,\\mathrm\{m\\,s^\{\-1\}\}trial, at unit subject mass\. Measured and inferred series are normalised separately\. Within each series, both legs are divided by the mean of that series’ left and right baseline values, so one divisor is used per source and the raw left–right asymmetry survives into the plotted curves\. Dividing each leg by its own baseline would remove part of that asymmetry\. Pearson correlations are invariant to such per\-series rescaling; the normalised RMS errors and the asymmetry statistics are not, and all quoted values use the shared\-divisor convention\.

Seeds\.A single trained model is used, on a single subject\. The error bars in the figures of this section are the spread over evaluation tiles of a trial and not over seeds; the tiled evaluation is deterministic\.

### 7\.2Evaluated baselines

No baseline is evaluated in this case\. The purpose is to compare the joint\-moment readout ofNewmark\-β\\beta\-DGNwith an independent inverse\-dynamics reference, not to rank forward simulators\. Forward prediction on articulated human motion is compared with the baselines in Section[5](https://arxiv.org/html/2609.30344#S5)\.

### 7\.3Hyperparameters

Settings not listed in Table[18](https://arxiv.org/html/2609.30344#S7.T18)follow Table[1](https://arxiv.org/html/2609.30344#S1.T1)\.

Supplementary Table 18:Biomechanics: hyperparameters\.At most20002000epochs with early stopping\. Prediction spandelta\_frame=30\\texttt\{delta\\\_frame\}=30frames\.

### 7\.4Additional results

The main text reports the knee\. Figures[6](https://arxiv.org/html/2609.30344#S7.F6)and[7](https://arxiv.org/html/2609.30344#S7.F7)give the hip and the ankle alongside it, for the cadence sweep of the main text and for the in\-distribution control\. The remaining sweeps behave comparably and are not reproduced here\.

Pooled across five evaluation experiments and both legs, the inferred profile agrees with the measured one at a Pearson correlation of0\.9430\.943at the hip,0\.8690\.869at the knee and0\.4920\.492at the ankle, with normalised RMS errors of0\.1710\.171,0\.2130\.213and0\.8870\.887\. The ankle is not recovered: it is anticorrelated with measurement in one condition, atr=−0\.322r=\-0\.322, and its magnitude is low by roughly an order of magnitude\. The left–right asymmetry of the gait is likewise absent\. The measured imbalance ranges from0\.0060\.006to0\.1700\.170across conditions while the inferred imbalance remains at0\.048±0\.0170\.048\\pm 0\.017, and the prediction error grows with the true asymmetry, atr=\+0\.76r=\+0\.76, reaching a mean error of17\.3%17\.3\\%in the condition where that asymmetry is greatest\.

A contributing limitation is that the joint\-moment readout is assembled from decoded internal segment forces and does not include the measured ground\-reaction contribution; this is particularly restrictive at the ankle, where only two distal markers contribute to the readout\.

Supplementary Figure 6:Inferred demand at all three joints across the cadence sweep\.Hip, knee and ankle, measured against inferred, for the conditions of the main\-text figure\. The hip and knee track the measurement; the ankle does not, for the reason given above\. Shaded groups lie within the training and validation data\.Supplementary Figure 7:In\-distribution control\.The same three joints evaluated on conditions seen during training\. Agreement at the knee \(r=0\.81r=0\.81\) is as good on these training conditions as on the withheld ones \(r=0\.88r=0\.88\)\. Readout accuracy is therefore limited by the moment model, not by distribution shift, and the inferred profiles are not artefacts of fitting\.

## 8Ablations

This section collects the ablations ofNewmark\-β\\beta\-DGN’s own components\. Each removes or replaces one ingredient of the update and reports its effect\. Four are established where they arise in the case sections and are summarised in Table[19](https://arxiv.org/html/2609.30344#S8.T19); the angular\-momentum channel, whose effect is on the inferred mechanics rather than on the rollout, is treated below\.

Supplementary Table 19:Ablations ofNewmark\-β\\beta\-DGN’s components\.Each row removes or replaces one ingredient of the update\. The first four are established in the sections cited; the last is treated in Section 8\.1\. Beam errors are whole\-body, mean over the reported configurations, seed 42\.### 8\.1The angular\-momentum channel

Newmark\-β\\beta\-DGNcarries a per\-node spin and a per\-edge angular\-momentum flux alongside the linear channel \(the Cosserat channel ofDynami\-CAL GraphNet\[[15](https://arxiv.org/html/2609.30344#bib.bib1)\], Section[10\.5](https://arxiv.org/html/2609.30344#S10.SS5)\)\. To isolate its contribution we train a variant that decodes only linear momentum and translational response operators, with no spin state and no angular flux; the rest of the update is unchanged\. The variant is SO\(3\)\-equivariant to3×10−63\\times 10^\{\-6\}on motion capture and exactly equivariant in double precision on the beam\.

Removing the channel does not cost rollout accuracy\. On motion capture the autoregressive error of the linear\-only variant is slightly lower than the full model, by66to13%13\\%across the five rollout steps \(three seeds, with a widening spread\), and on the biomechanics rollout the two are tied \(0\.4770\.477against0\.4800\.480\\,mm\)\. What changes is the inferred mechanics\. The joint moments of Section[7](https://arxiv.org/html/2609.30344#S7)are assembled from the decoded internal forces by𝐌=−∑i𝐫i×𝐅i\\mathbf\{M\}=\-\\sum\_\{i\}\\mathbf\{r\}\_\{i\}\\times\\mathbf\{F\}\_\{i\}, a construction identical in both models, yet without the angular channel the inferred demand flattens and no longer follows the measured trend across gait conditions \(Fig\.[8](https://arxiv.org/html/2609.30344#S8.F8)\)\. The channel therefore earns its place through the quality of the inferred mechanics, not through rollout error, and we retain it\.

Supplementary Figure 8:The angular\-momentum channel is needed for the inferred joint moments, not for the rollout\.Inferred against measured hip, knee and ankle demand across the cadence sweep, for the full model \(left\) and the variant with the angular\-momentum channel removed \(right\)\. Measured demand is black and grey; inferred is red\. With the channel the inferred hip and knee demand follows the measured trend; without it the inferred demand flattens toward unity and loses the trend\. The moment construction𝐌=−∑i𝐫i×𝐅i\\mathbf\{M\}=\-\\sum\_\{i\}\\mathbf\{r\}\_\{i\}\\times\\mathbf\{F\}\_\{i\}is identical in both, so the difference is in the decoded force\. The ankle is not recovered by either model, for the reason given in Section[7\.4](https://arxiv.org/html/2609.30344#S7.SS4)\.A mechanism is plausible but not established here\. The decoded angular flux is antisymmetric on each edge and sums to zero, so the channel exchanges angular momentum conservatively, whereas the linear\-only variant places no such constraint on its decoded force\. We hypothesise that this missing constraint is why the force\-based moment readout degrades\. We have not measured the angular\-momentum residual of either rollout, so the link is stated as a hypothesis rather than a demonstrated cause\.

## 9From implicit coupling to an operator\-weighted hub

The hub provides a tractable approximation of the non\-local coupling produced by a global implicit mechanical solve\. The derivation separates the two sides of the update\. On the left\-hand side, static condensation produces a dense effective interaction operator whose rank\-one approximation admits a star representation\. On the right\-hand side, the same condensation transfers the force drive from eliminated degrees of freedom to the retained degrees of freedom\.Newmark\-β\\beta\-DGNretains this star topology and approximates both effects without eliminating physical nodes or assembling a global Schur complement\.

### 9\.1Implicit time stepping produces non\-local coupling

Consider the semi\-discrete equations of motion over one linearised substep,

𝐌​𝐱¨\+𝐃glob​𝐱˙\+𝐊glob​𝐱=𝐟ext\.\\mathbf\{M\}\\ddot\{\\mathbf\{x\}\}\+\\mathbf\{D\}\_\{\\mathrm\{glob\}\}\\dot\{\\mathbf\{x\}\}\+\\mathbf\{K\}\_\{\\mathrm\{glob\}\}\\mathbf\{x\}=\\mathbf\{f\}\_\{\\mathrm\{ext\}\}\.\(3\)The mass matrix𝐌\\mathbf\{M\}is block diagonal, whereas the assembled stiffness𝐊glob\\mathbf\{K\}\_\{\\mathrm\{glob\}\}and damping𝐃glob\\mathbf\{D\}\_\{\\mathrm\{glob\}\}couple connected nodes\. Define

Δ​𝐯:=𝐯t\+1−𝐯t,Δ​𝐱:=𝐱t\+1−𝐱t\.\\Delta\\mathbf\{v\}:=\\mathbf\{v\}\_\{t\+1\}\-\\mathbf\{v\}\_\{t\},\\qquad\\Delta\\mathbf\{x\}:=\\mathbf\{x\}\_\{t\+1\}\-\\mathbf\{x\}\_\{t\}\.\(4\)For the average\-acceleration Newmark parametersγ=12\\gamma=\\tfrac\{1\}\{2\}andβ=14\\beta=\\tfrac\{1\}\{4\},

𝐌​Δ​𝐯=\[𝐟ext\+12​\(𝐟intt\+𝐟intt\+1\)\]​δ​t,Δ​𝐱=δ​t​\(𝐯t\+12​Δ​𝐯\)\.\\mathbf\{M\}\\Delta\\mathbf\{v\}=\\left\[\\mathbf\{f\}\_\{\\mathrm\{ext\}\}\+\\tfrac\{1\}\{2\}\\left\(\\mathbf\{f\}\_\{\\mathrm\{int\}\}^\{t\}\+\\mathbf\{f\}\_\{\\mathrm\{int\}\}^\{t\+1\}\\right\)\\right\]\\delta t,\\qquad\\Delta\\mathbf\{x\}=\\delta t\\left\(\\mathbf\{v\}\_\{t\}\+\\tfrac\{1\}\{2\}\\Delta\\mathbf\{v\}\\right\)\.\(5\)Within the substep, the linearised internal\-force change is

𝐟intt\+1=𝐟intt−𝐃glob​Δ​𝐯−𝐊glob​Δ​𝐱\.\\mathbf\{f\}\_\{\\mathrm\{int\}\}^\{t\+1\}=\\mathbf\{f\}\_\{\\mathrm\{int\}\}^\{t\}\-\\mathbf\{D\}\_\{\\mathrm\{glob\}\}\\Delta\\mathbf\{v\}\-\\mathbf\{K\}\_\{\\mathrm\{glob\}\}\\Delta\\mathbf\{x\}\.\(6\)Using the Newmark position increment in Equation \([6](https://arxiv.org/html/2609.30344#S9.E6)\) gives

𝐟intt\+1\\displaystyle\\mathbf\{f\}\_\{\\mathrm\{int\}\}^\{t\+1\}=𝐟intt−𝐃glob​Δ​𝐯−δ​t​𝐊glob​𝐯t−12​δ​t​𝐊glob​Δ​𝐯,\\displaystyle=\\mathbf\{f\}\_\{\\mathrm\{int\}\}^\{t\}\-\\mathbf\{D\}\_\{\\mathrm\{glob\}\}\\Delta\\mathbf\{v\}\-\\delta t\\,\\mathbf\{K\}\_\{\\mathrm\{glob\}\}\\mathbf\{v\}\_\{t\}\-\\tfrac\{1\}\{2\}\\delta t\\,\\mathbf\{K\}\_\{\\mathrm\{glob\}\}\\Delta\\mathbf\{v\},12​\(𝐟intt\+𝐟intt\+1\)\\displaystyle\\tfrac\{1\}\{2\}\\left\(\\mathbf\{f\}\_\{\\mathrm\{int\}\}^\{t\}\+\\mathbf\{f\}\_\{\\mathrm\{int\}\}^\{t\+1\}\\right\)=𝐟intt−12​𝐃glob​Δ​𝐯−12​δ​t​𝐊glob​𝐯t−14​δ​t​𝐊glob​Δ​𝐯\.\\displaystyle=\\mathbf\{f\}\_\{\\mathrm\{int\}\}^\{t\}\-\\tfrac\{1\}\{2\}\\mathbf\{D\}\_\{\\mathrm\{glob\}\}\\Delta\\mathbf\{v\}\-\\tfrac\{1\}\{2\}\\delta t\\,\\mathbf\{K\}\_\{\\mathrm\{glob\}\}\\mathbf\{v\}\_\{t\}\-\\tfrac\{1\}\{4\}\\delta t\\,\\mathbf\{K\}\_\{\\mathrm\{glob\}\}\\Delta\\mathbf\{v\}\.\(7\)Multiplying the second line byδ​t\\delta tand collecting all terms proportional toΔ​𝐯\\Delta\\mathbf\{v\}on the left gives

\(𝐌\+𝐂glob\)​Δ​𝐯⏟left\-hand side: implicit mechanical response=𝐫⏟right\-hand side: known force drive,\\underbrace\{\\left\(\\mathbf\{M\}\+\\mathbf\{C\}\_\{\\mathrm\{glob\}\}\\right\)\\Delta\\mathbf\{v\}\}\_\{\\text\{left\-hand side: implicit mechanical response\}\}=\\underbrace\{\\mathbf\{r\}\}\_\{\\text\{right\-hand side: known force drive\}\},\(8\)where

𝐂glob:=12​δ​t​𝐃glob\+14​δ​t2​𝐊glob,\\mathbf\{C\}\_\{\\mathrm\{glob\}\}:=\\tfrac\{1\}\{2\}\\delta t\\,\\mathbf\{D\}\_\{\\mathrm\{glob\}\}\+\\tfrac\{1\}\{4\}\\delta t^\{2\}\\mathbf\{K\}\_\{\\mathrm\{glob\}\},\(9\)and

𝐫:=δ​t​\(𝐟ext\+𝐟intt\)−12​δ​t2​𝐊glob​𝐯t\.\\mathbf\{r\}:=\\delta t\\left\(\\mathbf\{f\}\_\{\\mathrm\{ext\}\}\+\\mathbf\{f\}\_\{\\mathrm\{int\}\}^\{t\}\\right\)\-\\tfrac\{1\}\{2\}\\delta t^\{2\}\\mathbf\{K\}\_\{\\mathrm\{glob\}\}\\mathbf\{v\}\_\{t\}\.\(10\)Section[10](https://arxiv.org/html/2609.30344#S10)derives the complete update used byNewmark\-β\\beta\-DGN\.

The mass matrix is local, while𝐂glob\\mathbf\{C\}\_\{\\mathrm\{glob\}\}contains the inter\-node coupling\. Although the assembled system is sparse, eliminating some degrees of freedom generally produces dense effective coupling on the left\-hand side and dense force transfer on the right\-hand side\.

### 9\.2Static condensation of the implicit update

##### Retained and internal degrees of freedom\.

Static condensation divides the degrees of freedom into a retained setRRand an internal setII\. The retained set contains the nodes whose increments remain explicit in the reduced problem\. The internal set contains intermediate nodes that are eliminated algebraically\.

Eliminating an internal node does not discard its mechanical effect\. Its increment is solved in terms of the retained increments and the applied drive, and the result is substituted into the retained equations\. A path between two retained nodes can therefore be replaced by a direct effective coupling after the intermediate internal nodes have been eliminated\. This is how a sparse local operator can produce a dense reduced operator\.

The retained/internal partition is used here only as a thought experiment to reveal the required interaction topology\.Newmark\-β\\beta\-DGNdoes not eliminate physical nodes\.

##### The Schur complement\.

Consider the generic partitioned system

\[𝐀𝐁𝐄𝐃\]​\[𝐱𝐲\]=\[𝐟𝐠\],\\begin\{bmatrix\}\\mathbf\{A\}&\\mathbf\{B\}\\\\ \\mathbf\{E\}&\\mathbf\{D\}\\end\{bmatrix\}\\begin\{bmatrix\}\\mathbf\{x\}\\\\ \\mathbf\{y\}\\end\{bmatrix\}=\\begin\{bmatrix\}\\mathbf\{f\}\\\\ \\mathbf\{g\}\\end\{bmatrix\},\(11\)where𝐱\\mathbf\{x\}contains the retained unknowns and𝐲\\mathbf\{y\}contains the internal unknowns\. The internal block row gives

𝐲=𝐃−1​\(𝐠−𝐄𝐱\),\\mathbf\{y\}=\\mathbf\{D\}^\{\-1\}\\left\(\\mathbf\{g\}\-\\mathbf\{E\}\\mathbf\{x\}\\right\),\(12\)provided that𝐃\\mathbf\{D\}is invertible\. Substituting this result into the retained block row gives

\(𝐀−𝐁𝐃−1​𝐄\)​𝐱⏟condensed left\-hand side=𝐟−𝐁𝐃−1​𝐠⏟condensed right\-hand side\.\\underbrace\{\\left\(\\mathbf\{A\}\-\\mathbf\{B\}\\mathbf\{D\}^\{\-1\}\\mathbf\{E\}\\right\)\\mathbf\{x\}\}\_\{\\text\{condensed left\-hand side\}\}=\\underbrace\{\\mathbf\{f\}\-\\mathbf\{B\}\\mathbf\{D\}^\{\-1\}\\mathbf\{g\}\}\_\{\\text\{condensed right\-hand side\}\}\.\(13\)The matrix

𝐀−𝐁𝐃−1​𝐄\\mathbf\{A\}\-\\mathbf\{B\}\\mathbf\{D\}^\{\-1\}\\mathbf\{E\}\(14\)is the Schur complement of the internal block𝐃\\mathbf\{D\}\. It is the exact operator acting on the retained unknowns after the internal unknowns have been eliminated\. The same elimination transfers the internal drive𝐠\\mathbf\{g\}to the retained right\-hand side through−𝐁𝐃−1​𝐠\-\\mathbf\{B\}\\mathbf\{D\}^\{\-1\}\\mathbf\{g\}\.

Thus, static condensation has two simultaneous effects: it modifies the retained operator on the left\-hand side and transfers the internal drive to the retained right\-hand side\.

##### Application to the Newmark system\.

Partition the velocity increments and the known force drive in Equation \([8](https://arxiv.org/html/2609.30344#S9.E8)\) as

Δ​𝐯=\[Δ​𝐯RΔ​𝐯I\],𝐫=\[𝐫R𝐫I\]\.\\Delta\\mathbf\{v\}=\\begin\{bmatrix\}\\Delta\\mathbf\{v\}\_\{R\}\\\\ \\Delta\\mathbf\{v\}\_\{I\}\\end\{bmatrix\},\\qquad\\mathbf\{r\}=\\begin\{bmatrix\}\\mathbf\{r\}\_\{R\}\\\\ \\mathbf\{r\}\_\{I\}\\end\{bmatrix\}\.\(15\)Under the same partition,

𝐌=\[𝐌R𝟎𝟎𝐌I\],𝐂glob=\[𝐂R​R𝐂R​I𝐂I​R𝐂I​I\]\.\\mathbf\{M\}=\\begin\{bmatrix\}\\mathbf\{M\}\_\{R\}&\\mathbf\{0\}\\\\ \\mathbf\{0\}&\\mathbf\{M\}\_\{I\}\\end\{bmatrix\},\\qquad\\mathbf\{C\}\_\{\\mathrm\{glob\}\}=\\begin\{bmatrix\}\\mathbf\{C\}\_\{RR\}&\\mathbf\{C\}\_\{RI\}\\\\ \\mathbf\{C\}\_\{IR\}&\\mathbf\{C\}\_\{II\}\\end\{bmatrix\}\.\(16\)Here𝐂R​R\\mathbf\{C\}\_\{RR\}maps retained increments into the retained equations,𝐂R​I\\mathbf\{C\}\_\{RI\}maps internal increments into the retained equations,𝐂I​R\\mathbf\{C\}\_\{IR\}maps retained increments into the internal equations, and𝐂I​I\\mathbf\{C\}\_\{II\}maps internal increments into the internal equations\. These are precisely the four blocks of𝐂glob\\mathbf\{C\}\_\{\\mathrm\{glob\}\}defined in Equation \([9](https://arxiv.org/html/2609.30344#S9.E9)\)\.

Using these blocks, Equation \([8](https://arxiv.org/html/2609.30344#S9.E8)\) becomes

\[𝐌R\+𝐂R​R𝐂R​I𝐂I​R𝐌I\+𝐂I​I\]​\[Δ​𝐯RΔ​𝐯I\]=\[𝐫R𝐫I\]\.\\begin\{bmatrix\}\\mathbf\{M\}\_\{R\}\+\\mathbf\{C\}\_\{RR\}&\\mathbf\{C\}\_\{RI\}\\\\ \\mathbf\{C\}\_\{IR\}&\\mathbf\{M\}\_\{I\}\+\\mathbf\{C\}\_\{II\}\\end\{bmatrix\}\\begin\{bmatrix\}\\Delta\\mathbf\{v\}\_\{R\}\\\\ \\Delta\\mathbf\{v\}\_\{I\}\\end\{bmatrix\}=\\begin\{bmatrix\}\\mathbf\{r\}\_\{R\}\\\\ \\mathbf\{r\}\_\{I\}\\end\{bmatrix\}\.\(17\)The correspondence with Equation \([11](https://arxiv.org/html/2609.30344#S9.E11)\) is

𝐀=𝐌R\+𝐂R​R,𝐁=𝐂R​I,𝐄=𝐂I​R,𝐃=𝐌I\+𝐂I​I\.\\mathbf\{A\}=\\mathbf\{M\}\_\{R\}\+\\mathbf\{C\}\_\{RR\},\\qquad\\mathbf\{B\}=\\mathbf\{C\}\_\{RI\},\\qquad\\mathbf\{E\}=\\mathbf\{C\}\_\{IR\},\\qquad\\mathbf\{D\}=\\mathbf\{M\}\_\{I\}\+\\mathbf\{C\}\_\{II\}\.\(18\)The internal block row therefore gives

Δ​𝐯I=\(𝐌I\+𝐂I​I\)−1​\(𝐫I−𝐂I​R​Δ​𝐯R\)\.\\Delta\\mathbf\{v\}\_\{I\}=\\left\(\\mathbf\{M\}\_\{I\}\+\\mathbf\{C\}\_\{II\}\\right\)^\{\-1\}\\left\(\\mathbf\{r\}\_\{I\}\-\\mathbf\{C\}\_\{IR\}\\Delta\\mathbf\{v\}\_\{R\}\\right\)\.\(19\)Substituting this result into the retained block row gives

\[𝐌R\+𝐂R​R−𝐂R​I​\(𝐌I\+𝐂I​I\)−1​𝐂I​R\]​Δ​𝐯R⏟condensed left\-hand side=𝐫R−𝐂R​I​\(𝐌I\+𝐂I​I\)−1​𝐫I⏟condensed right\-hand side\.\\underbrace\{\\left\[\\mathbf\{M\}\_\{R\}\+\\mathbf\{C\}\_\{RR\}\-\\mathbf\{C\}\_\{RI\}\\left\(\\mathbf\{M\}\_\{I\}\+\\mathbf\{C\}\_\{II\}\\right\)^\{\-1\}\\mathbf\{C\}\_\{IR\}\\right\]\\Delta\\mathbf\{v\}\_\{R\}\}\_\{\\text\{condensed left\-hand side\}\}=\\underbrace\{\\mathbf\{r\}\_\{R\}\-\\mathbf\{C\}\_\{RI\}\\left\(\\mathbf\{M\}\_\{I\}\+\\mathbf\{C\}\_\{II\}\\right\)^\{\-1\}\\mathbf\{r\}\_\{I\}\}\_\{\\text\{condensed right\-hand side\}\}\.\(20\)Equation \([20](https://arxiv.org/html/2609.30344#S9.E20)\) is the exact condensed Newmark equation\. Its left\- and right\-hand\-side consequences are considered separately below\.

#### 9\.2\.1Left\-hand side: Schur\-complement coupling and the rank\-one virtual\-hub topology

To isolate the topology produced specifically by the stiffness–damping interaction, consider the static\-condensation limit

𝐌I=𝟎\.\\mathbf\{M\}\_\{I\}=\\mathbf\{0\}\.\(21\)The matrix𝐌I\\mathbf\{M\}\_\{I\}is the mass block associated with the internal degrees of freedom\. Setting𝐌I=𝟎\\mathbf\{M\}\_\{I\}=\\mathbf\{0\}means that the internal nodes are treated as massless and equilibrate without an inertial contribution\. This auxiliary limit is used only to expose the topology of𝐂glob\\mathbf\{C\}\_\{\\mathrm\{glob\}\}\. It does not mean that physical nodes are massless inNewmark\-β\\beta\-DGN\.

For this static interaction problem, the blocks in Equation \([11](https://arxiv.org/html/2609.30344#S9.E11)\) are

𝐀=𝐂R​R,𝐁=𝐂R​I,𝐄=𝐂I​R,𝐃=𝐂I​I\.\\mathbf\{A\}=\\mathbf\{C\}\_\{RR\},\\qquad\\mathbf\{B\}=\\mathbf\{C\}\_\{RI\},\\qquad\\mathbf\{E\}=\\mathbf\{C\}\_\{IR\},\\qquad\\mathbf\{D\}=\\mathbf\{C\}\_\{II\}\.\(22\)The condensed interaction operator is therefore

𝐂cond:=𝐂R​R−𝐂R​I​𝐂I​I−1​𝐂I​R\.\\mathbf\{C\}\_\{\\mathrm\{cond\}\}:=\\mathbf\{C\}\_\{RR\}\-\\mathbf\{C\}\_\{RI\}\\mathbf\{C\}\_\{II\}^\{\-1\}\\mathbf\{C\}\_\{IR\}\.\(23\)The Schur complement𝐂cond\\mathbf\{C\}\_\{\\mathrm\{cond\}\}is the exact interaction operator acting on the retained increments after the internal increments have been eliminated\. It combines the direct retained\-to\-retained coupling𝐂R​R\\mathbf\{C\}\_\{RR\}with the indirect coupling

−𝐂R​I​𝐂I​I−1​𝐂I​R\-\\mathbf\{C\}\_\{RI\}\\mathbf\{C\}\_\{II\}^\{\-1\}\\mathbf\{C\}\_\{IR\}\(24\)transmitted through the internal nodes\. Because𝐂I​I−1\\mathbf\{C\}\_\{II\}^\{\-1\}is generally dense,𝐂cond\\mathbf\{C\}\_\{\\mathrm\{cond\}\}can couple retained nodes that were not neighbours in the original graph\.

A uniform velocity increment,Δ​𝐯i=𝐚\\Delta\\mathbf\{v\}\_\{i\}=\\mathbf\{a\}at every node, produces no relative motion and therefore no stiffness or damping contribution\. Static condensation preserves this property, so each block row of the condensed interaction operator satisfies

\(𝐂cond\)i​i\+∑j∈Rj≠i\(𝐂cond\)i​j=𝟎,\(𝐂cond\)i​i=−∑j∈Rj≠i\(𝐂cond\)i​j\.\(\\mathbf\{C\}\_\{\\mathrm\{cond\}\}\)\_\{ii\}\+\\sum\_\{\\begin\{subarray\}\{c\}j\\in R\\\\ j\\neq i\\end\{subarray\}\}\(\\mathbf\{C\}\_\{\\mathrm\{cond\}\}\)\_\{ij\}=\\mathbf\{0\},\\qquad\(\\mathbf\{C\}\_\{\\mathrm\{cond\}\}\)\_\{ii\}=\-\\sum\_\{\\begin\{subarray\}\{c\}j\\in R\\\\ j\\neq i\\end\{subarray\}\}\(\\mathbf\{C\}\_\{\\mathrm\{cond\}\}\)\_\{ij\}\.\(25\)Consequently,

\(𝐂cond​Δ​𝐯R\)i=∑j∈Rj≠i\[−\(𝐂cond\)i​j\]​\(Δ​𝐯i−Δ​𝐯j\)\.\(\\mathbf\{C\}\_\{\\mathrm\{cond\}\}\\Delta\\mathbf\{v\}\_\{R\}\)\_\{i\}=\\sum\_\{\\begin\{subarray\}\{c\}j\\in R\\\\ j\\neq i\\end\{subarray\}\}\\left\[\-\(\\mathbf\{C\}\_\{\\mathrm\{cond\}\}\)\_\{ij\}\\right\]\\left\(\\Delta\\mathbf\{v\}\_\{i\}\-\\Delta\\mathbf\{v\}\_\{j\}\\right\)\.\(26\)
The collection of pairwise coefficients−\(𝐂cond\)i​j\-\(\\mathbf\{C\}\_\{\\mathrm\{cond\}\}\)\_\{ij\}can be dense\. In the isotropic scalar case, these coefficients can be approximated using one separable factor per retained node:

−\(𝐂cond\)i​j≈wi​wj​𝐈,wi≥0\.\-\(\\mathbf\{C\}\_\{\\mathrm\{cond\}\}\)\_\{ij\}\\approx w\_\{i\}w\_\{j\}\\mathbf\{I\},\\qquad w\_\{i\}\\geq 0\.\(27\)The scalar coefficient matrix with entrieswi​wjw\_\{i\}w\_\{j\}is the rank\-one outer product𝐰𝐰⊤\\mathbf\{w\}\\mathbf\{w\}^\{\\top\}\. Thus Equation \([27](https://arxiv.org/html/2609.30344#S9.E27)\) is a rank\-one approximation of the dense pairwise coupling coefficients\. It is not an assumption that the exact Schur complement𝐂cond\\mathbf\{C\}\_\{\\mathrm\{cond\}\}is rank one\.

Substituting Equation \([27](https://arxiv.org/html/2609.30344#S9.E27)\) into Equation \([26](https://arxiv.org/html/2609.30344#S9.E26)\) gives

\(𝐂cond​Δ​𝐯R\)i\(1\)\\displaystyle\(\\mathbf\{C\}\_\{\\mathrm\{cond\}\}\\Delta\\mathbf\{v\}\_\{R\}\)^\{\(1\)\}\_\{i\}=∑j≠iwi​wj​\(Δ​𝐯i−Δ​𝐯j\)\\displaystyle=\\sum\_\{j\\neq i\}w\_\{i\}w\_\{j\}\\left\(\\Delta\\mathbf\{v\}\_\{i\}\-\\Delta\\mathbf\{v\}\_\{j\}\\right\)=wi​\[\(∑jwj\)​Δ​𝐯i−∑jwj​Δ​𝐯j\]\.\\displaystyle=w\_\{i\}\\left\[\\left\(\\sum\_\{j\}w\_\{j\}\\right\)\\Delta\\mathbf\{v\}\_\{i\}\-\\sum\_\{j\}w\_\{j\}\\Delta\\mathbf\{v\}\_\{j\}\\right\]\.\(28\)Define the shared weighted increment

Δ​𝐯¯:=∑jwj​Δ​𝐯j∑jwj\.\\overline\{\\Delta\\mathbf\{v\}\}:=\\frac\{\\sum\_\{j\}w\_\{j\}\\Delta\\mathbf\{v\}\_\{j\}\}\{\\sum\_\{j\}w\_\{j\}\}\.\(29\)Equation \([28](https://arxiv.org/html/2609.30344#S9.E28)\) becomes

\(𝐂cond​Δ​𝐯R\)i\(1\)=ci​H​\(Δ​𝐯i−Δ​𝐯¯\),ci​H:=wi​∑jwj\.\(\\mathbf\{C\}\_\{\\mathrm\{cond\}\}\\Delta\\mathbf\{v\}\_\{R\}\)^\{\(1\)\}\_\{i\}=c\_\{iH\}\\left\(\\Delta\\mathbf\{v\}\_\{i\}\-\\overline\{\\Delta\\mathbf\{v\}\}\\right\),\\qquad c\_\{iH\}:=w\_\{i\}\\sum\_\{j\}w\_\{j\}\.\(30\)This is the interaction form of a star\. Each retained node couples only to the same shared quantityΔ​𝐯¯\\overline\{\\Delta\\mathbf\{v\}\}, represented by a virtual hub\. This establishes the rank\-one virtual\-hub topology on the left\-hand side\.

#### 9\.2\.2Right\-hand side: Transfer of the internal force drive

The exact condensed right\-hand side in Equation \([20](https://arxiv.org/html/2609.30344#S9.E20)\) is

𝐫cond:=𝐫R−𝐂R​I​\(𝐌I\+𝐂I​I\)−1​𝐫I\.\\mathbf\{r\}\_\{\\mathrm\{cond\}\}:=\\mathbf\{r\}\_\{R\}\-\\mathbf\{C\}\_\{RI\}\\left\(\\mathbf\{M\}\_\{I\}\+\\mathbf\{C\}\_\{II\}\\right\)^\{\-1\}\\mathbf\{r\}\_\{I\}\.\(31\)In the static\-condensation limit, this becomes

𝐫cond=𝐫R−𝐂R​I​𝐂I​I−1​𝐫I\.\\mathbf\{r\}\_\{\\mathrm\{cond\}\}=\\mathbf\{r\}\_\{R\}\-\\mathbf\{C\}\_\{RI\}\\mathbf\{C\}\_\{II\}^\{\-1\}\\mathbf\{r\}\_\{I\}\.\(32\)The term𝐫R\\mathbf\{r\}\_\{R\}is the force drive already acting on the retained degrees of freedom\. The second term transfers the force drive from the eliminated internal degrees of freedom into the retained equations\. Because𝐂I​I−1\\mathbf\{C\}\_\{II\}^\{\-1\}is generally dense, a force acting at one internal node can contribute to the effective right\-hand side of many retained nodes\.

The corresponding transfer can be written explicitly for the rank\-one star\. Letci​Hc\_\{iH\}denote the coupling between retained nodeiiand the hub, and let𝐫H\\mathbf\{r\}\_\{H\}denote the force drive associated with the condensed internal degrees of freedom\. The hub equation is

∑ici​H​\(Δ​𝐯H−Δ​𝐯i\)=𝐫H\.\\sum\_\{i\}c\_\{iH\}\\left\(\\Delta\\mathbf\{v\}\_\{H\}\-\\Delta\\mathbf\{v\}\_\{i\}\\right\)=\\mathbf\{r\}\_\{H\}\.\(33\)Solving for the hub increment gives

Δ​𝐯H=∑ici​H​Δ​𝐯i\+𝐫H∑ici​H\.\\Delta\\mathbf\{v\}\_\{H\}=\\frac\{\\sum\_\{i\}c\_\{iH\}\\Delta\\mathbf\{v\}\_\{i\}\+\\mathbf\{r\}\_\{H\}\}\{\\sum\_\{i\}c\_\{iH\}\}\.\(34\)Substitution into the retained\-node equation gives

ci​H​\[Δ​𝐯i−∑jcj​H​Δ​𝐯j∑jcj​H\]=𝐫i\+ci​H∑jcj​H​𝐫H\.c\_\{iH\}\\left\[\\Delta\\mathbf\{v\}\_\{i\}\-\\frac\{\\sum\_\{j\}c\_\{jH\}\\Delta\\mathbf\{v\}\_\{j\}\}\{\\sum\_\{j\}c\_\{jH\}\}\\right\]=\\mathbf\{r\}\_\{i\}\+\\frac\{c\_\{iH\}\}\{\\sum\_\{j\}c\_\{jH\}\}\\mathbf\{r\}\_\{H\}\.\(35\)The left\-hand side is the star interaction derived previously\. The right\-hand side shows that the force drive associated with the condensed internal degrees of freedom is distributed among the retained nodes according to their hub couplings\. Nodeiireceives the fraction

ci​H∑jcj​H​𝐫H\.\\frac\{c\_\{iH\}\}\{\\sum\_\{j\}c\_\{jH\}\}\\mathbf\{r\}\_\{H\}\.\(36\)These fractions preserve the complete internal drive:

∑ici​H∑jcj​H​𝐫H=𝐫H\.\\sum\_\{i\}\\frac\{c\_\{iH\}\}\{\\sum\_\{j\}c\_\{jH\}\}\\mathbf\{r\}\_\{H\}=\\mathbf\{r\}\_\{H\}\.\(37\)
If the hub represents only internal redistribution and carries no external drive, then𝐫H=𝟎\\mathbf\{r\}\_\{H\}=\\mathbf\{0\}\. In this case, the hub introduces no additional net force or linear impulse into the retained system\. This zero\-resultant condition does not by itself establish an energy balance, which would additionally require a condition on the total power of the transferred forces\.

The left\-hand\-side derivation therefore establishes the effective star operator, while the right\-hand\-side derivation shows how the same star transfers a force drive associated with condensed internal degrees of freedom\.

### 9\.3Implementation of the Virtual Hub inNewmark\-β\\beta\-DGN: A Rank\-One\-Motivated Approximation of Schur\-Complement\-Induced Non\-Local Coupling

The preceding derivation is a condensation thought experiment\.Newmark\-β\\beta\-DGNdoes not eliminate physical nodes, assemble𝐂cond\\mathbf\{C\}\_\{\\mathrm\{cond\}\}, compute𝐫cond\\mathbf\{r\}\_\{\\mathrm\{cond\}\}or recover the scalar factorswiw\_\{i\}\. Instead, it implements the two conclusions of the derivation directly\.

On the left\-hand side, the scalar star interaction

ci​H​\(Δ​𝐯i−Δ​𝐯H\)c\_\{iH\}\\left\(\\Delta\\mathbf\{v\}\_\{i\}\-\\Delta\\mathbf\{v\}\_\{H\}\\right\)\(38\)motivates an operator\-weighted virtual hub\.Newmark\-β\\beta\-DGNreplaces the scalar couplingci​Hc\_\{iH\}with learned stiffness and damping operators and computes the shared hub state from operator\-weighted equilibrium\.

On the right\-hand side, the force\-transfer result motivates learned force and angular\-flux messages carried by the virtual edges\. These messages make the hub an internal relay rather than an external actuator\.

#### 9\.3\.1Left\-hand\-side implementation: Operator\-weighted hub coordinates

The stiffness and damping channels act on different mechanical quantities\. At the level of the linearised edge contribution,

δ𝐟iK=−∑j∈𝒩⁡\(i\)𝐊i​j\(δ𝐱i−δ𝐱j\),δ𝐟iD=−∑j∈𝒩⁡\(i\)𝐃i​j\(δ𝐯i−δ𝐯j\)\.\\delta\\mathbf\{f\}^\{K\}\_\{i\}=\-\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\mathbf\{K\}\_\{ij\}\\left\(\\delta\\mathbf\{x\}\_\{i\}\-\\delta\\mathbf\{x\}\_\{j\}\\right\),\\qquad\\delta\\mathbf\{f\}^\{D\}\_\{i\}=\-\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\mathbf\{D\}\_\{ij\}\\left\(\\delta\\mathbf\{v\}\_\{i\}\-\\delta\\mathbf\{v\}\_\{j\}\\right\)\.\(39\)Define the nodal operators

𝐊i=∑j∈𝒩⁡\(i\)𝐊i​j,𝐃i=∑j∈𝒩⁡\(i\)𝐃i​j,𝐃rot,i=∑j∈𝒩⁡\(i\)𝐃rot,i​j\.\\mathbf\{K\}\_\{i\}=\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\mathbf\{K\}\_\{ij\},\\qquad\\mathbf\{D\}\_\{i\}=\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\mathbf\{D\}\_\{ij\},\\qquad\\mathbf\{D\}\_\{\\mathrm\{rot\},i\}=\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\mathbf\{D\}\_\{\\mathrm\{rot\},ij\}\.\(40\)
The virtual hub is massless and has no independently advanced state\. Its position, velocity and spin are obtained from equilibrium in the corresponding operator channels:

𝐱H\\displaystyle\\mathbf\{x\}\_\{H\}=\(∑i𝐊i\)−1​∑i𝐊i​𝐱i,\\displaystyle=\\left\(\\sum\_\{i\}\\mathbf\{K\}\_\{i\}\\right\)^\{\-1\}\\sum\_\{i\}\\mathbf\{K\}\_\{i\}\\mathbf\{x\}\_\{i\},𝐯H\\displaystyle\\mathbf\{v\}\_\{H\}=\(∑i𝐃i\)−1​∑i𝐃i​𝐯i,\\displaystyle=\\left\(\\sum\_\{i\}\\mathbf\{D\}\_\{i\}\\right\)^\{\-1\}\\sum\_\{i\}\\mathbf\{D\}\_\{i\}\\mathbf\{v\}\_\{i\},𝝎H\\displaystyle\\bm\{\\omega\}\_\{H\}=\(∑i𝐃rot,i\)−1​∑i𝐃rot,i​𝝎i\.\\displaystyle=\\left\(\\sum\_\{i\}\\mathbf\{D\}\_\{\\mathrm\{rot\},i\}\\right\)^\{\-1\}\\sum\_\{i\}\\mathbf\{D\}\_\{\\mathrm\{rot\},i\}\\bm\{\\omega\}\_\{i\}\.\(41\)These expressions are the operator\-valued counterparts of the shared weighted quantity derived for the rank\-one star\.

For example, if𝐊i=ki​𝐈\\mathbf\{K\}\_\{i\}=k\_\{i\}\\mathbf\{I\}, then

𝐱H=∑iki​𝐱i∑iki\.\\mathbf\{x\}\_\{H\}=\\frac\{\\sum\_\{i\}k\_\{i\}\\mathbf\{x\}\_\{i\}\}\{\\sum\_\{i\}k\_\{i\}\}\.\(42\)Only when allkik\_\{i\}are equal does this expression reduce to the arithmetic centroid\. The hub is therefore an operator\-weighted mechanical state, not a geometric centroid\.

The Newmark factors multiplying the stiffness and damping channels cancel from their corresponding equilibrium equations\. The hub position, velocity and spin are recomputed at every substep and carry no state between substeps\. The implementation symmetrises each summed operator and adds10−6​𝐈10^\{\-6\}\\mathbf\{I\}before the three3×33\\times 3solves\.

#### 9\.3\.2Right\-hand\-side implementation: Balanced virtual\-edge interactions

The right\-hand\-side derivation shows that the star acts as an internal relay: it transfers a force drive associated with condensed internal degrees of freedom to the retained degrees of freedom\. The virtual hub inNewmark\-β\\beta\-DGNplays the same relay role\. It is not an external actuator and carries no prescribed external load\.

The model does not explicitly compute𝐫H\\mathbf\{r\}\_\{H\}or the exact Schur\-complement force transfer\. Instead, it decodes force and angular\-flux contributions on the virtual edges\. Because these quantities represent internal transfer through an unloaded hub, they are projected to have zero sum:

𝐟~H→j=𝐟H→j−1nH​∑k𝐟H→k,𝐀~H→j=𝐀H→j−1nH​∑k𝐀H→k\.\\widetilde\{\\mathbf\{f\}\}\_\{H\\to j\}=\\mathbf\{f\}\_\{H\\to j\}\-\\frac\{1\}\{n\_\{H\}\}\\sum\_\{k\}\\mathbf\{f\}\_\{H\\to k\},\\qquad\\widetilde\{\\mathbf\{A\}\}\_\{H\\to j\}=\\mathbf\{A\}\_\{H\\to j\}\-\\frac\{1\}\{n\_\{H\}\}\\sum\_\{k\}\\mathbf\{A\}\_\{H\\to k\}\.\(43\)Therefore,

∑j𝐟~H→j=𝟎,∑j𝐀~H→j=𝟎\.\\sum\_\{j\}\\widetilde\{\\mathbf\{f\}\}\_\{H\\to j\}=\\mathbf\{0\},\\qquad\\sum\_\{j\}\\widetilde\{\\mathbf\{A\}\}\_\{H\\to j\}=\\mathbf\{0\}\.\(44\)The first identity ensures that the hub redistributes force among the physical nodes without adding a net linear resultant\. This is the model counterpart of the unloaded theoretical hub,𝐫H=𝟎\\mathbf\{r\}\_\{H\}=\\mathbf\{0\}\. The second identity imposes the corresponding zero\-sum condition on the angular\-flux coordinates\.

Because different virtual edges can use different reference points, zero net force does not necessarily imply zero net moment\. The virtual\-edge forces can retain the residual couple

𝐓H=∑j𝐫H​j0×𝐟~H→j\.\\mathbf\{T\}\_\{H\}=\\sum\_\{j\}\\mathbf\{r\}^\{0\}\_\{Hj\}\\times\\widetilde\{\\mathbf\{f\}\}\_\{H\\to j\}\.\(45\)The corresponding linear\-momentum, angular\-momentum and energy properties must therefore be stated separately\.

#### 9\.3\.3Computational cost

The hub adds2​N2Ndirected virtual edges and three3×33\\times 3operator\-weighted solves per substep\. Its additional cost is linear in the number of physical nodes\. The augmented graph has diameter two, while the physical\-edge message passing and independent nodal Newmark solves remain unchanged\.

## 10Virtual\-Hub Coupling and the Semi\-Implicit Nodal Newmark Update inNewmark\-β\\beta\-DGN

The preceding section established the virtual hub as a tractable representation of non\-local operator coupling and internal force transfer\. We now describe howNewmark\-β\\beta\-DGNcombines physical\-edge and virtual\-edge messages with independent nodal Newmark solves to advance the physical state\. The connection to classical Newmark integration identifies the numerical template; it does not make the learned update a classical global implicit solve\.

### 10\.1One learned mechanical substep

One observed transition of durationΔ​t\\Delta tis divided intoSSsubsteps of sizeδ​t=Δ​t/S\\delta t=\\Delta t/S\. At the beginning of each substep, the network decodes the net force acting on physical nodeii,

𝐛i=𝐟iext\+∑j∈𝒩⁡\(i\)𝐟i​j,\\mathbf\{b\}\_\{i\}=\\mathbf\{f\}^\{\\mathrm\{ext\}\}\_\{i\}\+\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\mathbf\{f\}\_\{ij\},\(46\)where𝒩⁡\(i\)\\mathcal\{N\}\(i\)includes both physical neighbours and the virtual hub\. Thus physical\-edge and virtual\-edge forces contribute directly to the right\-hand side of the nodal update\.

Each edge also carries a position\-response operator𝐊i​j\\mathbf\{K\}\_\{ij\}and a velocity\-response operator𝐃i​j\\mathbf\{D\}\_\{ij\}\. Their nodal sums are

𝐊i=∑j∈𝒩⁡\(i\)𝐊i​j,𝐃i=∑j∈𝒩⁡\(i\)𝐃i​j\.\\mathbf\{K\}\_\{i\}=\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\mathbf\{K\}\_\{ij\},\\qquad\\mathbf\{D\}\_\{i\}=\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\mathbf\{D\}\_\{ij\}\.\(47\)The virtual\-edge operators therefore contribute to the implicit mechanical response on the left\-hand side, while the virtual\-edge forces contribute to the drive on the right\-hand side\. This mirrors the operator and force\-transfer roles derived for the hub in Section[9](https://arxiv.org/html/2609.30344#S9)\.

The matrices𝐊i\\mathbf\{K\}\_\{i\}and𝐃i\\mathbf\{D\}\_\{i\}are learned, state\-dependent response coefficients\. They are not obtained by differentiating the decoded forces and are not assumed to be physical material tangents\. The inverse mass is likewise learned and parameterised as the positive isotropic operator

𝐌i−1=mi−1​𝐈,mi\>0\.\\mathbf\{M\}\_\{i\}^\{\-1\}=m\_\{i\}^\{\-1\}\\mathbf\{I\},\\qquad m\_\{i\}\>0\.\(48\)
Define the nodal increments

Δ​𝐯i:=𝐯i,t\+1−𝐯i,t,Δ​𝐱i:=𝐱i,t\+1−𝐱i,t\.\\Delta\\mathbf\{v\}\_\{i\}:=\\mathbf\{v\}\_\{i,t\+1\}\-\\mathbf\{v\}\_\{i,t\},\\qquad\\Delta\\mathbf\{x\}\_\{i\}:=\\mathbf\{x\}\_\{i,t\+1\}\-\\mathbf\{x\}\_\{i,t\}\.\(49\)The decoded operators represent the change in mechanical response over the substep as

Δ​𝐟iresp=−𝐃i​Δ​𝐯i−𝐊i​Δ​𝐱i\.\\Delta\\mathbf\{f\}^\{\\mathrm\{resp\}\}\_\{i\}=\-\\mathbf\{D\}\_\{i\}\\Delta\\mathbf\{v\}\_\{i\}\-\\mathbf\{K\}\_\{i\}\\Delta\\mathbf\{x\}\_\{i\}\.\(50\)Equation \([50](https://arxiv.org/html/2609.30344#S10.E50)\) defines how the learned operators enter the update\. It does not identify them as Jacobians of𝐛i\\mathbf\{b\}\_\{i\}\.

The velocity increment is determined from the trapezoidal response balance

𝐌i​Δ​𝐯i=𝐛i​δ​t−12​δ​t​𝐃i​Δ​𝐯i−12​δ​t​𝐊i​Δ​𝐱i\.\\mathbf\{M\}\_\{i\}\\Delta\\mathbf\{v\}\_\{i\}=\\mathbf\{b\}\_\{i\}\\delta t\-\\tfrac\{1\}\{2\}\\delta t\\,\\mathbf\{D\}\_\{i\}\\Delta\\mathbf\{v\}\_\{i\}\-\\tfrac\{1\}\{2\}\\delta t\\,\\mathbf\{K\}\_\{i\}\\Delta\\mathbf\{x\}\_\{i\}\.\(51\)The position increment is evaluated from the mean velocity,

Δ​𝐱i=12​δ​t​\(𝐯i,t\+𝐯i,t\+1\)=δ​t​\(𝐯i,t\+12​Δ​𝐯i\)\.\\Delta\\mathbf\{x\}\_\{i\}=\\tfrac\{1\}\{2\}\\delta t\\left\(\\mathbf\{v\}\_\{i,t\}\+\\mathbf\{v\}\_\{i,t\+1\}\\right\)=\\delta t\\left\(\\mathbf\{v\}\_\{i,t\}\+\\tfrac\{1\}\{2\}\\Delta\\mathbf\{v\}\_\{i\}\\right\)\.\(52\)Substituting Equation \([52](https://arxiv.org/html/2609.30344#S10.E52)\) into Equation \([51](https://arxiv.org/html/2609.30344#S10.E51)\), collecting the terms inΔ​𝐯i\\Delta\\mathbf\{v\}\_\{i\}, and multiplying by𝐌i−1\\mathbf\{M\}\_\{i\}^\{\-1\}gives the update implemented byNewmark\-β\\beta\-DGN:

\[𝐈\+12​δ​t​𝐌i−1​𝐃i\+14​δ​t2​𝐌i−1​𝐊i\]​Δ​𝐯i=𝐌i−1​𝐛i​δ​t−12​δ​t2​𝐌i−1​𝐊i​𝐯i,t\.\\boxed\{\\left\[\\mathbf\{I\}\+\\tfrac\{1\}\{2\}\\delta t\\,\\mathbf\{M\}\_\{i\}^\{\-1\}\\mathbf\{D\}\_\{i\}\+\\tfrac\{1\}\{4\}\\delta t^\{2\}\\mathbf\{M\}\_\{i\}^\{\-1\}\\mathbf\{K\}\_\{i\}\\right\]\\Delta\\mathbf\{v\}\_\{i\}=\\mathbf\{M\}\_\{i\}^\{\-1\}\\mathbf\{b\}\_\{i\}\\delta t\-\\tfrac\{1\}\{2\}\\delta t^\{2\}\\mathbf\{M\}\_\{i\}^\{\-1\}\\mathbf\{K\}\_\{i\}\\mathbf\{v\}\_\{i,t\}\}\.\(53\)The state is then advanced by

𝐯i,t\+1=𝐯i,t\+Δ​𝐯i,𝐱i,t\+1=𝐱i,t\+Δ​𝐱i\.\\mathbf\{v\}\_\{i,t\+1\}=\\mathbf\{v\}\_\{i,t\}\+\\Delta\\mathbf\{v\}\_\{i\},\\qquad\\mathbf\{x\}\_\{i,t\+1\}=\\mathbf\{x\}\_\{i,t\}\+\\Delta\\mathbf\{x\}\_\{i\}\.\(54\)
The physical and virtual messages make𝐛i\\mathbf\{b\}\_\{i\},𝐊i\\mathbf\{K\}\_\{i\}and𝐃i\\mathbf\{D\}\_\{i\}dependent on the current state of the complete graph\. Once these quantities have been decoded, Equation \([53](https://arxiv.org/html/2609.30344#S10.E53)\) is solved independently at each free physical node as a3×33\\times 3system\.

### 10\.2Connection to average\-acceleration Newmark integration and the meaning of semi\-implicit

For the average\-acceleration Newmark parametersγ=12\\gamma=\\tfrac\{1\}\{2\}andβ=14\\beta=\\tfrac\{1\}\{4\}, the classical velocity and position increments satisfy

Δ​𝐯=12​δ​t​\(𝐚t\+𝐚t\+1\),Δ​𝐱=δ​t​𝐯t\+12​δ​t​Δ​𝐯\.\\Delta\\mathbf\{v\}=\\tfrac\{1\}\{2\}\\delta t\\left\(\\mathbf\{a\}\_\{t\}\+\\mathbf\{a\}\_\{t\+1\}\\right\),\\qquad\\Delta\\mathbf\{x\}=\\delta t\\mathbf\{v\}\_\{t\}\+\\tfrac\{1\}\{2\}\\delta t\\Delta\\mathbf\{v\}\.\(55\)The second relation is exactly the mean\-velocity position rule in Equation \([52](https://arxiv.org/html/2609.30344#S10.E52)\)\.

For the classical linear system

𝐌​𝐱¨\+𝐃​𝐱˙\+𝐊𝐱=𝐟ext,\\mathbf\{M\}\\ddot\{\\mathbf\{x\}\}\+\\mathbf\{D\}\\dot\{\\mathbf\{x\}\}\+\\mathbf\{K\}\\mathbf\{x\}=\\mathbf\{f\}\_\{\\mathrm\{ext\}\},\(56\)the internal\-force change over one step is

−𝐃​Δ​𝐯−𝐊​Δ​𝐱\.\-\\mathbf\{D\}\\Delta\\mathbf\{v\}\-\\mathbf\{K\}\\Delta\\mathbf\{x\}\.\(57\)The average\-acceleration rule therefore produces the coefficientsδ​t/2\\delta t/2andδ​t2/4\\delta t^\{2\}/4that appear in Equation \([53](https://arxiv.org/html/2609.30344#S10.E53)\)\. This correspondence explains the numerical form of the learned update\. It does not identify the learned operators𝐊i\\mathbf\{K\}\_\{i\}and𝐃i\\mathbf\{D\}\_\{i\}with the physical matrices𝐊\\mathbf\{K\}and𝐃\\mathbf\{D\}\.

The learned update is implicit in each nodal increment becauseΔ​𝐯i\\Delta\\mathbf\{v\}\_\{i\}and the associated position increment enter the response being solved\. It is nevertheless not a classical global implicit Newmark step\. InNewmark\-β\\beta\-DGN,

1. 1\.the response operators are decoded from the current graph state rather than obtained by differentiating a prescribed mechanical residual;
2. 2\.one learned solve is performed per substep, without an inner Newton iteration; and
3. 3\.the global3​N×3​N3N\\times 3Nsolve is replaced by independent nodal solves, while message passing and the operator\-weighted hub supply non\-local information\.

The term*semi\-implicit*refers to this separation: each nodal response is treated implicitly, but the complete off\-diagonal mechanical system is not assembled or solved exactly\.

### 10\.3Response\-operator construction and nodal solvability

Each edge decoder emits six scalars and forms the lower\-triangular matrix

𝐋=\[d000s1d10s3s4d2\],\(d0,d1,d2\)=softplus⁡\(s0,s2,s5\)\+10−4\.\\mathbf\{L\}=\\begin\{bmatrix\}d\_\{0\}&0&0\\\\ s\_\{1\}&d\_\{1\}&0\\\\ s\_\{3\}&s\_\{4\}&d\_\{2\}\\end\{bmatrix\},\\qquad\(d\_\{0\},d\_\{1\},d\_\{2\}\)=\\operatorname\{softplus\}\(s\_\{0\},s\_\{2\},s\_\{5\}\)\+10^\{\-4\}\.\(58\)If𝐑i​j\\mathbf\{R\}\_\{ij\}is the edge\-local orthonormal frame, the corresponding operator in global coordinates is

𝐖i​j=𝐑i​j​𝐋𝐋⊤​𝐑i​j⊤\.\\mathbf\{W\}\_\{ij\}=\\mathbf\{R\}\_\{ij\}\\mathbf\{L\}\\mathbf\{L\}^\{\\top\}\\mathbf\{R\}\_\{ij\}^\{\\top\}\.\(59\)The positive diagonal of𝐋\\mathbf\{L\}makes𝐋𝐋⊤\\mathbf\{L\}\\mathbf\{L\}^\{\\top\}positive definite\. Under a rotation𝐐\\mathbf\{Q\},

𝐑i​j↦𝐐𝐑i​j,𝐖i​j↦𝐐𝐖i​j​𝐐⊤\.\\mathbf\{R\}\_\{ij\}\\mapsto\\mathbf\{Q\}\\mathbf\{R\}\_\{ij\},\\qquad\\mathbf\{W\}\_\{ij\}\\mapsto\\mathbf\{Q\}\\mathbf\{W\}\_\{ij\}\\mathbf\{Q\}^\{\\top\}\.\(60\)The edge operators, their nodal sums and the solution of Equation \([53](https://arxiv.org/html/2609.30344#S10.E53)\) are therefore rotation\-equivariant\.

Positive definiteness alone does not imply symmetry between the two directed traversals of an edge\. The equality𝐖j​i=𝐖i​j\\mathbf\{W\}\_\{ji\}=\\mathbf\{W\}\_\{ij\}follows from the interaction construction: the decoder coefficients are shared between the two traversals, and reversal of the local frame leaves the quadratic operator unchanged\.

For the isotropic inverse mass𝐌i−1=mi−1​𝐈\\mathbf\{M\}\_\{i\}^\{\-1\}=m\_\{i\}^\{\-1\}\\mathbf\{I\}, the coefficient matrix in the nodal solve is

𝐀i=𝐈\+12​δ​t​mi−1​𝐃i\+14​δ​t2​mi−1​𝐊i\.\\mathbf\{A\}\_\{i\}=\\mathbf\{I\}\+\\tfrac\{1\}\{2\}\\delta t\\,m\_\{i\}^\{\-1\}\\mathbf\{D\}\_\{i\}\+\\tfrac\{1\}\{4\}\\delta t^\{2\}m\_\{i\}^\{\-1\}\\mathbf\{K\}\_\{i\}\.\(61\)Because𝐊i\\mathbf\{K\}\_\{i\}and𝐃i\\mathbf\{D\}\_\{i\}are symmetric positive semidefinite,𝐀i\\mathbf\{A\}\_\{i\}is symmetric positive definite for everyδ​t\>0\\delta t\>0, with

λmin​\(𝐀i\)≥1,‖𝐀i−1‖2≤1\.\\lambda\_\{\\min\}\(\\mathbf\{A\}\_\{i\}\)\\geq 1,\\qquad\\left\\lVert\\mathbf\{A\}\_\{i\}^\{\-1\}\\right\\rVert\_\{2\}\\leq 1\.\(62\)Thus every nodal3×33\\times 3solve remains invertible as the step size or operator magnitude increases\. This guarantees solvability of an individual nodal system, not stability of the complete learned rollout\. The implementation adds10−8​𝐈10^\{\-8\}\\mathbf\{I\}as a floating\-point safeguard\.

The update is invariant to a uniform translation of the positions because the graph construction uses relative coordinates\. It is not invariant to a Galilean velocity shift: the term

−12​δ​t2​𝐌i−1​𝐊i​𝐯i,t\-\\tfrac\{1\}\{2\}\\delta t^\{2\}\\mathbf\{M\}\_\{i\}^\{\-1\}\\mathbf\{K\}\_\{i\}\\mathbf\{v\}\_\{i,t\}\(63\)changes under𝐯i,t↦𝐯i,t\+𝐜\\mathbf\{v\}\_\{i,t\}\\mapsto\\mathbf\{v\}\_\{i,t\}\+\\mathbf\{c\}\.

### 10\.4Classical stability correspondence and its scope

The stability statement concerns only the classical linear test problem\. Setm=1m=1,D=0D=0andK=ω2K=\\omega^\{2\}in Equation \([56](https://arxiv.org/html/2609.30344#S10.E56)\), and set the learned coefficients in Equation \([53](https://arxiv.org/html/2609.30344#S10.E53)\) equal to these exact linear coefficients\. Under this controlled substitution, the nodal update reduces to the classical average\-acceleration Newmark method\.

Let

Ω=δ​t​ω,a=4​Ω24\+Ω2\.\\Omega=\\delta t\\omega,\\qquad a=\\frac\{4\\Omega^\{2\}\}\{4\+\\Omega^\{2\}\}\.\(64\)The amplification matrix carrying\(xt,δ​t​vt\)\(x\_\{t\},\\delta tv\_\{t\}\)to\(xt\+1,δ​t​vt\+1\)\(x\_\{t\+1\},\\delta tv\_\{t\+1\}\)is

𝐆⁡\(Ω\)=\[1−a/21−a/4−a1−a/2\]\.\\mathbf\{G\}\(\\Omega\)=\\begin\{bmatrix\}1\-a/2&1\-a/4\\\\ \-a&1\-a/2\\end\{bmatrix\}\.\(65\)For every finiteΩ\>0\\Omega\>0,

det𝐆⁡\(Ω\)=1,ρ⁡\(𝐆⁡\(Ω\)\)=1\.\\det\\mathbf\{G\}\(\\Omega\)=1,\\qquad\\rho\\\!\\left\(\\mathbf\{G\}\(\\Omega\)\\right\)=1\.\(66\)The classical undamped update is therefore unconditionally stable and introduces no algorithmic decay\. It remains subject to phase error:

θ=Ω⁡\(1−112​Ω2\+𝒪⁡\(Ω4\)\),TnumT=1\+112​Ω2\+𝒪⁡\(Ω4\)\.\\theta=\\Omega\\left\(1\-\\tfrac\{1\}\{12\}\\Omega^\{2\}\+\\mathcal\{O\}\(\\Omega^\{4\}\)\\right\),\\qquad\\frac\{T\_\{\\mathrm\{num\}\}\}\{T\}=1\+\\tfrac\{1\}\{12\}\\Omega^\{2\}\+\\mathcal\{O\}\(\\Omega^\{4\}\)\.\(67\)Thus a large step can remain bounded while accumulating phase error\. The classical method is second\-order accurate for this linear problem\.

These standard properties do not establish unconditional stability ofNewmark\-β\\beta\-DGN\. In the trained model, the forces and response operators depend on the predicted state, the operators are not force derivatives, and the global off\-diagonal system is replaced by local solves, message passing and the virtual hub\. Stability of the complete rollout is therefore evaluated empirically through the beam rollouts, the explicit ablation and the stiffness\-dependent comparisons in Section[4\.7](https://arxiv.org/html/2609.30344#S4.SS7a)\.

### 10\.5Angular update

The angular state uses the same algebraic update:

\[𝐈\+12​δ​t​𝐈i−1​𝐃rot,i\+14​δ​t2​𝐈i−1​𝐊rot,i\]​Δ​𝝎i=𝐈i−1​𝝉i​δ​t−12​δ​t2​𝐈i−1​𝐊rot,i​𝝎i,t\.\\left\[\\mathbf\{I\}\+\\tfrac\{1\}\{2\}\\delta t\\,\\mathbf\{I\}\_\{i\}^\{\-1\}\\mathbf\{D\}\_\{\\mathrm\{rot\},i\}\+\\tfrac\{1\}\{4\}\\delta t^\{2\}\\mathbf\{I\}\_\{i\}^\{\-1\}\\mathbf\{K\}\_\{\\mathrm\{rot\},i\}\\right\]\\Delta\\bm\{\\omega\}\_\{i\}=\\mathbf\{I\}\_\{i\}^\{\-1\}\\bm\{\\tau\}\_\{i\}\\delta t\-\\tfrac\{1\}\{2\}\\delta t^\{2\}\\mathbf\{I\}\_\{i\}^\{\-1\}\\mathbf\{K\}\_\{\\mathrm\{rot\},i\}\\bm\{\\omega\}\_\{i,t\}\.\(68\)Here𝐈i−1\\mathbf\{I\}\_\{i\}^\{\-1\}is a learned positive isotropic inverse inertia,𝝉i\\bm\{\\tau\}\_\{i\}is the aggregated decoded torque, and𝐊rot,i\\mathbf\{K\}\_\{\\mathrm\{rot\},i\}and𝐃rot,i\\mathbf\{D\}\_\{\\mathrm\{rot\},i\}are learned rotational response operators\. The angular state is latent and receives no direct supervision\.

### 10\.6External forcing and constrained nodes

For the beam experiments, the applied load is observed\. Its known values across the outer interval are interpolated over theSSsubsteps and held constant within each substep\. For motion capture, protein dynamics and walking biomechanics, external forcing is unobserved and decoded from the current predicted state at the beginning of each substep\. The decoded value is held constant during the corresponding nodal solve, and no future forcing is supplied\.

Equation \([53](https://arxiv.org/html/2609.30344#S10.E53)\) is applied only to free physical nodes\. Clamped degrees of freedom are prescribed\. The virtual hub is not integrated as a physical degree of freedom; its position, velocity and spin are recomputed from the current physical state and learned operators at every substep\.

### 10\.7Explicit limit

As

𝐊i,𝐃i→𝟎,\\mathbf\{K\}\_\{i\},\\mathbf\{D\}\_\{i\}\\to\\mathbf\{0\},\(69\)Equation \([53](https://arxiv.org/html/2609.30344#S10.E53)\) reduces to

Δ​𝐯i=𝐌i−1​𝐛i​δ​t,Δ​𝐱i=12​δ​t​\(𝐯i,t\+𝐯i,t\+1\)\.\\Delta\\mathbf\{v\}\_\{i\}=\\mathbf\{M\}\_\{i\}^\{\-1\}\\mathbf\{b\}\_\{i\}\\delta t,\\qquad\\Delta\\mathbf\{x\}\_\{i\}=\\tfrac\{1\}\{2\}\\delta t\\left\(\\mathbf\{v\}\_\{i,t\}\+\\mathbf\{v\}\_\{i,t\+1\}\\right\)\.\(70\)This limit integrates the decoded force explicitly while retaining the trapezoidal position rule\. Its momentum properties depend on whether the contributing fluxes arise from antisymmetric physical\-edge exchanges or from the projected virtual\-hub interactions\. Section[11](https://arxiv.org/html/2609.30344#S11)treats these cases separately\.

## 11Momentum Balance of the Decoded Interactions and the Residual of the Independent Nodal Update

The decoded interactions and the numerical update have distinct conservation properties\. Physical\-edge interactions exchange linear and angular momentum pairwise, while the virtual hub enforces collective balance over its incident virtual edges\. These properties hold before time integration\.

The independent semi\-implicit nodal update preserves this interaction\-level balance only up to a finite\-substep linear\-momentum residual\. The residual is not an additional force produced by the decoder or the virtual hub\. It arises when the balanced decoded drive is processed by the learned response operators independently at each physical node\. Under the boundedness assumptions stated below, the absolute residual over one substep is𝒪⁡\(δ​t2\)\\mathcal\{O\}\(\\delta t^\{2\}\)\.

The interaction\-level analysis concerns linear and angular momentum\. The numerical\-update analysis concerns total linear momentum\. No energy\-conservation claim is made, because zero net force does not by itself imply zero mechanical power\.

For the dynamically integrated physical nodes𝒱dyn\\mathcal\{V\}\_\{\\mathrm\{dyn\}\}, define the linear\-momentum increment

Δ​𝐏:=𝐏t\+1−𝐏t=∑i∈𝒱dyn𝐌i​Δ​𝐯i\.\\Delta\\mathbf\{P\}:=\\mathbf\{P\}\_\{t\+1\}\-\\mathbf\{P\}\_\{t\}=\\sum\_\{i\\in\\mathcal\{V\}\_\{\\mathrm\{dyn\}\}\}\\mathbf\{M\}\_\{i\}\\Delta\\mathbf\{v\}\_\{i\}\.\(71\)The virtual hub is massless and is not included in this sum\. The total external force acting on the integrated nodes is

𝐅ext:=∑i∈𝒱dyn𝐟iext\.\\mathbf\{F\}\_\{\\mathrm\{ext\}\}:=\\sum\_\{i\\in\\mathcal\{V\}\_\{\\mathrm\{dyn\}\}\}\\mathbf\{f\}^\{\\mathrm\{ext\}\}\_\{i\}\.\(72\)Reactions associated with prescribed degrees of freedom must be included in the external impulse\.

### 11\.1Balance of the Decoded Momentum Exchanges

#### 11\.1\.1Pairwise Balance of Physical\-Edge Interactions

Consider a physical edge\(i,j\)\(i,j\)\. Its two directed traversals satisfy

𝐟i​j=−𝐟j​i,𝐀i​j=−𝐀j​i,𝐫i​j0=𝐫j​i0\.\\mathbf\{f\}\_\{ij\}=\-\\mathbf\{f\}\_\{ji\},\\qquad\\mathbf\{A\}\_\{ij\}=\-\\mathbf\{A\}\_\{ji\},\\qquad\\mathbf\{r\}^\{0\}\_\{ij\}=\\mathbf\{r\}^\{0\}\_\{ji\}\.\(73\)Here𝐟i​j\\mathbf\{f\}\_\{ij\}is the force transmitted fromjjtoii,𝐀i​j\\mathbf\{A\}\_\{ij\}is the decoded angular\-momentum flux, and𝐫i​j0\\mathbf\{r\}^\{0\}\_\{ij\}is the application point shared by the two traversals\. The spin torque delivered to nodeiiis

𝝉i​j=𝐀i​j−\(𝐫i−𝐫i​j0\)×𝐟i​j\.\\bm\{\\tau\}\_\{ij\}=\\mathbf\{A\}\_\{ij\}\-\\left\(\\mathbf\{r\}\_\{i\}\-\\mathbf\{r\}^\{0\}\_\{ij\}\\right\)\\times\\mathbf\{f\}\_\{ij\}\.\(74\)
The two forces have zero resultant:

𝐟i​j\+𝐟j​i=𝟎\.\\mathbf\{f\}\_\{ij\}\+\\mathbf\{f\}\_\{ji\}=\\mathbf\{0\}\.\(75\)
Angular momentum contains an orbital contribution from the force and an intrinsic contribution from the spin torque\. About an arbitrary origin𝐨\\mathbf\{o\}, the combined contribution of the two traversals is

\(𝐫i−𝐨\)×𝐟i​j\+𝝉i​j\+\(𝐫j−𝐨\)×𝐟j​i\+𝝉j​i\\displaystyle\\left\(\\mathbf\{r\}\_\{i\}\-\\mathbf\{o\}\\right\)\\times\\mathbf\{f\}\_\{ij\}\+\\bm\{\\tau\}\_\{ij\}\+\\left\(\\mathbf\{r\}\_\{j\}\-\\mathbf\{o\}\\right\)\\times\\mathbf\{f\}\_\{ji\}\+\\bm\{\\tau\}\_\{ji\}=\(𝐫i​j0−𝐨\)×\(𝐟i​j\+𝐟j​i\)\+𝐀i​j\+𝐀j​i=𝟎\.\\displaystyle\\qquad=\\left\(\\mathbf\{r\}^\{0\}\_\{ij\}\-\\mathbf\{o\}\\right\)\\times\\left\(\\mathbf\{f\}\_\{ij\}\+\\mathbf\{f\}\_\{ji\}\\right\)\+\\mathbf\{A\}\_\{ij\}\+\\mathbf\{A\}\_\{ji\}=\\mathbf\{0\}\.\(76\)Each physical edge therefore exchanges linear and angular momentum without creating either\.

#### 11\.1\.2Collective Balance of Virtual\-Hub Interactions

The virtual hub is a massless internal relay\. It redistributes interactions among the physical nodes, but it is not an independently integrated mechanical body and carries no prescribed external load\.

Physical\-edge balance is imposed pairwise\. Virtual\-hub balance is imposed collectively over all hub\-to\-body edges\. The projection derived in Section[9\.3\.2](https://arxiv.org/html/2609.30344#S9.SS3.SSS2)gives

∑j𝐟~H→j=𝟎,∑j𝐀~H→j=𝟎\.\\sum\_\{j\}\\widetilde\{\\mathbf\{f\}\}\_\{H\\to j\}=\\mathbf\{0\},\\qquad\\sum\_\{j\}\\widetilde\{\\mathbf\{A\}\}\_\{H\\to j\}=\\mathbf\{0\}\.\(77\)Thus the hub adds no net linear resultant, and the decoded angular\-momentum\-flux channel is also collectively balanced\.

On each virtual edge, the force and spin torque are referred to the same learned application point:

𝝉~H→j=𝐀~H→j−\(𝐫j−𝐫H​j0\)×𝐟~H→j\.\\widetilde\{\\bm\{\\tau\}\}\_\{H\\to j\}=\\widetilde\{\\mathbf\{A\}\}\_\{H\\to j\}\-\\left\(\\mathbf\{r\}\_\{j\}\-\\mathbf\{r\}^\{0\}\_\{Hj\}\\right\)\\times\\widetilde\{\\mathbf\{f\}\}\_\{H\\to j\}\.\(78\)The total angular contribution delivered to the physical nodes about an arbitrary origin𝐨\\mathbf\{o\}is therefore

ℳH​\(𝐨\)\\displaystyle\\mathcal\{M\}\_\{H\}\(\\mathbf\{o\}\):=∑j\[\(𝐫j−𝐨\)×𝐟~H→j\+𝝉~H→j\]\\displaystyle:=\\sum\_\{j\}\\left\[\\left\(\\mathbf\{r\}\_\{j\}\-\\mathbf\{o\}\\right\)\\times\\widetilde\{\\mathbf\{f\}\}\_\{H\\to j\}\+\\widetilde\{\\bm\{\\tau\}\}\_\{H\\to j\}\\right\]=∑j\[\(𝐫H​j0−𝐨\)×𝐟~H→j\+𝐀~H→j\]\.\\displaystyle=\\sum\_\{j\}\\left\[\\left\(\\mathbf\{r\}^\{0\}\_\{Hj\}\-\\mathbf\{o\}\\right\)\\times\\widetilde\{\\mathbf\{f\}\}\_\{H\\to j\}\+\\widetilde\{\\mathbf\{A\}\}\_\{H\\to j\}\\right\]\.\(79\)Using Equation \([77](https://arxiv.org/html/2609.30344#S11.E77)\), this can be expressed relative to any common hub point𝐫H\\mathbf\{r\}\_\{H\}as

ℳH​\(𝐨\)=∑j\(𝐫H​j0−𝐫H\)×𝐟~H→j\.\\mathcal\{M\}\_\{H\}\(\\mathbf\{o\}\)=\\sum\_\{j\}\\left\(\\mathbf\{r\}^\{0\}\_\{Hj\}\-\\mathbf\{r\}\_\{H\}\\right\)\\times\\widetilde\{\\mathbf\{f\}\}\_\{H\\to j\}\.\(80\)The arbitrary origin has disappeared because the projected forces have zero resultant\.

Equation \([80](https://arxiv.org/html/2609.30344#S11.E80)\) vanishes when all virtual\-edge application points coincide at𝐫H\\mathbf\{r\}\_\{H\}, or more generally when their force\-weighted moment about𝐫H\\mathbf\{r\}\_\{H\}is zero\. The virtual hub therefore enforces exact collective balance of the decoded force and angular\-momentum\-flux channels\. Exact total angular\-momentum balance about a common spatial origin additionally depends on the learned application points\. Their locations relative to the operator\-weighted hub have not been measured, so this geometric contribution is not quantified here\.

### 11\.2Linear\-Momentum Residual of the Independent Nodal Update

For clarity, we reproduce the quantities entering the nodal update from Section[10\.1](https://arxiv.org/html/2609.30344#S10.SS1)\. At physical nodeii, the decoded drive and the nodal response operators are

𝐛i=𝐟iext\+∑j∈𝒩⁡\(i\)𝐟i​j,𝐊i=∑j∈𝒩⁡\(i\)𝐊i​j,𝐃i=∑j∈𝒩⁡\(i\)𝐃i​j\.\\mathbf\{b\}\_\{i\}=\\mathbf\{f\}^\{\\mathrm\{ext\}\}\_\{i\}\+\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\mathbf\{f\}\_\{ij\},\\qquad\\mathbf\{K\}\_\{i\}=\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\mathbf\{K\}\_\{ij\},\\qquad\\mathbf\{D\}\_\{i\}=\\sum\_\{j\\in\\mathcal\{N\}\(i\)\}\\mathbf\{D\}\_\{ij\}\.\(81\)The physical\-edge antisymmetry and virtual\-hub projection established above give

∑i∈𝒱dyn𝐛i=𝐅ext\.\\sum\_\{i\\in\\mathcal\{V\}\_\{\\mathrm\{dyn\}\}\}\\mathbf\{b\}\_\{i\}=\\mathbf\{F\}\_\{\\mathrm\{ext\}\}\.\(82\)
For the isotropic mass𝐌i=mi​𝐈\\mathbf\{M\}\_\{i\}=m\_\{i\}\\mathbf\{I\}, the coefficient matrix is

𝐀i:=𝐈\+12​δ​t​mi−1​𝐃i\+14​δ​t2​mi−1​𝐊i\.\\mathbf\{A\}\_\{i\}:=\\mathbf\{I\}\+\\tfrac\{1\}\{2\}\\delta t\\,m\_\{i\}^\{\-1\}\\mathbf\{D\}\_\{i\}\+\\tfrac\{1\}\{4\}\\delta t^\{2\}m\_\{i\}^\{\-1\}\\mathbf\{K\}\_\{i\}\.\(83\)The independent nodal update from Equation \([53](https://arxiv.org/html/2609.30344#S10.E53)\) is reproduced here as

𝐀i​Δ​𝐯i=mi−1​𝐛i​δ​t−12​δ​t2​mi−1​𝐊i​𝐯i,t\.\\mathbf\{A\}\_\{i\}\\Delta\\mathbf\{v\}\_\{i\}=m\_\{i\}^\{\-1\}\\mathbf\{b\}\_\{i\}\\delta t\-\\tfrac\{1\}\{2\}\\delta t^\{2\}m\_\{i\}^\{\-1\}\\mathbf\{K\}\_\{i\}\\mathbf\{v\}\_\{i,t\}\.\(84\)Multiplying bymim\_\{i\}and solving for the momentum increment gives

Δ​𝐩i=𝐀i−1​\[𝐛i​δ​t−12​δ​t2​𝐊i​𝐯i,t\]\.\\Delta\\mathbf\{p\}\_\{i\}=\\mathbf\{A\}\_\{i\}^\{\-1\}\\left\[\\mathbf\{b\}\_\{i\}\\delta t\-\\tfrac\{1\}\{2\}\\delta t^\{2\}\\mathbf\{K\}\_\{i\}\\mathbf\{v\}\_\{i,t\}\\right\]\.\(85\)Equation \([85](https://arxiv.org/html/2609.30344#S11.E85)\) contains two response\-dependent contributions to the total momentum increment\.

##### Response of a balanced force pair\.

Consider an antisymmetric physical\-edge force𝐟j​i=−𝐟i​j\\mathbf\{f\}\_\{ji\}=\-\\mathbf\{f\}\_\{ij\}\. Its direct contribution to the total momentum increment is

Δ​𝐩pair​\(i​j\)=\(𝐀i−1−𝐀j−1\)​𝐟i​j​δ​t\.\\Delta\\mathbf\{p\}\_\{\\mathrm\{pair\}\}\(ij\)=\\left\(\\mathbf\{A\}\_\{i\}^\{\-1\}\-\\mathbf\{A\}\_\{j\}^\{\-1\}\\right\)\\mathbf\{f\}\_\{ij\}\\delta t\.\(86\)The decoded forces remain equal and opposite, but their momentum increments need not cancel when the two nodal coefficient matrices differ\. These matrices depend on the learned masses and on the nodal position\- and velocity\-response operators\.

The same reasoning applies to the collectively balanced virtual\-hub forces\. Their sum is zero before integration, but multiplication by different𝐀i−1\\mathbf\{A\}\_\{i\}^\{\-1\}matrices need not preserve that cancellation\. The virtual\-hub projection remains balanced; the finite\-substep residual is introduced afterwards by the independent nodal response\.

##### Known\-velocity position response\.

Equation \([85](https://arxiv.org/html/2609.30344#S11.E85)\) also contains the total contribution

Δ𝐏pos=−12δt2∑i𝐀i−1𝐊i𝐯i,t\.\\Delta\\mathbf\{P\}^\{\\,\\mathrm\{pos\}\}=\-\\tfrac\{1\}\{2\}\\delta t^\{2\}\\sum\_\{i\}\\mathbf\{A\}\_\{i\}^\{\-1\}\\mathbf\{K\}\_\{i\}\\mathbf\{v\}\_\{i,t\}\.\(87\)This term is the known\-velocity position response on the right\-hand side of Equation \([84](https://arxiv.org/html/2609.30344#S11.E84)\)\. Because it is evaluated independently at each physical node, its global sum is not constrained to vanish\.

Summing Equation \([85](https://arxiv.org/html/2609.30344#S11.E85)\), using Equation \([82](https://arxiv.org/html/2609.30344#S11.E82)\), and subtracting the external impulse gives the complete one\-substep residual:

𝐑P:=Δ​𝐏−δ​t​𝐅ext=δ​t​∑i\(𝐀i−1−𝐈\)​𝐛i−12​δ​t2​∑i𝐀i−1​𝐊i​𝐯i,t\.\\boxed\{\\mathbf\{R\}\_\{P\}:=\\Delta\\mathbf\{P\}\-\\delta t\\mathbf\{F\}\_\{\\mathrm\{ext\}\}=\\delta t\\sum\_\{i\}\\left\(\\mathbf\{A\}\_\{i\}^\{\-1\}\-\\mathbf\{I\}\\right\)\\mathbf\{b\}\_\{i\}\-\\tfrac\{1\}\{2\}\\delta t^\{2\}\\sum\_\{i\}\\mathbf\{A\}\_\{i\}^\{\-1\}\\mathbf\{K\}\_\{i\}\\mathbf\{v\}\_\{i,t\}\}\.\(88\)The first term appears because each component of the balanced drive is acted on by a different nodal matrix𝐀i−1\\mathbf\{A\}\_\{i\}^\{\-1\}\. The second is the known\-velocity position\-response contribution\. Both terms are specific to the independent nodal solves\.

To show why the same residual is absent from a fully assembled Newmark solve, we reproduce the corresponding global update:

\[𝐌\+12​δ​t​𝐃glob\+14​δ​t2​𝐊glob\]​Δ​𝐯=𝐛​δ​t−12​δ​t2​𝐊glob​𝐯t\.\\left\[\\mathbf\{M\}\+\\tfrac\{1\}\{2\}\\delta t\\,\\mathbf\{D\}\_\{\\mathrm\{glob\}\}\+\\tfrac\{1\}\{4\}\\delta t^\{2\}\\mathbf\{K\}\_\{\\mathrm\{glob\}\}\\right\]\\Delta\\mathbf\{v\}=\\mathbf\{b\}\\delta t\-\\tfrac\{1\}\{2\}\\delta t^\{2\}\\mathbf\{K\}\_\{\\mathrm\{glob\}\}\\mathbf\{v\}\_\{t\}\.\(89\)Let

𝐞⊤:=\[𝐈⋯𝐈\]\\mathbf\{e\}^\{\\top\}:=\\begin\{bmatrix\}\\mathbf\{I\}&\\cdots&\\mathbf\{I\}\\end\{bmatrix\}\(90\)denote the operator that sums the nodal blocks of a stacked vector\. The assembled stiffness and damping operators contain both nodal and cross\-node response blocks\. Their internal block sums vanish:

𝐞⊤​𝐊glob=𝟎,𝐞⊤​𝐃glob=𝟎\.\\mathbf\{e\}^\{\\top\}\\mathbf\{K\}\_\{\\mathrm\{glob\}\}=\\mathbf\{0\},\\qquad\\mathbf\{e\}^\{\\top\}\\mathbf\{D\}\_\{\\mathrm\{glob\}\}=\\mathbf\{0\}\.\(91\)Left\-multiplying Equation \([89](https://arxiv.org/html/2609.30344#S11.E89)\) by𝐞⊤\\mathbf\{e\}^\{\\top\}therefore eliminates the stiffness and damping terms on both sides:

𝐞⊤​𝐌​Δ​𝐯=δ​t​𝐞⊤​𝐛\.\\mathbf\{e\}^\{\\top\}\\mathbf\{M\}\\Delta\\mathbf\{v\}=\\delta t\\,\\mathbf\{e\}^\{\\top\}\\mathbf\{b\}\.\(92\)Using

𝐞⊤​𝐌​Δ​𝐯=Δ​𝐏,𝐞⊤​𝐛=𝐅ext,\\mathbf\{e\}^\{\\top\}\\mathbf\{M\}\\Delta\\mathbf\{v\}=\\Delta\\mathbf\{P\},\\qquad\\mathbf\{e\}^\{\\top\}\\mathbf\{b\}=\\mathbf\{F\}\_\{\\mathrm\{ext\}\},\(93\)gives

Δ​𝐏=δ​t​𝐅ext,𝐑P=𝟎\.\\Delta\\mathbf\{P\}=\\delta t\\mathbf\{F\}\_\{\\mathrm\{ext\}\},\\qquad\\mathbf\{R\}\_\{P\}=\\mathbf\{0\}\.\(94\)
The distinction is structural\. In the assembled system, the internal stiffness and damping contributions cancel when the coupled equations are summed, before the global system is inverted\. InNewmark\-β\\beta\-DGN, the global system is replaced by independent nodal solves, so the balanced drives are acted on separately by the matrices𝐀i−1\\mathbf\{A\}\_\{i\}^\{\-1\}and the cross\-node cancellation is not imposed\. The virtual hub and message passing provide a tractable approximation of the resulting non\-local coupling, while Equation \([88](https://arxiv.org/html/2609.30344#S11.E88)\) quantifies the remaining one\-substep linear\-momentum residual\.

### 11\.3Asymptotic Order of the Residual and Conservative Explicit\-Response Limit

For bounded masses and response operators, Equation \([83](https://arxiv.org/html/2609.30344#S11.E83)\) gives

𝐀i−1=𝐈\+𝒪⁡\(δ​t\)\.\\mathbf\{A\}\_\{i\}^\{\-1\}=\\mathbf\{I\}\+\\mathcal\{O\}\(\\delta t\)\.\(95\)The first term in Equation \([88](https://arxiv.org/html/2609.30344#S11.E88)\) is therefore𝒪⁡\(δ​t2\)\\mathcal\{O\}\(\\delta t^\{2\}\)\. The second already contains the factorδ​t2\\delta t^\{2\}and has the same order for bounded states and response operators\. Consequently,

𝐑P=𝒪⁡\(δ​t2\)\.\\mathbf\{R\}\_\{P\}=\\mathcal\{O\}\(\\delta t^\{2\}\)\.\(96\)The absolute one\-substep residual is second order inδ​t\\delta t\. Relative to a non\-zero impulse of orderδ​t\\delta t, it is first order inδ​t\\delta t\. This is a local asymptotic result, not a bound on accumulated conservation drift\.

The explicit\-response limit is reproduced here for completeness\. As𝐊i,𝐃i→𝟎\\mathbf\{K\}\_\{i\},\\mathbf\{D\}\_\{i\}\\to\\mathbf\{0\},

𝐀i⟶𝐈,Δ​𝐩i⟶𝐛i​δ​t\.\\mathbf\{A\}\_\{i\}\\longrightarrow\\mathbf\{I\},\\qquad\\Delta\\mathbf\{p\}\_\{i\}\\longrightarrow\\mathbf\{b\}\_\{i\}\\delta t\.\(97\)Using Equation \([82](https://arxiv.org/html/2609.30344#S11.E82)\) then gives

Δ​𝐏⟶δ​t​𝐅ext,𝐑P⟶𝟎\.\\Delta\\mathbf\{P\}\\longrightarrow\\delta t\\mathbf\{F\}\_\{\\mathrm\{ext\}\},\\qquad\\mathbf\{R\}\_\{P\}\\longrightarrow\\mathbf\{0\}\.\(98\)Thus balanced decoded forces and heterogeneous masses do not themselves produce the residual\. It enters through the learned semi\-implicit response terms\.

### 11\.4Numerical Verification

We verify the two contributions in Equation \([88](https://arxiv.org/html/2609.30344#S11.E88)\) using a1212\-node system with random positive\-definite position\- and velocity\-response operators, non\-uniform masses, exactly antisymmetric physical\-edge forces, and no external force\. The relative residual is

R⁡\(δ​t\):=‖∑imi​Δ​𝐯i‖∑i‖mi​Δ​𝐯i‖\.R\(\\delta t\):=\\frac\{\\left\\lVert\\sum\_\{i\}m\_\{i\}\\Delta\\mathbf\{v\}\_\{i\}\\right\\rVert\}\{\\sum\_\{i\}\\left\\lVert m\_\{i\}\\Delta\\mathbf\{v\}\_\{i\}\\right\\rVert\}\.\(99\)
Two controlled modifications isolate the two terms in Equation \([88](https://arxiv.org/html/2609.30344#S11.E88)\)\. Replacing the node\-dependent𝐀i−1\\mathbf\{A\}\_\{i\}^\{\-1\}matrices by one common matrix removes the first term because∑i𝐛i=𝟎\\sum\_\{i\}\\mathbf\{b\}\_\{i\}=\\mathbf\{0\}in this test\. Suppressing−12​δ​t2​𝐊i​𝐯i,t\-\\tfrac\{1\}\{2\}\\delta t^\{2\}\\mathbf\{K\}\_\{i\}\\mathbf\{v\}\_\{i,t\}removes the second term\. These are diagnostic checks, not alternative implementations ofNewmark\-β\\beta\-DGN\.

The observed order between the two reported substep sizes is

p:=log⁡\[R⁡\(0\.2\)/R⁡\(0\.0125\)\]log⁡\(16\)\.p:=\\frac\{\\log\\left\[R\(0\.2\)/R\(0\.0125\)\\right\]\}\{\\log\(16\)\}\.\(100\)A value near one corresponds to a first\-order relative residual and, equivalently, a second\-order absolute residual\.

Supplementary Table 20:Linear\-momentum residual of the independent nodal solve\.A common coefficient matrix removes the first term in Equation \([88](https://arxiv.org/html/2609.30344#S11.E88)\)\. Suppressing the known\-velocity position\-response term removes the second\.The original nodal update in row A exhibits the predicted approximately first\-order relative residual\. Retaining either contribution alone, rows B and D, gives the same asymptotic order\. Removing both contributions, row C, reduces the residual to single\-precision round\-off\. The pairwise expression in Equation \([86](https://arxiv.org/html/2609.30344#S11.E86)\) agrees component\-wise with the measured contribution to2\.4×10−72\.4\\times 10^\{\-7\}\.

These checks verify the derived algebra for prescribed response operators\. They do not measure accumulated conservation drift in a trained rollout\.

### 11\.5Scope of the Conservation Statements

The physical\-edge decoder conserves linear and angular momentum exactly under the antisymmetry and shared\-application\-point conditions in Equation \([73](https://arxiv.org/html/2609.30344#S11.E73)\)\. The virtual hub enforces exact collective balance of its decoded force and angular\-momentum\-flux channels\. Its total moment about a common spatial origin additionally depends on the learned application points through Equation \([80](https://arxiv.org/html/2609.30344#S11.E80)\); this geometric contribution has not been measured\.

The decoded interactions therefore retain the intended conservation structure before integration\. The independent semi\-implicit nodal update introduces the linear\-momentum residual in Equation \([88](https://arxiv.org/html/2609.30344#S11.E88)\)\. Under the stated boundedness assumptions, the residual is𝒪⁡\(δ​t2\)\\mathcal\{O\}\(\\delta t^\{2\}\)over one substep and vanishes in the explicit\-response limit\.

The numerical derivation concerns total linear momentum\. It does not establish a corresponding finite\-substep result for the complete angular update\. Conservation drift has also not been measured in the trained rollouts\. The walking systems exchange momentum with the ground, the beam with its clamp, and the protein with its represented environment\. Quantifying physical drift in these systems would require a complete accounting of external impulses and boundary reactions\.