Verifier-Guided Model Discovery for Physical Dynamical Systems with Pretrained Symbolic Transformers

# Verifier-Guided Model Discovery for Physical Dynamical Systems with Pretrained Symbolic Transformers
Source: [https://arxiv.org/html/2608.02662](https://arxiv.org/html/2608.02662)
Farbod Faraji Francesco Belardinelli Department of Computing, Huxley Building, Imperial College London 180 Queen’s Gate, South Kensington, London SW7 2RH, United Kingdom

###### Abstract

Reliable forecasting of nonlinear physical systems underpins scientific discovery and engineering decision\-making\. Yet high\-fidelity simulations are prohibitively costly, and machine\-learning surrogates can be opaque and encode assumptions about system dynamics, limiting generalizability\. Pretrained transformers mapping synthetic ODE trajectories to equations offer interpretable alternatives, promising transfer without system\-specific equation knowledge\. Transferring them reliably to high\-dimensional physical data, however, remains an open challenge\. We develop a verifier\-guided \(VG\) workflow around ODEFormer as a symbolic backbone, using dynamical and physical\-admissibility criteria to select from a multi\-trajectory candidate equation pool, enabling transfer\. On canonical Van der Pol oscillators, VG outperforms the original ODEFormer workflow across held\-out initial conditions\. We then address vortex shedding—a phenomenon occurring in atmospheric and plasma systems of societal relevance—through coordinate reduction and symbolic discovery at fixed and varying Reynolds numbers\. VG discovers fixed\-parameter reduced\-order equations that recover the fundamental shedding oscillator and higher harmonics without a wake\-specific candidate library or prescribed Navier–Stokes structure, while the cross\-parameter model generalizes to withheld regimes\. Reconstruction fidelity alone did not determine symbolic discoverability, highlighting the importance of compatibility between latent dynamics and the backbone’s pretraining distribution\. This work establishes a verifier\-guided neural\-to\-symbolic methodology for interpretable and physically auditable forecasting in the natural sciences\.

## 1Introduction

In order to predict the evolution of nonlinear physical systems, computational efficiency alone is insufficient when the predictions of a model cannot be interpreted or audited\. For AI in the natural sciences, forecasts should be accompanied by evidence that the inferred dynamics are physically credible and by clearly identified limits\. Symbolic reduced\-order models support such scrutiny: their explicit equations can be examined against dynamical and physical requirements\(Faraji and Reza[2025](https://arxiv.org/html/2608.02662#bib.bib37)\), while remaining efficient enough for parameter exploration, design optimization, and control\.

Pretrained symbolic transformers map trajectories to such explicit models without case\-specific candidate libraries\. Their reliable transfer from synthetic pretraining to physical systems, however, requires more than generating equations reproducing individual trajectories\. Pretraining uses synthetic equation–trajectory pairs drawn from a prescribed symbolic vocabulary, while reduced physical coordinates may omit unresolved dynamics needed for a closed autonomous description\. Generated equations may still be unstable, fail across initial conditions or parameters, or violate physical requirements\. We address this challenge with a verifier\-guided \(VG\) workflow that uses executable tests of dynamical behavior and physical admissibility to select equations across multiple trajectories\.

Testing this approach beyond directly observed ODEs requires a high\-dimensional physical problem\. Vortex shedding provides such a setting with broad scientific and societal relevance: atmospheric von Kármán vortex streets occur in island wakes\(Japan Meteorological Agency[2019](https://arxiv.org/html/2608.02662#bib.bib2); Etling[1990](https://arxiv.org/html/2608.02662#bib.bib3)\); related instabilities drive vortex\-induced loading\(Williamson and Govardhan[2004](https://arxiv.org/html/2608.02662#bib.bib4)\)and arise in plasma flows past obstacles\(Gruszeckiet al\.[2010](https://arxiv.org/html/2608.02662#bib.bib5); Bailunget al\.[2020](https://arxiv.org/html/2608.02662#bib.bib6)\)\. The controlled cylinder\-flow case examined in this work thus connects auditable reduced\-order discovery to atmospheric observation, aerospace and transportation, infrastructure safety, and space\-plasma environments\.

The principal contributions of this work are:

- •Verifier\-guided symbolic discovery\.We place executable tests of dynamical behavior and physical admissibility at the center of multi\-trajectory model selection using a pretrained symbolic transformer\.
- •Generalization across initial conditions\.On fixed Van der Pol oscillators, VG outperforms the baseline single\-trajectory transformer pipeline for every held\-out initial condition while closely recovering the governing structure\.
- •Transfer to high\-dimensional physical data\.At fixed and varying Reynolds numbers, VG recovers shedding dynamics and yields stable predictions at withheld interpolation and extrapolation values; representation fidelity alone does not determine symbolic discoverability\.

## 2Related Work

#### Predictive and physics\-informed dynamics learning\.

Data\-driven modeling of dynamical systems offers several routes to efficient forecasting\. Optimized dynamic mode decomposition \(OPT\-DMD\) identifies coherent spatiotemporal modes and their best\-fit linear evolution, reducing noise\-induced bias and, when constrained, enabling stable forecasting\(Askham and Kutz[2018](https://arxiv.org/html/2608.02662#bib.bib10); Farajiet al\.[2024](https://arxiv.org/html/2608.02662#bib.bib11)\)\. Neural ODEs parameterize continuous\-time vector fields, while neural operators, including the Fourier Neural Operator, learn mappings between function spaces for families of PDE solutions\(Chenet al\.[2018](https://arxiv.org/html/2608.02662#bib.bib12); Liet al\.[2021](https://arxiv.org/html/2608.02662#bib.bib13)\)\. These approaches do not generally recover explicit nonlinear governing laws: OPT\-DMD uses linear evolution in its coordinates, while neural ODEs and neural operators retain opaque parameterizations\. Physics\-informed neural networks incorporate known governing\-equation residuals and boundary or initial conditions, whereas structure\-preserving architectures encode properties such as Hamiltonian conservation\(Raissiet al\.[2019](https://arxiv.org/html/2608.02662#bib.bib14); Greydanuset al\.[2019](https://arxiv.org/html/2608.02662#bib.bib15)\)\. Such methods require the relevant physical structure to be known and supplied in advance\. Residual penalties promote but do not guarantee constraint satisfaction; hard constraints can enforce known, expressible requirements\(Luet al\.[2021](https://arxiv.org/html/2608.02662#bib.bib16)\)\.

#### Sparse\-regression approaches to symbolic discovery\.

Symbolic dynamics discovery instead seeks explicit governing equations\. Sparse Identification of Nonlinear Dynamics \(SINDy\) identifies parsimonious systems through sparse regression over a prescribed candidate library\(Bruntonet al\.[2016](https://arxiv.org/html/2608.02662#bib.bib9)\)\. Weak\-form, ensemble, and Bayesian variants reduce sensitivity to numerical differentiation, improve robustness to limited or noisy data, and quantify model uncertainty\(Messenger and Bortz[2021](https://arxiv.org/html/2608.02662#bib.bib17); Faselet al\.[2022](https://arxiv.org/html/2608.02662#bib.bib18); Funget al\.[2025](https://arxiv.org/html/2608.02662#bib.bib19)\)\. Physical structure can also enter regression: constrained sparse Galerkin regression combines POD reduction with energy\-preserving constraints in fluid reduced\-order models\(Loiseau and Brunton[2018](https://arxiv.org/html/2608.02662#bib.bib20)\)\. Phi Method identifies discretized evolution operators from a candidate library, avoiding explicit derivatives and continuous\-time integration\(Farajiet al\.[2025](https://arxiv.org/html/2608.02662#bib.bib21)\)\. AutoSINDy\(Basiri and Nicholson[2026](https://arxiv.org/html/2608.02662#bib.bib22)\)combines PySR candidate generation\(Cranmer[2023](https://arxiv.org/html/2608.02662#bib.bib36)\), library curation, and sparse identification to mitigate the difficulty of prescribing an expressive library\.

#### Search\-based and generative symbolic dynamics discovery\.

Many general\-purpose search\-based and neural generative symbolic regression methods were initially developed for static relations,y=f​\(𝐱\)y=f\(\\mathbf\{x\}\), with later extensions to dynamics\. DynAIFeynman applies AI Feynman’s separability, symmetry, and dimensional\-consistency tests to states paired with finite\-difference derivatives, while ProGED samples grammar\-defined structures and fits their parameters against estimated derivatives or by simulating candidate ODEs\(Udrescu and Tegmark[2020](https://arxiv.org/html/2608.02662#bib.bib23); Weilbachet al\.[2021](https://arxiv.org/html/2608.02662#bib.bib24); Brenceet al\.[2021](https://arxiv.org/html/2608.02662#bib.bib25); Omejcet al\.[2024](https://arxiv.org/html/2608.02662#bib.bib26)\)\. Neural generators learn distributions over expressions: Deep Symbolic Regression uses a recurrent generator with risk\-seeking policy gradients, while pretrained transformer approaches include NeSymReS, SymbolicGPT, and end\-to\-end symbolic regression; TPSR further adds Monte Carlo tree search with accuracy and complexity feedback\(Petersenet al\.[2021](https://arxiv.org/html/2608.02662#bib.bib27); Biggioet al\.[2021](https://arxiv.org/html/2608.02662#bib.bib41); Valipouret al\.[2021](https://arxiv.org/html/2608.02662#bib.bib28); Kamiennyet al\.[2022](https://arxiv.org/html/2608.02662#bib.bib42); Shojaeeet al\.[2023](https://arxiv.org/html/2608.02662#bib.bib29)\)\. ODEFormer specializes this formulation to dynamics discovery by mapping multivariate trajectories directly to coupled ODE systems\(d’Ascoliet al\.[2024](https://arxiv.org/html/2608.02662#bib.bib30)\)\. MIO trains a trajectory\-to\-equation transformer to infer shared dynamics from multiple trajectories of the same system\(Şahinet al\.[2025](https://arxiv.org/html/2608.02662#bib.bib35)\)\. Diffusion\-based symbolic regression provides a non\-autoregressive alternative through masked iterative denoising and dataset\-specific reinforcement learning\(Bastianiet al\.[2025](https://arxiv.org/html/2608.02662#bib.bib31)\)\. These approaches expand or accelerate symbolic search but retain representational choices through observed variables, primitives, grammar, objectives, or pretraining\.

#### Toward knowledge\-guided and verified discovery\.

Recent work increasingly makes scientific knowledge operational during equation discovery\. LLM\-assisted physics\-informed symbolic regression adds language\-model assessments to the search objective; prior\-guided methods use executable constraint programs to steer evolutionary search; and Latent Grammar Flow embeds stability constraints in grammar rules or conditions generation on them\(Taskinet al\.[2026](https://arxiv.org/html/2608.02662#bib.bib32); Xiaoet al\.[2026](https://arxiv.org/html/2608.02662#bib.bib33); Yuet al\.[2026](https://arxiv.org/html/2608.02662#bib.bib34)\)\. These approaches place scientific consistency within search or generation\. Our complementary objective is to transfer an existing pretrained symbolic transformer without modifying or retraining its generator\. We leave our chosen pretrained backbone, ODEFormer\(d’Ascoliet al\.[2024](https://arxiv.org/html/2608.02662#bib.bib30)\), unchanged and use executable dynamical and physical\-admissibility tests on pooled candidates generated independently across trajectories to select a common equation system\. This enables transfer from synthetic pretraining to reduced physical coordinates while retaining symbolic transparency and avoiding a case\-specific candidate library or prescribed governing structure\.

## 3Background and Problem Formulation

ConsiderMMobserved trajectories𝐱\(m\)​\(t\)\\mathbf\{x\}^\{\(m\)\}\(t\), wheremmindexes different initial conditions, time windows, or parameter values\. Symbolic dynamics discovery seeks an explicit model

𝐱˙=𝐟​\(𝐱\),𝐟:𝒳⊆ℝd→ℝd,\\dot\{\\mathbf\{x\}\}=\\mathbf\{f\}\(\\mathbf\{x\}\),\\qquad\\mathbf\{f\}:\\mathcal\{X\}\\subseteq\\mathbb\{R\}^\{d\}\\rightarrow\\mathbb\{R\}^\{d\},\(1\)whose initial\-value problems yield finite trajectories over the state domain and time horizons of interest\. Equation \([1](https://arxiv.org/html/2608.02662#S3.E1)\) treats the supplied coordinates as a sufficient state for a deterministic, autonomous, first\-order description\. We do not prescribe a case\-specific parametric form or candidate library for𝐟\\mathbf\{f\}in this work; the representable candidates are nevertheless shaped by the learned distribution of the symbolic backbone\.

We use ODEFormer as this backbone because it maps multivariate trajectories directly to coupled symbolic ODEs without numerical derivative targets, case\-specific retraining, or a manually constructed function library\(d’Ascoliet al\.[2024](https://arxiv.org/html/2608.02662#bib.bib30)\)\. Its encoder processes tokenized time–state observations, and its autoregressive decoder constructs the right\-hand sides of the ODE system as symbolic sequences in prefix notation\. ODEFormer thus induces an implicit structural prior through its synthetic pretraining distribution, including the symbolic vocabulary, expression\-tree complexity, coefficient and trajectory distributions, and training\-time filtering of unstable or uninformative trajectories\. At inference, beam sampling produces multiple candidate equation systems\(d’Ascoliet al\.[2024](https://arxiv.org/html/2608.02662#bib.bib30)\)\. The beam size controls how many candidate sequences are explored, while the temperature controls the concentration of the token distribution and thus candidate diversity\. We use a beam size of2020and temperature of0\.10\.1; their sensitivity is examined on the canonical Van der Pol problem in Technical Appendix[S1](https://arxiv.org/html/2608.02662#A1)\. Whereas the original ODEFormer workflow selects candidates by reconstruction of a single input trajectory, we reconsider how the generated candidates are evaluated and selected \(Section[4](https://arxiv.org/html/2608.02662#S4)\)\.

For a spatially distributed physical state𝝎​\(t\)∈ℝN\\boldsymbol\{\\omega\}\(t\)\\in\\mathbb\{R\}^\{N\}, our objective is a symbolic ODE describing the temporal evolution of its dominant coordinates, rather than direct discovery of the underlying PDE\. Let an encoderℰ\\mathcal\{E\}and reconstruction mapℛ\\mathcal\{R\}define

𝐳​\(t\)=ℰ​\(𝝎​\(t\)\),𝝎​\(t\)≈ℛ​\(𝐳​\(t\)\),𝐳˙=𝐠​\(𝐳\),\\mathbf\{z\}\(t\)=\\mathcal\{E\}\\\!\\left\(\\boldsymbol\{\\omega\}\(t\)\\right\),\\qquad\\boldsymbol\{\\omega\}\(t\)\\approx\\mathcal\{R\}\\\!\\left\(\\mathbf\{z\}\(t\)\\right\),\\qquad\\dot\{\\mathbf\{z\}\}=\\mathbf\{g\}\(\\mathbf\{z\}\),\(2\)where𝐳​\(t\)∈ℝd\\mathbf\{z\}\(t\)\\in\\mathbb\{R\}^\{d\}andd≪Nd\\ll N\. Eq\. \([2](https://arxiv.org/html/2608.02662#S3.E2)\) provides a low\-dimensional phase\-space model whose equilibria, limit cycles, stability, frequencies, and modal couplings can be examined using established dynamical\-systems tools\. The primary role of coordinate reduction here is to define an interpretable state for the dominant temporal dynamics\. The retained coordinates may omit higher\-order dynamics important for closure; in that case,𝐠\\mathbf\{g\}represents an effective autonomous approximation within the chosen coordinates and the backbone’s symbolic hypothesis class\.

For trajectories observed at several values of a governing parameterpp, we append the parameter to the state and constrain it to remain constant:

𝐳˙=𝐠​\(𝐳,p\),p˙=0\.\\dot\{\\mathbf\{z\}\}=\\mathbf\{g\}\(\\mathbf\{z\},p\),\\qquad\\dot\{p\}=0\.\(3\)The augmentation in Eq\. \([3](https://arxiv.org/html/2608.02662#S3.E3)\) preserves the autonomous form expected by the backbone while allowing a single symbolic system to represent parameter\-dependent dynamics\. The resulting problem is to select an equation system that remains dynamically and physically admissible across multiple trajectories and at parameter values excluded from symbolic discovery\.

## 4Methodology: Verifier\-Guided Symbolic Model Discovery

Figure[1](https://arxiv.org/html/2608.02662#S4.F1)illustrates the VG workflow\. Directly observed ODE states enter unchanged, whereas high\-dimensional fields are first reduced by POD or an autoencoder\. The frozen backbone processes each trajectory or window separately; generated equations are then pooled, ranked, verified, and coefficient\-refined\.

![Refer to caption](https://arxiv.org/html/2608.02662v1/x1.png)Figure 1:Verifier\-guided symbolic model discovery\. Independently decoded candidates are pooled, ranked, verified, and coefficient\-refined before an admissible symbolic ODE is returned\.#### Candidate generation, pooling, and ranking\.

Let𝝃\(m\)\\boldsymbol\{\\xi\}^\{\(m\)\}denote the trajectory or window supplied to the backbone:𝝃\(m\)=𝐱\(m\)\\boldsymbol\{\\xi\}^\{\(m\)\}=\\mathbf\{x\}^\{\(m\)\}for observed ODE states,𝝃\(m\)=𝐳\(m\)\\boldsymbol\{\\xi\}^\{\(m\)\}=\\mathbf\{z\}^\{\(m\)\}for reduced fields, and𝝃\(m\)=\[𝐳\(m\)⊤,p\(m\)\]⊤\\boldsymbol\{\\xi\}^\{\(m\)\}=\[\\mathbf\{z\}^\{\(m\)\\top\},p^\{\(m\)\}\]^\{\\top\}for cross\-parameter discovery\. For inputmm, the backbone retainsBmB\_\{m\}decoded candidate systems indexed bybb:

𝒞m=\{𝐡m,b\}b=1Bm,𝒞=⋃m=1M𝒞m\.\\mathcal\{C\}\_\{m\}=\\\{\\mathbf\{h\}\_\{m,b\}\\\}\_\{b=1\}^\{B\_\{m\}\},\\qquad\\mathcal\{C\}=\\bigcup\_\{m=1\}^\{M\}\\mathcal\{C\}\_\{m\}\.\(4\)In Eq\. \([4](https://arxiv.org/html/2608.02662#S4.E4)\), symbolic canonicalization removes duplicates before the union, which constitutes pooling\.

Each𝐡∈𝒞\\mathbf\{h\}\\in\\mathcal\{C\}is ranked across discovery trajectories by Eq\. \([5](https://arxiv.org/html/2608.02662#S4.E5)\):

ℒroll​\(𝐡\)=1M​∑m=1M1d​∑j=1dMSEi⁡\(ξ^i​j\(m\),ξi​j\(m\)\)Vari⁡\(ξi​j\(m\)\)\+ϵ,\\mathcal\{L\}\_\{\\mathrm\{roll\}\}\(\\mathbf\{h\}\)=\\frac\{1\}\{M\}\\sum\_\{m=1\}^\{M\}\\frac\{1\}\{d\}\\sum\_\{j=1\}^\{d\}\\frac\{\\operatorname\{MSE\}\_\{i\}\\left\(\\widehat\{\\xi\}\_\{ij\}^\{\(m\)\},\\xi\_\{ij\}^\{\(m\)\}\\right\)\}\{\\operatorname\{Var\}\_\{i\}\\left\(\\xi\_\{ij\}^\{\(m\)\}\\right\)\+\\epsilon\},\(5\)whereddis the number of modeled dynamical coordinates,iiindexes time samples, and𝝃^\(m\)\\widehat\{\\boldsymbol\{\\xi\}\}^\{\(m\)\}is the candidate rollout initialized from the first observed state of trajectorymm\. The appended parameter is excluded from theddcoordinates in cross\-parameter discovery\. Nonfinite rollouts receive infinite loss, and the lowest\-loss candidates form the shortlist𝒞K\\mathcal\{C\}\_\{K\}; for cross\-parameter discovery, the shortlist also retains the best candidates with explicit parameter dependence\.

#### Executable verification and selection\.

The central component of VG is a suite of executable verifiers applied to𝒞K\\mathcal\{C\}\_\{K\}\. Each candidate defines the vector field𝝃˙=𝐡​\(𝝃\)\\dot\{\\boldsymbol\{\\xi\}\}=\\mathbf\{h\}\(\\boldsymbol\{\\xi\}\)\. Equation \([6](https://arxiv.org/html/2608.02662#S4.E6)\) measures its local agreement with observed dynamics:

Vloc​\(𝐡\)=1d​∑j=1dMSEm,i⁡\(hj​\(𝝃i\(m\)\),ξ˙i​j\(m\)\)Varm,i⁡\(ξ˙i​j\(m\)\)\+ϵV\_\{\\mathrm\{loc\}\}\(\\mathbf\{h\}\)=\\frac\{1\}\{d\}\\sum\_\{j=1\}^\{d\}\\frac\{\\operatorname\{MSE\}\_\{m,i\}\\left\(h\_\{j\}\(\\boldsymbol\{\\xi\}\_\{i\}^\{\(m\)\}\),\\dot\{\\xi\}\_\{ij\}^\{\(m\)\}\\right\)\}\{\\operatorname\{Var\}\_\{m,i\}\\left\(\\dot\{\\xi\}\_\{ij\}^\{\(m\)\}\\right\)\+\\epsilon\}\(6\)whereξ˙i​j\(m\)\\dot\{\\xi\}\_\{ij\}^\{\(m\)\}is the numerically estimated rate of coordinatejjat observed sampleii, and the error and variance are evaluated across discovery trajectories and times\. A low score therefore indicates that the candidate vector field reproduces the locally observed directions and rates of state evolution\.

Complementing this local test, rollout\-based verifiers evaluate each candidate from prescribed starting states: observed states selected as independent initial conditions\. They require finite and bounded evolution, agreement with observed oscillation amplitude and frequency or period, and consistency of these features across starting states\. For cross\-parameter systems, they additionally require a zero parameter equation, negligible integrated parameter drift, and explicit parameter dependence in at least one state equation\.

Equation \([7](https://arxiv.org/html/2608.02662#S4.E7)\) defines the admissible set and selects within it:

𝒜\\displaystyle\\mathcal\{A\}=\{𝐡∈𝒞K:𝗉𝖺𝗌𝗌ℓ​\(𝐡\)=1​for every verifier​ℓ\},\\displaystyle=\\left\\\{\\mathbf\{h\}\\in\\mathcal\{C\}\_\{K\}:\\mathsf\{pass\}\_\{\\ell\}\(\\mathbf\{h\}\)=1\\ \\text\{for every verifier \}\\ell\\right\\\},\(7\)𝐡⋆\\displaystyle\\mathbf\{h\}^\{\\star\}=arg⁡min𝐡∈𝒜​ℒroll​\(𝐡\)\.\\displaystyle=\\underset\{\\mathbf\{h\}\\in\\mathcal\{A\}\}\{\\arg\\min\}\\,\\mathcal\{L\}\_\{\\mathrm\{roll\}\}\(\\mathbf\{h\}\)\.Verification thus acts as an executable model\-selection constraint rather than a post\-hoc diagnostic\. These tests establish empirical admissibility over prescribed domains and horizons\. Technical Appendix[S2](https://arxiv.org/html/2608.02662#A2)reports their detailed definitions, thresholds, case\-specific settings, and tolerance\-sensitivity evaluation\.

#### Identifiability and coefficient refinement\.

Viewed as an augmented observation map, the verifiers contract the set of trajectory\-consistent candidates whenever they reject an otherwise indistinguishable model while retaining an admissible one\. They can therefore improve identifiability within the generated candidate class; Technical Appendix[S2](https://arxiv.org/html/2608.02662#A2)formalizes this perspective and its local sensitivity\-rank interpretation\.

With the selected symbolic structure fixed, its decoded constants are jointly refined by Nelder–Mead minimization ofℒroll\\mathcal\{L\}\_\{\\mathrm\{roll\}\}across the discovery trajectories\. The refined equation is reverified; if it becomes inadmissible, VG returns the admissible unoptimized equation𝐡⋆\\mathbf\{h\}^\{\\star\}\.

## 5Empirical Evaluation

Evaluation considers the canonical Van der Pol \(VdP\) oscillator and high\-dimensional flow past a stationary cylinder\. VdP tests whether VG generalizes across unseen initial conditions, while the cylinder\-flow problem tests transfer from synthetic ODE pretraining to reduced vortex\-shedding dynamics at fixed and withheld values of Reynolds number\.

### 5\.1Canonical Dynamical System Case: Van der Pol Oscillator

The Van der Pol oscillator is a canonical self\-excited system whose dynamics range from nearly harmonic motion to increasingly pronounced relaxation oscillations\(van der Pol[1926](https://arxiv.org/html/2608.02662#bib.bib1)\)\. Its first\-order form is given in Eq\. \([8](https://arxiv.org/html/2608.02662#S5.E8)\):

x˙0\\displaystyle\\dot\{x\}\_\{0\}=x1,\\displaystyle=x\_\{1\},\(8\)x˙1\\displaystyle\\dot\{x\}\_\{1\}=μ​\(1−x02\)​x1−x0\\displaystyle=\\mu\(1\-x\_\{0\}^\{2\}\)x\_\{1\}\-x\_\{0\}=μ​x1−x0−μ​x02​x1\.\\displaystyle=\\mu x\_\{1\}\-x\_\{0\}\-\\mu x\_\{0\}^\{2\}x\_\{1\}\.
In Eq\. \([8](https://arxiv.org/html/2608.02662#S5.E8)\),x˙0=x1\\dot\{x\}\_\{0\}=x\_\{1\}makesx1x\_\{1\}velocity\-like; inx˙1\\dot\{x\}\_\{1\},−x0\-x\_\{0\}restores whileμ​x1\\mu x\_\{1\}injects energy near the origin and−μ​x02​x1\-\\mu x\_\{0\}^\{2\}x\_\{1\}dissipates it at larger amplitudes\. Their balance produces a stable limit cycle\. We considerμ=0\.5\\mu=0\.5and1\.51\.5as separate fixed systems, representing weak and more pronounced nonlinear relaxation, respectively\. This provides a controlled test of whether a discovered equation represents the surrounding vector field and transfers across initial conditions, rather than merely reproducing one observed orbit\.

#### Experimental setup\.

At each value ofμ\\mu, the original ODEFormer protocol\(d’Ascoliet al\.[2024](https://arxiv.org/html/2608.02662#bib.bib30)\)generated, ranked, and coefficient\-optimized candidates from a single trajectory\. We repeated this protocol for three preset initial conditions and report, for each test initial condition, the median across successful rollouts\. VG instead decoded eight preset windows, ranked the pooled candidates over 24 windows from 12 preset initial conditions, verified the ten highest\-ranked candidates, and applied the same coefficient optimization\. Both protocols were evaluated on the same eight unseen initial conditions\. Candidate\-pool and test initial conditions were sampled independently from\[−1,1\]2\[\-1,1\]^\{2\}using fixed seeds and were disjoint; Technical Appendix[S3](https://arxiv.org/html/2608.02662#A3)reports the exact assignments\.

![Refer to caption](https://arxiv.org/html/2608.02662v1/x2.png)Figure 2:Held\-out rollout errors across eight test initial conditions\. Circles denote the original ODEFormer protocol and squares the VG workflow; open and filled markers indicate pre\- and post\-optimization results\. Gray lines pair identical initial conditions, and black bars show medians\.
#### Results\.

Figure[2](https://arxiv.org/html/2608.02662#S5.F2)shows that the VG model outperforms the single\-trajectory protocol for all eight held\-out initial conditions at both values ofμ\\mu\. Forμ=0\.5\\mu=0\.5, the median post\-optimization error decreases from1\.211\.21to6\.99×10−46\.99\\times 10^\{\-4\}; forμ=1\.5\\mu=1\.5, it decreases from4\.5×10−24\.5\\times 10^\{\-2\}to3\.27×10−33\.27\\times 10^\{\-3\}\. The larger separation atμ=0\.5\\mu=0\.5is consistent with its longer transient: a single observed orbit constrains less of the surrounding vector field, whereas the VG workflow evaluates a common equation across multiple approaches to the limit cycle\. One repeat of the original ODEFormer protocol failed on the test set, and coefficient optimization degraded another\.

At eachμ\\mu, the candidate with the lowest multi\-trajectory rollout error was also the only member of the ten\-candidate shortlist to satisfy all verifier checks\. Verification therefore in this canonical case distinguished a generalizable closed dynamical model from alternatives that reproduced only individual features, such as the oscillation period or limit\-cycle geometry\.

After coefficient optimization, Eqs\. \([9](https://arxiv.org/html/2608.02662#S5.E9)\)–\([10](https://arxiv.org/html/2608.02662#S5.E10)\) were selected:

μ=0\.5:x˙0\\displaystyle\\mu=5:\\qquad\\dot\{x\}\_\{0\}=0\.997​x1−0\.012​x0,\\displaystyle=997x\_\{1\}\-012x\_\{0\},\(9\)x˙1\\displaystyle\\dot\{x\}\_\{1\}=0\.542​x1−0\.997​x0\\displaystyle=542x\_\{1\}\-997x\_\{0\}−0\.543​x02​x1−0\.0039​x0​x1,\\displaystyle\\quad\-543x\_\{0\}^\{2\}x\_\{1\}\-0039x\_\{0\}x\_\{1\},
and

μ=1\.5:x˙0\\displaystyle\\mu=5:\\quad\\dot\{x\}\_\{0\}=0\.996​x1\\displaystyle=996x\_\{1\}\(10\)\+0\.000676​\(0\.998\+1\.273​x0\)2,\\displaystyle\\quad\+000676\(998\+273x\_\{0\}\)^\{2\},x˙1\\displaystyle\\dot\{x\}\_\{1\}=1\.440​x1−1\.002​x0\\displaystyle=440x\_\{1\}\-002x\_\{0\}−1\.461​x02​x1−0\.0080​x0​x1\.\\displaystyle\\quad\-461x\_\{0\}^\{2\}x\_\{1\}\-0080x\_\{0\}x\_\{1\}\.
Both systems recover the defining Van der Pol structure: the kinematic relationx˙0≈x1\\dot\{x\}\_\{0\}\\approx x\_\{1\}, linear restoring dynamics through−x0\-x\_\{0\}, and the opposing linear and cubic damping terms inx˙1\\dot\{x\}\_\{1\}\. Their principal coefficients closely approximate the corresponding ground\-truth values, while the additionalx0x\_\{0\},x0​x1x\_\{0\}x\_\{1\}, and weak quadratic contributions provide small corrections\. The discovered equations are therefore dynamically close approximations rather than exact symbolic replicas of the ground\-truth system, supporting evaluation through behavior across initial conditions and physical admissibility rather than expression matching alone\.

### 5\.2High\-Dimensional Physical System Case: Flow Past a Stationary Cylinder at a Fixed Reynolds Number

The flow\-past\-a\-cylinder case provides the first test of whether the VG workflow can recover compact dynamics from high\-dimensional physical fields rather than directly observed ODE states\. From a physics standpoint, above the onset of wake instability in this test case configuration, vortices detach alternately from the separated shear layers and form a periodic von Kármán street\. For two\-dimensional incompressible flow, the nondimensional spanwise vorticityω\\omegasatisfies

∂ω∂t\+\(𝐮⋅∇\)​ω\\displaystyle\\frac\{\\partial\\omega\}\{\\partial t\}\+\(\\mathbf\{u\}\\cdot\\nabla\)\\omega=1R​e​∇2ω,\\displaystyle=\\frac\{1\}\{Re\}\\nabla^\{2\}\\omega,\(11\)∇⋅𝐮\\displaystyle\\nabla\\cdot\\mathbf\{u\}=0,R​e=U∞​Dν\.\\displaystyle=0,\\qquad Re=\\frac\{U\_\{\\infty\}D\}\{\\nu\}\.
In Eq\. \([11](https://arxiv.org/html/2608.02662#S5.E11)\),\(𝐮⋅∇\)​ω\(\\mathbf\{u\}\\cdot\\nabla\)\\omegarepresents advection of vorticity by the local flow, whereasR​e−1​∇2ωRe^\{\-1\}\\nabla^\{2\}\\omegarepresents viscous diffusion of vorticity gradients\. IncreasingR​eRereduces the relative contribution of this diffusion and permits the wake instability to develop\. AtR​e=300Re=300, the fixed value considered here, the two\-dimensional simulation exhibits established periodic shedding\. As the physical state is a spatial field containing 60,000 vorticity values at each time, dimensionality reduction is required before applying the VG workflow for symbolic model discovery\. This is performed using proper orthogonal decomposition \(POD\)\(Sirovich[1987](https://arxiv.org/html/2608.02662#bib.bib7)\)\.

#### Experimental setup\.

Vorticity snapshots were generated with the open\-sourceViscousFlow\.jlsolver\(Eldredge[2021](https://arxiv.org/html/2608.02662#bib.bib38)\), which implements the immersed\-layer method\(Eldredge[2022](https://arxiv.org/html/2608.02662#bib.bib8)\)\. A unit\-diameter cylinder was centered inx/D∈\[−1,5\]x/D\\in\[\-1,5\]andy/D∈\[−2,2\]y/D\\in\[\-2,2\], withU∞=1U\_\{\\infty\}=1, a2∘2^\{\\circ\}free\-stream incidence, and a no\-slip surface\. The outputs were interpolated onto a300×200300\\times 200grid and sampled atΔ​t∗=Δ​t​U∞/D=0\.1\\Delta t^\{\*\}=\\Delta tU\_\{\\infty\}/D=0\.1, yielding 1,000 snapshots\. Further numerical and preprocessing details are reported byFarajiet al\.\([2025](https://arxiv.org/html/2608.02662#bib.bib21)\)\.

The first 500 snapshots were discarded to isolate established shedding\. The remainder was divided chronologically into 250 development, 100 validation, and 150 final\-test snapshots\. The development interval defined the reduced coordinates and supported candidate generation, ranking, verification, and coefficient optimization; the validation set was used to select the POD representation and decoding configuration before testing on the remaining 150 snapshots\.

After subtracting the development\-interval mean, snapshot POD produced near\-equal mode pairs representing quadrature components of the periodic wake\. The first three pairs yielded six standardized coordinates and retained 94\.52% of the fluctuation energy\. Validation selected this rank\-6 representation and eight decoding windows from the tested ranks44and66and window counts88and2424\. Candidates decoded from eight five\-time\-unit windows were pooled and ranked over 24 development windows; the ten highest\-ranked were verified, and the selected structure was optimized over the same 24 windows\. Technical Appendix[S4](https://arxiv.org/html/2608.02662#A4)reports the selection matrix and POD spectrum\.

#### Results\.

The selected VG model integrated successfully from all eight predeclared starting points in the final\-test interval\. The window and full\-interval coordinate errors in Table[1](https://arxiv.org/html/2608.02662#S5.T1)show that the discovered symbolic system retains the recurrent shedding dynamics across different phases of the cycle\. At the field level, the six\-mode representation incurs 12\.53% error before symbolic forecasting; the 15\.56% end\-to\-end result therefore reflects the combined effects of coordinate reduction and dynamics prediction\. Figure[3](https://arxiv.org/html/2608.02662#S5.F3)shows that the wake\-averaged vorticity preserves the dominant phase and amplitude relative to both the simulation and POD reconstruction\. Representative vorticity\-field snapshot comparisons separating POD truncation from symbolic model forecasting are provided in Technical Appendix[S4](https://arxiv.org/html/2608.02662#A4)\.

Table 1:Fixed\-R​eRefinal\-test performance of the selected VG model\. Coordinate rollout errors are mean\-square errors normalized by the variance of each standardized POD coordinate\. Field errors are temporal means of snapshot\-wise relativeL2L\_\{2\}errors\.![Refer to caption](https://arxiv.org/html/2608.02662v1/x3.png)Figure 3:Wake\-averaged nondimensional vorticity over the final\-test interval atR​e=300Re=300, evaluated over0\.5≤x/D≤50\.5\\leq x/D\\leq 5and\|y/D\|≤1\.5\\lvert y/D\\rvert\\leq 1\.5\.Equation \([12](https://arxiv.org/html/2608.02662#S5.E12)\) gives the coefficient\-optimized VG model in standardized rank\-6 POD coordinatesz1,…,z6z\_\{1\},\\ldots,z\_\{6\}:

z˙1\\displaystyle\\dot\{z\}\_\{1\}=1\.001​z2\\displaystyle=001z\_\{2\}\(12\)−0\.0679​\(11\.746−0\.8668​z2\)−1\.052,\\displaystyle\\quad\-0679\(1746\-8668z\_\{2\}\)^\{\-1\.052\},z˙2\\displaystyle\\dot\{z\}\_\{2\}=−1\.123​z1,\\displaystyle=\-123z\_\{1\},z˙3\\displaystyle\\dot\{z\}\_\{3\}=−2\.099​z4,\\displaystyle=\-099z\_\{4\},z˙4\\displaystyle\\dot\{z\}\_\{4\}=2\.120​z3,\\displaystyle=120z\_\{3\},z˙5\\displaystyle\\dot\{z\}\_\{5\}=3\.121​z6\\displaystyle=121z\_\{6\}−0\.1164​sin⁡\(0\.1144\+11\.036​z2\),\\displaystyle\\quad\-1164\\sin\(1144\+1036z\_\{2\}\),z˙6\\displaystyle\\dot\{z\}\_\{6\}=−0\.1458−3\.233​z5\.\\displaystyle=\-1458\-233z\_\{5\}\.
The dominant linear parts of the\(z1,z2\)\(z\_\{1\},z\_\{2\}\),\(z3,z4\)\(z\_\{3\},z\_\{4\}\), and\(z5,z6\)\(z\_\{5\},z\_\{6\}\)subsystems define oscillators with angular frequencies 1\.060, 2\.109, and 3\.177, respectively\. Their ratio of1:1\.99:3\.001:1\.99:3\.00identifies the fundamental shedding cycle and its second and third harmonics\. The reciprocal contribution toz˙1\\dot\{z\}\_\{1\}deforms the fundamental oscillator from exact linear motion, while the sinusoidal dependence ofz˙5\\dot\{z\}\_\{5\}onz2z\_\{2\}couples the fundamental and third\-harmonic coordinate pairs\. The constant term inz˙6\\dot\{z\}\_\{6\}produces a small shift in the equilibrium of the third oscillator\. Because the POD coordinates are global projections of the vorticity field, these terms describe coupling within the reduced representation and cannot be assigned independently to localized flow mechanisms\.

Notably, the VG workflow obtained this system without a wake\-specific candidate library or prescribed Navier–Stokes structure: the pretrained backbone searched its learned vocabulary, and verification selected an admissible model reproducing the shedding dynamics\.

### 5\.3Cross\-Parameter Extension: Flow Past a Cylinder across Reynolds Numbers

The fixed\-R​eReresult presented in the previous subsection motivates a more demanding test: whether the VG workflow can derive a single symbolic reduced\-order model capable of representing vortex\-shedding dynamics as the governing flow parameter, the Reynolds number, varies\.

#### Experimental setup\.

Model discovery usedR​e∈\{150,200,250,300,350,400,450\}Re\\in\\\{150,200,250,300,350,400,450\\\}\. The withheld test values wereR​e∈\{175,275,425\}Re\\in\\\{175,275,425\\\}for interpolation andR​e=500Re=500for extrapolation\.

The simulations and field preprocessing followed the fixed\-R​eRecase\. After mean subtraction, a common shallow autoencoder \(AE\) with one 256\-neuron hidden layer compressed each 60,000\-component snapshot into three latent coordinates\. A Reynolds\-number coordinate, centered and scaled over the discovery range, was appended and constrained to remain constant during rollout\. Candidates generated independently from the seven full\-length discovery trajectories were pooled and ranked across all seven; the shortlist was verified, and the selected structure was coefficient\-optimized over the same trajectories\.

#### Results\.

The selected shallow three\-coordinate autoencoder had a mean field\-reconstruction error of 5\.83%, the lowest latent roughness \(0\.0224\), and the highest spectral concentration \(0\.926\) among several tested encodings\. Roughness measures temporal curvature relative to first\-order variation, while spectral concentration measures the fraction of nonzero\-frequency power in the three dominant Fourier components\. Although several deeper or four\-coordinate encodings reconstructed more accurately, among the controlled encodings advanced to VG, only a deeper three\-coordinate case yielded a raw admissible equation, with rollout loss 2\.17 compared with 1\.04 for the selected shallow encoding; coefficient refinement reduced the latter to 0\.200 while preserving admissibility\. Reconstruction fidelity therefore did not determine symbolic discoverability\. Technical Appendix[S5](https://arxiv.org/html/2608.02662#A5)reports more details on the outcomes of architecture and initialization\-seed sensitivity, latent diagnostics, and corresponding symbolic\-model discovery\.

In standardized latent coordinatesz1,z2,z3z\_\{1\},z\_\{2\},z\_\{3\}, withrrthe normalized Reynolds\-number coordinate, Eq\. \([13](https://arxiv.org/html/2608.02662#S5.E13)\) gives the coefficient\-optimized system that satisfied all verifier checks:

z˙1\\displaystyle\\dot\{z\}\_\{1\}=1\.143​z3\+0\.0422​r\+0\.201​z2​z3,\\displaystyle=143z\_\{3\}\+0422r\+201z\_\{2\}z\_\{3\},\(13\)z˙2\\displaystyle\\dot\{z\}\_\{2\}=0\.798​z1,\\displaystyle=798z\_\{1\},z˙3\\displaystyle\\dot\{z\}\_\{3\}=−0\.995​z1,\\displaystyle=\-995z\_\{1\},r˙\\displaystyle\\dot\{r\}=0\.\\displaystyle=0\.
The1\.143​z31\.143z\_\{3\}and−0\.995​z1\-0\.995z\_\{1\}terms define the leading linear oscillator in the\(z1,z3\)\(z\_\{1\},z\_\{3\}\)coordinates, whose linearized angular frequency is1\.143×0\.995=1\.066\\sqrt\{1\.143\\times 0\.995\}=1\.066\. The evolution ofz2z\_\{2\}is coupled to the same cycle, while the quadratic contribution0\.201​z2​z30\.201z\_\{2\}z\_\{3\}feeds this coordinate back intoz˙1\\dot\{z\}\_\{1\}, introducing a nonlinear deformation of the underlying oscillator\. The additive0\.0422​r0\.0422rterm shifts the latent vector field with Reynolds number, andr˙=0\\dot\{r\}=0preserves the selected parameter value throughout each rollout\.

Becausez˙3=−\(0\.995/0\.798\)​z˙2\\dot\{z\}\_\{3\}=\-\(0\.995/0\.798\)\\dot\{z\}\_\{2\}, the combinationz3\+1\.247​z2z\_\{3\}\+1\.247z\_\{2\}is conserved to the precision of the reported coefficients\. Consequently, for each fixed value ofrr, trajectories of the discovered model remain on a two\-dimensional invariant surface within the three\-dimensional latent space\. The Reynolds\-number coordinate therefore conditions a shared nonlinear oscillator rather than acting as an independently evolving input\.

Table 2:Performance at Reynolds numbers withheld from symbolic discovery\. Latent NMSE is the coordinate\-wise mean rollout error normalized by the variance of each reference latent coordinate\. Field quantities are temporal means of snapshot\-wise relative spatialL2L\_\{2\}errors\. I and E denote interpolation and extrapolation, respectively\.The selected VG model produced finite full\-length rollouts at all four withheld Reynolds numbers\. Referring to Table[2](https://arxiv.org/html/2608.02662#S5.T2), generalization did not deteriorate monotonically from interpolation to extrapolation: theR​e=500Re=500extrapolation case produced the lowest latent NMSE \(0\.036\) and end\-to\-end field error \(13\.8%\), whereas theR​e=275Re=275interpolation case was the most difficult, with corresponding errors of 0\.247 and 33\.0%\. Figure[4](https://arxiv.org/html/2608.02662#S5.F4)provides important physical context for the latter value\. Although displacement of the predicted vortices increases the spatialL2L\_\{2\}error, the VG field retains the alternating wake structure and principal vortex locations\. The error is therefore strongly influenced by phase and spatial displacement accumulated during rollout rather than by a loss of the vortex\-shedding regime\. AtR​e=500Re=500, the close agreement between the autoencoder reconstruction and VG prediction is consistent with its comparatively low VG\-to\-AE field error of 11\.9%\.

Complete latent\-coordinate rollouts, wake\-averaged vorticity histories, additional field comparisons, and extended error metrics for all four test Reynolds numbers are provided in Technical Appendix[S6](https://arxiv.org/html/2608.02662#A6)\.

![Refer to caption](https://arxiv.org/html/2608.02662v1/x4.png)Figure 4:Simulation, autoencoder \(AE\) reconstruction, and VG\-predicted vorticity att∗=75t^\{\*\}=75for withheldR​e=275Re=275\(interpolation, I\) andR​e=500Re=500\(extrapolation, E\)\. All panels share a symmetric color scale\.

## 6Conclusions

The principal takeaway is that a frozen symbolic transformer pretrained on synthetic ODEs can support physical model discovery beyond its pretraining distribution when paired with multi\-trajectory executable verification\. Across both the Van der Pol and cylinder\-flow cases, transfer depended not only on generating plausible equations but on selecting a common model that satisfied the prescribed dynamical and physical\-admissibility tests\. Model interpretation is, however, coordinate\-dependent: the recovered Van der Pol equations closely approximate the governing structure in directly observed states, whereas the equations discovered in POD and autoencoder coordinates are effective reduced\-order laws rather than unique reductions of the Navier–Stokes equations\.

The representation results also clarify a central limitation\. Symbolic discovery requires the supplied coordinates to admit an approximately closed autonomous description compatible with the backbone’s learned vocabulary and pretraining distribution; unresolved modes, memory, forcing, or parameter dependence can prevent this\. Pooling and verification enable transfer with a frozen pretrained backbone but cannot supply dynamics or structures absent from its candidate class\. These limits motivate tighter neuro\-symbolic integration of learned representations and symbolic model discovery\(Marraet al\.[2024](https://arxiv.org/html/2608.02662#bib.bib39); Cranmeret al\.[2020](https://arxiv.org/html/2608.02662#bib.bib40)\), specifically through joint verified coordinate–dynamics learning that directly integrates executable constraints encoding fundamental physical principles \(e\.g\., conservation, invariance, and symmetry\) into search or generation\. These directions are particularly important for higher\-dimensional, multiscale systems whose reduced coordinates are only approximately closed, including many plasma systems, and would strengthen the basis for interpretable and auditable AI across the natural sciences\.

## Code and Data Availability

The code and data bundle supporting the results reported in this work will be made publicly available upon acceptance for publication\.

## References

- Variable projection methods for an optimized dynamic mode decomposition\.SIAM Journal on Applied Dynamical Systems17\(1\),pp\. 380–416\.External Links:[Document](https://dx.doi.org/10.1137/M1124176)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px1.p1.1)\.
- Y\. Bailung, B\. Chutia, T\. Deka, A\. Boruah, S\. K\. Sharma, S\. Kumar, J\. Chutia, Y\. Nakamura, and H\. Bailung \(2020\)Vortex formation in a strongly coupled dusty plasma flow past an obstacle\.Physics of Plasmas27\(12\),pp\. 123702\.External Links:[Document](https://dx.doi.org/10.1063/5.0022356)Cited by:[§1](https://arxiv.org/html/2608.02662#S1.p3.1)\.
- M\. A\. Basiri and C\. Nicholson \(2026\)Discovery of nonlinear dynamics with automated basis function generation\.arXiv preprint arXiv:2605\.09696\.External Links:[Document](https://dx.doi.org/10.48550/arXiv.2605.09696)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px2.p1.1)\.
- Z\. Bastiani, R\. M\. Kirby, J\. Hochhalter, and S\. Zhe \(2025\)Diffusion\-based symbolic regression\.arXiv preprint arXiv:2505\.24776\.External Links:[Link](https://arxiv.org/abs/2505.24776)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px3.p1.1)\.
- L\. Biggio, T\. Bendinelli, A\. Neitz, A\. Lucchi, and G\. Parascandolo \(2021\)Neural symbolic regression that scales\.InProceedings of the 38th International Conference on Machine Learning,Proceedings of Machine Learning Research, Vol\.139,pp\. 936–945\.Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px3.p1.1)\.
- J\. Brence, L\. Todorovski, and S\. Džeroski \(2021\)Probabilistic grammars for equation discovery\.Knowledge\-Based Systems224,pp\. 107077\.External Links:[Document](https://dx.doi.org/10.1016/j.knosys.2021.107077)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px3.p1.1)\.
- S\. L\. Brunton, J\. L\. Proctor, and J\. N\. Kutz \(2016\)Discovering governing equations from data by sparse identification of nonlinear dynamical systems\.Proceedings of the National Academy of Sciences113\(15\),pp\. 3932–3937\.External Links:[Document](https://dx.doi.org/10.1073/pnas.1517384113)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px2.p1.1)\.
- R\. T\. Q\. Chen, Y\. Rubanova, J\. Bettencourt, and D\. K\. Duvenaud \(2018\)Neural ordinary differential equations\.InAdvances in Neural Information Processing Systems,Vol\.31,pp\. 6571–6583\.Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px1.p1.1)\.
- M\. Cranmer, A\. Sanchez\-Gonzalez, P\. Battaglia, R\. Xu, K\. Cranmer, D\. Spergel, and S\. Ho \(2020\)Discovering symbolic models from deep learning with inductive biases\.InAdvances in Neural Information Processing Systems,Vol\.33,pp\. 17429–17442\.Cited by:[§6](https://arxiv.org/html/2608.02662#S6.p2.1)\.
- M\. Cranmer \(2023\)Interpretable machine learning for science with PySR and SymbolicRegression\.jl\.arXiv preprint arXiv:2305\.01582\.Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px2.p1.1)\.
- S\. d’Ascoli, S\. Becker, A\. Mathis, P\. Schwaller, and N\. Kilbertus \(2024\)ODEFormer: symbolic regression of dynamical systems with transformers\.InThe Twelfth International Conference on Learning Representations,External Links:[Link](https://openreview.net/forum?id=TzoHLiGVMo)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px3.p1.1),[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px4.p1.1),[§3](https://arxiv.org/html/2608.02662#S3.p2.2),[§5\.1](https://arxiv.org/html/2608.02662#S5.SS1.SSS0.Px1.p1.2)\.
- J\. D\. Eldredge \(2021\)ViscousFlow\.jl: a framework for simulating viscous incompressible flows\.Note:ZenodoExternal Links:[Document](https://dx.doi.org/10.5281/zenodo.4776967),[Link](https://doi.org/10.5281/zenodo.4776967)Cited by:[§5\.2](https://arxiv.org/html/2608.02662#S5.SS2.SSS0.Px1.p1.6)\.
- J\. D\. Eldredge \(2022\)A method of immersed layers on cartesian grids, with application to incompressible flows\.Journal of Computational Physics448,pp\. 110716\.External Links:[Document](https://dx.doi.org/10.1016/j.jcp.2021.110716)Cited by:[§5\.2](https://arxiv.org/html/2608.02662#S5.SS2.SSS0.Px1.p1.6)\.
- D\. Etling \(1990\)Mesoscale vortex shedding from large islands: a comparison with laboratory experiments of rotating stratified flows\.Meteorology and Atmospheric Physics43,pp\. 145–151\.External Links:[Document](https://dx.doi.org/10.1007/BF01028117)Cited by:[§1](https://arxiv.org/html/2608.02662#S1.p3.1)\.
- F\. Faraji, M\. Reza, A\. Knoll, and J\. N\. Kutz \(2024\)Dynamic mode decomposition for data\-driven analysis and reduced\-order modelling ofE×BE\\times Bplasmas: ii\. dynamics forecasting\.Journal of Physics D: Applied Physics57\(6\),pp\. 065202\.External Links:[Document](https://dx.doi.org/10.1088/1361-6463/ad0911)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px1.p1.1)\.
- F\. Faraji, M\. Reza, and A\. Knoll \(2025\)Discovery of discretized differential equations from data: benchmarking and application to a plasma system\.Journal of Applied Physics137\(12\),pp\. 123301\.External Links:[Document](https://dx.doi.org/10.1063/5.0254956)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px2.p1.1),[§5\.2](https://arxiv.org/html/2608.02662#S5.SS2.SSS0.Px1.p1.6)\.
- F\. Faraji and M\. Reza \(2025\)Machine learning applications to computational plasma physics and reduced\-order plasma modeling: a perspective\.Journal of Physics D: Applied Physics58\(10\),pp\. 102002\.External Links:[Document](https://dx.doi.org/10.1088/1361-6463/ada167)Cited by:[§1](https://arxiv.org/html/2608.02662#S1.p1.1)\.
- U\. Fasel, J\. N\. Kutz, B\. W\. Brunton, and S\. L\. Brunton \(2022\)Ensemble\-SINDy: robust sparse model discovery in the low\-data, high\-noise limit, with active learning and control\.Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences478\(2260\),pp\. 20210904\.External Links:[Document](https://dx.doi.org/10.1098/rspa.2021.0904)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px2.p1.1)\.
- L\. Fung, U\. Fasel, and M\. P\. Juniper \(2025\)Rapid bayesian identification of sparse nonlinear dynamics from scarce and noisy data\.Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences481\(2307\),pp\. 20240200\.External Links:[Document](https://dx.doi.org/10.1098/rspa.2024.0200)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px2.p1.1)\.
- S\. Greydanus, M\. Dzamba, and J\. Yosinski \(2019\)Hamiltonian neural networks\.InAdvances in Neural Information Processing Systems,Vol\.32,pp\. 15379–15389\.Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px1.p1.1)\.
- M\. Gruszecki, V\. M\. Nakariakov, T\. Van Doorsselaere, and T\. D\. Arber \(2010\)Phenomenon of alfvénic vortex shedding\.Physical Review Letters105\(5\),pp\. 055004\.External Links:[Document](https://dx.doi.org/10.1103/PhysRevLett.105.055004)Cited by:[§1](https://arxiv.org/html/2608.02662#S1.p3.1)\.
- Japan Meteorological Agency \(2019\)Collection of images captured by himawari\-8/9: karman vortex\.Note:[https://www\.jma\.go\.jp/jma/jma\-eng/satellite/introduction/image\.html](https://www.jma.go.jp/jma/jma-eng/satellite/introduction/image.html)Accessed: 2026\-07\-22Cited by:[§1](https://arxiv.org/html/2608.02662#S1.p3.1)\.
- P\. Kamienny, S\. d’Ascoli, G\. Lample, and F\. Charton \(2022\)End\-to\-end symbolic regression with transformers\.InAdvances in Neural Information Processing Systems,Vol\.35,pp\. 10269–10281\.Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px3.p1.1)\.
- Z\. Li, N\. Kovachki, K\. Azizzadenesheli, B\. Liu, K\. Bhattacharya, A\. Stuart, and A\. Anandkumar \(2021\)Fourier neural operator for parametric partial differential equations\.InInternational Conference on Learning Representations,External Links:[Link](https://openreview.net/forum?id=c8P9NQVtmnO)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px1.p1.1)\.
- J\. Loiseau and S\. L\. Brunton \(2018\)Constrained sparse galerkin regression\.Journal of Fluid Mechanics838,pp\. 42–67\.External Links:[Document](https://dx.doi.org/10.1017/jfm.2017.823)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px2.p1.1)\.
- L\. Lu, R\. Pestourie, W\. Yao, Z\. Wang, F\. Verdugo, and S\. G\. Johnson \(2021\)Physics\-informed neural networks with hard constraints for inverse design\.SIAM Journal on Scientific Computing43\(6\),pp\. B1105–B1132\.External Links:[Document](https://dx.doi.org/10.1137/21M1397908)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px1.p1.1)\.
- G\. Marra, S\. Dumančić, R\. Manhaeve, and L\. De Raedt \(2024\)From statistical relational to neurosymbolic artificial intelligence: a survey\.Artificial Intelligence328,pp\. 104062\.External Links:[Document](https://dx.doi.org/10.1016/j.artint.2023.104062)Cited by:[§6](https://arxiv.org/html/2608.02662#S6.p2.1)\.
- D\. A\. Messenger and D\. M\. Bortz \(2021\)Weak SINDy for partial differential equations\.Journal of Computational Physics443,pp\. 110525\.External Links:[Document](https://dx.doi.org/10.1016/j.jcp.2021.110525)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px2.p1.1)\.
- N\. Omejc, B\. Gec, J\. Brence, L\. Todorovski, and S\. Džeroski \(2024\)Probabilistic grammars for modeling dynamical systems from coarse, noisy, and partial data\.Machine Learning113,pp\. 7689–7721\.External Links:[Document](https://dx.doi.org/10.1007/s10994-024-06522-1)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px3.p1.1)\.
- B\. K\. Petersen, M\. Landajuela, T\. N\. Mundhenk, C\. P\. Santiago, S\. K\. Kim, and J\. T\. Kim \(2021\)Deep symbolic regression: recovering mathematical expressions from data via risk\-seeking policy gradients\.InInternational Conference on Learning Representations,External Links:[Link](https://openreview.net/forum?id=m5Qsh0kBQG)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px3.p1.1)\.
- M\. Raissi, P\. Perdikaris, and G\. E\. Karniadakis \(2019\)Physics\-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations\.Journal of Computational Physics378,pp\. 686–707\.External Links:[Document](https://dx.doi.org/10.1016/j.jcp.2018.10.045)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px1.p1.1)\.
- Y\. E\. Şahin, N\. Kilbertus, and S\. Becker \(2025\)Predicting symbolic ODEs from multiple trajectories\.InNeurIPS 2025 Workshop on Machine Learning and the Physical Sciences,External Links:[Document](https://dx.doi.org/10.48550/arXiv.2510.23295),[Link](https://arxiv.org/abs/2510.23295)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px3.p1.1)\.
- P\. Shojaee, K\. Meidani, A\. Barati Farimani, and C\. K\. Reddy \(2023\)Transformer\-based planning for symbolic regression\.InAdvances in Neural Information Processing Systems,Vol\.36,pp\. 45907–45919\.External Links:[Link](https://proceedings.neurips.cc/paper_files/paper/2023/hash/8ffb4e3118280a66b192b6f06e0e2596-Abstract-Conference.html)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px3.p1.1)\.
- L\. Sirovich \(1987\)Turbulence and the dynamics of coherent structures\. part i: coherent structures\.Quarterly of Applied Mathematics45\(3\),pp\. 561–571\.External Links:[Document](https://dx.doi.org/10.1090/qam/910462)Cited by:[§5\.2](https://arxiv.org/html/2608.02662#S5.SS2.p3.4)\.
- B\. Taskin, W\. Xie, and T\. Lazebnik \(2026\)Knowledge integration for physics\-informed symbolic regression using pre\-trained large language models\.Scientific Reports16\(1\),pp\. 1614\.External Links:[Document](https://dx.doi.org/10.1038/s41598-026-35327-6),[Link](https://doi.org/10.1038/s41598-026-35327-6)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px4.p1.1)\.
- S\. Udrescu and M\. Tegmark \(2020\)AI Feynman: a physics\-inspired method for symbolic regression\.Science Advances6\(16\),pp\. eaay2631\.External Links:[Document](https://dx.doi.org/10.1126/sciadv.aay2631)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px3.p1.1)\.
- M\. Valipour, B\. You, M\. Panju, and A\. Ghodsi \(2021\)SymbolicGPT: a generative transformer model for symbolic regression\.arXiv preprint arXiv:2106\.14131\.External Links:[Link](https://arxiv.org/abs/2106.14131)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px3.p1.1)\.
- Balth\. van der Pol \(1926\)On “relaxation\-oscillations”\.The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science2\(11\),pp\. 978–992\.External Links:[Document](https://dx.doi.org/10.1080/14786442608564127)Cited by:[§5\.1](https://arxiv.org/html/2608.02662#S5.SS1.p1.1)\.
- J\. Weilbach, S\. Gerwinn, C\. Weilbach, and M\. Kandemir \(2021\)Inferring the structure of ordinary differential equations\.arXiv preprint arXiv:2107\.07345\.External Links:[Link](https://arxiv.org/abs/2107.07345)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px3.p1.1)\.
- C\. H\. K\. Williamson and R\. Govardhan \(2004\)Vortex\-induced vibrations\.Annual Review of Fluid Mechanics36,pp\. 413–455\.External Links:[Document](https://dx.doi.org/10.1146/annurev.fluid.36.050802.122128)Cited by:[§1](https://arxiv.org/html/2608.02662#S1.p3.1)\.
- J\. Xiao, X\. Chen, J\. Peng, Q\. Wang, M\. Jia, Z\. Lai, G\. Yu, D\. Li, T\. Li, and J\. Liu \(2026\)Prior\-guided symbolic regression: towards scientific consistency in equation discovery\.arXiv preprint arXiv:2602\.13021\.External Links:[Document](https://dx.doi.org/10.48550/arXiv.2602.13021),[Link](https://arxiv.org/abs/2602.13021)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px4.p1.1)\.
- K\. Yu, E\. Chatzi, and G\. Kissas \(2026\)Neuro\-symbolic ODE discovery with latent grammar flow\.arXiv preprint arXiv:2604\.16232\.External Links:[Document](https://dx.doi.org/10.48550/arXiv.2604.16232),[Link](https://arxiv.org/abs/2604.16232)Cited by:[§2](https://arxiv.org/html/2608.02662#S2.SS0.SSS0.Px4.p1.1)\.

## Technical Appendices

These technical appendices provide the supplementary results and analyses referenced in the paper\. They report beam\-size–sampling\-temperature candidate\-generation sensitivity \(Section[S1](https://arxiv.org/html/2608.02662#A1)\), detailed verifier definitions and their identifiability interpretation \(Section[S2](https://arxiv.org/html/2608.02662#A2)\), Van der Pol test\-case details \(Section[S3](https://arxiv.org/html/2608.02662#A3)\), the fixed\-Reynolds\-number representation study \(Section[S4](https://arxiv.org/html/2608.02662#A4)\), the autoencoder sensitivity analysis \(Section[S5](https://arxiv.org/html/2608.02662#A5)\), and extended results from the multi\-Reynolds\-number test case \(Section[S6](https://arxiv.org/html/2608.02662#A6)\)\.

## Appendix S1Sensitivity to Beam Size and Sampling Temperature

The beam size and sampling temperature are candidate\-symbolic\-equation generation \(decoding\) settings inherited from the symbolic backbone\. The beam size determines how many candidate sequences are retained from each input window, while the temperature controls the concentration of the sampled token distribution and hence the diversity of the resulting equations\. Because both settings affect the composition and size of the candidate pool available to VG, their sensitivity was assessed before the final experiments\. The Van der Pol oscillator provides a suitable diagnostic because its states are observed directly and its governing structure is known, avoiding the additional influence of coordinate reduction\. This section concerns selection of the decoding settings only; Section[S3](https://arxiv.org/html/2608.02662#A3)reports the complete Van der Pol trajectory sets, selected equations, and final\-test behavior\.

### S1\.1Comparison protocol

Beam sizes 20 and 50 were crossed with temperatures 0\.1 and 0\.2, giving four pre\-specified decoding configurations\. For each of the two oscillator parameters,μ=0\.5\\mu=0\.5and1\.51\.5, the same 12 model\-construction trajectories were sampled over0≤t≤200\\leq t\\leq 20\. Each trajectory contributed an early window over0≤t≤100\\leq t\\leq 10and a late window over10≤t≤2010\\leq t\\leq 20, giving the 24\-window model\-construction bank\. Eight preset members of this bank, termed*candidate\-generation windows*below, were processed individually by VG to generate equation candidates\. These candidates were pooled and deduplicated before rollout ranking across all 24 model\-construction windows and evaluation by the verifier suite\. The same 24\-window bank was used for coefficient optimization of the selected system\. Performance was measured on eight additional validation trajectories that were disjoint from both the model\-construction trajectories and the final\-test trajectories\. Their exact initial conditions are reported in Section[S3](https://arxiv.org/html/2608.02662#A3)\.

For each candidate\-generation window, VG retained a number of sequences equal to the beam size\. The eight candidate\-generation windows therefore produced at most 160 candidates for beam size 20 and 400 for beam size 50 before canonical deduplication; the resulting pools contained 160 and 398–400 distinct systems, respectively\. Each pool was ranked by its mean rollout error across the 24 model\-construction windows\. Section[S2](https://arxiv.org/html/2608.02662#A2)defines the verifier metrics and acceptance decisions, and Table[S2](https://arxiv.org/html/2608.02662#A2.T2)summarizes their numerical settings for the Van der Pol and cylinder\-flow cases\. The ten highest\-ranked systems were evaluated using these verifiers\. The admissible system with the lowest rollout error was coefficient\-optimized over the same 24\-window bank and then verified again\. Random seeds, trajectories, window selections, verifier settings, and optimization procedure were identical across the four configurations\. When no member of a shortlist passed the pre\-optimization verifiers, that outcome was retained; the rollout\-ranked candidate was optimized only to complete the diagnostic comparison\.

### S1\.2Validation performance and admissibility

Table[S1](https://arxiv.org/html/2608.02662#A1.T1)reports the validation error of the lowest\-ranked pre\-optimization candidate satisfying the verifier suite and the corresponding post\-optimization result\. The error is the coordinate\-wise mean rollout error normalized by the variance of each reference state\. The configuration with beam size 20 and temperature 0\.1 produced an admissible system before and after coefficient optimization at both values ofμ\\mu, with the lowest post\-optimization validation median in each case\. Increasing the temperature at beam size 20 increased the remaining error, particularly forμ=1\.5\\mu=1\.5\. The larger beam did not provide a consistent benefit: at temperature 0\.1 no pre\-optimization candidate was admissible forμ=0\.5\\mu=0\.5, while the post\-optimization system forμ=1\.5\\mu=1\.5no longer satisfied the local vector\-field verifier\. Beam size 50 also expanded the candidate pool by approximately a factor of 2\.5 and correspondingly increased the rollout ranking burden\.

Table S1:Sensitivity of validation rollout error and verifier admissibility to beam size and sampling temperature\. “Pre\-opt\. candidate” indicates whether at least one member of the ten\-system shortlist passed every applicable verifier\. Bold rows identify the configuration carried forward\.∗No shortlist member passed all pre\-optimization verifiers\. The reported errors and post\-optimization decision describe the rollout\-ranked candidate retained for diagnostic purposes; this candidate was not treated as a verifier\-guided pre\-optimization selection\.

These results support beam size 20 and temperature 0\.1 as a reasonable operating point rather than indicating that the decoding settings are immaterial\. This configuration provided the most accurate admissible post\-optimization system for both oscillator parameters while using the smaller candidate pool\. It was consequently fixed before the final\-test comparison reported in Section[5\.1](https://arxiv.org/html/2608.02662#S5.SS1)and was carried forward to the cylinder\-flow experiments\.

### S1\.3Early\- and late\-window diagnostic

The selected decoding configuration was also examined separately on the 12 early and 12 late model\-construction windows to determine whether its aggregate rollout score concealed a marked difference between transient\-rich and established\-cycle behavior\. This comparison did not introduce an additional selection rule\. Forμ=0\.5\\mu=0\.5, the median pre\-optimization NMSE was 0\.114 on the early windows and 0\.0982 on the late windows; coefficient optimization reduced these values to7\.39×10−47\.39\\times 10^\{\-4\}and3\.19×10−43\.19\\times 10^\{\-4\}, respectively\. Forμ=1\.5\\mu=1\.5, the corresponding medians were 0\.00537 and 0\.00266 before optimization and1\.26×10−31\.26\\times 10^\{\-3\}and8\.07×10−48\.07\\times 10^\{\-4\}afterwards\. The higher errors on the early windows are consistent with their inclusion of trajectories approaching the limit cycle, which sample a broader portion of the surrounding vector field\. The low post\-optimization errors in both regimes show that the selected configuration was supported by accurate rollouts across both the transient\-rich early windows and the more regularly oscillatory late windows\.

## Appendix S2Verifier Definitions, Identifiability, and Sensitivity Considerations

### S2\.1Constrained model selection

The verifier\-guided workflow described in Section[4](https://arxiv.org/html/2608.02662#S4)of the paper separates dynamical admissibility from rollout\-based ranking\. Let𝒞K\\mathcal\{C\}\_\{K\}denote the shortlist obtained after pooling, deduplicating, and ranking candidates across the discovery trajectories\. For a candidate𝐡∈𝒞K\\mathbf\{h\}\\in\\mathcal\{C\}\_\{K\}, verifierℓ\\ellevaluates a metricmℓ​\(𝐡\)m\_\{\\ell\}\(\\mathbf\{h\}\)against a toleranceτℓ\\tau\_\{\\ell\}\. An inequality\-based decision can be written compactly as

𝗉𝖺𝗌𝗌ℓ​\(𝐡\)=𝕀​\[mℓ​\(𝐡\)≤τℓ\],\\mathsf\{pass\}\_\{\\ell\}\(\\mathbf\{h\}\)=\\mathbb\{I\}\\\!\\left\[m\_\{\\ell\}\(\\mathbf\{h\}\)\\leq\\tau\_\{\\ell\}\\right\],\(S1\)where𝕀​\[⋅\]\\mathbb\{I\}\[\\cdot\]is the indicator function andτℓ\\tau\_\{\\ell\}is the largest accepted value of the metric\. This form applies because the numerical metrics used here measure discrepancy, variability, drift, or growth, for which smaller values indicate closer agreement or greater admissibility\. Conditions on symbolic structure, such as an identically zero parameter equation, are evaluated directly\. The admissible set and selected model are then

𝒜=\{𝐡∈𝒞K:𝗉𝖺𝗌𝗌ℓ​\(𝐡\)=1​for every applicable verifier​ℓ\},𝐡⋆=arg⁡min𝐡∈𝒜​ℒroll​\(𝐡\)\.\\mathcal\{A\}=\\left\\\{\\mathbf\{h\}\\in\\mathcal\{C\}\_\{K\}:\\mathsf\{pass\}\_\{\\ell\}\(\\mathbf\{h\}\)=1\\ \\text\{for every applicable verifier \}\\ell\\right\\\},\\qquad\\mathbf\{h\}^\{\\star\}=\\underset\{\\mathbf\{h\}\\in\\mathcal\{A\}\}\{\\arg\\min\}\\,\\mathcal\{L\}\_\{\\mathrm\{roll\}\}\(\\mathbf\{h\}\)\.\(S2\)Equations \([S1](https://arxiv.org/html/2608.02662#A2.E1)\) and \([S2](https://arxiv.org/html/2608.02662#A2.E2)\) therefore define admissibility separately fromℒroll\\mathcal\{L\}\_\{\\mathrm\{roll\}\}\. Rollout error orders only the candidates that satisfy every applicable requirement\. When𝒜\\mathcal\{A\}is empty, no admissible symbolic model is selected\.

### S2\.2Verifier definitions and interpretation

#### Local vector\-field agreement\.

Let𝝃i\(m\)∈ℝdξ\\boldsymbol\{\\xi\}\_\{i\}^\{\(m\)\}\\in\\mathbb\{R\}^\{d\_\{\\xi\}\}be the full state at sampleiiof discovery trajectorymm\. Its numerically estimated rate in evaluated coordinatej=1,…,dj=1,\\ldots,dis the scalarξ˙i​j\(m\)\\dot\{\\xi\}\_\{ij\}^\{\(m\)\}, whilehj​\(𝝃i\(m\)\)h\_\{j\}\(\\boldsymbol\{\\xi\}\_\{i\}^\{\(m\)\}\)is the corresponding scalar component of the candidate vector field evaluated at the full state\. The local discrepancy and its admissibility condition are given by

Vloc​\(𝐡\)\\displaystyle V\_\{\\mathrm\{loc\}\}\(\\mathbf\{h\}\)=1d​∑j=1dMSEm,i⁡\(hj​\(𝝃i\(m\)\),ξ˙i​j\(m\)\)Varm,i⁡\(ξ˙i​j\(m\)\)\+ϵ,\\displaystyle=\\frac\{1\}\{d\}\\sum\_\{j=1\}^\{d\}\\frac\{\\operatorname\{MSE\}\_\{m,i\}\\left\(h\_\{j\}\(\\boldsymbol\{\\xi\}\_\{i\}^\{\(m\)\}\),\\dot\{\\xi\}\_\{ij\}^\{\(m\)\}\\right\)\}\{\\operatorname\{Var\}\_\{m,i\}\\left\(\\dot\{\\xi\}\_\{ij\}^\{\(m\)\}\\right\)\+\\epsilon\},\(S3\)𝗉𝖺𝗌𝗌loc​\(𝐡\)\\displaystyle\\mathsf\{pass\}\_\{\\mathrm\{loc\}\}\(\\mathbf\{h\}\)=𝕀​\[Vloc​\(𝐡\)≤2​min𝐡~∈𝒞K⁡Vloc​\(𝐡~\)\],\\displaystyle=\\mathbb\{I\}\\\!\\left\[V\_\{\\mathrm\{loc\}\}\(\\mathbf\{h\}\)\\leq 2\\min\_\{\\widetilde\{\\mathbf\{h\}\}\\in\\mathcal\{C\}\_\{K\}\}V\_\{\\mathrm\{loc\}\}\(\\widetilde\{\\mathbf\{h\}\}\)\\right\],whereϵ=10−12\\epsilon=10^\{\-12\}prevents division by zero\. For the directly observed and fixed\-parameter systems,dξ=dd\_\{\\xi\}=d; for the cross\-Re\\mathrm\{Re\}system,dξ=d\+1d\_\{\\xi\}=d\+1because the input also contains the normalized Reynolds coordinater=\(Re−300\)/150r=\(\\mathrm\{Re\}\-300\)/150, while the score is evaluated over theddlatent\-state equations\. The variance normalization in Eq\. \([S3](https://arxiv.org/html/2608.02662#A2.E3)\) gives each evaluated coordinate comparable influence despite different derivative scales\. Because derivative scale and numerical differentiation error also differ between the observed Van der Pol states, POD coordinates, and autoencoder coordinates, the tolerance is defined relative to the lowest local discrepancy within the same shortlist\. A candidate passes when its instantaneous direction and rate remain within a factor of two of that case\-specific reference\.

#### Long\-horizon boundedness\.

Let𝝃^𝐡\(s\)​\(t\)\\widehat\{\\boldsymbol\{\\xi\}\}\_\{\\mathbf\{h\}\}^\{\(s\)\}\(t\)denote the trajectory obtained by integrating candidate𝐡\\mathbf\{h\}from observed starting statess, and letBobs=maxm,i⁡‖𝝃i\(m\)‖2B\_\{\\mathrm\{obs\}\}=\\max\_\{m,i\}\\\|\\boldsymbol\{\\xi\}\_\{i\}^\{\(m\)\}\\\|\_\{2\}denote the largest state norm in the corresponding reference data\. Boundedness requires

𝗉𝖺𝗌𝗌bound​\(𝐡\)=𝕀​\[𝝃^𝐡\(s\)​\(t\)​is finite for all​s,t∧maxs,t⁡‖𝝃^𝐡\(s\)​\(t\)‖2≤cbound​Bobs\]\.\\mathsf\{pass\}\_\{\\mathrm\{bound\}\}\(\\mathbf\{h\}\)=\\mathbb\{I\}\\\!\\left\[\\widehat\{\\boldsymbol\{\\xi\}\}\_\{\\mathbf\{h\}\}^\{\(s\)\}\(t\)\\text\{ is finite for all \}s,t\\ \\land\\ \\max\_\{s,t\}\\left\\\|\\widehat\{\\boldsymbol\{\\xi\}\}\_\{\\mathbf\{h\}\}^\{\(s\)\}\(t\)\\right\\\|\_\{2\}\\leq c\_\{\\mathrm\{bound\}\}B\_\{\\mathrm\{obs\}\}\\right\]\.\(S4\)In Eq\. \([S4](https://arxiv.org/html/2608.02662#A2.E4)\), the norm and its reference scale are evaluated only over the latent coordinates and separately at each Reynolds number for the cross\-Re\\mathrm\{Re\}case\. The factor iscbound=1\.5c\_\{\\mathrm\{bound\}\}=1\.5for Van der Pol andcbound=3c\_\{\\mathrm\{bound\}\}=3for the reduced flow coordinates\. This check rejects numerically divergent solutions and trajectories whose state magnitude grows well beyond that represented by the observations\.

#### Oscillation amplitude\.

After removal of the initial verification interval, the amplitude of coordinatejjis estimated by the robust percentile half\-rangeAj​\(𝝃\)=12​\[Q0\.99​\(ξj\)−Q0\.01​\(ξj\)\]A\_\{j\}\(\\boldsymbol\{\\xi\}\)=\\tfrac\{1\}\{2\}\[Q\_\{0\.99\}\(\\xi\_\{j\}\)\-Q\_\{0\.01\}\(\\xi\_\{j\}\)\]\. LetAj\(m\)A\_\{j\}^\{\(m\)\}be the amplitude of reference trajectorymm, and letA^s,j\(m\)\\widehat\{A\}\_\{s,j\}^\{\(m\)\}be the corresponding amplitude of the candidate rollout from starting statess\. Ifℳv\\mathcal\{M\}\_\{v\}indexes the reference cases used for verification and𝒥A\\mathcal\{J\}\_\{A\}the coordinates whose amplitudes are evaluated, then

mamp​\(𝐡\)\\displaystyle m\_\{\\mathrm\{amp\}\}\(\\mathbf\{h\}\)=1\|ℳv\|​\|𝒥A\|​∑m∈ℳv∑j∈𝒥A\|medians⁡A^s,j\(m\)−Aj\(m\)\|\|Aj\(m\)\|\+ϵ,\\displaystyle=\\frac\{1\}\{\|\\mathcal\{M\}\_\{v\}\|\\,\|\\mathcal\{J\}\_\{A\}\|\}\\sum\_\{m\\in\\mathcal\{M\}\_\{v\}\}\\sum\_\{j\\in\\mathcal\{J\}\_\{A\}\}\\frac\{\\left\|\\operatorname\{median\}\_\{s\}\\widehat\{A\}\_\{s,j\}^\{\(m\)\}\-A\_\{j\}^\{\(m\)\}\\right\|\}\{\|A\_\{j\}^\{\(m\)\}\|\+\\epsilon\},\(S5\)𝗉𝖺𝗌𝗌amp​\(𝐡\)\\displaystyle\\mathsf\{pass\}\_\{\\mathrm\{amp\}\}\(\\mathbf\{h\}\)=𝕀​\[mamp​\(𝐡\)≤0\.15\]\.\\displaystyle=\\mathbb\{I\}\\\!\\left\[m\_\{\\mathrm\{amp\}\}\(\\mathbf\{h\}\)\\leq 0\.15\\right\]\.For Van der Pol,𝒥A\\mathcal\{J\}\_\{A\}contains the observed oscillatory coordinate, andℳv\\mathcal\{M\}\_\{v\}contains one pooled reference case whose target is formed from the discovery trajectories\. For the fixed\-Re\\mathrm\{Re\}flow case,ℳv\\mathcal\{M\}\_\{v\}contains the singleRe=300\\mathrm\{Re\}=300reference trajectory; for the cross\-Re\\mathrm\{Re\}case, it contains the three reference trajectories atRe=150,300,\\mathrm\{Re\}=150,300,and450450\. All retained coordinates belong to𝒥A\\mathcal\{J\}\_\{A\}in both flow cases\. The verifier therefore tests whether the candidate reproduces the size of the recurrent motion, not only its pointwise path, as expressed by Eq\. \([S5](https://arxiv.org/html/2608.02662#A2.E5)\)\.

#### Period or dominant frequency\.

Letqj\(m\)q\_\{j\}^\{\(m\)\}denote the characteristic temporal scale of coordinatejjin reference casemm, and letq^s,j\(m\)\\widehat\{q\}\_\{s,j\}^\{\(m\)\}denote the same quantity from the candidate rollout\. Its relative discrepancy is

mtime​\(𝐡\)\\displaystyle m\_\{\\mathrm\{time\}\}\(\\mathbf\{h\}\)=1\|ℳv\|​\|𝒥q\|​∑m∈ℳv∑j∈𝒥q\|medians⁡q^s,j\(m\)−qj\(m\)\|\|qj\(m\)\|\+ϵ,\\displaystyle=\\frac\{1\}\{\|\\mathcal\{M\}\_\{v\}\|\\,\|\\mathcal\{J\}\_\{q\}\|\}\\sum\_\{m\\in\\mathcal\{M\}\_\{v\}\}\\sum\_\{j\\in\\mathcal\{J\}\_\{q\}\}\\frac\{\\left\|\\operatorname\{median\}\_\{s\}\\widehat\{q\}\_\{s,j\}^\{\(m\)\}\-q\_\{j\}^\{\(m\)\}\\right\|\}\{\|q\_\{j\}^\{\(m\)\}\|\+\\epsilon\},\(S6\)𝗉𝖺𝗌𝗌time​\(𝐡\)\\displaystyle\\mathsf\{pass\}\_\{\\mathrm\{time\}\}\(\\mathbf\{h\}\)=𝕀​\[mtime​\(𝐡\)≤0\.10\]\.\\displaystyle=\\mathbb\{I\}\\\!\\left\[m\_\{\\mathrm\{time\}\}\(\\mathbf\{h\}\)\\leq 0\.10\\right\]\.For Van der Pol,qqis the median interval between successive peaks after the transient and𝒥q\\mathcal\{J\}\_\{q\}contains the observed oscillatory coordinate\. For each reduced flow coordinate,qjq\_\{j\}is instead the dominant nonzero Fourier frequency,qj=arg​maxν\>0⁡\|ℱ​\[ξj−ξ¯j\]​\(ν\)\|2q\_\{j\}=\\operatorname\*\{arg\\,max\}\_\{\\nu\>0\}\|\\mathcal\{F\}\[\\xi\_\{j\}\-\\overline\{\\xi\}\_\{j\}\]\(\\nu\)\|^\{2\}\. Peak spacing measures the period of the directly observed oscillator, whereas the spectral definition treats all reduced flow coordinates uniformly\. In both cases, Eq\. \([S6](https://arxiv.org/html/2608.02662#A2.E6)\) tests preservation of the dominant oscillation rate\.

#### Starting\-state consistency\.

Each candidate is integrated from three observed starting states\. For a rollout propertyu∈\{A,q\}u\\in\\\{A,q\\\}, whereAAis amplitude andqqis period or dominant frequency, define

Cu​\(𝐡\)\\displaystyle C\_\{u\}\(\\mathbf\{h\}\)=1\|ℳv\|​\|𝒥u\|​∑m∈ℳv∑j∈𝒥ustds⁡\(u^s,j\(m\)\)means⁡\|u^s,j\(m\)\|\+ϵ,\\displaystyle=\\frac\{1\}\{\|\\mathcal\{M\}\_\{v\}\|\\,\|\\mathcal\{J\}\_\{u\}\|\}\\sum\_\{m\\in\\mathcal\{M\}\_\{v\}\}\\sum\_\{j\\in\\mathcal\{J\}\_\{u\}\}\\frac\{\\operatorname\{std\}\_\{s\}\\left\(\\widehat\{u\}\_\{s,j\}^\{\(m\)\}\\right\)\}\{\\operatorname\{mean\}\_\{s\}\\left\|\\widehat\{u\}\_\{s,j\}^\{\(m\)\}\\right\|\+\\epsilon\},\(S7\)𝗉𝖺𝗌𝗌cons​\(𝐡\)\\displaystyle\\mathsf\{pass\}\_\{\\mathrm\{cons\}\}\(\\mathbf\{h\}\)=𝕀​\[CA​\(𝐡\)≤0\.10∧Cq​\(𝐡\)≤0\.05\]\.\\displaystyle=\\mathbb\{I\}\\\!\\left\[C\_\{A\}\(\\mathbf\{h\}\)\\leq 0\.10\\ \\land\\ C\_\{q\}\(\\mathbf\{h\}\)\\leq 0\.05\\right\]\.Here the hat denotes a quantity measured from an integrated candidate trajectory, and the standard deviation and mean are taken across the three starting states\. A stable recurrent model should approach comparable amplitudes and time scales when initiated at different observed phases\. Equation \([S7](https://arxiv.org/html/2608.02662#A2.E7)\) excludes candidates for which those long\-time properties depend strongly on the selected starting state\.

#### Parameter conservation and participation\.

For the cross\-Re\\mathrm\{Re\}model, lethrh\_\{r\}denote the candidate equation for the normalized Reynolds coordinate defined above and lethjh\_\{j\},j=1,…,dzj=1,\\ldots,d\_\{z\}, denote the equations for thedzd\_\{z\}latent coordinates\. The parameter check is

𝗉𝖺𝗌𝗌par\(𝐡\)=𝕀\[hr≡0∧maxs,t\|r^𝐡\(s\)\(t\)−r\(s\)\(0\)\|≤10−8∧∃j≤dz:∂hj∂r≢0\]\.\\mathsf\{pass\}\_\{\\mathrm\{par\}\}\(\\mathbf\{h\}\)=\\mathbb\{I\}\\\!\\left\[h\_\{r\}\\equiv 0\\ \\land\\ \\max\_\{s,t\}\|\\,\\widehat\{r\}\_\{\\mathbf\{h\}\}^\{\(s\)\}\(t\)\-r^\{\(s\)\}\(0\)\\,\|\\leq 10^\{\-8\}\\ \\land\\ \\exists\\,j\\leq d\_\{z\}:\\frac\{\\partial h\_\{j\}\}\{\\partial r\}\\not\\equiv 0\\right\]\.\(S8\)The first condition preserves the autonomous augmentationr˙=0\\dot\{r\}=0symbolically, and the second confirms the same conservation in numerical rollouts\. The third requires the latent dynamics to depend explicitly onrr; carrying an unused constant coordinate is not sufficient for a shared cross\-parameter model\.

Table[S2](https://arxiv.org/html/2608.02662#A2.T2)summarizes the settings associated with Eqs\. \([S3](https://arxiv.org/html/2608.02662#A2.E3)\)–\([S8](https://arxiv.org/html/2608.02662#A2.E8)\)\.

The cross\-Re\\mathrm\{Re\}shortlist is the deduplicated union of the ten candidates with the lowest multi\-trajectory rollout errors and the ten lowest\-error candidates containing explicit Reynolds\-number dependence in at least one latent\-state equation\. This two\-branch construction retains hypotheses capable of representing a shared parametric law when rollout ranking favors a parameter\-independent approximation over the discovery trajectories\. Overlap between the branches produced the 17 distinct candidates reported in Table[S2](https://arxiv.org/html/2608.02662#A2.T2)\.

For the cross\-Re\\mathrm\{Re\}integration, one model\-time unit corresponds to0\.10\.1physical\-time units\. The numerical horizon and discard interval of 300 and 50 model\-time units therefore correspond to the reported values of 30 and 5\.

The numerical thresholds in Table[S2](https://arxiv.org/html/2608.02662#A2.T2)were held constant within each experimental setting\. Their sensitivity is examined for the representative fixed\-Re\\mathrm\{Re\}case in Section[S2\.7](https://arxiv.org/html/2608.02662#A2.SS7)\.

Table S2:Verifier settings used in the three experimental settings\.∗The two cross\-Re\\mathrm\{Re\}shortlist branches yielded 17 distinct candidates\.

∗∗Cross\-Re\\mathrm\{Re\}horizon and discard values are reported in physical\-time units\.

### S2\.3Identifiability interpretation

Identifiability concerns whether the observations and modeling assumptions distinguish one structure or coefficient vector from plausible alternatives\. This issue is central to trajectory\-based symbolic discovery because different equations can reproduce a finite observed path while defining different vector fields away from it, different long\-time attractors, or different responses to a changed initial condition or parameter\. A low rollout error on the supplied trajectories is therefore not, by itself, evidence that the dynamics have been uniquely identified\.

The claim made in this work is intentionally restricted\. The verifier suite does not establish global structural identifiability over an unrestricted class of differential equations\. It improves discrimination within the finite hypothesis class generated by the symbolic backbone and over the specified state domains, initial conditions, parameter values, and verification horizons\. To express this distinction, letℋG\\mathcal\{H\}\_\{G\}be the generated, deduplicated hypothesis class,R​\(𝐟\)R\(\\mathbf\{f\}\)its multi\-trajectory rollout discrepancy, andε\\varepsilona chosen rollout tolerance\. The set of candidates that are indistinguishable at this rollout resolution is

ℰR​\(ε\)=\{𝐟∈ℋG:R​\(𝐟\)≤ε\}\.\\mathcal\{E\}\_\{R\}\(\\varepsilon\)=\\left\\\{\\mathbf\{f\}\\in\\mathcal\{H\}\_\{G\}:R\(\\mathbf\{f\}\)\\leq\\varepsilon\\right\\\}\.\(S9\)Equation \([S9](https://arxiv.org/html/2608.02662#A2.E9)\) describes equivalence only with respect to the sampled rollout discrepancy\. LetℋA⊆ℋG\\mathcal\{H\}\_\{A\}\\subseteq\\mathcal\{H\}\_\{G\}denote the candidates satisfying all applicable dynamical and physical\-admissibility conditions\. Verification then restricts the rollout\-consistent set to

ℰV​G​\(ε\)=ℰR​\(ε\)∩ℋA\.\\mathcal\{E\}\_\{VG\}\(\\varepsilon\)=\\mathcal\{E\}\_\{R\}\(\\varepsilon\)\\cap\\mathcal\{H\}\_\{A\}\.\(S10\)The intersection in Eq\. \([S10](https://arxiv.org/html/2608.02662#A2.E10)\) improves practical identifiability when it removes at least one rollout\-consistent alternative while retaining an admissible candidate\. The additional discrimination comes from properties not represented by a scalar rollout ordering alone, including local vector\-field agreement, long\-horizon boundedness, recurrent behavior across starting states, and parameter conservation and participation\.

### S2\.4Local sensitivity\-rank interpretation

The preceding set\-based view compares distinct candidates in a finite hypothesis class\. A complementary local view asks whether the rollout and verification diagnostics can distinguish small coefficient changes within one symbolic structure\. Let𝜽∈ℝnθ\\boldsymbol\{\\theta\}\\in\\mathbb\{R\}^\{n\_\{\\theta\}\}contain itsnθn\_\{\\theta\}free coefficients,𝐲​\(𝜽\)\\mathbf\{y\}\(\\boldsymbol\{\\theta\}\)collect the sampled candidate rollouts, and𝐦​\(𝜽\)\\mathbf\{m\}\(\\boldsymbol\{\\theta\}\)collect the real\-valued verifier metrics before thresholding\. The mapping from coefficients to the combined vector of rollout samples and verifier metrics has the local sensitivity matrix

𝐉V​G​\(𝜽\)=\[∂𝐲/∂𝜽∂𝐦/∂𝜽\]\.\\mathbf\{J\}\_\{VG\}\(\\boldsymbol\{\\theta\}\)=\\begin\{bmatrix\}\\partial\\mathbf\{y\}/\\partial\\boldsymbol\{\\theta\}\\\\ \\partial\\mathbf\{m\}/\\partial\\boldsymbol\{\\theta\}\\end\{bmatrix\}\.\(S11\)Each column of𝐉V​G\\mathbf\{J\}\_\{VG\}in Eq\. \([S11](https://arxiv.org/html/2608.02662#A2.E11)\) describes how the measured quantities change under a perturbation of one coefficient\. Appending verifier sensitivities cannot reduce the rank available from the rollout samples:

rank⁡𝐉V​G≥rank⁡\(∂𝐲/∂𝜽\)\.\\operatorname\{rank\}\\mathbf\{J\}\_\{VG\}\\geq\\operatorname\{rank\}\\left\(\\partial\\mathbf\{y\}/\\partial\\boldsymbol\{\\theta\}\\right\)\.\(S12\)A strict inequality in Eq\. \([S12](https://arxiv.org/html/2608.02662#A2.E12)\) means that at least one coefficient direction unresolved by the sampled rollouts affects an additional diagnostic\. Full column rank is a sufficient local condition for coefficient identifiability within the fixed structure, subject to numerical conditioning and symbolic symmetries\. This rank analysis provides an interpretation of the information contributed by verification; it was not an additional computational step in the experiments\. The binary decisions themselves are not differentiated, and metrics based on peak locations, percentiles, or discrete spectral maxima are interpreted locally only while the selected features remain unchanged\.

### S2\.5Implications for joint coordinate and equation discovery

In the present workflow, POD or autoencoder coordinates are selected before symbolic discovery, so identifiability is assessed conditionally on that representation\. A future joint formulation could instead allow verification information to influence the coordinates themselves\. Letℰϕ\\mathcal\{E\}\_\{\\boldsymbol\{\\phi\}\}be an encoder with parametersϕ\\boldsymbol\{\\phi\},𝐳ϕ=ℰϕ​\(𝝎\)\\mathbf\{z\}\_\{\\boldsymbol\{\\phi\}\}=\\mathcal\{E\}\_\{\\boldsymbol\{\\phi\}\}\(\\boldsymbol\{\\omega\}\)its latent state, and𝐠𝜽\\mathbf\{g\}\_\{\\boldsymbol\{\\theta\}\}a symbolic latent vector field with coefficients𝜽\\boldsymbol\{\\theta\}\. The sensitivity of reconstruction, latent rollout, and verifier quantities to both parameter sets could be represented by

𝐉joint=∂∂\(ϕ,𝜽\)​\[ℛϕ​\(𝐳ϕ\)𝐳^ϕ,𝜽​\(t\)𝐦​\(ϕ,𝜽\)\],\\mathbf\{J\}\_\{\\mathrm\{joint\}\}=\\frac\{\\partial\}\{\\partial\(\\boldsymbol\{\\phi\},\\boldsymbol\{\\theta\}\)\}\\begin\{bmatrix\}\\mathcal\{R\}\_\{\\boldsymbol\{\\phi\}\}\(\\mathbf\{z\}\_\{\\boldsymbol\{\\phi\}\}\)\\\\ \\widehat\{\\mathbf\{z\}\}\_\{\\boldsymbol\{\\phi\},\\boldsymbol\{\\theta\}\}\(t\)\\\\ \\mathbf\{m\}\(\\boldsymbol\{\\phi\},\\boldsymbol\{\\theta\}\)\\end\{bmatrix\},\(S13\)whereℛϕ\\mathcal\{R\}\_\{\\boldsymbol\{\\phi\}\}is the reconstruction map,𝐳^ϕ,𝜽​\(t\)\\widehat\{\\mathbf\{z\}\}\_\{\\boldsymbol\{\\phi\},\\boldsymbol\{\\theta\}\}\(t\)is the integrated latent trajectory, and𝐦\\mathbf\{m\}contains differentiable surrogates of the verifier properties\. Equation \([S13](https://arxiv.org/html/2608.02662#A2.E13)\) shows how additional dynamical information could distinguish embeddings that reconstruct the field similarly but differ in approximate closure, boundedness, recurrence, or parameter dependence\. Such a formulation could improve the identifiability of latent coordinates compatible with symbolic discovery\. It remains a prospective extension; the autoencoder used here was optimized without verifier feedback\.

### S2\.6Parameter dependence in the cross\-Re\\mathrm\{Re\}model

In the cross\-Re\\mathrm\{Re\}experiment, each decoding is performed for one trajectory whose appended normalized Reynolds coordinaterris constant\. A single trajectory therefore cannot identify how the latent dynamics vary with Reynolds number\. For example, on a trajectory withr=r0≠0r=r\_\{0\}\\neq 0, the terms

α​zj=αr0​r​zjwhen​r=r0\\alpha z\_\{j\}=\\frac\{\\alpha\}\{r\_\{0\}\}\\,rz\_\{j\}\\qquad\\text\{when \}r=r\_\{0\}\(S14\)are observationally indistinguishable\. The occurrence ofrrin an individually decoded expression is therefore a candidate structural hypothesis, rather than evidence that one decoding has resolved the cross\-trajectory parameter dependence\. The ambiguity in Eq\. \([S14](https://arxiv.org/html/2608.02662#A2.E14)\) is resolved only by evaluating a shared candidate across different values ofrr\.

Such terms can nevertheless arise becauserris supplied as a state variable and the symbolic pretraining distribution includes cross\-variable products and related dependencies\. Variations in latent amplitude, frequency, and trajectory geometry at different Reynolds numbers can thus be represented by expressions involving the constant parameter channel\. The cross\-trajectory evidence enters after decoding: candidates are pooled and evaluated as shared equations over all seven discovery Reynolds numbers\. The parameter verifier in Eq\. \([S8](https://arxiv.org/html/2608.02662#A2.E8)\) additionally requires conservation ofrrand its participation in at least one latent equation\. Parameter dependence in the selected model consequently emerges from candidate generation followed by cross\-Re\\mathrm\{Re\}pooling, evaluation, and verification, rather than from simultaneous observation of all trajectories during an individual decode\.

### S2\.7Sensitivity to verifier tolerances

The numerical tolerances in Table[S2](https://arxiv.org/html/2608.02662#A2.T2)determine the boundary of the admissible set and therefore warrant a sensitivity check\. We use the fixed\-Re=300\\mathrm\{Re\}=300cylinder case as a representative example because it combines the local, boundedness, amplitude, frequency, and starting\-state verifiers without the additional parameter constraints of the cross\-Re\\mathrm\{Re\}model\. The same perturbation procedure can be applied to the Van der Pol and cross\-Re\\mathrm\{Re\}cases, although stability in this representative case does not imply invariance of their admissible sets\.

The sensitivity matrix crosses POD ranks four and six with candidate generation from either 8 evenly distributed development windows or all 24 development windows\. The two even POD ranks retain complete modal pairs while testing a more compact and a more energetic representation; the two window counts compare selective and exhaustive candidate generation from the same window bank\. Candidate ranking, verification, and coefficient refinement use all 24 windows in every configuration\. Section[S4\.2](https://arxiv.org/html/2608.02662#A4.SS2)provides the full representation and selection comparison\.

Letγ∈\{0\.75,1\.0,1\.5\}\\gamma\\in\\\{0\.75,1\.0,1\.5\\\}jointly scale the local vector\-field agreement factor in Eq\. \([S3](https://arxiv.org/html/2608.02662#A2.E3)\) and the amplitude, frequency, amplitude\-CV, and frequency\-CV tolerances\. Thusγ<1\\gamma<1imposes stricter quantitative agreement,γ=1\\gamma=1recovers the nominal settings in Table[S2](https://arxiv.org/html/2608.02662#A2.T2), andγ\>1\\gamma\>1relaxes them\. Finite integration and boundedness remain at their nominal settings because they serve a different role: finite integration is a Boolean requirement for a valid rollout, while the boundedness limit defines the admissible state envelope relative to the observed dynamics\. Changing either would alter the minimum dynamical\-admissibility requirement rather than isolate sensitivity to the agreement tolerances\.

The number of admissible candidates, ordered byγ=\(0\.75,1\.0,1\.5\)\\gamma=\(0\.75,1\.0,1\.5\), was:

- •POD rank four, 8 decoding windows:\(10,10,10\)\(10,10,10\);
- •POD rank four, 24 decoding windows:\(10,10,10\)\(10,10,10\);
- •POD rank six, 8 decoding windows:\(5,6,9\)\(5,6,9\); and
- •POD rank six, 24 decoding windows:\(9,10,10\)\(9,10,10\)\.

The admissible set therefore expands under looser tolerances for the rank\-six cases, whereas every shortlisted rank\-four candidate satisfies even the stricter setting\. In all four configurations, however, the lowest\-rollout admissible equation selected at the nominal tolerance remained the selected equation throughout the sensitivity range\. The nominally selected pre\-optimization equation and the next\-lowest\-rollout admissible alternative for the selected rank\-six, eight\-window configuration are reported in Section[S4\.2](https://arxiv.org/html/2608.02662#A4.SS2), together with the complete selection matrix\.

These results show that the verifier thresholds influence which secondary candidates enter the admissible set, as expected for any thresholded selection rule, but do not materially change the final outcome in this representative fixed\-Re\\mathrm\{Re\}case\. The nominal values provide a reasonable balance: they reject candidates with substantial local or oscillatory disagreement without requiring near\-exact agreement from an effective reduced\-order model\. Tightening or relaxing them by 25% and 50%, respectively, preserves the selected equation in every configuration\.

## Appendix S3Van der Pol Test\-Case Details and Extended Results

Section[S1](https://arxiv.org/html/2608.02662#A1)used the Van der Pol oscillator to select the beam size and sampling temperature inherited from the symbolic backbone\. The present section documents the final experiment conducted after those settings had been fixed\. It provides the exact trajectory sets and their evaluation roles, the selected pre\- and post\-optimization equations, their verifier outcomes, and representative final\-test rollouts\. The candidate\-generation \(decoding\) comparison is not repeated here: the results below use only the selected beam size of 20 and temperature of 0\.1\.

### S3\.1Trajectory Sets and Evaluation Roles

For bothμ=0\.5\\mu=0\.5and1\.51\.5, trajectories were integrated over0≤t≤200\\leq t\\leq 20at 301 equally spaced times\. The same initial\-condition sets were used at both parameter values\. The complete fixed\-seed sets and their distinct roles are shown below: the model\-construction bank is on the left, the validation set used in the beam–temperature comparison is in the center, and the final\-test set is on the right\.

Model\-construction bank

Early window:0≤t≤100\\leq t\\leq 10\(151 samples\)\. Late window:10≤t≤2010\\leq t\\leq 20\(151 samples\)\.

Beam–temperature validation

Final test

As shown in the left\-most table, the 12 model\-construction trajectories each contributed one early and one late window, giving a 24\-window bank\. The eight windows used as inputs for candidate generation were fixed before the beam–temperature comparison: trajectory 1 contributed its late window; trajectory 2 contributed both windows; trajectories 5 and 6 contributed their late windows; trajectory 8 contributed both windows; and trajectory 9 contributed its early window\. The saved bank orders the windows by trajectory, with each trajectory’s early window immediately followed by its late window\. Under this storage convention, the eight candidate\-generation inputs have global indices\(3,4,5,11,13,16,17,18\)\(3,4,5,11,13,16,17,18\)\. Candidate pooling, rollout ranking, verification, and coefficient optimization then used evidence from all 24 model\-construction windows, rather than only the eight inputs from which candidates were generated\.

The long\-horizon verifier rollouts began from model\-construction initial conditions 0, 4, and 8 in the left\-most table\. This reuse was intentional: the verifiers assess dynamical admissibility from multiple locations in the model\-construction domain and do not provide a held\-out performance estimate\. The eight initial conditions in the center table were used only for the beam–temperature comparison in Section[S1](https://arxiv.org/html/2608.02662#A1)\. The eight initial conditions in the right\-most table were reserved for final testing after the candidate\-generation settings and equations had been fixed\. These two evaluation sets are disjoint from the model\-construction bank and from one another\.

### S3\.2Selected equations and verifier outcomes

The canonical system hasx˙0=x1\\dot\{x\}\_\{0\}=x\_\{1\}andx˙1=μ​x1−x0−μ​x02​x1\\dot\{x\}\_\{1\}=\\mu x\_\{1\}\-x\_\{0\}\-\\mu x\_\{0\}^\{2\}x\_\{1\}\. The selected VG equations are written below in the factorized form returned by the symbolic workflow\. Coefficients are rounded only for presentation; verification and rollout evaluation used the full\-precision values\.

𝝁=0\.5\\boldsymbol\{\\mu=0\.5\}\.

Pre\-optimization:x˙0\\displaystyle\\text\{Pre\-optimization:\}\\qquad\\dot\{x\}\_\{0\}=1\.0239​x1−0\.0548​x0,\\displaystyle=0239x\_\{1\}\-0548x\_\{0\},\(S15\)x˙1\\displaystyle\\dot\{x\}\_\{1\}=0\.2763​x1−1\.0463​x0−0\.0769​x1​\(0\.0736​x0\+3\.6009​x02\)\.\\displaystyle=2763x\_\{1\}\-0463x\_\{0\}\-0769x\_\{1\}\\\!\\left\(0\.0736x\_\{0\}\+3\.6009x\_\{0\}^\{2\}\\right\)\.Post\-optimization:x˙0\\displaystyle\\text\{Post\-optimization:\}\\qquad\\dot\{x\}\_\{0\}=0\.99726​x1−0\.012400​x0,\\displaystyle=99726x\_\{1\}\-012400x\_\{0\},\(S16\)x˙1\\displaystyle\\dot\{x\}\_\{1\}=0\.54166​x1−0\.99708​x0−0\.13251​x1​\(0\.029451​x0\+4\.0969​x02\)\.\\displaystyle=54166x\_\{1\}\-99708x\_\{0\}\-13251x\_\{1\}\\\!\\left\(0\.029451x\_\{0\}\+4\.0969x\_\{0\}^\{2\}\\right\)\.
𝝁=1\.5\\boldsymbol\{\\mu=1\.5\}\.

Pre\-optimization:x˙0\\displaystyle\\text\{Pre\-optimization:\}\\qquad\\dot\{x\}\_\{0\}=0\.0007​\(1\+1\.2715​x0\)2\+0\.9889​x1,\\displaystyle=0007\\left\(1\+1\.2715x\_\{0\}\\right\)^\{2\}\+9889x\_\{1\},\(S17\)x˙1\\displaystyle\\dot\{x\}\_\{1\}=1\.4616​x1−0\.9972​x0−0\.1104​x1​\(0\.0690​x0\+12\.5443​x02\)\.\\displaystyle=4616x\_\{1\}\-9972x\_\{0\}\-1104x\_\{1\}\\\!\\left\(0\.0690x\_\{0\}\+12\.5443x\_\{0\}^\{2\}\\right\)\.Post\-optimization:x˙0\\displaystyle\\text\{Post\-optimization:\}\\qquad\\dot\{x\}\_\{0\}=0\.00067573​\(0\.99818\+1\.2726​x0\)2\+0\.99637​x1,\\displaystyle=00067573\\left\(0\.99818\+1\.2726x\_\{0\}\\right\)^\{2\}\+99637x\_\{1\},\(S18\)x˙1\\displaystyle\\dot\{x\}\_\{1\}=1\.43998​x1−1\.00151​x0−0\.11351​x1​\(0\.070367​x0\+12\.8741​x02\)\.\\displaystyle=43998x\_\{1\}\-00151x\_\{0\}\-11351x\_\{1\}\\\!\\left\(0\.070367x\_\{0\}\+12\.8741x\_\{0\}^\{2\}\\right\)\.
Comparison of each pre\- and post\-optimization pair shows that coefficient optimization preserves the generated symbolic structure and changes only its numerical constants; it does not add or remove terms\. For each value ofμ\\mu, only one of the ten rollout\-ranked candidates passed the complete pre\-optimization verifier suite\. That candidate also had the lowest aggregate rollout error\. The verifier suite therefore provided an independent admissibility requirement without changing the selected candidate in these two cases\.

The selected equations were verified again after coefficient optimization\. Table[S3](https://arxiv.org/html/2608.02662#A3.T3)places each measured outcome beside the corresponding Van der Pol acceptance requirement\. The local vector\-field limit is shortlist\-relative and therefore differs between the two parameter values; the remaining limits are the fixed settings defined in Section[S2](https://arxiv.org/html/2608.02662#A2)\.

Table S3:Post\-optimization verifier outcomes for the selected Van der Pol equations\. Percentages are reported relative to the corresponding ground\-truth oscillation statistic\.The local vector\-field scores are 2\.17% and 4\.59% of their respective limits forμ=0\.5\\mu=0\.5and1\.51\.5\. Across both systems, the amplitude errors remain below 0\.75%, the period errors below 0\.27%, and the cross\-start coefficients of variation below 0\.38%\. The post\-optimization equations thus satisfy the numerical\-integration, boundedness, local\-dynamics, oscillation, and attractor\-consistency requirements from all three verifier initial conditions\. This reverification establishes dynamical admissibility after the constants have changed; Section[S3\.3](https://arxiv.org/html/2608.02662#A3.SS3)separately assesses predictive generalization from initial conditions excluded from model construction\.

### S3\.3Representative held\-out rollouts

Across the eight final\-test initial conditions, coefficient optimization reduced the median VG rollout NMSE from 0\.287 to6\.99×10−46\.99\\times 10^\{\-4\}forμ=0\.5\\mu=0\.5, and from 0\.0139 to3\.27×10−33\.27\\times 10^\{\-3\}forμ=1\.5\\mu=1\.5\. These values differ from the validation errors in Table[S1](https://arxiv.org/html/2608.02662#A1.T1): the validation errors selected the candidate\-generation settings, whereas the present values were obtained from the previously unused final\-test trajectories after those settings and the equations had been fixed\.

Figure[S1](https://arxiv.org/html/2608.02662#A3.F1)shows one representative final\-test case for each parameter\. The representative\-case rule was defined over the eight final\-test trajectories before inspecting the plots\. An eligible trajectory required complete rollouts from the original ODEFormer equation before coefficient optimization and from the VG equation both before and after optimization\. Among these eligible trajectories, we selected the case whose post\-optimization VG NMSE was closest to the median over all eight final\-test trajectories\. When the two observations bracketing the median were equally close, the lower trajectory index was used\. This rule gives test trajectory 2, with𝐱​\(0\)=\(−0\.811743,−0\.133746\)\\mathbf\{x\}\(0\)=\(\-0\.811743,\-0\.133746\), forμ=0\.5\\mu=0\.5, and test trajectory 3, with𝐱​\(0\)=\(−0\.041897,−0\.680522\)\\mathbf\{x\}\(0\)=\(\-0\.041897,\-0\.680522\), forμ=1\.5\\mu=1\.5\. The post\-optimization original ODEFormer rollout is included in the figure to complete the before–after comparison but did not influence representative case selection\.

![Refer to caption](https://arxiv.org/html/2608.02662v1/x5.png)Figure S1:Representative final\-test Van der Pol rollouts\. Ground\-truth trajectories are compared with equations from the original single\-trajectory ODEFormer workflow and the VG workflow, in each case before and after coefficient optimization\.Forμ=0\.5\\mu=0\.5, both original ODEFormer equations begin to separate markedly from the ground truth after approximatelyt=10t=10\. Their oscillation amplitudes then grow rapidly, and neither preserves the bounded ground\-truth cycle over the full displayed interval\. The pre\-optimization VG equation remains bounded and oscillatory but accumulates a visible phase discrepancy\. Coefficient optimization brings the VG trajectory and phase portrait into close agreement with the ground truth\. Forμ=1\.5\\mu=1\.5, all four learned rollouts remain oscillatory\. The original ODEFormer equations before and after optimization trace visibly different cycles, whereas the post\-optimization VG rollout follows the ground\-truth cycle throughout the interval; its smaller residual phase and shape differences are consistent with the final\-test NMSE\.

The verifiers therefore contribute more than a check on numerical solvability: before selection, they require finite and bounded integration, local vector\-field agreement, accurate oscillation amplitude and period, and consistent attractor statistics across several initial conditions\. The final\-test rollouts then provide separate evidence that the selected, reverified equations preserve the relevant oscillatory dynamics from initial conditions excluded from model construction\. Together with the aggregate comparison in Figure[2](https://arxiv.org/html/2608.02662#S5.F2), these results show that the VG advantage is supported by explicit multi\-start dynamical requirements as well as by held\-out predictive accuracy, rather than by low average rollout error alone\.

## Appendix S4Fixed\-Reynolds\-Number Representation and Selection Details

This section expands the fixed\-Re=300\\mathrm\{Re\}=300experiment reported in Section[5\.2](https://arxiv.org/html/2608.02662#S5.SS2)\. After the initial 500 snapshots were discarded, the remaining time series was divided chronologically into a 250\-snapshot development interval, a 100\-snapshot validation interval, and a 150\-snapshot final\-test interval\. The development interval supplied the POD representation and supported candidate generation, multi\-window ranking and verification, and coefficient optimization\. The validation interval was used to choose the POD rank and decoding\-window configuration; the final\-test interval remained held out until those choices were complete\. Sections[S4\.1](https://arxiv.org/html/2608.02662#A4.SS1),[S4\.2](https://arxiv.org/html/2608.02662#A4.SS2), and[S4\.3](https://arxiv.org/html/2608.02662#A4.SS3)respectively examine the representation, configuration selection, and final\-test behavior\.

### S4\.1POD spectrum and mode\-pair diagnostics

POD was applied to fluctuations about the development\-interval mean, rather than to the total vorticity fields\. If𝝎​\(ti\)\\boldsymbol\{\\omega\}\(t\_\{i\}\)is the vectorized field at development timetit\_\{i\}and𝝎¯D\\overline\{\\boldsymbol\{\\omega\}\}\_\{\\\!D\}is their temporal mean, the fluctuation snapshot matrix and its singular\-value decomposition are

𝐗′\\displaystyle\\mathbf\{X\}^\{\\prime\}=\[𝝎​\(t1\)−𝝎¯D⋯𝝎​\(tND\)−𝝎¯D\]=𝐔​𝚺​𝐕𝖳,\\displaystyle=\\begin\{bmatrix\}\\boldsymbol\{\\omega\}\(t\_\{1\}\)\-\\overline\{\\boldsymbol\{\\omega\}\}\_\{\\\!D\}&\\cdots&\\boldsymbol\{\\omega\}\(t\_\{N\_\{D\}\}\)\-\\overline\{\\boldsymbol\{\\omega\}\}\_\{\\\!D\}\\end\{bmatrix\}=\\mathbf\{U\}\\boldsymbol\{\\Sigma\}\\mathbf\{V\}^\{\\mathsf\{T\}\},\(S19\)Er\\displaystyle E\_\{r\}=∑j=1rσj2∑j=1Nsσj2\.\\displaystyle=\\frac\{\\sum\_\{j=1\}^\{r\}\\sigma\_\{j\}^\{2\}\}\{\\sum\_\{j=1\}^\{N\_\{s\}\}\\sigma\_\{j\}^\{2\}\}\.HereND=250N\_\{D\}=250,σj\\sigma\_\{j\}is thejjth singular value,NsN\_\{s\}is the number of nonzero singular values, andErE\_\{r\}is the fraction of fluctuation energy represented by the firstrrmodes\. The mean field is restored only when the reduced coordinates are reconstructed in physical space\.

The near\-equal singular values in Figure[S2](https://arxiv.org/html/2608.02662#A4.F2)form three leading pairs, as expected for quadrature components of an oscillatory wake\. We therefore compared two even ranks that preserve complete pairs: rank four gives a compact representation of the first two pairs, whereas rank six adds the third pair and retains a larger fraction of the wake fluctuations\. This comparison tests whether the additional oscillatory content improves field reconstruction sufficiently to justify a higher\-dimensional symbolic system\. Ranks four and six retain 85\.61% and 94\.52% of the fluctuation energy, respectively\.

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

Figure S2:POD singular\-value spectrum and cumulative fluctuation energy at fixedRe=300\\mathrm\{Re\}=300\.For a consecutive pair\(2​k−1,2​k\)\(2k\-1,2k\), we quantify its relative energy imbalance by

δk=\|σ2​k−12−σ2​k2\|σ2​k−12\+σ2​k2\.\\delta\_\{k\}=\\frac\{\|\\sigma\_\{2k\-1\}^\{2\}\-\\sigma\_\{2k\}^\{2\}\|\}\{\\sigma\_\{2k\-1\}^\{2\}\+\\sigma\_\{2k\}^\{2\}\}\.\(S20\)The pair\-energy column in Table[S4](https://arxiv.org/html/2608.02662#A4.T4)is the percentage of the total fluctuation energy carried by the two modes in that row; the cumulative column gives the energy retained through that pair\. The first pair therefore accounts for 73\.44% of the wake fluctuations\. The second adds 12\.17 percentage points, bringing the rank\-four representation to 85\.61%, and the third adds a further 8\.91 percentage points, bringing rank six to 94\.52%\. Their coefficient trajectories vary on successively faster time scales, corresponding to the fundamental shedding oscillation and its higher harmonics\. The imbalanceδk\\delta\_\{k\}equals zero for an exactly energy\-matched pair and increases as one member becomes dominant\. Its small and decreasing values show that each retained pair uses both quadrature components, most closely for modes 5–6\. Table[S4](https://arxiv.org/html/2608.02662#A4.T4)therefore complements the cumulative spectrum by separating the energy contributed by each oscillatory pair\.

Table S4:POD mode\-pair energy and imbalance diagnostics\.
### S4\.2POD\-rank and decoding\-window selection matrix

The representation and decoding choices motivated in Section[S4\.1](https://arxiv.org/html/2608.02662#A4.SS1)were compared without using the final\-test interval\. For each rank, the VG workflow generated candidates independently from either eight evenly distributed windows or all 24 windows in the development interval\. The eight\-window setting tests whether broad temporal coverage is sufficient without decoding every available window; the 24\-window setting instead maximizes the diversity of generated candidates\. The resulting pools contained 160 and 480 equation systems, respectively\. In every configuration, candidate ranking, verification, and coefficient optimization used all 24 windows\. Only the number of windows supplied independently to the symbolic backbone during candidate generation was varied\.

Table[S5](https://arxiv.org/html/2608.02662#A4.T5)reports the validation comparison\. The rank\-four models have the smallest coordinate\-space errors, but their POD truncation error exceeds 20%, and their end\-to\-end field errors consequently remain above 21%\. Rank six reduces the representation error to 12\.70%\. Within this rank, decoding from eight windows gives the smaller continuous rollout and end\-to\-end field errors\. It was therefore chosen before the final\-test interval was evaluated\.

Table S5:Fixed\-Re\\mathrm\{Re\}validation matrix\. “Window” and “Continuous” are median coordinate errors; POD, VG–POD, and end\-to\-end are mean relative field errors\. All errors are percentages, and “Adm\.” gives the number of verifier\-admissible candidates in the ten\-candidate shortlist\.Table[S5](https://arxiv.org/html/2608.02662#A4.T5)reports the admissible counts obtained at the nominal verifier settings\. Section[S2\.7](https://arxiv.org/html/2608.02662#A2.SS7)examines how these counts change when the quantitative agreement tolerances are tightened or relaxed\. For the selected rank\-six, eight\-window configuration, Eq\. \([S21](https://arxiv.org/html/2608.02662#A4.E21)\) below is the pre\-optimization equation selected at the nominal settings, and Eq\. \([S22](https://arxiv.org/html/2608.02662#A4.E22)\) is the next\-lowest\-rollout admissible alternative\. These are the two equations referenced in Section[S2\.7](https://arxiv.org/html/2608.02662#A2.SS7)\. Although the number of admissible rank\-six candidates changes across the tested tolerance factors, Eq\. \([S21](https://arxiv.org/html/2608.02662#A4.E21)\) remains selected in every case\.

Let𝐳=\(z1,…,z6\)𝖳\\mathbf\{z\}=\(z\_\{1\},\\ldots,z\_\{6\}\)^\{\\mathsf\{T\}\}denote the standardized POD coordinates, consistent with the notation used for the post\-optimization system in Eq\. \([12](https://arxiv.org/html/2608.02662#S5.E12)\)\. The pre\-optimization system with the lowest development\- window rollout error, whose structure was advanced to coefficient optimization, was

z˙1\\displaystyle\\dot\{z\}\_\{1\}=1\.0044​z2−0\.0700​\(12\.3400−0\.9594​z2\)−1,\\displaystyle=0044z\_\{2\}\-0700\(23400\-9594z\_\{2\}\)^\{\-1\},z˙2\\displaystyle\\dot\{z\}\_\{2\}=−1\.1105​z1,\\displaystyle=\-1105z\_\{1\},\(S21\)z˙3\\displaystyle\\dot\{z\}\_\{3\}=−2\.1602​z4,\\displaystyle=\-1602z\_\{4\},z˙4\\displaystyle\\dot\{z\}\_\{4\}=2\.1867​z3,\\displaystyle=1867z\_\{3\},z˙5\\displaystyle\\dot\{z\}\_\{5\}=3\.4432​z6−0\.1382​sin⁡\(0\.1052\+9\.1998​z2\),\\displaystyle=4432z\_\{6\}\-1382\\sin\(1052\+1998z\_\{2\}\),z˙6\\displaystyle\\dot\{z\}\_\{6\}=−0\.1399−2\.8694​z5\.\\displaystyle=\-1399\-8694z\_\{5\}\.The next\-lowest\-rollout admissible system provides a representative alternative:

z˙1\\displaystyle\\dot\{z\}\_\{1\}=1\.1208​z2,\\displaystyle=1208z\_\{2\},z˙2\\displaystyle\\dot\{z\}\_\{2\}=−1\.0604​z1,\\displaystyle=\-0604z\_\{1\},\(S22\)z˙3\\displaystyle\\dot\{z\}\_\{3\}=−2\.2995​z4,\\displaystyle=\-2995z\_\{4\},z˙4\\displaystyle\\dot\{z\}\_\{4\}=2\.0672​z3,\\displaystyle=0672z\_\{3\},z˙5\\displaystyle\\dot\{z\}\_\{5\}=3\.1395​z6−0\.2202​sin⁡\(11\.8700\+148\.5203​z6\),\\displaystyle=1395z\_\{6\}\-2202\\sin\(18700\+485203z\_\{6\}\),z˙6\\displaystyle\\dot\{z\}\_\{6\}=−0\.1463−3\.1228​z5\.\\displaystyle=\-1463\-1228z\_\{5\}\.Both systems preserve three oscillator pairs and satisfy all applicable verifiers, but their nonlinear structures differ\. In Eq\. \([S21](https://arxiv.org/html/2608.02662#A4.E21)\),−0\.0700​\(12\.3400−0\.9594​z2\)−1\-0\.0700\(12\.3400\-0\.9594z\_\{2\}\)^\{\-1\}adds a reciprocal correction toz˙1\\dot\{z\}\_\{1\}in the first pair, while−0\.1382​sin⁡\(0\.1052\+9\.1998​z2\)\-0\.1382\\sin\(0\.1052\+9\.1998z\_\{2\}\)couples the first pair toz˙5\\dot\{z\}\_\{5\}in the third pair\. Equation \([S22](https://arxiv.org/html/2608.02662#A4.E22)\) leaves the first pair linear and instead introduces−0\.2202​sin⁡\(11\.8700\+148\.5203​z6\)\-0\.2202\\sin\(11\.8700\+148\.5203z\_\{6\}\)as a nonlinear correction within the third pair\. Both systems also contain a small constant offset inz˙6\\dot\{z\}\_\{6\}\.

At the nominal settings, six of the ten shortlisted systems for this configuration satisfy every verifier\. The coexistence of Eqs\. \([S21](https://arxiv.org/html/2608.02662#A4.E21)\) and \([S22](https://arxiv.org/html/2608.02662#A4.E22)\) therefore illustrates the practical non\-uniqueness discussed in Section[S2\.3](https://arxiv.org/html/2608.02662#A2.SS3): verification excludes four rollout\-competitive candidates but does not imply a unique symbolic structure within the admissible set\. Rollout error provides the remaining selection criterion and chooses Eq\. \([S21](https://arxiv.org/html/2608.02662#A4.E21)\)\. This supports improved discrimination within the generated hypothesis class, rather than global structural identifiability\. Coefficient optimization of the selected equation produced the post\-optimization system reported in Eq\. \([12](https://arxiv.org/html/2608.02662#S5.E12)\)\.

### S4\.3Extended coordinate and field comparisons

The final\-test interval was evaluated only after the rank\-six, eight\-window configuration had been chosen\. In Figures[S3](https://arxiv.org/html/2608.02662#A4.F3)and[S4](https://arxiv.org/html/2608.02662#A4.F4),aja\_\{j\}denotes the plotted standardized POD coefficient and corresponds tozjz\_\{j\}in Eqs\. \([S21](https://arxiv.org/html/2608.02662#A4.E21)\) and \([S22](https://arxiv.org/html/2608.02662#A4.E22)\)\. The first oscillator pair remains closely aligned throughout the interval\. The second pair develops the most visible phase difference, while the third retains the higher\-frequency oscillation with moderate amplitude and phase differences\. The phase portraits in Figure[S4](https://arxiv.org/html/2608.02662#A4.F4)\(a\)–\(c\) show that all three predicted pairs remain bounded and continue to trace closed oscillatory paths\. The first predicted orbit remains close to its POD counterpart; the higher\-harmonic pairs show greater differences in orbit geometry\. Because these pairs carry substantially less fluctuation energy than the first pair, their coordinate\-level differences have a smaller influence on the reconstructed vorticity field\.

![Refer to caption](https://arxiv.org/html/2608.02662v1/x7.png)Figure S3:Final\-test fixed\-Re\\mathrm\{Re\}rollouts in the six standardized POD coordinates\.Field error was decomposed into three contributions\. For the simulated vorticity field𝝎​\(t\)\\boldsymbol\{\\omega\}\(t\), its rank\-six POD projection𝝎POD​\(t\)\\boldsymbol\{\\omega\}\_\{\\mathrm\{POD\}\}\(t\), and the field reconstructed from the symbolic rollout𝝎^VG​\(t\)\\widehat\{\\boldsymbol\{\\omega\}\}\_\{\\mathrm\{VG\}\}\(t\), the instantaneous relative errors are

ePOD​\(t\)\\displaystyle e\_\{\\mathrm\{POD\}\}\(t\)=‖𝝎POD​\(t\)−𝝎​\(t\)‖2‖𝝎​\(t\)‖2,\\displaystyle=\\frac\{\\\|\\boldsymbol\{\\omega\}\_\{\\mathrm\{POD\}\}\(t\)\-\\boldsymbol\{\\omega\}\(t\)\\\|\_\{2\}\}\{\\\|\\boldsymbol\{\\omega\}\(t\)\\\|\_\{2\}\},\(S23\)eVG​\-​POD​\(t\)\\displaystyle e\_\{\\mathrm\{VG\\text\{\-\}POD\}\}\(t\)=‖𝝎^VG​\(t\)−𝝎POD​\(t\)‖2‖𝝎POD​\(t\)‖2,\\displaystyle=\\frac\{\\\|\\widehat\{\\boldsymbol\{\\omega\}\}\_\{\\mathrm\{VG\}\}\(t\)\-\\boldsymbol\{\\omega\}\_\{\\mathrm\{POD\}\}\(t\)\\\|\_\{2\}\}\{\\\|\\boldsymbol\{\\omega\}\_\{\\mathrm\{POD\}\}\(t\)\\\|\_\{2\}\},eend​\(t\)\\displaystyle e\_\{\\mathrm\{end\}\}\(t\)=‖𝝎^VG​\(t\)−𝝎​\(t\)‖2‖𝝎​\(t\)‖2\.\\displaystyle=\\frac\{\\\|\\widehat\{\\boldsymbol\{\\omega\}\}\_\{\\mathrm\{VG\}\}\(t\)\-\\boldsymbol\{\\omega\}\(t\)\\\|\_\{2\}\}\{\\\|\\boldsymbol\{\\omega\}\(t\)\\\|\_\{2\}\}\.Their time\-averaged values are 12\.53%, 8\.77%, and 15\.56%, respectively\. Figure[S4](https://arxiv.org/html/2608.02662#A4.F4)\(d\) shows their variation over the final\-test interval\. All three errors remain bounded without systematic growth\. The smaller VG\-to\-POD error indicates that the symbolic rollout remains closer to the retained six\-mode trajectory than the POD reconstruction is to the simulation\. The end\-to\-end error reflects both the POD truncation and the symbolic\-rollout discrepancy, but it is not their arithmetic sum because the error fields are neither collinear nor normalized by the same reference\.

![Refer to caption](https://arxiv.org/html/2608.02662v1/x8.png)Figure S4:Final\-test fixed\-Re\\mathrm\{Re\}comparisons: phase portraits for \(a\)a1a\_\{1\}–a2a\_\{2\}, \(b\)a3a\_\{3\}–a4a\_\{4\}, and \(c\)a5a\_\{5\}–a6a\_\{6\}; \(d\) POD truncation, VG\-to\-POD, and end\-to\-end field errors\.Representative vorticity fields are shown in Figure[S5](https://arxiv.org/html/2608.02662#A4.F5)\. The POD reconstruction preserves the alternating vortex street but smooths the smaller\-scale vorticity gradients omitted by the six\-mode representation\. At all three times, the VG reconstruction closely follows the shedding phase and large\-scale organization of the POD field\. Its remaining deviation from the simulation reflects the POD truncation visible in the middle column together with the smaller additional differences introduced by the symbolic coordinate rollout\.

![Refer to caption](https://arxiv.org/html/2608.02662v1/x9.png)Figure S5:Simulation, rank\-six POD reconstruction, and VG prediction at three times in the final\-test fixed\-Re\\mathrm\{Re\}interval\.

## Appendix S5Autoencoder Representation–Discoverability Sensitivity

### S5\.1Controlled architectures, seeds, and optimization protocol

Unlike POD coordinates, which follow directly from a covariance decomposition of a specified snapshot set up to sign and rotations within degenerate subspaces, nonlinear autoencoders do not define a unique latent embedding\. Different architectures, initializations, and optimization paths may yield similar field reconstructions while producing latent trajectories with different temporal geometry\. This non\-uniqueness motivates examining representation fidelity together with properties relevant to symbolic discoverability\.

The controlled representation study crossed latent dimensions three and four with shallow and deep fully connected autoencoders\. The shallow encoder used one hidden layer of width 256, whereas the deep encoder used hidden widths 4096 and 256; both used symmetric decoders, as shown in Figure[S6](https://arxiv.org/html/2608.02662#A5.F6)\. The independently prepared autoencoder used for the cross\-Re\\mathrm\{Re\}result in Section[5\.3](https://arxiv.org/html/2608.02662#S5.SS3)has the same architecture as the shallow three\-coordinate autoencoders tested here\. It was trained before the controlled study without a prespecified, recorded initialization seed and is therefore retained as a separate comparison rather than treated as an additional controlled seed\.

![Refer to caption](https://arxiv.org/html/2608.02662v1/x10.png)Figure S6:Fully connected shallow and deep autoencoder architectures used in the representation study\. Hidden layers use ReLU activation; latent and output layers are linear\.Each configuration was initialized from seeds 17, 29, and 43 and optimized using the seven Reynolds numbers employed for symbolic model construction and the same 500 post\-transient snapshots per Reynolds number\. All comparisons in Table[S6](https://arxiv.org/html/2608.02662#A5.T6)use these post\-transient fields\. One predeclared seed per architecture was carried forward to the symbolic discovery comparison\.

### S5\.2Definition of the latent\-representation diagnostics

Three diagnostics characterize the temporal properties of each embedding before symbolic model construction\. Let𝐙\(m\)∈ℝT×dz\\mathbf\{Z\}^\{\(m\)\}\\in\\mathbb\{R\}^\{T\\times d\_\{z\}\}contain the latent trajectory at themmth Reynolds number, whereTTis the number of snapshots anddzd\_\{z\}the latent dimension\. Each coordinate is standardized using its meanμj\\mu\_\{j\}and standard deviationσj\\sigma\_\{j\}on the concatenated trajectories:Z~t,j\(m\)=\(Zt,j\(m\)−μj\)/\(σj\+ϵ\)\\widetilde\{Z\}\_\{t,j\}^\{\(m\)\}=\(Z\_\{t,j\}^\{\(m\)\}\-\\mu\_\{j\}\)/\(\\sigma\_\{j\}\+\\epsilon\)\. WithΔ\\DeltaandΔ2\\Delta^\{2\}denoting first and second temporal differences, respectively, the latent roughness is

ℛz=1MR​∑m=1MRmeant,j⁡\[\(Δ2​Z~t,j\(m\)\)2\]meant,j⁡\[\(Δ​Z~t,j\(m\)\)2\]\+ϵ,\\mathcal\{R\}\_\{z\}=\\frac\{1\}\{M\_\{R\}\}\\sum\_\{m=1\}^\{M\_\{R\}\}\\frac\{\\operatorname\{mean\}\_\{t,j\}\\left\[\(\\Delta^\{2\}\\widetilde\{Z\}\_\{t,j\}^\{\(m\)\}\)^\{2\}\\right\]\}\{\\operatorname\{mean\}\_\{t,j\}\\left\[\(\\Delta\\widetilde\{Z\}\_\{t,j\}^\{\(m\)\}\)^\{2\}\\right\]\+\\epsilon\},\(S24\)whereMR=7M\_\{R\}=7is the number of Reynolds numbers\. The denominator normalizes by the overall first\-difference scale, so Eq\. \([S24](https://arxiv.org/html/2608.02662#A5.E24)\) measures rapid changes in temporal slope rather than latent amplitude\. Smaller values indicate smoother trajectories\.

For spectral concentration, letPk,j\(m\)P\_\{k,j\}^\{\(m\)\}be the discrete Fourier power of the centered standardized coordinate at positive, nonzero frequency binkk, and let𝒦3\(m,j\)\\mathcal\{K\}\_\{3\}^\{\(m,j\)\}contain its three largest\-power bins\. The reported diagnostic is

𝒮z=1MR​dz​∑m=1MR∑j=1dz∑k∈𝒦3\(m,j\)Pk,j\(m\)∑k\>0Pk,j\(m\)\+ϵ\.\\mathcal\{S\}\_\{z\}=\\frac\{1\}\{M\_\{R\}d\_\{z\}\}\\sum\_\{m=1\}^\{M\_\{R\}\}\\sum\_\{j=1\}^\{d\_\{z\}\}\\frac\{\\sum\_\{k\\in\\mathcal\{K\}\_\{3\}^\{\(m,j\)\}\}P\_\{k,j\}^\{\(m\)\}\}\{\\sum\_\{k\>0\}P\_\{k,j\}^\{\(m\)\}\+\\epsilon\}\.\(S25\)Values of𝒮z\\mathcal\{S\}\_\{z\}in Eq\. \([S25](https://arxiv.org/html/2608.02662#A5.E25)\) closer to one indicate that most temporal variation is organized around a small number of frequencies, as expected for a dominant shedding oscillator and its harmonics\.

Finally, concatenate the unstandardized trajectories into𝐙all\\mathbf\{Z\}\_\{\\mathrm\{all\}\}and letλ1,…,λdz≥0\\lambda\_\{1\},\\ldots,\\lambda\_\{d\_\{z\}\}\\geq 0be the eigenvalues of its sample covariance matrix\. The covariance participation ratio is

dPR=\(∑j=1dzλj\)2∑j=1dzλj2\+ϵ\.d\_\{\\mathrm\{PR\}\}=\\frac\{\\left\(\\sum\_\{j=1\}^\{d\_\{z\}\}\\lambda\_\{j\}\\right\)^\{2\}\}\{\\sum\_\{j=1\}^\{d\_\{z\}\}\\lambda\_\{j\}^\{2\}\+\\epsilon\}\.\(S26\)This effective dimension equals one when nearly all covariance is concentrated in one direction and increases as variance is distributed more evenly among latent coordinates\. Equations \([S24](https://arxiv.org/html/2608.02662#A5.E24)\)–\([S26](https://arxiv.org/html/2608.02662#A5.E26)\) describe complementary properties: temporal regularity, spectral organization, and covariance use of the available latent dimensions\.

### S5\.3Latent\-representation diagnostics evaluation

Table[S6](https://arxiv.org/html/2608.02662#A5.T6)compares reconstruction fidelity with the three latent diagnostics\. The smallest reconstruction error, 4\.494%, was attained by the deep three\-coordinate model initialized with seed 43\. This representation did not have the smoothest latent trajectories: its roughness was 0\.0465, compared with 0\.0224 for the independent shallow reference\. The latter also had the largest spectral concentration, 0\.9264, despite its higher reconstruction error of 5\.826%\. The comparison therefore separates field reconstruction from temporal properties that affect compatibility with the symbolic backbone\. As the symbolic\-discovery outcomes in Table[S7](https://arxiv.org/html/2608.02662#A5.T7)subsequently show, an accurate snapshot reconstruction need not yield latent dynamics that are readily expressed within the backbone’s learned symbolic distribution\.

Table S6:Autoencoder representation sensitivity across architecture and initialization seed\. Field error is the temporal mean of snapshot\-wise relative spatialL2L\_\{2\}error; the latent diagnostics are defined in Eqs\. \([S24](https://arxiv.org/html/2608.02662#A5.E24)\)–\([S26](https://arxiv.org/html/2608.02662#A5.E26)\)\.∗Seed carried forward for the corresponding controlled architecture\.

∗∗The independent shallow reference is the separately prepared autoencoder used for the cross\-parameter result reported in Section[5\.3](https://arxiv.org/html/2608.02662#S5.SS3)\. It has the same architecture as the shallow three\-coordinate autoencoders but no prespecified, recorded initialization seed and is not one of the 12 controlled runs\.

### S5\.4Corresponding symbolic\-discovery outcomes

The same beam size of 20 and sampling temperature of 0\.1 were used for the controlled symbolic\-discovery comparison and for the independent shallow reference\. Among the four controlled encodings identified by asterisks in Table[S6](https://arxiv.org/html/2608.02662#A5.T6), only AE3D produced a pre\-optimization candidate satisfying every verifier\. Its multi\-trajectory rollout loss was 2\.169; the post\-optimization system was not admissible\. By contrast, the independently prepared shallow reference used for the cross\-parameter result in Section[5\.3](https://arxiv.org/html/2608.02662#S5.SS3)yielded an admissible pre\-optimization candidate with loss 1\.042, and post\-optimization reduced the loss to 0\.200 while preserving admissibility\. These outcomes are summarized in Table[S7](https://arxiv.org/html/2608.02662#A5.T7)\.

Table S7:Symbolic\-discovery outcomes for the representations advanced from the controlled study and for the independent shallow reference\.The comparison supports a specific conclusion: reconstruction fidelity alone did not determine whether VG could identify an admissible cross\-parameter equation\. It does not imply that roughness or spectral concentration is by itself sufficient for discoverability; these diagnostics describe relevant representation properties rather than a complete selection criterion\. This result further supports the joint coordinate\-and\-equation discovery direction discussed in Section[S2\.5](https://arxiv.org/html/2608.02662#A2.SS5), in which dynamical verification can help distinguish latent embeddings with similar reconstruction fidelity but different compatibility with symbolic discovery\.

## Appendix S6Extended Results for the Cross\-Parameter Cylinder\-Flow Case

### S6\.1Representation and evaluation protocol

The cross\-Re\\mathrm\{Re\}study used the same shallow autoencoder as the result reported in Section[5\.3](https://arxiv.org/html/2608.02662#S5.SS3); Section[S5](https://arxiv.org/html/2608.02662#A5)examined this representation further\. The autoencoder operates on mean\-subtracted \(fluctuation\) vorticity fields\. Its three latent coordinates were centered and scaled using the seven Reynolds numbers employed for symbolic model construction,

Re∈\{150,200,250,300,350,400,450\}\.\\mathrm\{Re\}\\in\\\{150,200,250,300,350,400,450\\\}\.The Reynolds number was appended asr=\(Re−300\)/150r=\(\\mathrm\{Re\}\-300\)/150\. The symbolic backbone was given the uniformly sampled trajectories at unit time increments, so its model time isτ=t∗/0\.1\\tau=t^\{\*\}/0\.1, where consecutive physical snapshots are separated byΔ​t∗=0\.1\\Delta t^\{\*\}=0\.1\. This rescaling changes the numerical values of the equation coefficients but not the represented trajectories\. At each Reynolds number, the model state was therefore𝝃=\(z1,z2,z3,r\)𝖳\\boldsymbol\{\\xi\}=\(z\_\{1\},z\_\{2\},z\_\{3\},r\)^\{\\mathsf\{T\}\}\.

For every Reynolds number, the first 500 simulation snapshots were removed as the wake\-development transient\. Model construction and evaluation use the subsequent 500 post\-transient snapshots, spanning approximately50≤t∗≤10050\\leq t^\{\*\}\\leq 100\. The withheld cases comprise interpolation atRe=175,275,\\mathrm\{Re\}=175,275,and425425and extrapolation atRe=500\\mathrm\{Re\}=500\. None of these four latent trajectories contributed to generation, ranking, verification, or coefficient optimization of the symbolic model\. The three interpolation values were also excluded when fitting the autoencoder\. The autoencoder was fitted using fields atRe=500\\mathrm\{Re\}=500, but the corresponding latent trajectory was excluded from symbolic model construction\. Thus, “extrapolation” atRe=500\\mathrm\{Re\}=500refers specifically to evaluating the Reynolds\-conditioned equation atr=4/3r=4/3, beyond its construction interval−1≤r≤1\-1\\leq r\\leq 1; it does not refer to extrapolation of the field representation\.

Candidates generated from these seven trajectories were verified using the cross\-Re\\mathrm\{Re\}settings in Table[S2](https://arxiv.org/html/2608.02662#A2.T2)\. The two\-branch shortlist contained 17 distinct systems, of which four satisfied all applicable verifiers\. Coefficient optimization was then applied to the admissible candidates, followed by the same verifier suite\. The resulting system was chosen before the withheld Reynolds numbers were examined\. To retain the notation used elsewhere in the supplement, dots in Eqs\. \([S27](https://arxiv.org/html/2608.02662#A6.E27)\) and \([S28](https://arxiv.org/html/2608.02662#A6.E28)\) denote differentiation with respect to the model timeτ\\tau\. The pre\-optimization and post\-optimization systems are

z˙1\\displaystyle\\dot\{z\}\_\{1\}=0\.1070​z3−0\.0013​r\+0\.0089​z2​z3,\\displaystyle=0\.1070z\_\{3\}\-0\.0013r\+0\.0089z\_\{2\}z\_\{3\},z˙2\\displaystyle\\dot\{z\}\_\{2\}=0\.1003​z1,\\displaystyle=0\.1003z\_\{1\},z˙3\\displaystyle\\dot\{z\}\_\{3\}=−0\.1098​z1,\\displaystyle=\-0\.1098z\_\{1\},r˙\\displaystyle\\dot\{r\}=0,\\displaystyle=0,\(S27\)z˙1\\displaystyle\\dot\{z\}\_\{1\}=0\.114317​z3\+0\.004217​r\+0\.020148​z2​z3,\\displaystyle=0\.114317z\_\{3\}\+0\.004217r\+0\.020148z\_\{2\}z\_\{3\},z˙2\\displaystyle\\dot\{z\}\_\{2\}=0\.079842​z1,\\displaystyle=0\.079842z\_\{1\},z˙3\\displaystyle\\dot\{z\}\_\{3\}=−0\.099544​z1,\\displaystyle=\-0\.099544z\_\{1\},r˙\\displaystyle\\dot\{r\}=0\.\\displaystyle=0\.\(S28\)Since one unit ofτ\\taucorresponds to0\.10\.1units oft∗t^\{\*\}, multiplying the right\-hand sides of Eq\. \([S28](https://arxiv.org/html/2608.02662#A6.E28)\) by ten gives the physical\-time equation reported in Eq\. \([13](https://arxiv.org/html/2608.02662#S5.E13)\)\.

### S6\.2Latent\-coordinate rollouts at withheld Reynolds numbers

Figure[S7](https://arxiv.org/html/2608.02662#A6.F7)compares the latent trajectories obtained by encoding the simulated fields with the pre\-optimization and post\-optimization VG rollouts\. We retainzjz\_\{j\}for these nonlinear autoencoder coordinates, consistent with the notation in Section[5\.3](https://arxiv.org/html/2608.02662#S5.SS3)and distinct from the POD coefficientsaja\_\{j\}in Figures[S3](https://arxiv.org/html/2608.02662#A4.F3)and[S4](https://arxiv.org/html/2608.02662#A4.F4)\. For an encoded trajectory𝐳​\(tn\)\\mathbf\{z\}\(t\_\{n\}\)and its prediction𝐳^​\(tn\)\\widehat\{\\mathbf\{z\}\}\(t\_\{n\}\), the latent normalized mean\-squared error is

ℰz=13​∑j=13N−1​∑n=1N\(z^j​\(tn\)−zj​\(tn\)\)2Varn⁡\[zj​\(tn\)\]\+ϵ,\\mathcal\{E\}\_\{z\}=\\frac\{1\}\{3\}\\sum\_\{j=1\}^\{3\}\\frac\{N^\{\-1\}\\sum\_\{n=1\}^\{N\}\\left\(\\widehat\{z\}\_\{j\}\(t\_\{n\}\)\-z\_\{j\}\(t\_\{n\}\)\\right\)^\{2\}\}\{\\operatorname\{Var\}\_\{n\}\\\!\\left\[z\_\{j\}\(t\_\{n\}\)\\right\]\+\\epsilon\},\(S29\)whereN=500N=500andϵ\\epsilonis a small numerical constant\. This metric weights each coordinate relative to its own variation; the conserved parameter coordinate is not included\.

Coefficient optimization reducesℰz\\mathcal\{E\}\_\{z\}at every withheld Reynolds number:

- •Re=175\\mathrm\{Re\}=175:2\.312→0\.1422\.312\\rightarrow 0\.142, a 93\.9% reduction;
- •Re=275\\mathrm\{Re\}=275:0\.797→0\.2470\.797\\rightarrow 0\.247, a 69\.0% reduction;
- •Re=425\\mathrm\{Re\}=425:0\.188→0\.0430\.188\\rightarrow 0\.043, a 77\.2% reduction;
- •Re=500\\mathrm\{Re\}=500:0\.714→0\.0360\.714\\rightarrow 0\.036, a 94\.9% reduction\.

All post\-optimization rollouts remain finite and oscillatory over the complete interval\. The largest remaining error occurs atRe=275\\mathrm\{Re\}=275, where a progressive phase difference is visible in all three coordinates without a loss of the oscillatory regime\. The error ordering is not monotonic with distance from the Reynolds numbers used for model construction\. It also depends on how accurately the selected equation represents the variation of latent frequency, phase, and amplitude across the parameter range, and on the accumulation of small frequency differences over a long rollout\. In Eq\. \([S28](https://arxiv.org/html/2608.02662#A6.E28)\),rrenters explicitly through an additive term inz˙1\\dot\{z\}\_\{1\}, while the oscillator couplings are shared across Reynolds number\. The model can therefore represent a Reynolds\-dependent shift of the common oscillator, but it does not allow every coupling or frequency coefficient to vary independently withRe\\mathrm\{Re\}\. The strong result atRe=500\\mathrm\{Re\}=500indicates that the post\-optimization equation continues the same periodic shedding dynamics beyond the parameter range used to construct it\.

![Refer to caption](https://arxiv.org/html/2608.02662v1/x11.png)Figure S7:Autoencoder\-encoded trajectories and the corresponding pre\-optimization and post\-optimization VG rollouts at the withheld Reynolds numbers\. I denotes interpolation within the symbolic model\-construction range, and E denotes extrapolation beyond that range\.
### S6\.3Spatially averaged wake evolution

The latent comparison can be related to an observable quantity in the physical field by averaging the nondimensional vorticity over the wake region

𝒲=\{\(x/D,y/D\):0\.5≤x/D≤5,\|y/D\|≤1\.5\}\.\\mathcal\{W\}=\\left\\\{\(x/D,y/D\):0\.5\\leq x/D\\leq 5,\\;\|y/D\|\\leq 1\.5\\right\\\}\.For a fieldω∗​\(𝐱,t∗\)=ω​D/U∞\\omega^\{\*\}\(\\mathbf\{x\},t^\{\*\}\)=\\omega D/U\_\{\\infty\}, its spatial average is

⟨ω∗⟩𝒲​\(t∗\)=1\|𝒲\|​∫𝒲ω∗​\(𝐱,t∗\)​d𝐱\.\\left\\langle\\omega^\{\*\}\\right\\rangle\_\{\\mathcal\{W\}\}\(t^\{\*\}\)=\\frac\{1\}\{\|\\mathcal\{W\}\|\}\\int\_\{\\mathcal\{W\}\}\\omega^\{\*\}\(\\mathbf\{x\},t^\{\*\}\)\\,\\mathrm\{d\}\\mathbf\{x\}\.\(S30\)This integral measures the net balance of positive and negative vorticity within the specified region and provides a scalar record of shedding phase and amplitude\.

![Refer to caption](https://arxiv.org/html/2608.02662v1/x12.png)Figure S8:Spatially averaged wake vorticity from the simulation, autoencoder reconstruction, and post\-optimization VG prediction\. I denotes interpolation within the symbolic model\-construction range; E denotes extrapolation beyond that range\.The autoencoder traces in Figure[S8](https://arxiv.org/html/2608.02662#A6.F8)closely follow the simulation at all four Reynolds numbers, indicating that the common representation preserves this integral wake quantity\. The VG traces retain the oscillation frequency and amplitude range throughout each withheld interval\. Agreement is strongest atRe=425\\mathrm\{Re\}=425and500500\. AtRe=175\\mathrm\{Re\}=175the predicted phase gradually separates from the reconstructed field, whileRe=275\\mathrm\{Re\}=275exhibits the largest late\-time phase difference, consistent with its larger latent error\. All four cases remain within the periodic vortex\-shedding regime, whose dominant oscillator provides a common dynamical structure across Reynolds number\. The residual differences are consistent with the selected Reynolds\-conditioned equation capturing this shared oscillator more accurately than the detailed variation of its phase and amplitude\.

### S6\.4Field reconstruction and error decomposition

The field comparison uses the same three\-part error decomposition as Eq\. \([S23](https://arxiv.org/html/2608.02662#A4.E23)\)\. Here the POD reconstruction𝝎POD\\boldsymbol\{\\omega\}\_\{\\mathrm\{POD\}\}is replaced by the autoencoder reconstruction𝝎AE\\boldsymbol\{\\omega\}\_\{\\mathrm\{AE\}\}\. Accordingly,eAEe\_\{\\mathrm\{AE\}\}andeVG​\-​AEe\_\{\\mathrm\{VG\\text\{\-\}AE\}\}replaceePODe\_\{\\mathrm\{POD\}\}andeVG​\-​PODe\_\{\\mathrm\{VG\\text\{\-\}POD\}\}, respectively, while the end\-to\-end definition is unchanged\.

Table[S8](https://arxiv.org/html/2608.02662#A6.T8)summarizes the latent and field errors over all 500 snapshots\. The mean autoencoder error varies only from 5\.90% to 6\.96% across the four cases, whereas the VG\-to\-autoencoder error varies more substantially\. The difference in end\-to\-end accuracy is therefore associated mainly with the symbolic rollout rather than a deterioration of the common field representation\. As in the fixed\-Re\\mathrm\{Re\}decomposition, the end\-to\-end error is not the arithmetic sum of its two components\.

Table S8:Extended cross\-Re\\mathrm\{Re\}latent and field errors\. Field entries give the temporal mean with the median in parentheses, in percent\.Figure[S9](https://arxiv.org/html/2608.02662#A6.F9)consolidates all four withheld Reynolds numbers att∗=75t^\{\*\}=75, including theRe=275\\mathrm\{Re\}=275and500500cases presented in the paper, so that the field predictions can be compared directly\. The autoencoder consistently smooths smaller spatial features while preserving the alternating wake\. The VG prediction reconstructs the dominant structures available in the three\-coordinate representation at every Reynolds number\. Its greater displacement from the encoded field atRe=275\\mathrm\{Re\}=275is consistent with the accumulated phase difference in Figures[S7](https://arxiv.org/html/2608.02662#A6.F7)and[S8](https://arxiv.org/html/2608.02662#A6.F8); the closer latent rollouts atRe=425\\mathrm\{Re\}=425and500500correspond to more closely aligned large\-scale fields\.

As discussed in Section[S6\.2](https://arxiv.org/html/2608.02662#A6.SS2), the recovered equation represents a Reynolds\-conditioned common oscillator rather than a complete symbolic description of all parameter\-dependent wake dynamics\. Even with this qualification, obtaining one explicit system that remains bounded and predictive at four Reynolds numbers excluded from symbolic model construction—including one outside the construction range—is nontrivial for reduced coordinates of a high\-dimensional flow, particularly when using a synthetically pretrained symbolic transformer as the backbone\. Parametric symbolic modeling of physical flow systems remains an open challenge; these results demonstrate the value of combining shared oscillatory structure, explicit parameter dependence, and cross\-trajectory verification\.

![Refer to caption](https://arxiv.org/html/2608.02662v1/x13.png)Figure S9:Simulation, autoencoder reconstruction, and VG prediction att∗=75t^\{\*\}=75for all Reynolds numbers excluded from symbolic model construction\. The first three rows are interpolation cases; the final row is extrapolation of the symbolic dynamics\.