CoSynFlow: Conformal Symplectic Neural Flows for Cross-System Prediction of Dissipative Hamiltonian Dynamics

# Conformal Symplectic Neural Flows for Cross-System Prediction of Dissipative Hamiltonian Dynamics
Source: [https://arxiv.org/html/2608.00571](https://arxiv.org/html/2608.00571)
Baige Xu RIKEN AIP 744 Motooka Nishi\-ku, Fukuoka, Japan baige\.xu@riken\.jp &Takaharu Yaguchi IMI, Kyushu University RIKEN AIP 744 Motooka Nishi\-ku, Fukuoka, Japan yaguchi@imi\.kyushu\-u\.ac\.jp

###### Abstract

Learning solution operators for differential equations is a central problem in scientific machine learning\. However, many neural operator methods optimize prediction accuracy without explicitly enforcing the geometric structure of the dynamics\. Structure\-preserving models such as SympNets and Symplectic Neural Flows address this issue for conservative Hamiltonian systems by preserving the symplectic form\. In dissipative Hamiltonian systems with conformal symplectic structure, however, the symplectic form evolves according to a conformal factor determined by the dissipation\. We proposeCoSynFlow, a conformal symplectic neural flow for learning continuous\-time solution maps of dissipative Hamiltonian dynamics\. CoSynFlow composes symplectic shear maps with explicit conformal scaling, preserving the conformal symplectic structure by construction\. By conditioning it on a finite\-dimensional Hamiltonian descriptor and the dissipation parameter, a single trained model predicts solution maps for unseen systems without retraining\. CoSynFlow keeps the structure error at machine precision, attains the lowest long\-horizon error, and admits physics\-informed training\.

## 1Introduction

Learning solution maps of differential equations is a central goal of scientific machine learning\. Neural operators are a natural framework for this task, since one trained model represents solutions across varying initial conditions, coefficients, or governing equations\[[1](https://arxiv.org/html/2608.00571#bib.bib1),[2](https://arxiv.org/html/2608.00571#bib.bib2),[3](https://arxiv.org/html/2608.00571#bib.bib3)\]\. However, they treat the solution map as an unconstrained function approximation, with no guarantee that the predicted dynamics respect the geometric structure of the underlying equations\.

Geometric structure is fundamental to Hamiltonian dynamics\. The flow of a conservative Hamiltonian system preserves the symplectic form, and therefore also the phase\-space volume\. Hamiltonian Neural Networks, SympNets, and Symplectic Neural Flows encode this geometry at the level of learned vector fields, maps, or flow architectures\[[4](https://arxiv.org/html/2608.00571#bib.bib4),[5](https://arxiv.org/html/2608.00571#bib.bib5),[6](https://arxiv.org/html/2608.00571#bib.bib6)\]\. However, existing map\- and flow\-based architectures are designed around ordinary symplecticity and are commonly trained for a fixed Hamiltonian system\.

The setting considered here raises three challenges\. First, ordinary symplectic preservation is not the appropriate constraint for dissipative Hamiltonian dynamics\. A broad class of such systems can be formulated as conformal Hamiltonian systems\[[7](https://arxiv.org/html/2608.00571#bib.bib7)\], whose vector fieldXXand flowϕt\\phi\_\{t\}, the solution map at timett, satisfy

ℒX​ω=γ​ω,ϕt∗​ω=eγ​t​ω,\\displaystyle\\mathcal\{L\}\_\{X\}\\omega=\\gamma\\omega,\\qquad\\phi\_\{t\}^\{\*\}\\omega=e^\{\\gamma t\}\\omega,whereγ<0\\gamma<0is the dissipation parameter\. The symplectic form is scaled by a prescribed conformal factor rather than preserved exactly\. Second, conventional integrators and learned fixed\-step maps reach a target time through repeated local updates, which motivates a single model that represents the solution maps over a continuous range of query times\. Third, a structure\-preserving model trained for one Hamiltonian does not provide a shared solution model for other Hamiltonians, and applying such models to a system family requires separate training for each system\.

We proposeCoSynFlow, a system\-conditioned conformal symplectic neural flow for learning continuous\-time solution maps of dissipative Hamiltonian dynamics\. Each time\-ttmap is a composition of symmetric blocks built from gradient shear maps and explicit conformal scaling\. The shears are symplectic because they are generated by scalar potentials, and the layer\-wise scaling factors are normalized to multiply toeγ​te^\{\\gamma t\}, so the map is exactly conformal symplectic for arbitrary trainable parameters\. It takes the query time as a continuous input rather than a fixed step\. An operator\-learning conditioner maps a finite\-dimensional Hamiltonian descriptor and the dissipation parameter to the block parameters, so one trained model predicts the solution map of an unseen system without retraining\. We train CoSynFlow on randomly generated Hamiltonian systems and evaluate it on held\-out systems from the training distribution and on four analytic benchmark systems outside it, across several dissipation rates, and beyond the training interval by composing the learned maps\. We compare it with an ordinary symplectic neural flow, unstructured and softly constrained direct\-flow models, and representative neural operator architectures, reporting predictive accuracy, energy behavior, and the conformal symplectic structure error\.

The main contributions of this work are as follows:

- •Exact conformal symplecticity\.For any trainable parameters, the resulting time\-ttmap is exactly conformal symplectic by construction, matching the geometry of dissipative Hamiltonian dynamics\.
- •Continuous\-time flow learning\.The model produces the time\-ttmap for anyttin a single evaluation and composes to predict beyond the training horizon with the exact conformal factor\. Differentiability inttalso admits physics\-informed training on the residual of the governing equation, which is undefined for a fixed\-step map\.
- •Cross\-system prediction\.One shared model predicts the solution map of an unseen dissipative Hamiltonian system without retraining, and conditioning leaves the exact structure intact\.

## 2Related Work

#### Neural Operators for Solution\-Map Learning\.

Deep Operator Networks \(DeepONets\), Fourier Neural Operators \(FNOs\), Physics\-Informed Neural Operators \(PINOs\) and Transolver learn maps between function spaces and predict across varying input functions and system parameters\[[1](https://arxiv.org/html/2608.00571#bib.bib1),[2](https://arxiv.org/html/2608.00571#bib.bib2),[8](https://arxiv.org/html/2608.00571#bib.bib8),[9](https://arxiv.org/html/2608.00571#bib.bib9)\]\. None of them imposes Hamiltonian geometry on the learned operator by construction\. CoSynFlow instead uses operator learning to condition a finite\-time flow whose geometry is fixed by the architecture\.

#### Structure\-Preserving Learning for Hamiltonian Dynamics\.

Hamiltonian Neural Networks \(HNNs\) and Lagrangian Neural Networks \(LNNs\) learn scalar Hamiltonian or Lagrangian functions and derive the dynamics from them, while Symplectic ODE\-Net \(SymODEN\) incorporates control into a Hamiltonian vector\-field model\[[4](https://arxiv.org/html/2608.00571#bib.bib4),[10](https://arxiv.org/html/2608.00571#bib.bib10),[11](https://arxiv.org/html/2608.00571#bib.bib11)\]\. SympNets parameterize exactly symplectic evolution maps, and Symplectic Neural Flows \(SympFlow\) construct time\-dependent symplectic flow maps\[[5](https://arxiv.org/html/2608.00571#bib.bib5),[6](https://arxiv.org/html/2608.00571#bib.bib6)\]\. Poisson Neural Networks and Neural Symplectic Form extend structure\-preserving learning to Poisson systems and noncanonical coordinates\[[12](https://arxiv.org/html/2608.00571#bib.bib12),[13](https://arxiv.org/html/2608.00571#bib.bib13)\]\. These methods address conservative dynamics\. CoSynFlow instead targets dissipative systems, whose flows scale rather than preserve the symplectic form\.

#### Learning Dissipative Hamiltonian Dynamics\.

Dissipative Hamiltonian Neural Networks \(D\-HNNs\) separate Hamiltonian and Rayleigh\-dissipative components, while Dissipative SymODEN and port\-Hamiltonian neural networks encode dissipation and external inputs in structured continuous\-time dynamics\[[14](https://arxiv.org/html/2608.00571#bib.bib14),[15](https://arxiv.org/html/2608.00571#bib.bib15),[16](https://arxiv.org/html/2608.00571#bib.bib16)\]\. Discrete\-gradient models and energy\-behavior\-preserving integrators enforce conservation or dissipation laws in discrete time, and GENERIC\-informed networks encode thermodynamic consistency\[[17](https://arxiv.org/html/2608.00571#bib.bib17),[18](https://arxiv.org/html/2608.00571#bib.bib18),[19](https://arxiv.org/html/2608.00571#bib.bib19)\]\. These approaches identify unknown dynamics from trajectory data, whereas CoSynFlow learns a shared family of solution maps for systems that are specified, and predicts unseen systems without retraining\.

## 3Preliminaries and Problem Setting

### 3\.1Conformal Symplectic Structure

Let\(ℳ,ω\)\(\\mathcal\{M\},\\omega\)be an exact symplectic manifold, so thatω=−d​θ\\omega=\-\\mathrm\{d\}\\thetafor a one\-formθ\\theta\. For a sufficiently smooth HamiltonianH:ℳ→ℝH:\\mathcal\{M\}\\rightarrow\\mathbb\{R\}, the Hamiltonian vector fieldXHX\_\{H\}is defined byιXH​ω=d​H\\iota\_\{X\_\{H\}\}\\omega=\\mathrm\{d\}H, and Cartan’s formula together withd​ω=0\\mathrm\{d\}\\omega=0givesℒXH​ω=0\\mathcal\{L\}\_\{X\_\{H\}\}\\omega=0\. The Liouville vector fieldZZis defined byιZ​ω=−θ\\iota\_\{Z\}\\omega=\-\\thetaand satisfiesℒZ​ω=ω\\mathcal\{L\}\_\{Z\}\\omega=\\omegafor the same reason\. We consider vector fields of the form

X\\displaystyle X=XH\+γ​Z,γ<0,\\displaystyle=X\_\{H\}\+\\gamma Z,\\qquad\\gamma<0,\(1\)which by linearity of the Lie derivative satisfyℒX​ω=γ​ω\\mathcal\{L\}\_\{X\}\\omega=\\gamma\\omegaand are called conformal symplectic\. Letϕt\\phi\_\{t\}denote the flow ofXX, so thatz​\(t\)=ϕt​\(z0\)z\(t\)=\\phi\_\{t\}\(z\_\{0\}\)withz​\(0\)=z0z\(0\)=z\_\{0\}\. Sincedd​t​ϕt∗​ω=ϕt∗​ℒX​ω=γ​ϕt∗​ω\\frac\{\\mathrm\{d\}\}\{\\mathrm\{d\}t\}\\phi\_\{t\}^\{\*\}\\omega=\\phi\_\{t\}^\{\*\}\\mathcal\{L\}\_\{X\}\\omega=\\gamma\\,\\phi\_\{t\}^\{\*\}\\omegaandϕ0∗​ω=ω\\phi\_\{0\}^\{\*\}\\omega=\\omega,

ϕt∗​ω\\displaystyle\\phi\_\{t\}^\{\*\}\\omega=eγ​t​ω\.\\displaystyle=e^\{\\gamma t\}\\omega\.\(2\)Forγ<0\\gamma<0the symplectic form contracts by the prescribed factoreγ​te^\{\\gamma t\}\. This is the geometric relation imposed on the learned maps\.

Throughout this work we takeℳ=ℝ2​d\\mathcal\{M\}=\\mathbb\{R\}^\{2d\}with its canonical structure, writingz=\(q,p\)z=\(q,p\)withq,p∈ℝdq,p\\in\\mathbb\{R\}^\{d\}the generalized positions and momenta, andJ=\(0Id−Id0\)J=\\left\(\\begin\{smallmatrix\}0&I\_\{d\}\\\\ \-I\_\{d\}&0\\end\{smallmatrix\}\\right\)\. Thenθ=∑i=1dpi​d​qi\\theta=\\sum\_\{i=1\}^\{d\}p\_\{i\}\\,\\mathrm\{d\}q\_\{i\},ω=∑i=1dd​qi∧d​pi\\omega=\\sum\_\{i=1\}^\{d\}\\mathrm\{d\}q\_\{i\}\\wedge\\mathrm\{d\}p\_\{i\},XH=J​∇HX\_\{H\}=J\\nabla HandZ​\(q,p\)=\(0,p\)Z\(q,p\)=\(0,p\), and \([2](https://arxiv.org/html/2608.00571#S3.E2)\) readsD​ϕt​\(z\)⊤​J​D​ϕt​\(z\)=eγ​t​JD\\phi\_\{t\}\(z\)^\{\\top\}J\\,D\\phi\_\{t\}\(z\)=e^\{\\gamma t\}J\.

### 3\.2Dissipative Hamiltonian Systems

In canonical coordinates, \([1](https://arxiv.org/html/2608.00571#S3.E1)\) takes the form

z˙\\displaystyle\\dot\{z\}=J​∇H​\(z\)\+γ​\(000Id\)​z,\\displaystyle=J\\nabla H\(z\)\+\\gamma\\begin\{pmatrix\}0&0\\\\ 0&I\_\{d\}\\end\{pmatrix\}z,\(3\)or equivalentlyq˙=∇pH​\(q,p\)\\dot\{q\}=\\nabla\_\{p\}H\(q,p\)andp˙=−∇qH​\(q,p\)\+γ​p\\dot\{p\}=\-\\nabla\_\{q\}H\(q,p\)\+\\gamma p\. Along solutionsd​H/d​t=γ​p⊤​∇pH​\(q,p\)\\mathrm\{d\}H/\\mathrm\{d\}t=\\gamma\\,p^\{\\top\}\\nabla\_\{p\}H\(q,p\), so the energy is no longer conserved and its instantaneous rate of change is proportional toγ\\gamma\. The dissipation acts on the momentum variables alone, andγ=0\\gamma=0recovers the conservative case\. The solution flow of \([3](https://arxiv.org/html/2608.00571#S3.E3)\) isϕt\\phi\_\{t\}and satisfies \([2](https://arxiv.org/html/2608.00571#S3.E2)\)\.

### 3\.3Cross\-System Solution\-Map Learning

Letℋ\\mathcal\{H\}be a family of sufficiently smooth Hamiltonians onℝ2​d\\mathbb\{R\}^\{2d\}and letΓ⊂\(−∞,0\)\\Gamma\\subset\(\-\\infty,0\)be a set of dissipation parameters\. Each pair\(H,γ\)∈ℋ×Γ\(H,\\gamma\)\\in\\mathcal\{H\}\\times\\Gammadefines a system of the form \([3](https://arxiv.org/html/2608.00571#S3.E3)\), whose solution flowϕH,γt\\phi\_\{H,\\gamma\}^\{t\}satisfies\(ϕH,γt\)∗​ω=eγ​t​ω\(\\phi\_\{H,\\gamma\}^\{t\}\)^\{\*\}\\omega=e^\{\\gamma t\}\\omega\.

In the cross\-system setting the Hamiltonian is not fixed, and we assume thatHHandγ\\gammaare known for each system\. To condition a shared model on the Hamiltonian, letr1,…,rm∈ℝ2​dr\_\{1\},\\ldots,r\_\{m\}\\in\\mathbb\{R\}^\{2d\}be fixed sensor points and define the finite\-dimensional descriptorhH=ℰ​\(H\):=\(H​\(r1\),…,H​\(rm\)\)∈ℝmh\_\{H\}=\\mathcal\{E\}\(H\):=\(H\(r\_\{1\}\),\\ldots,H\(r\_\{m\}\)\)\\in\\mathbb\{R\}^\{m\}and the system conditionξ:=\(hH,γ\)\\xi:=\(h\_\{H\},\\gamma\)\. The known Hamiltonian enters the structure\-preserving construction directly, whilehHh\_\{H\}conditions its shared trainable components\.

The objective is to learn one model for the entire family\. Given an initial statez0z\_\{0\}, a query timet∈\[0,T\]t\\in\[0,T\], a known HamiltonianHHand the system conditionξ\\xi, the model approximates

Φ^θt​\(z0;H,ξ\)≈ϕH,γt​\(z0\),\\displaystyle\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\(z\_\{0\};H,\\xi\)\\approx\\phi\_\{H,\\gamma\}^\{t\}\(z\_\{0\}\),so that a single set of trainable parameters represents finite\-time solution maps for different Hamiltonians and dissipation parameters, rather than a separate model for each system\. In addition to predictive accuracy, every conditioned map is required to satisfy\(Φ^θt​\(⋅;H,ξ\)\)∗​ω=eγ​t​ω\(\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\(\\,\\cdot\\,;H,\\xi\)\)^\{\*\}\\omega=e^\{\\gamma t\}\\omega\.

![Refer to caption](https://arxiv.org/html/2608.00571v1/x1.png)Figure 1:Overview of CoSynFlow\. The system conditioner𝒪θ\\mathcal\{O\}\_\{\\theta\}maps a system to the parameters of every flow block\. The query time is read by every block, so one trained model is evaluated at anyttinstead of at a fixed step\. The shears preserveω\\omegaand only the scalings change it, so∑ℓαℓ=1\\sum\_\{\\ell\}\\alpha\_\{\\ell\}=1makes the map conformal symplectic for arbitrary trainable parameters\.

## 4Methodology: CoSynFlow

We first define the flow blocks for a fixed system and then condition their trainable components on the system specification\. Figure[1](https://arxiv.org/html/2608.00571#S3.F1)summarizes the resulting model\.

### 4\.1Conformal Symplectic Flow Blocks

For a fixed HamiltonianHHand dissipation parameterγ\\gamma, we approximate the solution flowϕt\\phi\_\{t\}by a neural flowΦ^θt\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}satisfying \([2](https://arxiv.org/html/2608.00571#S3.E2)\)\. For each layerℓ=1,…,L\\ell=1,\\ldots,L, letA~ℓ\\widetilde\{A\}\_\{\\ell\}andB~ℓ\\widetilde\{B\}\_\{\\ell\}be scalar potentials and define

Pℓt/2​\(q,p\)\\displaystyle P\_\{\\ell\}^\{t/2\}\(q,p\)=\(q,p−t2​∇qA~ℓ​\(q;t\)\),\\displaystyle=\\left\(q,\\;p\-\\tfrac\{t\}\{2\}\\nabla\_\{q\}\\widetilde\{A\}\_\{\\ell\}\(q;t\)\\right\),\(4\)Qℓt​\(q,p\)\\displaystyle Q\_\{\\ell\}^\{t\}\(q,p\)=\(q\+t​∇pB~ℓ​\(p;t\),p\),\\displaystyle=\\left\(q\+t\\nabla\_\{p\}\\widetilde\{B\}\_\{\\ell\}\(p;t\),\\;p\\right\),\(5\)Cℓt/2​\(q,p\)\\displaystyle C\_\{\\ell\}^\{t/2\}\(q,p\)=\(q,eγ​t​αℓ/2​p\),\\displaystyle=\\left\(q,\\;e^\{\\gamma t\\alpha\_\{\\ell\}/2\}p\\right\),\(6\)where∑ℓ=1Lαℓ=1\\sum\_\{\\ell=1\}^\{L\}\\alpha\_\{\\ell\}=1\. The shearsPℓt/2P\_\{\\ell\}^\{t/2\}andQℓtQ\_\{\\ell\}^\{t\}are symplectic because they are generated by scalar potentials, while\(Cℓt/2\)∗​ω=eγ​t​αℓ/2​ω\(C\_\{\\ell\}^\{t/2\}\)^\{\*\}\\omega=e^\{\\gamma t\\alpha\_\{\\ell\}/2\}\\omega\. We use the symmetric block

Fℓt\\displaystyle F\_\{\\ell\}^\{t\}=Cℓt/2∘Pℓt/2∘Qℓt∘Pℓt/2∘Cℓt/2,\\displaystyle=C\_\{\\ell\}^\{t/2\}\\circ P\_\{\\ell\}^\{t/2\}\\circ Q\_\{\\ell\}^\{t\}\\circ P\_\{\\ell\}^\{t/2\}\\circ C\_\{\\ell\}^\{t/2\},\(7\)which satisfies\(Fℓt\)∗​ω=eγ​t​αℓ​ω\(F\_\{\\ell\}^\{t\}\)^\{\*\}\\omega=e^\{\\gamma t\\alpha\_\{\\ell\}\}\\omega, and the complete mapΦ^θt=FLt∘⋯∘F1t\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}=F\_\{L\}^\{t\}\\circ\\cdots\\circ F\_\{1\}^\{t\}therefore satisfies\(Φ^θt\)∗​ω=eγ​t​ω\(\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\)^\{\*\}\\omega=e^\{\\gamma t\}\\omega\. The explicit factors ofttalso giveΦ^θ0=id\\widehat\{\\Phi\}\_\{\\theta\}^\{0\}=\\mathrm\{id\}and allow the model to be evaluated at continuously varying query times\.

### 4\.2System\-Conditioned CoSynFlow

The blocks above are defined for a single system\. To share one model across systems, we learn an operator that maps a system to every system\-dependent component of its blocks,

𝒪θ:\(H,γ\)⟼\(\{αℓ\}ℓ=1L,\{A~ℓ,B~ℓ\}ℓ=1L\)\.\\displaystyle\\mathcal\{O\}\_\{\\theta\}:\\ \(H,\\gamma\)\\longmapsto\\big\(\\\{\\alpha\_\{\\ell\}\\\}\_\{\\ell=1\}^\{L\},\\;\\\{\\widetilde\{A\}\_\{\\ell\},\\widetilde\{B\}\_\{\\ell\}\\\}\_\{\\ell=1\}^\{L\}\\big\)\.Its trainable components read the system only through the descriptorhHh\_\{H\}and the dissipation parameter, so they are conditioned onξ=\(hH,γ\)\\xi=\(h\_\{H\},\\gamma\)\. How this information reaches the potential networks is not essential to the construction\. Supplying it as an additional input and using it to generate the network parameters both yield potentials that vary with the system\. The conditioned flow is

Φ^θt​\(⋅;H,ξ\)=FL,H,ξt∘⋯∘F1,H,ξt,\\displaystyle\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\(\\,\\cdot\\,;H,\\xi\)=F\_\{L,H,\\xi\}^\{t\}\\circ\\cdots\\circ F\_\{1,H,\\xi\}^\{t\},\(8\)whereHHenters through the block potentials andξ\\xithrough the trainable components\. The condition specifies the system and the query time selects which map of the flow is evaluated, so𝒪θ\\mathcal\{O\}\_\{\\theta\}is evaluated once per system and reused across query times\.

For every fixed\(H,ξ,t\)\(H,\\xi,t\), each shear potential still depends only onqqor only onpp, the scaling allocations remain independent of the phase\-space state, and∑ℓαℓ​\(ξ\)=1\\sum\_\{\\ell\}\\alpha\_\{\\ell\}\(\\xi\)=1holds for everyξ\\xi\. The preceding geometric argument therefore applies to every conditioned system\.

### 4\.3Model Instantiation

Any component of the Hamiltonian that depends on the positions alone or on the momenta alone can be evaluated directly in the corresponding shear\. We denote such components byHqH\_\{q\}andHpH\_\{p\}, while the descriptor enters only through𝒪θ\\mathcal\{O\}\_\{\\theta\}\. The scalar block potentials are

A~ℓ​\(q;t,ξ\)\\displaystyle\\widetilde\{A\}\_\{\\ell\}\(q;t,\\xi\)=aℓ​\(ξ\)​Hq​\(q\)\+δ​Aℓ​\(q;t,ξ\),\\displaystyle=a\_\{\\ell\}\(\\xi\)H\_\{q\}\(q\)\+\\delta A\_\{\\ell\}\(q;t,\\xi\),\(9\)B~ℓ​\(p;t,ξ\)\\displaystyle\\widetilde\{B\}\_\{\\ell\}\(p;t,\\xi\)=bℓ​\(ξ\)​Hp​\(p\)\+δ​Bℓ​\(p;t,ξ\),\\displaystyle=b\_\{\\ell\}\(\\xi\)H\_\{p\}\(p\)\+\\delta B\_\{\\ell\}\(p;t,\\xi\),\(10\)whereaℓa\_\{\\ell\}andbℓb\_\{\\ell\}are scalar coefficients andδ​Aℓ\\delta A\_\{\\ell\}andδ​Bℓ\\delta B\_\{\\ell\}are trainable scalar potentials, one pair per block\. The first terms supply the dynamics that is known in closed form and the trainable terms enrich the finite\-depth flow\. Since both enter only as scalar potentials, the shears remain gradient maps, and their gradients are taken by automatic differentiation\.

A Fourier neural operator\[[2](https://arxiv.org/html/2608.00571#bib.bib2)\]encodes\(hH,γ\)\(h\_\{H\},\\gamma\)into a system representation\. A linear head maps it toaℓa\_\{\\ell\},bℓb\_\{\\ell\}and the logits whose softmax givesαℓ\\alpha\_\{\\ell\}, which enforces∑ℓαℓ​\(ξ\)=1\\sum\_\{\\ell\}\\alpha\_\{\\ell\}\(\\xi\)=1\. The same representation conditionsδ​Aℓ\\delta A\_\{\\ell\}andδ​Bℓ\\delta B\_\{\\ell\}through FiLM layers\. Each network takes its own phase\-space variable and features of the query time as input\. Thus, althoughγ\\gammaalso conditions the learned components, the total conformal factor remains exactlyeγ​te^\{\\gamma t\}by construction\.

The trainable potentials are initialized to zero and the coefficients toaℓ=bℓ=αℓ=1/La\_\{\\ell\}=b\_\{\\ell\}=\\alpha\_\{\\ell\}=1/L, so that at initialization CoSynFlow is exactly the classical conformal Strang splitting ofHq\+HpH\_\{q\}\+H\_\{p\}with stept/Lt/L\. Training then adapts the coefficients and the potentials while retaining exact conformal symplecticity\.

## 5Theoretical Guarantees

We give two guarantees for this construction\. Theorem[1](https://arxiv.org/html/2608.00571#Thmtheorem1)states that the learned map is exactly conformal symplectic for every value of the trainable parameters\. Theorem[2](https://arxiv.org/html/2608.00571#Thmtheorem2)states that this exact constraint does not limit the approximation capability of the flow\. Full proofs are given in the appendix\.

###### Theorem 1\(Exact conformal symplecticity\)

Fix a HamiltonianHHand a system conditionξ=\(hH,γ\)\\xi=\(h\_\{H\},\\gamma\)\. Suppose that for everyt∈\[0,T\]t\\in\[0,T\]and everyℓ=1,…,L\\ell=1,\\ldots,Lthe scalar potentialsA~ℓ​\(⋅;t,ξ\)\\widetilde\{A\}\_\{\\ell\}\(\\cdot\\,;t,\\xi\)andB~ℓ​\(⋅;t,ξ\)\\widetilde\{B\}\_\{\\ell\}\(\\cdot\\,;t,\\xi\)areC2C^\{2\}onℝd\\mathbb\{R\}^\{d\}, and that the allocationsαℓ​\(ξ\)\\alpha\_\{\\ell\}\(\\xi\)are independent of the phase\-space variable and satisfy∑ℓ=1Lαℓ​\(ξ\)=1\\sum\_\{\\ell=1\}^\{L\}\\alpha\_\{\\ell\}\(\\xi\)=1\. LetΦ^θt​\(⋅;H,ξ\)\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\(\\cdot\\,;H,\\xi\)be the conditioned composition \([8](https://arxiv.org/html/2608.00571#S4.E8)\) of the symmetric blocks \([7](https://arxiv.org/html/2608.00571#S4.E7)\)\. Then, for arbitrary trainable parameters and everyt∈\[0,T\]t\\in\[0,T\], the mapΦ^θt​\(⋅;H,ξ\)\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\(\\cdot\\,;H,\\xi\)is aC1C^\{1\}diffeomorphism ofℝ2​d\\mathbb\{R\}^\{2d\}with an explicit inverse and satisfies

\(Φ^θt​\(⋅;H,ξ\)\)∗​ω\\displaystyle\\left\(\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\(\\cdot\\,;H,\\xi\)\\right\)^\{\*\}\\omega=eγ​t​ω,\\displaystyle=e^\{\\gamma t\}\\omega,\(11\)equivalentlyD​Φ^θt​\(z\)⊤​J​D​Φ^θt​\(z\)=eγ​t​JD\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\(z\)^\{\\top\}J\\,D\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\(z\)=e^\{\\gamma t\}Jfor allz∈ℝ2​dz\\in\\mathbb\{R\}^\{2d\}\. MoreoverΦ^θ0​\(⋅;H,ξ\)=id\\widehat\{\\Phi\}\_\{\\theta\}^\{0\}\(\\cdot\\,;H,\\xi\)=\\mathrm\{id\}anddetD​Φ^θt​\(z;H,ξ\)=ed​γ​t\\det D\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\(z;H,\\xi\)=e^\{d\\gamma t\}\.

The smoothness assumption holds for the model instantiation, since the Hamiltonian components are smooth and the trainable potentials use smooth activations\. The system representation is a constant when differentiating with respect to the phase\-space variables, so each shear remains the gradient of a scalar potential\.

###### Corollary 1\(Exact structure under composition\)

Under the assumptions of Theorem[1](https://arxiv.org/html/2608.00571#Thmtheorem1), for any query timest1,…,tk∈\[0,T\]t\_\{1\},\\ldots,t\_\{k\}\\in\[0,T\],

\(Φ^θtk∘⋯∘Φ^θt1\)∗​ω=eγ​\(t1\+⋯\+tk\)​ω,\\displaystyle\\left\(\\widehat\{\\Phi\}\_\{\\theta\}^\{t\_\{k\}\}\\circ\\cdots\\circ\\widehat\{\\Phi\}\_\{\\theta\}^\{t\_\{1\}\}\\right\)^\{\*\}\\omega=e^\{\\gamma\(t\_\{1\}\+\\cdots\+t\_\{k\}\)\}\\omega,where all maps are conditioned on the same\(H,ξ\)\(H,\\xi\)\. The relation therefore remains exact under any finite number of compositions, even whent1\+⋯\+tkt\_\{1\}\+\\cdots\+t\_\{k\}exceedsTTand trajectory errors accumulate\.

*Proof sketch\.*Each gradient shear has a block\-triangular Jacobian whose off\-diagonal block is a symmetric Hessian, so every shear is symplectic, and the two half\-scalings in blockℓ\\ellcontribute the factoreγ​t​αℓ​\(ξ\)e^\{\\gamma t\\alpha\_\{\\ell\}\(\\xi\)\}\. Conformal factors multiply under composition, so∑ℓαℓ​\(ξ\)=1\\sum\_\{\\ell\}\\alpha\_\{\\ell\}\(\\xi\)=1gives the total factoreγ​te^\{\\gamma t\}\. Every factor is invertible by reversing its own update, every factor is the identity att=0t=0because of the explicit factors oftt, and the shears have determinant one while the scalings giveed​γ​te^\{d\\gamma t\}\. Details are in the appendix\.

We next show that this exact constraint costs no approximation capability\. The result is stated for a single system, andξ\\xiis suppressed from the notation\.

###### Theorem 2\(Uniform approximation\)

LetK⊂ℝ2​dK\\subset\\mathbb\{R\}^\{2d\}be compact, letT\>0T\>0andγ<0\\gamma<0, letH∈C2​\(ℝ2​d\)H\\in C^\{2\}\(\\mathbb\{R\}^\{2d\}\), and letφt\\varphi^\{t\}be the flow of the conformal Hamiltonian system \([3](https://arxiv.org/html/2608.00571#S3.E3)\)\. Assume that there exists a compact setK′⊂ℝ2​dK^\{\\prime\}\\subset\\mathbb\{R\}^\{2d\}such that

φs​\(K\)⊂int⁡K′\\displaystyle\\varphi^\{s\}\(K\)\\subset\\operatorname\{int\}K^\{\\prime\}for everys∈\[0,T\]s\\in\[0,T\]\. Then, for everyε\>0\\varepsilon\>0, there exist a depthLLand parametersθL\\theta\_\{L\}, independent oftt, such that

supz∈Ksupt∈\[0,T\]‖Φ^θLt​\(z\)−φt​\(z\)‖<ε\.\\displaystyle\\sup\_\{z\\in K\}\\sup\_\{t\\in\[0,T\]\}\\left\\\|\\widehat\{\\Phi\}\_\{\\theta\_\{L\}\}^\{t\}\(z\)\-\\varphi^\{t\}\(z\)\\right\\\|<\\varepsilon\.Moreover, for everyt∈\[0,T\]t\\in\[0,T\],

\(Φ^θLt\)∗​ω=eγ​t​ω\.\\displaystyle\\left\(\\widehat\{\\Phi\}\_\{\\theta\_\{L\}\}^\{t\}\\right\)^\{\*\}\\omega=e^\{\\gamma t\}\\omega\.

*Proof sketch\.*The change of variablesy=\(q,e−γ​t​p\)y=\(q,e^\{\-\\gamma t\}p\)turns the conformal flow into the symplectic flow of the time\-dependent HamiltonianH~​\(t,y\)=e−γ​t​H​\(yq,eγ​t​yp\)\\widetilde\{H\}\(t,y\)=e^\{\-\\gamma t\}H\(y\_\{q\},e^\{\\gamma t\}y\_\{p\}\), which is againC2C^\{2\}inyy\. On a compact setH~\\widetilde\{H\}is approximated by a Hamiltonian that is polynomial in the phase\-space variable, and the resulting flow is split into substeps, each of which is approximated by the flow of a Hamiltonian whose two parts depend on the positions and on the momenta separately\[[20](https://arxiv.org/html/2608.00571#bib.bib20)\]\. Such flows are realized by the gradient shears, whose potentials take the query time as an input, and the scalings are moved between the blocks usingCλ∘PA=Pλ​A∘CλC\_\{\\lambda\}\\circ P\_\{A\}=P\_\{\\lambda A\}\\circ C\_\{\\lambda\}andCλ∘QB=QBλ∘CλC\_\{\\lambda\}\\circ Q\_\{B\}=Q\_\{B\_\{\\lambda\}\}\\circ C\_\{\\lambda\}withBλ​\(p\)=λ​B​\(p/λ\)B\_\{\\lambda\}\(p\)=\\lambda B\(p/\\lambda\)\. A single parameter vector therefore covers the whole interval, and Theorem[1](https://arxiv.org/html/2608.00571#Thmtheorem1)supplies the exact conformal factor\. The argument adapts the universality result for symplectic neural flows\[[6](https://arxiv.org/html/2608.00571#bib.bib6)\]to the conformal setting, the additional ingredients being the change of variables and the exact reinsertion of the scalings into the block structure\. Details are in the appendix\.

## 6Numerical Experiments

### 6\.1Experimental Setup

#### Data\.

We instantiate the construction on two\-degree\-of\-freedom natural HamiltoniansH​\(q,p\)=12​∥p∥2\+V​\(q\)H\(q,p\)=\\tfrac\{1\}\{2\}\\lVert p\\rVert^\{2\}\+V\(q\)on the position domain𝒬=\[−3,3\]2\\mathcal\{Q\}=\[\-3,3\]^\{2\}, so thatHq=VH\_\{q\}=VandHp=12​∥p∥2H\_\{p\}=\\tfrac\{1\}\{2\}\\lVert p\\rVert^\{2\}in \([9](https://arxiv.org/html/2608.00571#S4.E9)\) and \([10](https://arxiv.org/html/2608.00571#S4.E10)\)\. Each potential is a random quartic backbone plus a Gaussian process perturbation of amplitudeδ∈\[0\.2,1\.5\]\\delta\\in\[0\.2,1\.5\], which gives single well, double well and multiple well landscapes, and the dissipation parameter is drawn log\-uniformly with\|γ\|∈\[0\.1,0\.4\]\\lvert\\gamma\\rvert\\in\[0\.1,0\.4\]\. The descriptor consists of the values ofVVon a uniform32×3232\\times 32grid over𝒬\\mathcal\{Q\}\. Reference trajectories are integrated in double precision on\[0,T\]\[0,T\]withT=4T=4, from initial conditions confined to the domain by an energy budget\. The dataset contains38,00038\{,\}000training,1,0001\{,\}000validation and1,0001\{,\}000test systems with3232trajectories each, and training pairs\(z,τ,z\+\)\(z,\\tau,z^\{\+\}\)are sampled online with stratified time gaps\. Sampling ranges and integration settings are given in the appendix\.

#### Evaluation systems\.

Beyond the held\-out random systems we evaluate on four analytic benchmark systems not seen during training\. These are an anisotropic oscillator \(B1\), a coupled Duffing system \(B2\), a Mexican hat potential \(B3\) and a quartic coupled oscillator \(B4\)\. They share the quartic form of the backbone and carry no Gaussian process perturbation, so they lie outside the training distribution, in which every system hasδ≥0\.2\\delta\\geq 0\.2\. Each is evaluated from6464initial conditions drawn from its own energy sublevel set\. Their coefficients are given in the appendix\.

#### Metrics\.

On the held\-out random systems we report the normalized mean squared training lossℒpred\\mathcal\{L\}\_\{\\mathrm\{pred\}\}\. On the benchmark systems the state errorℰz​\(t\)\\mathcal\{E\}\_\{z\}\(t\)is∥z^i​\(t\)−zi​\(t\)∥/σ\\lVert\\widehat\{z\}\_\{i\}\(t\)\-z\_\{i\}\(t\)\\rVert/\\sigmaand the energy errorℰH​\(t\)\\mathcal\{E\}\_\{H\}\(t\)is\|H​\(z^i​\(t\)\)−H​\(zi​\(t\)\)\|/\(Emax−Vmin\)\\lvert H\(\\widehat\{z\}\_\{i\}\(t\)\)\-H\(z\_\{i\}\(t\)\)\\rvert/\(E\_\{\\max\}\-V\_\{\\min\}\), both averaged overN=64N=64initial conditions, withσ=Emax−Vmin\\sigma=\\sqrt\{E\_\{\\max\}\-V\_\{\\min\}\}\. Geometric consistency is measured by the conformal symplectic structure error

ℰω​\(z,t\)\\displaystyle\\mathcal\{E\}\_\{\\omega\}\(z,t\)=∥D​Φ^θt​\(z\)⊤​J​D​Φ^θt​\(z\)−eγ​t​J∥F∥eγ​t​J∥F,\\displaystyle=\\frac\{\\big\\lVert D\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\(z\)^\{\\top\}J\\,D\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\(z\)\-e^\{\\gamma t\}J\\big\\rVert\_\{F\}\}\{\\big\\lVert e^\{\\gamma t\}J\\big\\rVert\_\{F\}\},evaluated in double precision\.

#### Long\-horizon protocol\.

Predictions beyond the trained horizon are produced by composition rather than by a single large query\. We partition\[0,10​T\]\[0,10T\]into ten windows of lengthTT, advance the anchors byzk\+1=Φ^θT​\(zk\)z\_\{k\+1\}=\\widehat\{\\Phi\}\_\{\\theta\}^\{T\}\(z\_\{k\}\), and fill each window with direct queriesΦ^θτ​\(zk\)\\widehat\{\\Phi\}\_\{\\theta\}^\{\\tau\}\(z\_\{k\}\)forτ∈\(0,T\]\\tau\\in\(0,T\]\. Every query time therefore lies inside the trained range and crossing a window boundary is pure composition, so Corollary[1](https://arxiv.org/html/2608.00571#Thmcorollary1)applies at every step\. These results use\|γ\|=0\.1\\lvert\\gamma\\rvert=0\.1, the least damped end of the training range, since at larger dissipation the state has decayed to the fixed point well beforet=10​Tt=10T\.

#### Models and training\.

We compare CoSynFlow, which usesL=12L=12blocks, with five models\.Symplectic Flowuses the same conditioned blocks with the conformal scaling removed, so it represents ordinary symplectic flow models and is also the direct ablation of the conformal component\.MLP Flowis an unstructured direct\-flow model with the same conditioner, andMLP Flow with soft constraintadds a conformal symplectic Jacobian penalty to it\.DeepONetreads the system in its branch network and\(z0,t\)\(z\_\{0\},t\)in its trunk network, andTransolver Flowreplaces the conditioner of MLP Flow with a Physics\-Attention encoder\. We further consider two ablations of CoSynFlow\.Fixed Allocationreplaces the learned scaling allocations byαℓ=1/L\\alpha\_\{\\ell\}=1/L, andNo Correctionremoves the trainable potentials\. All models share the same training data, normalization and optimization budget, and every baseline matches the parameter count of CoSynFlow within3%3\\%, with DeepONet and Transolver using their native encoders\. Results are means over three seeds, each run on a single H100 GPU partition\. Remaining settings are in the appendix\.

### 6\.2Accuracy and Geometric Consistency

Table 1:Accuracy on held\-out random systems and geometric consistency of the learned map att=2t=2\. Lower is better in both columns\.In Table[1](https://arxiv.org/html/2608.00571#S6.T1), CoSynFlow attains the lowest in\-distribution MSE, with a margin over the best baseline far larger than the seed variation, and both ablations outperform every baseline by more than an order of magnitude\. The structure error separates the models far more sharply\. CoSynFlow and the two ablations keep it at double precision round\-off, fourteen orders of magnitude below the unstructured and softly constrained models, as predicted by Theorem[1](https://arxiv.org/html/2608.00571#Thmtheorem1)\. The differences among the three are round\-off noise\. After sixteen compositions, a total time of7\.5​T7\.5T, the three are still of order10−1310^\{\-13\}while every other model exceeds2020, which is the quantitative form of Corollary[1](https://arxiv.org/html/2608.00571#Thmcorollary1)\.

Fixing the allocations toαℓ=1/L\\alpha\_\{\\ell\}=1/Lpreserves the exact structure, since only the normalization enters the guarantee, and costs15%15\\%in accuracy\. Removing the trainable potentials leaves the classical conformal Strang splitting and costs a factor of1\.81\.8\. Removing the conformal scaling instead gives Symplectic Flow, whose structure error is exactlye−γ​t−1e^\{\-\\gamma t\}\-1and is identical across seeds, since it is fixed by the geometry the model enforces rather than by what it learns\. The three ablations therefore separate the roles of the conformal factor, the allocation and the learned potentials\.

### 6\.3Long\-Horizon Prediction

Table 2:Long\-horizon error on the four benchmark systems, averaged overt∈\(0,10​T\]t\\in\(0,10T\]and over6464initial conditions at\|γ\|=0\.1\\lvert\\gamma\\rvert=0\.1\. Predictions beyondTTare obtained by composing the learned flow over ten windows\. The last column is the energy error averaged over the four systems\. DeepONet diverges on B3\. Per\-system energy errors are in the appendix\.![Refer to caption](https://arxiv.org/html/2608.00571v1/x2.png)Figure 2:Trajectories in the position plane\(q1,q2\)\(q\_\{1\},q\_\{2\}\)for the four benchmark systems over0≤t≤10​T0\\leq t\\leq 10Tat\|γ\|=0\.1\\lvert\\gamma\\rvert=0\.1, one initial condition each, with the reference in thick gray\. Predictions are obtained by composing the learned flow over ten windows of lengthTT\.![Refer to caption](https://arxiv.org/html/2608.00571v1/x3.png)Figure 3:First position componentq1q\_\{1\}versus time for the same trajectories as Figure[2](https://arxiv.org/html/2608.00571#S6.F2), with the reference in thick gray\. Faint vertical lines mark the ten composition window boundaries, inside which the model is queried directly at times in\(0,T\]\(0,T\]\.Table[2](https://arxiv.org/html/2608.00571#S6.T2)reports accuracy over ten times the trained horizon\. CoSynFlow attains the lowest state error on three of the four benchmark systems and the lowest average error in both metrics, ahead of the best baseline by nearly a factor of five in state error and twelve in energy error\. Symplectic Flow has the largest state error of all models, three times that of the unconstrained MLP Flow, so ordinary symplecticity is not a weaker constraint here but the wrong one, and enforcing it exactly is worse than enforcing nothing at all\.

Figures[2](https://arxiv.org/html/2608.00571#S6.F2)and[3](https://arxiv.org/html/2608.00571#S6.F3)show why the gap widens with time\. The reference trajectories spiral toward a potential minimum\. CoSynFlow contracts phase\-space volume at exactly the rateed​γ​te^\{d\\gamma t\}and tracks the reference on all four systems over the full horizon\. Symplectic Flow preserves volume and cannot contract at all, so its amplitude never decays, and on the coupled Duffing system it leaves the well in which the true solution settles\. The unstructured models contract at a learned rather than prescribed rate, lose the oscillation phase within a few windows, and on the Mexican\-hat potential settle in a different well from the reference, so their error there is qualitative rather than numerical, whereas the residual deviation of CoSynFlow on that system stays within the correct well\. Theorem[1](https://arxiv.org/html/2608.00571#Thmtheorem1)constrains how the map contracts phase space, not where a trajectory goes, so the amplitude and the phase are both learned\. The predictions also cross the ten window boundaries without a visible break, so composition introduces no artefact of its own\.

When a single system is learned in isolation rather than across the family, the same architecture reaches an MSE as low as3×10−73\\times 10^\{\-7\}, so the errors reported here reflect the cost of amortizing one model over a family rather than a limitation of the construction\.

### 6\.4Extrapolation in Dissipation

![Refer to caption](https://arxiv.org/html/2608.00571v1/x4.png)Figure 4:State errorℰz\\mathcal\{E\}\_\{z\}averaged overt∈\(0,T\]t\\in\(0,T\], on a logarithmic scale\. The shaded region lies outside the training range\|γ\|∈\[0\.1,0\.4\]\\lvert\\gamma\\rvert\\in\[0\.1,0\.4\]and extends to three times its upper end\.Figure[4](https://arxiv.org/html/2608.00571#S6.F4)varies the dissipation parameter over an interval that reaches three times the largest value seen during training\. Every model has a shallow minimum inside the training range and degrades in both directions, but the error of CoSynFlow varies by less than a factor of three over the whole sweep and stays below every baseline at every value of\|γ\|\\lvert\\gamma\\rvert\. The margin grows from about ten times at the edge of the training range to about twenty times at the far end\. The scaling factoreγ​t​αℓe^\{\\gamma t\\alpha\_\{\\ell\}\}is applied analytically, so the conformal factor is correct at anyγ\\gamma\. Accuracy outside the training range remains empirical, since the coefficients and potentials are also conditioned onγ\\gamma, which the baselines must learn in full\. Symplectic Flow degrades fastest and monotonically, because the factor it fails to represent grows with\|γ\|\\lvert\\gamma\\rvert\.

### 6\.5Physics\-Informed Training

BecauseΦ^θt\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}is differentiable in the query time, the residual of the governing equation can be evaluated directly on the learned flow,

ℒPI​\(θ\)=1\|𝒞\|​∑\(z,t\)∈𝒞∥∂tΦ^θt​\(z\)−X​\(Φ^θt​\(z\)\)∥24​\(σ/T\)2,\\displaystyle\\mathcal\{L\}\_\{\\mathrm\{PI\}\}\(\\theta\)=\\frac\{1\}\{\\lvert\\mathcal\{C\}\\rvert\}\\sum\_\{\(z,t\)\\in\\mathcal\{C\}\}\\frac\{\\big\\lVert\\partial\_\{t\}\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\(z\)\-X\\big\(\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\(z\)\\big\)\\big\\rVert^\{2\}\}\{4\\,\(\\sigma/T\)^\{2\}\},whereXXis the vector field \([3](https://arxiv.org/html/2608.00571#S3.E3)\) andσ/T\\sigma/Tis a velocity scale that makes the term dimensionless\. The collocation set𝒞\\mathcal\{C\}consists of states drawn from the energy sublevel set and query times drawn from\[0,T\]\[0,T\], resampled at every step, so the term supplies supervision at no data cost\. A model of a fixed\-step map cannot form this residual, since its output is defined only at multiples of the step and carries no derivative in the query time\.

We test the term on a single fixed system, the coupled Duffing benchmark B2 at\|γ\|=0\.2\\lvert\\gamma\\rvert=0\.2, with only two reference trajectories, the regime in which supervision that needs no data is expected to matter\. The same architecture is trained twice under an identical budget, once withℒpred\\mathcal\{L\}\_\{\\mathrm\{pred\}\}alone and once withℒpred\+λ​ℒPI\\mathcal\{L\}\_\{\\mathrm\{pred\}\}\+\\lambda\\mathcal\{L\}\_\{\\mathrm\{PI\}\}atλ=1\\lambda=1, with three seeds each, and errors are measured on6464held\-out trajectories\.

Table 3:Physics\-informed training of CoSynFlow on a single system with two training trajectories\. The two rows differ only in the training objective\.In Table[3](https://arxiv.org/html/2608.00571#S6.T3), the residual term lowers the state error by a factor of four and the energy error by six, and a similar gain is observed after composition to10​T10T\. Every seed with the term beats every seed without it, and the spread across seeds falls by almost an order of magnitude\. The structure error is unaffected by the additional term, as reported in the appendix\. The term therefore changes what the model learns and never whether the geometry holds\.

## 7Conclusion

We presented CoSynFlow, a system\-conditioned neural flow for the solution maps of dissipative Hamiltonian dynamics\. Each time\-ttmap composes symmetric blocks of gradient shears and conformal scalings\. The shears are symplectic because they are generated by scalar potentials, and the scaling factors are normalized to multiply toeγ​te^\{\\gamma t\}, so the map is exactly conformal symplectic for arbitrary trainable parameters, at every query time and under any number of compositions\. Because the query time is a continuous differentiable input, one trained model is evaluated at any time, composed beyond the training horizon, and trained with a physics\-informed residual\. Conditioning on a Hamiltonian descriptor and the dissipation parameter extends the construction to a family of systems without weakening the guarantee\.

Across four unseen benchmark systems the model attains the lowest state and energy error among structured and unstructured baselines at a matched parameter count, while holding the structure error at double precision round\-off, both at a single query time and under repeated composition\. Accuracy is retained well outside the range of dissipation rates seen during training, where the conformal factor is applied analytically\. The comparison also shows that a structure can be enforced exactly and still be the wrong one, since ordinary symplecticity performs worse here than no constraint at all\.

Being differentiable and exactly structured, the learned map can stand in for a solver in loops where many solutions of one system family are required\. The physics\-informed residual needs no reference trajectory, which points to building such surrogates from far less data\. More broadly, factoring a flow into structure preserving parts and one analytic part that carries the prescribed geometric weight is not specific to conformal symplecticity, and we expect it to transfer to other structured dynamics\.

## Appendix AProofs

### A\.1Proof of Theorem 1 and Corollary 1

###### Proof 1\(Proof of Theorem 1\)

Writeωz​\(u,v\)=u⊤​J​v\\omega\_\{z\}\(u,v\)=u^\{\\top\}Jv, so that for aC1C^\{1\}mapΨ\\Psiwe have

\(Ψ∗​ω\)z​\(u,v\)=u⊤​D​Ψ​\(z\)⊤​J​D​Ψ​\(z\)​v,\\displaystyle\(\\Psi^\{\*\}\\omega\)\_\{z\}\(u,v\)=u^\{\\top\}D\\Psi\(z\)^\{\\top\}J\\,D\\Psi\(z\)\\,v,andΨ∗​ω=c​ω\\Psi^\{\*\}\\omega=c\\,\\omegaholds precisely whenD​Ψ​\(z\)⊤​J​D​Ψ​\(z\)=c​JD\\Psi\(z\)^\{\\top\}J\\,D\\Psi\(z\)=cJfor everyzz\. We call such a map conformal symplectic with factorcc\. IfΨ1\\Psi\_\{1\}andΨ2\\Psi\_\{2\}are conformal symplectic with factorsc1c\_\{1\}andc2c\_\{2\}, the chain rule gives

D​\(Ψ2∘Ψ1\)⊤​J​D​\(Ψ2∘Ψ1\)=D​Ψ1⊤​\[D​Ψ2⊤​J​D​Ψ2\]​D​Ψ1=c1​c2​J,D\(\\Psi\_\{2\}\\circ\\Psi\_\{1\}\)^\{\\top\}J\\,D\(\\Psi\_\{2\}\\circ\\Psi\_\{1\}\)=D\\Psi\_\{1\}^\{\\top\}\\big\[D\\Psi\_\{2\}^\{\\top\}J\\,D\\Psi\_\{2\}\\big\]D\\Psi\_\{1\}=c\_\{1\}c\_\{2\}\\,J,\(12\)where the bracket is evaluated atΨ1​\(z\)\\Psi\_\{1\}\(z\)\. Determinants also multiply, sincedet\(M​N\)=detM​detN\\det\(MN\)=\\det M\\det N\.

SetAℓ:=t2​A~ℓ​\(⋅;t,ξ\)A\_\{\\ell\}:=\\tfrac\{t\}\{2\}\\widetilde\{A\}\_\{\\ell\}\(\\,\\cdot\\,;t,\\xi\),Bℓ:=t​B~ℓ​\(⋅;t,ξ\)B\_\{\\ell\}:=t\\,\\widetilde\{B\}\_\{\\ell\}\(\\,\\cdot\\,;t,\\xi\)andλℓ:=eγ​t​αℓ​\(ξ\)/2\\lambda\_\{\\ell\}:=e^\{\\gamma t\\alpha\_\{\\ell\}\(\\xi\)/2\}\. By hypothesisAℓA\_\{\\ell\}andBℓB\_\{\\ell\}areC2C^\{2\}onℝd\\mathbb\{R\}^\{d\}, andλℓ\>0\\lambda\_\{\\ell\}\>0\. Becauseαℓ​\(ξ\)\\alpha\_\{\\ell\}\(\\xi\)does not depend on the phase\-space variable, eachλℓ\\lambda\_\{\\ell\}is a constant, and because the system representation is held fixed while differentiating in\(q,p\)\(q,p\), the factors ofFℓtF\_\{\\ell\}^\{t\}are the three maps

PA​\(q,p\)=\(q,p−∇A​\(q\)\),QB​\(q,p\)=\(q\+∇B​\(p\),p\),Cλ​\(q,p\)=\(q,λ​p\)P\_\{A\}\(q,p\)=\\big\(q,\\;p\-\\nabla A\(q\)\\big\),\\qquad Q\_\{B\}\(q,p\)=\\big\(q\+\\nabla B\(p\),\\;p\\big\),\\qquad C\_\{\\lambda\}\(q,p\)=\(q,\\;\\lambda p\)withA=AℓA=A\_\{\\ell\},B=BℓB=B\_\{\\ell\}andλ=λℓ\\lambda=\\lambda\_\{\\ell\}\.

For the shears, writeS:=∇2A​\(q\)S:=\\nabla^\{2\}A\(q\)andR:=∇2B​\(p\)R:=\\nabla^\{2\}B\(p\)\. Then

D​PA⊤​J​D​PA=\(I−S⊤0I\)​\(−SI−I0\)=\(S⊤−SI−I0\)=J,D​QB⊤​J​D​QB=\(0I−IR⊤−R\)=J,DP\_\{A\}^\{\\top\}J\\,DP\_\{A\}=\\begin\{pmatrix\}I&\-S^\{\\top\}\\\\ 0&I\\end\{pmatrix\}\\begin\{pmatrix\}\-S&I\\\\ \-I&0\\end\{pmatrix\}=\\begin\{pmatrix\}S^\{\\top\}\-S&I\\\\ \-I&0\\end\{pmatrix\}=J,\\qquad DQ\_\{B\}^\{\\top\}J\\,DQ\_\{B\}=\\begin\{pmatrix\}0&I\\\\ \-I&R^\{\\top\}\-R\\end\{pmatrix\}=J,since the Hessian of aC2C^\{2\}scalar function is symmetric\. Both Jacobians are block unitriangular, sodetD​PA=detD​QB=1\\det DP\_\{A\}=\\det DQ\_\{B\}=1\. For the scaling,D​Cλ=diag⁡\(I,λ​I\)DC\_\{\\lambda\}=\\operatorname\{diag\}\(I,\\lambda I\)gives

D​Cλ⊤​J​D​Cλ=\(0λ​I−λ​I0\)=λ​J,detD​Cλ=λd\.DC\_\{\\lambda\}^\{\\top\}J\\,DC\_\{\\lambda\}=\\begin\{pmatrix\}0&\\lambda I\\\\ \-\\lambda I&0\\end\{pmatrix\}=\\lambda J,\\qquad\\det DC\_\{\\lambda\}=\\lambda^\{d\}\.
Each blockFℓtF\_\{\\ell\}^\{t\}consists of three shears of factor one and two scalings of factorλℓ\\lambda\_\{\\ell\}, so by \([12](https://arxiv.org/html/2608.00571#A1.E12)\)

\(Fℓt\)∗​ω=λℓ2​ω=eγ​t​αℓ​\(ξ\)​ω,detD​Fℓt=ed​γ​t​αℓ​\(ξ\),\\big\(F\_\{\\ell\}^\{t\}\\big\)^\{\*\}\\omega=\\lambda\_\{\\ell\}^\{2\}\\,\\omega=e^\{\\gamma t\\alpha\_\{\\ell\}\(\\xi\)\}\\,\\omega,\\qquad\\det DF\_\{\\ell\}^\{t\}=e^\{d\\gamma t\\alpha\_\{\\ell\}\(\\xi\)\},and composing theLLblocks with∑ℓ=1Lαℓ​\(ξ\)=1\\sum\_\{\\ell=1\}^\{L\}\\alpha\_\{\\ell\}\(\\xi\)=1yields

\(Φ^θt\)∗​ω=∏ℓ=1Leγ​t​αℓ​\(ξ\)​ω=eγ​t​ω,detD​Φ^θt=∏ℓ=1Led​γ​t​αℓ​\(ξ\)=ed​γ​t\.\\big\(\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\\big\)^\{\*\}\\omega=\\prod\_\{\\ell=1\}^\{L\}e^\{\\gamma t\\alpha\_\{\\ell\}\(\\xi\)\}\\,\\omega=e^\{\\gamma t\}\\omega,\\qquad\\det D\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}=\\prod\_\{\\ell=1\}^\{L\}e^\{d\\gamma t\\alpha\_\{\\ell\}\(\\xi\)\}=e^\{d\\gamma t\}\.
Each factor is a bijection ofℝ2​d\\mathbb\{R\}^\{2d\}whose inverse is a map of the same type,

PA−1​\(q,p\)=\(q,p\+∇A​\(q\)\),QB−1​\(q,p\)=\(q−∇B​\(p\),p\),Cλ−1=C1/λ,P\_\{A\}^\{\-1\}\(q,p\)=\\big\(q,\\;p\+\\nabla A\(q\)\\big\),\\qquad Q\_\{B\}^\{\-1\}\(q,p\)=\\big\(q\-\\nabla B\(p\),\\;p\\big\),\\qquad C\_\{\\lambda\}^\{\-1\}=C\_\{1/\\lambda\},well defined becausePAP\_\{A\}leavesqqunchanged,QBQ\_\{B\}leavesppunchanged andλ\>0\\lambda\>0\. All of these areC1C^\{1\}sinceA,B∈C2A,B\\in C^\{2\}, so every factor is a globalC1C^\{1\}diffeomorphism ofℝ2​d\\mathbb\{R\}^\{2d\}and so isΦ^θt\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}, with

\(Φ^θt\)−1=\(F1t\)−1∘⋯∘\(FLt\)−1,\(Fℓt\)−1=Cλℓ−1∘PAℓ−1∘QBℓ−1∘PAℓ−1∘Cλℓ−1\.\\big\(\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\\big\)^\{\-1\}=\\big\(F\_\{1\}^\{t\}\\big\)^\{\-1\}\\circ\\cdots\\circ\\big\(F\_\{L\}^\{t\}\\big\)^\{\-1\},\\qquad\\big\(F\_\{\\ell\}^\{t\}\\big\)^\{\-1\}=C\_\{\\lambda\_\{\\ell\}\}^\{\-1\}\\circ P\_\{A\_\{\\ell\}\}^\{\-1\}\\circ Q\_\{B\_\{\\ell\}\}^\{\-1\}\\circ P\_\{A\_\{\\ell\}\}^\{\-1\}\\circ C\_\{\\lambda\_\{\\ell\}\}^\{\-1\}\.Finally, the prefactors inAℓA\_\{\\ell\}andBℓB\_\{\\ell\}do not depend on the phase\-space variable, so att=0t=0both gradients vanish andλℓ=1\\lambda\_\{\\ell\}=1\. Every factor is then the identity andΦ^θ0=id\\widehat\{\\Phi\}\_\{\\theta\}^\{0\}=\\mathrm\{id\}\.

###### Proof 2\(Proof of Corollary 1\)

The hypotheses of Theorem 1 hold for everyt∈\[0,T\]t\\in\[0,T\], so\(Φ^θti\)∗​ω=eγ​ti​ω\\big\(\\widehat\{\\Phi\}\_\{\\theta\}^\{t\_\{i\}\}\\big\)^\{\*\}\\omega=e^\{\\gamma t\_\{i\}\}\\omegafor eachii, and applying \([12](https://arxiv.org/html/2608.00571#A1.E12)\) repeatedly gives

\(Φ^θtk∘⋯∘Φ^θt1\)∗​ω=∏i=1keγ​ti​ω=eγ​\(t1\+⋯\+tk\)​ω\.\\Big\(\\widehat\{\\Phi\}\_\{\\theta\}^\{t\_\{k\}\}\\circ\\cdots\\circ\\widehat\{\\Phi\}\_\{\\theta\}^\{t\_\{1\}\}\\Big\)^\{\*\}\\omega=\\prod\_\{i=1\}^\{k\}e^\{\\gamma t\_\{i\}\}\\,\\omega=e^\{\\gamma\(t\_\{1\}\+\\cdots\+t\_\{k\}\)\}\\omega\.The identity holds for anyt1,…,tk∈\[0,T\]t\_\{1\},\\ldots,t\_\{k\}\\in\[0,T\], including sums that exceedTT\.

### A\.2Proof of Theorem 2

###### Proof 3\(Proof of Theorem 2\)

*Change of variables\.*Fors∈ℝs\\in\\mathbb\{R\}letΛs:=Ce−γ​s\\Lambda\_\{s\}:=C\_\{e^\{\-\\gamma s\}\}, that isΛs​\(q,p\)=\(q,e−γ​s​p\)\\Lambda\_\{s\}\(q,p\)=\(q,e^\{\-\\gamma s\}p\), soΛ0=id\\Lambda\_\{0\}=\\mathrm\{id\}\. Letz​\(t\)=φt​\(z0\)z\(t\)=\\varphi^\{t\}\(z\_\{0\}\)and puty​\(t\):=Λt​\(z​\(t\)\)y\(t\):=\\Lambda\_\{t\}\(z\(t\)\), whose components are writteny=\(yq,yp\)y=\(y\_\{q\},y\_\{p\}\), so thatyq=qy\_\{q\}=qandyp=e−γ​t​py\_\{p\}=e^\{\-\\gamma t\}p\. Then

y˙q=∇pH​\(z\),y˙p=−γ​e−γ​t​zp\+e−γ​t​z˙p=−e−γ​t​∇qH​\(z\),\\dot\{y\}\_\{q\}=\\nabla\_\{p\}H\(z\),\\qquad\\dot\{y\}\_\{p\}=\-\\gamma e^\{\-\\gamma t\}z\_\{p\}\+e^\{\-\\gamma t\}\\dot\{z\}\_\{p\}=\-e^\{\-\\gamma t\}\\nabla\_\{q\}H\(z\),and with

H~​\(t,y\):=e−γ​t​H​\(yq,eγ​t​yp\)\\widetilde\{H\}\(t,y\):=e^\{\-\\gamma t\}H\\big\(y\_\{q\},\\;e^\{\\gamma t\}y\_\{p\}\\big\)\(13\)one has∇yqH~​\(t,y\)=e−γ​t​∇qH​\(z\)\\nabla\_\{y\_\{q\}\}\\widetilde\{H\}\(t,y\)=e^\{\-\\gamma t\}\\nabla\_\{q\}H\(z\)and∇ypH~​\(t,y\)=∇pH​\(z\)\\nabla\_\{y\_\{p\}\}\\widetilde\{H\}\(t,y\)=\\nabla\_\{p\}H\(z\)forz=\(yq,eγ​t​yp\)z=\(y\_\{q\},e^\{\\gamma t\}y\_\{p\}\), soy˙=J​∇yH~​\(t,y\)\\dot\{y\}=J\\nabla\_\{y\}\\widetilde\{H\}\(t,y\)\. Writingψt\\psi^\{t\}for the flow of this equation and usingΛ0=id\\Lambda\_\{0\}=\\mathrm\{id\},

φt=Ceγ​t∘ψt\.\\varphi^\{t\}=C\_\{e^\{\\gamma t\}\}\\circ\\psi^\{t\}\.\(14\)By \([13](https://arxiv.org/html/2608.00571#A1.E13)\) andH∈C2​\(ℝ2​d\)H\\in C^\{2\}\(\\mathbb\{R\}^\{2d\}\), the functionH~​\(t,⋅\)\\widetilde\{H\}\(t,\\cdot\)isC2C^\{2\}onℝ2​d\\mathbb\{R\}^\{2d\}for eachttand depends smoothly ontt\.

*Transfer of the containment assumption\.*SetKy′:=⋃s∈\[0,T\]Λs​\(K′\)K^\{\\prime\}\_\{y\}:=\\bigcup\_\{s\\in\[0,T\]\}\\Lambda\_\{s\}\(K^\{\\prime\}\), the image of the compact set\[0,T\]×K′\[0,T\]\\times K^\{\\prime\}under the continuous map\(s,z\)↦Λs​\(z\)\(s,z\)\\mapsto\\Lambda\_\{s\}\(z\), hence compact\. Forz∈Kz\\in Kands∈\[0,T\]s\\in\[0,T\], \([14](https://arxiv.org/html/2608.00571#A1.E14)\) givesψs​\(z\)=Λs​\(φs​\(z\)\)\\psi^\{s\}\(z\)=\\Lambda\_\{s\}\(\\varphi^\{s\}\(z\)\), and sinceΛs\\Lambda\_\{s\}is a linear isomorphism,

ψs​\(K\)⊂Λs​\(int⁡K′\)=int⁡Λs​\(K′\)⊆int⁡Ky′\.\\psi^\{s\}\(K\)\\subset\\Lambda\_\{s\}\\big\(\\operatorname\{int\}K^\{\\prime\}\\big\)=\\operatorname\{int\}\\Lambda\_\{s\}\(K^\{\\prime\}\)\\subseteq\\operatorname\{int\}K^\{\\prime\}\_\{y\}\.Every trajectory ofH~\\widetilde\{H\}issued fromKKtherefore stays in the fixed compact setKy′K^\{\\prime\}\_\{y\}on\[0,T\]\[0,T\]\. The statement of\[[6](https://arxiv.org/html/2608.00571#bib.bib6), Theorem 1\]assumes forward invariance, but its proof uses containment only in this form, through a Gronwall estimate together with a stopping\-time argument on a compact neighbourhood\.

*Approximation ofψt\\psi^\{t\}\.*We apply the argument of\[[6](https://arxiv.org/html/2608.00571#bib.bib6), Appendix F\]toH~\\widetilde\{H\}onKy′K^\{\\prime\}\_\{y\}and recall its five steps\. First, by the Weierstrass approximation theoremH~\\widetilde\{H\}is replaced by a Hamiltonian polynomial in the phase\-space variable with coefficients continuous intt, whose flow is uniformly close toψt\\psi^\{t\}onK×\[0,T\]K\\times\[0,T\]\. Second, that flow is split intoNNsubsteps of lengtht/Nt/N\. Third, each substep is approximated toO​\(t2/N2\)O\(t^\{2\}/N^\{2\}\)by the time\-one flow of a Hamiltonian whose two parts depend on the positions and on the momenta separately\[[20](https://arxiv.org/html/2608.00571#bib.bib20), Lemma 1\], and the two parts are then replaced by single\-hidden\-layer networks, which are dense inC1C^\{1\}on compact sets\[[21](https://arxiv.org/html/2608.00571#bib.bib21)\]\. Fourth, a Lie–Trotter splitting turns each such time\-one flow into a composition of maps

\(q,p\)↦\(q\+∇p\[V~​\(t,p\)−V~​\(0,p\)\],p\),\(q,p\)↦\(q,p−∇q\[V~​\(t,q\)−V~​\(0,q\)\]\),\\big\(q,\\,p\\big\)\\mapsto\\big\(q\+\\nabla\_\{p\}\[\\,\\widetilde\{V\}\(t,p\)\-\\widetilde\{V\}\(0,p\)\\,\],\\;p\\big\),\\qquad\\big\(q,\\,p\\big\)\\mapsto\\big\(q,\\;p\-\\nabla\_\{q\}\[\\,\\widetilde\{V\}\(t,q\)\-\\widetilde\{V\}\(0,q\)\\,\]\\big\),\(15\)each pair contributing an errorO​\(1/R2\)O\(1/R^\{2\}\)overRRsubsteps\. Fifth, the approximations are combined, the composition of such maps being again of the same form\. The conclusion is that for everyε1\>0\\varepsilon\_\{1\}\>0there existM∈ℕM\\in\\mathbb\{N\}and potentialsV~1q,…,V~Mq\\widetilde\{V\}^\{q\}\_\{1\},\\ldots,\\widetilde\{V\}^\{q\}\_\{M\},V~1p,…,V~Mp\\widetilde\{V\}^\{p\}\_\{1\},\\ldots,\\widetilde\{V\}^\{p\}\_\{M\}, eachC2C^\{2\}in the phase\-space variable andC1C^\{1\}intt, such that

St:=QFMt∘PGMt∘⋯∘QF1t∘PG1t,Gmt:=V~mq​\(t,⋅\)−V~mq​\(0,⋅\),Fmt:=V~mp​\(t,⋅\)−V~mp​\(0,⋅\),S^\{t\}:=Q\_\{F^\{t\}\_\{M\}\}\\circ P\_\{G^\{t\}\_\{M\}\}\\circ\\cdots\\circ Q\_\{F^\{t\}\_\{1\}\}\\circ P\_\{G^\{t\}\_\{1\}\},\\qquad\\begin\{aligned\} G^\{t\}\_\{m\}&:=\\widetilde\{V\}^\{q\}\_\{m\}\(t,\\cdot\)\-\\widetilde\{V\}^\{q\}\_\{m\}\(0,\\cdot\),\\\\ F^\{t\}\_\{m\}&:=\\widetilde\{V\}^\{p\}\_\{m\}\(t,\\cdot\)\-\\widetilde\{V\}^\{p\}\_\{m\}\(0,\\cdot\),\\end\{aligned\}\(16\)satisfiessupz∈Ksupt∈\[0,T\]‖St​\(z\)−ψt​\(z\)‖<ε1\\sup\_\{z\\in K\}\\sup\_\{t\\in\[0,T\]\}\\left\\lVert S^\{t\}\(z\)\-\\psi^\{t\}\(z\)\\right\\rVert<\\varepsilon\_\{1\}, withPAP\_\{A\}andQBQ\_\{B\}the gradient shears of the proof of Theorem 1\.

*Normal form of the architecture\.*BesidesPA∘PA′=PA\+A′P\_\{A\}\\circ P\_\{A^\{\\prime\}\}=P\_\{A\+A^\{\\prime\}\}andCλ∘Cμ=Cλ​μC\_\{\\lambda\}\\circ C\_\{\\mu\}=C\_\{\\lambda\\mu\},

PA∘Cλ=Cλ∘PA/λ,QB∘Cλ=Cλ∘QBλ,Bλ​\(p\):=1λ​B​\(λ​p\),P\_\{A\}\\circ C\_\{\\lambda\}=C\_\{\\lambda\}\\circ P\_\{A/\\lambda\},\\qquad Q\_\{B\}\\circ C\_\{\\lambda\}=C\_\{\\lambda\}\\circ Q\_\{B^\{\\lambda\}\},\\qquad B^\{\\lambda\}\(p\):=\\tfrac\{1\}\{\\lambda\}B\(\\lambda p\),\(17\)the second because∇Bλ​\(p\)=∇B​\(λ​p\)\\nabla B^\{\\lambda\}\(p\)=\\nabla B\(\\lambda p\)\. Applying \([17](https://arxiv.org/html/2608.00571#A1.E17)\) inside one block givesFℓt=Cλℓ2∘PAℓ/λℓ∘QBℓλℓ∘PAℓ/λℓF\_\{\\ell\}^\{t\}=C\_\{\\lambda\_\{\\ell\}^\{2\}\}\\circ P\_\{A\_\{\\ell\}/\\lambda\_\{\\ell\}\}\\circ Q\_\{B\_\{\\ell\}^\{\\lambda\_\{\\ell\}\}\}\\circ P\_\{A\_\{\\ell\}/\\lambda\_\{\\ell\}\}, and moving every scaling to the left across the remaining blocks yields

Φ^θt=Ceγ​t∘PG^L\+1∘QF^L∘PG^L∘⋯∘QF^1∘PG^1,\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}=C\_\{e^\{\\gamma t\}\}\\circ P\_\{\\widehat\{G\}\_\{L\+1\}\}\\circ Q\_\{\\widehat\{F\}\_\{L\}\}\\circ P\_\{\\widehat\{G\}\_\{L\}\}\\circ\\cdots\\circ Q\_\{\\widehat\{F\}\_\{1\}\}\\circ P\_\{\\widehat\{G\}\_\{1\}\},\(18\)the total scaling being∏ℓλℓ2=eγ​t\\prod\_\{\\ell\}\\lambda\_\{\\ell\}^\{2\}=e^\{\\gamma t\}by∑ℓαℓ​\(ξ\)=1\\sum\_\{\\ell\}\\alpha\_\{\\ell\}\(\\xi\)=1\. Writingcℓ:=\(λℓ​∏j<ℓλj2\)−1\>0c\_\{\\ell\}:=\\big\(\\lambda\_\{\\ell\}\\prod\_\{j<\\ell\}\\lambda\_\{j\}^\{2\}\\big\)^\{\-1\}\>0andκℓ:=λℓ​∏j<ℓλj2\\kappa\_\{\\ell\}:=\\lambda\_\{\\ell\}\\prod\_\{j<\\ell\}\\lambda\_\{j\}^\{2\}, the merged potentials are

F^ℓ=\(Bℓ\)κℓ,G^1=c1​A1,G^ℓ=cℓ−1​Aℓ−1\+cℓ​Aℓ\(2≤ℓ≤L\),G^L\+1=cL​AL\.\\widehat\{F\}\_\{\\ell\}=\\big\(B\_\{\\ell\}\\big\)^\{\\kappa\_\{\\ell\}\},\\qquad\\widehat\{G\}\_\{1\}=c\_\{1\}A\_\{1\},\\qquad\\widehat\{G\}\_\{\\ell\}=c\_\{\\ell\-1\}A\_\{\\ell\-1\}\+c\_\{\\ell\}A\_\{\\ell\}\\ \\ \(2\\leq\\ell\\leq L\),\\qquad\\widehat\{G\}\_\{L\+1\}=c\_\{L\}A\_\{L\}\.
*Matching\.*TakeL=M\+1L=M\+1\. PrescribingG^ℓ=Gℓt\\widehat\{G\}\_\{\\ell\}=G^\{t\}\_\{\\ell\}forℓ≤M\\ell\\leq MdeterminesA1,…,AMA\_\{1\},\\ldots,A\_\{M\}recursively, since the relations are triangular with diagonal entriescℓ≠0c\_\{\\ell\}\\neq 0\. PrescribingF^ℓ=Fℓt\\widehat\{F\}\_\{\\ell\}=F^\{t\}\_\{\\ell\}forℓ≤M\\ell\\leq MdeterminesBℓB\_\{\\ell\}, because\(Bλ\)1/λ=B\(B^\{\\lambda\}\)^\{1/\\lambda\}=B; setBM\+1=0B\_\{M\+1\}=0, soQF^M\+1=idQ\_\{\\widehat\{F\}\_\{M\+1\}\}=\\mathrm\{id\}and the last two kicks of \([18](https://arxiv.org/html/2608.00571#A1.E18)\) merge intoPG^M\+1\+G^M\+2P\_\{\\widehat\{G\}\_\{M\+1\}\+\\widehat\{G\}\_\{M\+2\}\}\. Choosing

AM\+1=−cM​AMcM\+1\+cM\+2A\_\{M\+1\}=\-\\frac\{c\_\{M\}A\_\{M\}\}\{c\_\{M\+1\}\+c\_\{M\+2\}\}makes this merged kick the identity, and \([18](https://arxiv.org/html/2608.00571#A1.E18)\) becomesCeγ​t∘StC\_\{e^\{\\gamma t\}\}\\circ S^\{t\}\.

*Admissibility of the potentials\.*By the abbreviations introduced in the proof of Theorem 1, the block potentials areA~ℓ​\(⋅;t,ξ\)=2t​Aℓ\\widetilde\{A\}\_\{\\ell\}\(\\,\\cdot\\,;t,\\xi\)=\\tfrac\{2\}\{t\}A\_\{\\ell\}andB~ℓ​\(⋅;t,ξ\)=1t​Bℓ\\widetilde\{B\}\_\{\\ell\}\(\\,\\cdot\\,;t,\\xi\)=\\tfrac\{1\}\{t\}B\_\{\\ell\}\. EachAℓA\_\{\\ell\}andBℓB\_\{\\ell\}is a fixed linear combination of the differences \([16](https://arxiv.org/html/2608.00571#A1.E16)\), and forV~\\widetilde\{V\}of classC1C^\{1\}intt

V~​\(t,⋅\)−V~​\(0,⋅\)t→t→0∂tV~​\(0,⋅\),\\frac\{\\widetilde\{V\}\(t,\\cdot\)\-\\widetilde\{V\}\(0,\\cdot\)\}\{t\}\\xrightarrow\[t\\to 0\]\{\}\\partial\_\{t\}\\widetilde\{V\}\(0,\\cdot\),soA~ℓ\\widetilde\{A\}\_\{\\ell\}andB~ℓ\\widetilde\{B\}\_\{\\ell\}extend continuously tot=0t=0and areC2C^\{2\}in the phase\-space variable for everyt∈\[0,T\]t\\in\[0,T\]\. Setaℓ=bℓ=0a\_\{\\ell\}=b\_\{\\ell\}=0, so thatA~ℓ=δ​Aℓ\\widetilde\{A\}\_\{\\ell\}=\\delta A\_\{\\ell\}andB~ℓ=δ​Bℓ\\widetilde\{B\}\_\{\\ell\}=\\delta B\_\{\\ell\}\. These networks are dense inC1C^\{1\}on compact subsets ofℝd×\[0,T\]\\mathbb\{R\}^\{d\}\\times\[0,T\]\[[22](https://arxiv.org/html/2608.00571#bib.bib22),[21](https://arxiv.org/html/2608.00571#bib.bib21)\], so for everyη\>0\\eta\>0there are parameters with‖∇δ​Aℓ−∇Aℓ‖∞≤η\\left\\lVert\\nabla\\delta A\_\{\\ell\}\-\\nabla A\_\{\\ell\}\\right\\rVert\_\{\\infty\}\\leq\\etaand‖∇δ​Bℓ−∇Bℓ‖∞≤η\\left\\lVert\\nabla\\delta B\_\{\\ell\}\-\\nabla B\_\{\\ell\}\\right\\rVert\_\{\\infty\}\\leq\\etaon the relevant compact sets\. Each factor of \([18](https://arxiv.org/html/2608.00571#A1.E18)\) is Lipschitz there, so composing3​L\+13L\+1of them gives a constantc\>0c\>0, depending only on those Lipschitz constants, with

supz∈Ksupt∈\[0,T\]‖Φ^θt​\(z\)−Ceγ​t​\(St​\(z\)\)‖≤c​η\.\\sup\_\{z\\in K\}\\sup\_\{t\\in\[0,T\]\}\\left\\lVert\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\(z\)\-C\_\{e^\{\\gamma t\}\}\\big\(S^\{t\}\(z\)\\big\)\\right\\rVert\\leq c\\,\\eta\.
*Conclusion\.*Sinceγ<0\\gamma<0, the mapCeγ​tC\_\{e^\{\\gamma t\}\}is Lipschitz with constanteγ​t≤1e^\{\\gamma t\}\\leq 1, so by \([14](https://arxiv.org/html/2608.00571#A1.E14)\)

‖Φ^θt​\(z\)−φt​\(z\)‖≤‖Φ^θt​\(z\)−Ceγ​t​\(St​\(z\)\)‖\+‖Ceγ​t​\(St​\(z\)\)−Ceγ​t​\(ψt​\(z\)\)‖≤c​η\+ε1\.\\left\\lVert\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\(z\)\-\\varphi^\{t\}\(z\)\\right\\rVert\\leq\\left\\lVert\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\(z\)\-C\_\{e^\{\\gamma t\}\}\\big\(S^\{t\}\(z\)\\big\)\\right\\rVert\+\\left\\lVert C\_\{e^\{\\gamma t\}\}\\big\(S^\{t\}\(z\)\\big\)\-C\_\{e^\{\\gamma t\}\}\\big\(\\psi^\{t\}\(z\)\\big\)\\right\\rVert\\leq c\\,\\eta\+\\varepsilon\_\{1\}\.Takingε1=ε/2\\varepsilon\_\{1\}=\\varepsilon/2andη=ε/\(2​c\)\\eta=\\varepsilon/\(2c\)gives the stated bound withL=M\+1L=M\+1and parametersθL\\theta\_\{L\}that do not depend ontt\. The potentials areC2C^\{2\}in the phase\-space variable for everyt∈\[0,T\]t\\in\[0,T\]and the allocations satisfy∑ℓαℓ​\(ξ\)=1\\sum\_\{\\ell\}\\alpha\_\{\\ell\}\(\\xi\)=1, so Theorem 1 gives\(Φ^θt\)∗​ω=eγ​t​ω\\big\(\\widehat\{\\Phi\}\_\{\\theta\}^\{t\}\\big\)^\{\*\}\\omega=e^\{\\gamma t\}\\omega\.

## Appendix BExperimental Details

### B\.1Data Generation

The experiments use two\-degree\-of\-freedom natural HamiltoniansH​\(q,p\)=12​∥p∥2\+V​\(q\)H\(q,p\)=\\tfrac\{1\}\{2\}\\lVert p\\rVert^\{2\}\+V\(q\)on the position domain𝒬=\[−3,3\]2\\mathcal\{Q\}=\[\-3,3\]^\{2\}, so thatHq=VH\_\{q\}=VandHp=12​∥p∥2H\_\{p\}=\\tfrac\{1\}\{2\}\\lVert p\\rVert^\{2\}in the block potentials of the main paper\. Each potential isV=Vq\+VGPV=V\_\{\\mathrm\{q\}\}\+V\_\{\\mathrm\{GP\}\}with a quartic backbone

Vq​\(q\)=12​\(b1​q12\+b2​q22\)\+c​q1​q2\+14​\(a1​q14\+a2​q24\+2​a3​q12​q22\)\.V\_\{\\mathrm\{q\}\}\(q\)=\\tfrac\{1\}\{2\}\\big\(b\_\{1\}q\_\{1\}^\{2\}\+b\_\{2\}q\_\{2\}^\{2\}\\big\)\+c\\,q\_\{1\}q\_\{2\}\+\\tfrac\{1\}\{4\}\\big\(a\_\{1\}q\_\{1\}^\{4\}\+a\_\{2\}q\_\{2\}^\{4\}\+2a\_\{3\}q\_\{1\}^\{2\}q\_\{2\}^\{2\}\\big\)\.The coefficients are drawn independently withb1,b2∼U​\[−1,2\.25\]b\_\{1\},b\_\{2\}\\sim U\[\-1,2\.25\],c∼U​\[−0\.5,0\.5\]c\\sim U\[\-0\.5,0\.5\]andai=0\.5​ui2a\_\{i\}=0\.5\\,u\_\{i\}^\{2\}withui∼U​\[0,1\]u\_\{i\}\\sim U\[0,1\], so the quartic coefficients are concentrated near zero\. Negativebib\_\{i\}are admitted, so the family contains single well, double well and multiple well landscapes\.

The perturbation is a Gaussian process with a squared exponential kernel, realized withD=128D=128random Fourier features,

VGP​\(q\)=β​2D​∑k=1Dak​cos⁡\(wk⊤​q\+ϕk\),wk∼𝒩​\(0,I2\)/ℓ,ϕk∼U​\[0,2​π\],ak∼𝒩​\(0,1\),V\_\{\\mathrm\{GP\}\}\(q\)=\\beta\\sqrt\{\\tfrac\{2\}\{D\}\}\\sum\_\{k=1\}^\{D\}a\_\{k\}\\cos\\big\(w\_\{k\}^\{\\top\}q\+\\phi\_\{k\}\\big\),\\qquad w\_\{k\}\\sim\\mathcal\{N\}\(0,I\_\{2\}\)/\\ell,\\quad\\phi\_\{k\}\\sim U\[0,2\\pi\],\\quad a\_\{k\}\\sim\\mathcal\{N\}\(0,1\),with length scaleℓ∼U​\[0\.6,1\.8\]\\ell\\sim U\[0\.6,1\.8\]\. The amplitudeβ\\betais set so that

β​maxq∈𝒢⁡∥∇VGP0​\(q\)∥=δ⋅medianq∈𝒢⁡∥∇Vq​\(q\)∥,\\beta\\,\\max\_\{q\\in\\mathcal\{G\}\}\\lVert\\nabla V\_\{\\mathrm\{GP\}\}^\{\\,0\}\(q\)\\rVert=\\delta\\cdot\\operatorname\{median\}\_\{q\\in\\mathcal\{G\}\}\\lVert\\nabla V\_\{\\mathrm\{q\}\}\(q\)\\rVert,whereVGP0V\_\{\\mathrm\{GP\}\}^\{\\,0\}is the unscaled feature sum and𝒢\\mathcal\{G\}the descriptor grid\. The parameterδ\\deltais therefore the ratio of the largest perturbation force to the median backbone force, and it is drawn log\-uniformly withδ∈\[0\.2,1\.5\]\\delta\\in\[0\.2,1\.5\]\. Since every training system hasδ≥0\.2\\delta\\geq 0\.2, the benchmark systems, which haveδ=0\\delta=0, lie outside the training distribution\. The length scale satisfiesℓ≥0\.6\\ell\\geq 0\.6, which is3\.23\.2times the descriptor spacing6/32=0\.18756/32=0\.1875, so the sampled potentials are resolved by the descriptor\. The dissipation parameter is drawn log\-uniformly with\|γ\|∈\[0\.1,0\.4\]\\lvert\\gamma\\rvert\\in\[0\.1,0\.4\]\.

The descriptor is the vector of values ofVVon a uniform32×3232\\times 32grid over𝒬\\mathcal\{Q\}\. In the notation of the main paper this corresponds to sensor pointsri=\(qi,0\)r\_\{i\}=\(q\_\{i\},0\)withqiq\_\{i\}on that grid, for whichH​\(ri\)=V​\(qi\)H\(r\_\{i\}\)=V\(q\_\{i\}\)\.

A proposed system is accepted only ifV∂−Vmin≥0\.5V\_\{\\partial\}\-V\_\{\\min\}\\geq 0\.5, whereVminV\_\{\\min\}is the minimum ofVVover the grid andV∂V\_\{\\partial\}its minimum on the grid boundary\. This rejection removes systems whose well is too shallow to confine a trajectory\. The energy budget is thenEmax=Vmin\+0\.9​\(V∂−Vmin\)E\_\{\\max\}=V\_\{\\min\}\+0\.9\\,\(V\_\{\\partial\}\-V\_\{\\min\}\)\.

Initial conditions are drawn by rejection\. Positions are uniform on𝒬\\mathcal\{Q\}, momenta arep=r​u​pmaxp=r\\,u\\,p\_\{\\max\}withuuuniform on the unit circle,r=v1/2r=v^\{1/2\}forv∼U​\[0,1\]v\\sim U\[0,1\]andpmax=2​\(Emax−Vmin\)p\_\{\\max\}=\\sqrt\{2\(E\_\{\\max\}\-V\_\{\\min\}\)\}, and a draw is kept whenH​\(q,p\)≤EmaxH\(q,p\)\\leq E\_\{\\max\}\. The pool is oversampled by a factor of3232and doubled up to five times if a system does not reach its quota\. Along solutions the energy is nonincreasing, so a trajectory starting belowEmaxE\_\{\\max\}cannot reach a part of the domain whereVVexceedsEmaxE\_\{\\max\}\. SinceV∂V\_\{\\partial\}is evaluated on the boundary of the descriptor grid rather than on∂𝒬\\partial\\mathcal\{Q\}, this was checked numerically on4040validation systems, comprising12801280trajectories: no trajectory left𝒬\\mathcal\{Q\}, none exceeded its energy budget by more than10−610^\{\-6\}, and the energy was nonincreasing along every stored trajectory to the same tolerance\.

Reference trajectories are integrated with a fourth\-order Runge–Kutta scheme in double precision, using10241024substeps on the fixed interval\[0,T\]\[0,T\]withT=4T=4, and each trajectory is stored at6565equally spaced states\. Initial conditions are rounded to single precision before integration, so the stored initial state coincides with the first stored trajectory state\.

The generator draws40,00040\{,\}000systems from a single seed, of which1,0001\{,\}000are assigned to validation and1,0001\{,\}000to test, leaving38,00038\{,\}000for training, with3232trajectories per system\. The system\-defining parametersw,ϕ,a,βw,\\phi,a,\\betaand the backbone coefficients are stored in double precision, so the analytic potential is reconstructed exactly at training and evaluation time; trajectories are stored in single precision\.

Training pairs\(z,τ,z\+\)\(z,\\tau,z^\{\+\}\)are formed online, wherezzandz\+z^\{\+\}are two stored states of the same trajectory andτ\\tauis the elapsed time between them\. Each step draws512512systems and1616pairs per system, for a batch of81928192\. The time gaps are stratified: the index range\[1,65\)\[1,65\)is divided into1616contiguous bins and one gap is drawn from each bin per system, so every batch contains the full range of query times\.

This construction is used in place of an existing collection because the setting requires the Hamiltonian and the dissipation parameter of every system to be known exactly, both entering the architecture directly; a family in which the distance from the training distribution is controlled, which is what places the benchmark systems outside it; and reference solutions accurate enough that the structure error can be measured at the level of double precision round\-off\. Published Hamiltonian trajectory datasets provide trajectories for a small number of fixed systems and do not supply these three properties\.

### B\.2Sampling Ranges and Settings

Table 4:Sampling ranges, integration settings and hyperparameters\.QuantityValueDomain and descriptorPosition domain𝒬\\mathcal\{Q\}\[−3,3\]2\[\-3,3\]^\{2\}Descriptor grid32×3232\\times 32Grid spacing0\.18750\.1875System familyb1,b2b\_\{1\},b\_\{2\}U​\[−1,2\.25\]U\[\-1,2\.25\]ccU​\[−0\.5,0\.5\]U\[\-0\.5,0\.5\]a1,a2,a3a\_\{1\},a\_\{2\},a\_\{3\}0\.5​U​\[0,1\]20\.5\\,U\[0,1\]^\{2\}Random Fourier featuresDD128128Length scaleℓ\\ellU​\[0\.6,1\.8\]U\[0\.6,1\.8\]Perturbation ratioδ\\deltalog⁡U​\[0\.2,1\.5\]\\log U\[0\.2,1\.5\]Dissipation\|γ\|\\lvert\\gamma\\rvertlog⁡U​\[0\.1,0\.4\]\\log U\[0\.1,0\.4\]Acceptance thresholdV∂−VminV\_\{\\partial\}\-V\_\{\\min\}≥0\.5\\geq 0\.5Energy budget fraction0\.90\.9IntegrationSchemeRK4, double precisionHorizonTT44Substeps10241024Step sizehh3\.906×10−33\.906\\times 10^\{\-3\}Stored states per trajectory6565DatasetSystems \(train / val / test\)38,00038\{,\}000/1,0001\{,\}000/1,0001\{,\}000Trajectories per system3232Generator seed2026072020260720OptimizationSteps80,00080\{,\}000Systems per step512512Pairs per system1616Batch size81928192OptimizerAdamWLearning rate10−410^\{\-4\}, cosine to10−610^\{\-6\}Warmup steps30003000Weight decay10−410^\{\-4\}Gradient clipping0\.50\.5Validation interval10001000stepsValidation systems512512Training seeds0,11,22ArchitectureBlocksLL1212FNO width / modes / layers3232/88/44Condition dimension128128Correction network width6464Softmax logit bound88Soft\-penalty weightλsoft\\lambda\_\{\\mathrm\{soft\}\}11Long\-horizon evaluation\|γ\|\\lvert\\gamma\\rvert0\.10\.1Windows1010Query times per window3232Initial conditions6464Reference RK4 steps per unit time512512The final checkpoint is the one with the lowest validation loss\. No hyperparameter search was carried out; the values above were fixed before the reported runs\.

### B\.3Benchmark Systems

The four benchmark systems haveV=VqV=V\_\{\\mathrm\{q\}\}with the coefficients of Table[5](https://arxiv.org/html/2608.00571#A2.T5), that is, no Gaussian process perturbation\. Initial conditions are drawn from each system’s own energy sublevel set by the same rejection scheme as the training data, withEmaxE\_\{\\max\}computed on the native32×3232\\times 32grid\. Reference trajectories use RK4 in double precision with step size approximately10−210^\{\-2\}, which agrees with an adaptive high\-order reference to about10−910^\{\-9\}on all four systems\.

Table 5:Backbone coefficients of the four benchmark systems\.Figure[5](https://arxiv.org/html/2608.00571#A2.F5)shows the resulting potentials\.

![Refer to caption](https://arxiv.org/html/2608.00571v1/x5.png)Figure 5:PotentialsVVof the four benchmark systems on\[−3,3\]2\[\-3,3\]^\{2\}\.
### B\.4Computing Environment

All models were trained and evaluated on one MIG partition of profile1g\.12gbof an NVIDIA H100, providing approximately11\.511\.5GB of device memory\. The software stack is PyTorch 2\.3\.1 with CUDA 12\.1 on Red Hat Enterprise Linux 8\. Latency figures reported elsewhere in this document were measured on the same partition and are not comparable with full\-card measurements\.

## Appendix CAdditional Results

### C\.1Per\-System Energy Error

Table[6](https://arxiv.org/html/2608.00571#A3.T6)gives the per\-system energy errors summarized by the last column of the main paper\. DeepONet diverges on B3, which is why its mean is dominated by that system\.

Table 6:Long\-horizon energy errorℰH\\mathcal\{E\}\_\{H\}on the four benchmark systems, averaged overt∈\(0,10​T\]t\\in\(0,10T\]and over6464initial conditions at\|γ\|=0\.1\\lvert\\gamma\\rvert=0\.1, under the protocol of the main paper\.
### C\.2Energy Error over the Long Horizon

The energy errors of Table[6](https://arxiv.org/html/2608.00571#A3.T6)are averages overt∈\(0,10​T\]t\\in\(0,10T\]\. Figure[6](https://arxiv.org/html/2608.00571#A3.F6)and Figure[7](https://arxiv.org/html/2608.00571#A3.F7)resolve them in time, under the same protocol and at\|γ\|=0\.1\\lvert\\gamma\\rvert=0\.1\.

![Refer to caption](https://arxiv.org/html/2608.00571v1/x6.png)\(a\)B1, anisotropic oscillator
![Refer to caption](https://arxiv.org/html/2608.00571v1/x7.png)\(b\)B2, coupled Duffing
![Refer to caption](https://arxiv.org/html/2608.00571v1/x8.png)\(c\)B3, Mexican hat
![Refer to caption](https://arxiv.org/html/2608.00571v1/x9.png)\(d\)B4, quartic coupled

Figure 6:Energy errorℰH​\(t\)\\mathcal\{E\}\_\{H\}\(t\)for CoSynFlow and the five baselines over0≤t≤10​T0\\leq t\\leq 10Tat\|γ\|=0\.1\\lvert\\gamma\\rvert=0\.1, averaged over6464initial conditions, on a logarithmic scale\.![Refer to caption](https://arxiv.org/html/2608.00571v1/x10.png)\(a\)B1, anisotropic oscillator
![Refer to caption](https://arxiv.org/html/2608.00571v1/x11.png)\(b\)B2, coupled Duffing
![Refer to caption](https://arxiv.org/html/2608.00571v1/x12.png)\(c\)B3, Mexican hat
![Refer to caption](https://arxiv.org/html/2608.00571v1/x13.png)\(d\)B4, quartic coupled

Figure 7:Energy errorℰH​\(t\)\\mathcal\{E\}\_\{H\}\(t\)for CoSynFlow, Symplectic Flow and the two ablations, under the protocol of Figure[6](https://arxiv.org/html/2608.00571#A3.F6)\.In the position plane and in the time series ofq1q\_\{1\}, the three structured models are visually indistinguishable over the whole horizon, so those plots are not reproduced here for the ablation group\. Their differences appear only in the aggregate errors of Table[6](https://arxiv.org/html/2608.00571#A3.T6)and of the main paper\.

### C\.3Structure Error under Composition

The main paper reports the structure error att=2t=2and after composing the learned map sixteen times, a total time of7\.5​T7\.5T\. Table[7](https://arxiv.org/html/2608.00571#A3.T7)gives both\.

Table 7:Conformal symplectic structure error, evaluated in double precision atγ=−0\.3\\gamma=\-0\.3\.
### C\.4Physics\-Informed Training

Table[8](https://arxiv.org/html/2608.00571#A3.T8)gives the full results of the physics\-informed experiment\. The model is CoSynFlow instantiated for the single system B2 at\|γ\|=0\.2\\lvert\\gamma\\rvert=0\.2, trained on two reference trajectories for40004000steps withλ=1\\lambda=1, with three seeds\. Errors are measured on6464held\-out trajectories, within the trained horizon and after composition to10​T10T\.

Table 8:Physics\-informed training of CoSynFlow on a single system with two training trajectories, three seeds\. The last column is the ratio of the two means\.For each of the three metrics, every seed trained with the residual term achieves a lower error than every seed trained without it\. The relative standard deviation across seeds falls from45%45\\%to5\.8%5\.8\\%for the state error within the trained horizon\.

## References

- \[1\]Lu Lu, Pengzhan Jin, and George Em Karniadakis\.Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators\.arXiv preprint arXiv:1910\.03193, 2019\.
- \[2\]Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar\.Fourier neural operator for parametric partial differential equations\.arXiv preprint arXiv:2010\.08895, 2020\.
- \[3\]Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar\.Neural operator: Learning maps between function spaces with applications to pdes\.Journal of Machine Learning Research, 24\(89\):1–97, 2023\.
- \[4\]Samuel Greydanus, Misko Dzamba, and Jason Yosinski\.Hamiltonian neural networks\.Advances in neural information processing systems, 32, 2019\.
- \[5\]Pengzhan Jin, Zhen Zhang, Aiqing Zhu, Yifa Tang, and George Em Karniadakis\.Sympnets: Intrinsic structure\-preserving symplectic networks for identifying hamiltonian systems\.Neural Networks, 132:166–179, 2020\.
- \[6\]Priscilla Canizares, Davide Murari, Carola\-Bibiane Schönlieb, Ferdia Sherry, and Zakhar Shumaylov\.Symplectic neural flows for modeling and discovery\.arXiv preprint arXiv:2412\.16787, 2024\.
- \[7\]Robert McLachlan and Matthew Perlmutter\.Conformal hamiltonian systems\.Journal of Geometry and Physics, 39\(4\):276–300, 2001\.
- \[8\]Zongyi Li, Hongkai Zheng, Nikola Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzadenesheli, and Anima Anandkumar\.Physics\-informed neural operator for learning partial differential equations\.ACM/IMS Journal of Data Science, 1\(3\):1–27, 2024\.
- \[9\]Haixu Wu, Huakun Luo, Haowen Wang, Jianmin Wang, and Mingsheng Long\.Transolver: A fast transformer solver for pdes on general geometries\.arXiv preprint arXiv:2402\.02366, 2024\.
- \[10\]Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho\.Lagrangian neural networks\.arXiv preprint arXiv:2003\.04630, 2020\.
- \[11\]Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty\.Symplectic ode\-net: Learning hamiltonian dynamics with control\.arXiv preprint arXiv:1909\.12077, 2019\.
- \[12\]Pengzhan Jin, Zhen Zhang, Ioannis G Kevrekidis, and George Em Karniadakis\.Learning poisson systems and trajectories of autonomous systems via poisson neural networks\.IEEE Transactions on Neural Networks and Learning Systems, 34\(11\):8271–8283, 2022\.
- \[13\]Yuhan Chen, Takashi Matsubara, and Takaharu Yaguchi\.Neural symplectic form: Learning hamiltonian equations on general coordinate systems\.Advances in Neural Information Processing Systems, 34:16659–16670, 2021\.
- \[14\]Andrew Sosanya and Sam Greydanus\.Dissipative hamiltonian neural networks: Learning dissipative and conservative dynamics separately\.arXiv preprint arXiv:2201\.10085, 2022\.
- \[15\]Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty\.Dissipative symoden: Encoding hamiltonian dynamics with dissipation and control into deep learning\.arXiv preprint arXiv:2002\.08860, 2020\.
- \[16\]Shaan A Desai, Marios Mattheakis, David Sondak, Pavlos Protopapas, and Stephen J Roberts\.Port\-hamiltonian neural networks for learning explicit time\-dependent dynamical systems\.Physical Review E, 104\(3\):034312, 2021\.
- \[17\]Ge\-Peng Ji, Deng\-Ping Fan, Yu\-Cheng Chou, Dengxin Dai, Alexander Liniger, and Luc Van Gool\.Deep gradient learning for efficient camouflaged object detection\.Machine Intelligence Research, 20\(1\):92–108, 2023\.
- \[18\]Elena Celledoni, Brynjulf Owren, Chong Shen, Baige Xu, and Takaharu Yaguchi\.Uepi: Universal energy\-behavior\-preserving integrators for energy conservative/dissipative differential equations\.Advances in Neural Information Processing Systems, 38:43717–43751, 2026\.
- \[19\]Zhen Zhang, Yeonjong Shin, and George Em Karniadakis\.Gfinns: Generic formalism informed neural networks for deterministic and stochastic dynamical systems\.Philosophical Transactions of the Royal Society A, 380\(2229\):20210207, 2022\.
- \[20\]Dmitry Turaev\.Polynomial approximations of symplectic dynamics and richness of chaos in non\-hyperbolic area\-preserving maps\.Nonlinearity, 16\(1\):123–135, 2003\.
- \[21\]Allan Pinkus\.Approximation theory of the mlp model in neural networks\.Acta numerica, 8:143–195, 1999\.
- \[22\]Kurt Hornik\.Approximation capabilities of multilayer feedforward networks\.Neural networks, 4\(2\):251–257, 1991\.