Paper deep dive
Verifier-Guided Model Discovery for Physical Dynamical Systems with Pretrained Symbolic Transformers
Farbod Faraji, Francesco Belardinelli
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%
Last extracted: 8/5/2026, 3:54:36 AM
Summary
The paper introduces a Verifier-Guided (VG) workflow for symbolic model discovery using pretrained transformers (ODEFormer). It addresses the challenge of transferring symbolic models from synthetic pretraining to high-dimensional physical systems by using dynamical and physical-admissibility criteria to select equations from a multi-trajectory candidate pool. The method is evaluated on Van der Pol oscillators and vortex shedding phenomena, demonstrating superior generalization across initial conditions and Reynolds numbers compared to baseline methods, while maintaining interpretability and physical consistency.
Entities (7)
Relation Signals (5)
Verifier-Guided Workflow → uses → ODEFormer
confidence 98% · We develop a verifier-guided (VG) workflow around ODEFormer as a symbolic backbone
Verifier-Guided Workflow → applies → Physical Admissibility
confidence 95% · using dynamical and physical-admissibility criteria to select from a multi-trajectory candidate equation pool
Van der Pol Oscillator → testedwith → Verifier-Guided Workflow
confidence 95% · On canonical Van der Pol oscillators, VG outperforms the original ODEFormer workflow
Verifier-Guided Workflow → outperforms → ODEFormer
confidence 92% · On canonical Van der Pol oscillators, VG outperforms the original ODEFormer workflow across held-out initial conditions
Verifier-Guided Workflow → discovers → Vortex Shedding
confidence 90% · We then address vortex shedding... through coordinate reduction and symbolic discovery
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
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.
Tags
Links
- Source: https://arxiv.org/abs/2608.02662v1
- Canonical: https://arxiv.org/abs/2608.02662v1
Trouble viewing inline? Open PDF directly →
Full Text
121,365 characters extracted from source content.
Expand or collapse full text
Verifier-Guided Model Discovery for Physical Dynamical Systems with Pretrained Symbolic Transformers 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. 1 Introduction 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), 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; Etling 1990); related instabilities drive vortex-induced loading (Williamson and Govardhan 2004) and arise in plasma flows past obstacles (Gruszecki et al. 2010; Bailung et al. 2020). 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. 2 Related 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; Faraji et al. 2024). 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 (Chen et al. 2018; Li et al. 2021). 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 (Raissi et al. 2019; Greydanus et al. 2019). 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 (Lu et al. 2021). 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 (Brunton et al. 2016). 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; Fasel et al. 2022; Fung et al. 2025). 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). Phi Method identifies discretized evolution operators from a candidate library, avoiding explicit derivatives and continuous-time integration (Faraji et al. 2025). AutoSINDy (Basiri and Nicholson 2026) combines PySR candidate generation (Cranmer 2023), 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(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; Weilbach et al. 2021; Brence et al. 2021; Omejc et al. 2024). 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 (Petersen et al. 2021; Biggio et al. 2021; Valipour et al. 2021; Kamienny et al. 2022; Shojaee et al. 2023). ODEFormer specializes this formulation to dynamics discovery by mapping multivariate trajectories directly to coupled ODE systems (d’Ascoli et al. 2024). MIO trains a trajectory-to-equation transformer to infer shared dynamics from multiple trajectories of the same system (Şahin et al. 2025). Diffusion-based symbolic regression provides a non-autoregressive alternative through masked iterative denoising and dataset-specific reinforcement learning (Bastiani et al. 2025). 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 (Taskin et al. 2026; Xiao et al. 2026; Yu et al. 2026). 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’Ascoli et al. 2024), 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. 3 Background and Problem Formulation Consider M observed trajectories (m)(t)x^(m)(t), where m indexes different initial conditions, time windows, or parameter values. Symbolic dynamics discovery seeks an explicit model ˙=(),:⊆ℝd→ℝd, x=f(x), :X ^d ^d, (1) whose initial-value problems yield finite trajectories over the state domain and time horizons of interest. Equation (1) 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 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’Ascoli et al. 2024). 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’Ascoli et al. 2024). 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 of 2020 and temperature of 0.10.1; their sensitivity is examined on the canonical Van der Pol problem in Technical Appendix S1. 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). For a spatially distributed physical state (t)∈ℝN ω(t) ^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 ℰE and reconstruction map ℛR define (t)=ℰ((t)),(t)≈ℛ((t)),˙=(),z(t)=E\! ( ω(t) ), ω(t) \! (z(t) ), z=g(z), (2) where (t)∈ℝdz(t) ^d and d≪Nd N. Eq. (2) 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, 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 parameter p, we append the parameter to the state and constrain it to remain constant: ˙=(,p),p˙=0. z=g(z,p), p=0. (3) The augmentation in Eq. (3) 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. 4 Methodology: Verifier-Guided Symbolic Model Discovery Figure 1 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. 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) ξ^(m) denote the trajectory or window supplied to the backbone: (m)=(m) ξ^(m)=x^(m) for observed ODE states, (m)=(m) ξ^(m)=z^(m) for reduced fields, and (m)=[(m)⊤,p(m)]⊤ ξ^(m)=[z^(m) ,p^(m)] for cross-parameter discovery. For input m, the backbone retains BmB_m decoded candidate systems indexed by b: m=m,bb=1Bm,=⋃m=1Mm.C_m=\h_m,b\_b=1^B_m, = _m=1^MC_m. (4) In Eq. (4), symbolic canonicalization removes duplicates before the union, which constitutes pooling. Each ∈h is ranked across discovery trajectories by Eq. (5): ℒroll()=1M∑m=1M1d∑j=1dMSEi(ξ^ij(m),ξij(m))Vari(ξij(m))+ϵ,L_roll(h)= 1M _m=1^M 1d _j=1^d MSE_i ( ξ_ij^(m), _ij^(m) )Var_i ( _ij^(m) )+ε, (5) where d is the number of modeled dynamical coordinates, i indexes time samples, and ^(m) ξ^(m) is the candidate rollout initialized from the first observed state of trajectory m. The appended parameter is excluded from the d coordinates in cross-parameter discovery. Nonfinite rollouts receive infinite loss, and the lowest-loss candidates form the shortlist KC_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 KC_K. Each candidate defines the vector field ˙=() ξ=h( ξ). Equation (6) measures its local agreement with observed dynamics: Vloc()=1d∑j=1dMSEm,i(hj(i(m)),ξ˙ij(m))Varm,i(ξ˙ij(m))+ϵV_loc(h)= 1d _j=1^d MSE_m,i (h_j( ξ_i^(m)), ξ_ij^(m) )Var_m,i ( ξ_ij^(m) )+ε (6) where ξ˙ij(m) ξ_ij^(m) is the numerically estimated rate of coordinate j at observed sample i, 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) defines the admissible set and selects within it: =∈K:ℓ()=1for every verifier ℓ, = \h _K: pass_ (h)=1\ for every verifier \, (7) ⋆ =argmin∈ℒroll(). = h \,L_roll(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 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 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 ℒrollL_roll across the discovery trajectories. The refined equation is reverified; if it becomes inadmissible, VG returns the admissible unoptimized equation ⋆h . 5 Empirical 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.1 Canonical 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). Its first-order form is given in Eq. (8): x˙0 x_0 =x1, =x_1, (8) x˙1 x_1 =μ(1−x02)x1−x0 =μ(1-x_0^2)x_1-x_0 =μx1−x0−μx02x1. =μ x_1-x_0-μ x_0^2x_1. In Eq. (8), x˙0=x1 x_0=x_1 makes x1x_1 velocity-like; in x˙1 x_1, −x0-x_0 restores while μx1μ x_1 injects energy near the origin and −μx02x1-μ x_0^2x_1 dissipates it at larger amplitudes. Their balance produces a stable limit cycle. We consider μ=0.5μ=0.5 and 1.51.5 as 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 μ, the original ODEFormer protocol (d’Ascoli et al. 2024) 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 reports the exact assignments. 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 shows that the VG model outperforms the single-trajectory protocol for all eight held-out initial conditions at both values of μ. For μ=0.5μ=0.5, the median post-optimization error decreases from 1.211.21 to 6.99×10−46.99× 10^-4; for μ=1.5μ=1.5, it decreases from 4.5×10−24.5× 10^-2 to 3.27×10−33.27× 10^-3. The larger separation at μ=0.5μ=0.5 is 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 μ, 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)–(10) were selected: μ=0.5:x˙0 μ=5: x_0 =0.997x1−0.012x0, =997x_1-012x_0, (9) x˙1 x_1 =0.542x1−0.997x0 =542x_1-997x_0 −0.543x02x1−0.0039x0x1, -543x_0^2x_1-0039x_0x_1, and μ=1.5:x˙0 μ=5: x_0 =0.996x1 =996x_1 (10) +0.000676(0.998+1.273x0)2, +000676(998+273x_0)^2, x˙1 x_1 =1.440x1−1.002x0 =440x_1-002x_0 −1.461x02x1−0.0080x0x1. -461x_0^2x_1-0080x_0x_1. Both systems recover the defining Van der Pol structure: the kinematic relation x˙0≈x1 x_0≈ x_1, linear restoring dynamics through −x0-x_0, and the opposing linear and cubic damping terms in x˙1 x_1. Their principal coefficients closely approximate the corresponding ground-truth values, while the additional x0x_0, x0x1x_0x_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.2 High-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 ω satisfies ∂ω∂t+(⋅∇)ω ∂ω∂ t+(u·∇)ω =1Re∇2ω, = 1Re∇^2ω, (11) ∇⋅ ∇·u =0,Re=U∞Dν. =0, Re= U_∞Dν. In Eq. (11), (⋅∇)ω(u·∇)ω represents advection of vorticity by the local flow, whereas Re−1∇2ωRe^-1∇^2ω represents viscous diffusion of vorticity gradients. Increasing ReRe reduces the relative contribution of this diffusion and permits the wake instability to develop. At Re=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). Experimental setup. Vorticity snapshots were generated with the open-source ViscousFlow.jl solver (Eldredge 2021), which implements the immersed-layer method (Eldredge 2022). A unit-diameter cylinder was centered in x/D∈[−1,5]x/D∈[-1,5] and y/D∈[−2,2]y/D∈[-2,2], with U∞=1U_∞=1, a 2∘2 free-stream incidence, and a no-slip surface. The outputs were interpolated onto a 300×200300× 200 grid and sampled at Δt∗=ΔtU∞/D=0.1 t^*= tU_∞/D=0.1, yielding 1,000 snapshots. Further numerical and preprocessing details are reported by Faraji et al. (2025). 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 ranks 44 and 66 and window counts 88 and 2424. 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 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 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 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. Quantity Error (%) Median five-time-unit (single-window) coordinate rollout 3.49 Full-interval coordinate rollout 6.89 POD field reconstruction 12.53 VG field relative to POD reconstruction 8.77 End-to-end field prediction 15.56 Table 1: Fixed-ReRe final-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 relative L2L_2 errors. Figure 3: Wake-averaged nondimensional vorticity over the final-test interval at Re=300Re=300, evaluated over 0.5≤x/D≤50.5≤ x/D≤ 5 and |y/D|≤1.5 y/D ≤ 1.5. Equation (12) gives the coefficient-optimized VG model in standardized rank-6 POD coordinates z1,…,z6z_1,…,z_6: z˙1 z_1 =1.001z2 =001z_2 (12) −0.0679(11.746−0.8668z2)−1.052, -0679(1746-8668z_2)^-1.052, z˙2 z_2 =−1.123z1, =-123z_1, z˙3 z_3 =−2.099z4, =-099z_4, z˙4 z_4 =2.120z3, =120z_3, z˙5 z_5 =3.121z6 =121z_6 −0.1164sin(0.1144+11.036z2), -1164 (1144+1036z_2), z˙6 z_6 =−0.1458−3.233z5. =-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 of 1:1.99:3.001:1.99:3.00 identifies the fundamental shedding cycle and its second and third harmonics. The reciprocal contribution to z˙1 z_1 deforms the fundamental oscillator from exact linear motion, while the sinusoidal dependence of z˙5 z_5 on z2z_2 couples the fundamental and third-harmonic coordinate pairs. The constant term in z˙6 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.3 Cross-Parameter Extension: Flow Past a Cylinder across Reynolds Numbers The fixed-ReRe result 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 used Re∈150,200,250,300,350,400,450Re∈\150,200,250,300,350,400,450\. The withheld test values were Re∈175,275,425Re∈\175,275,425\ for interpolation and Re=500Re=500 for extrapolation. The simulations and field preprocessing followed the fixed-ReRe case. 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 reports more details on the outcomes of architecture and initialization-seed sensitivity, latent diagnostics, and corresponding symbolic-model discovery. In standardized latent coordinates z1,z2,z3z_1,z_2,z_3, with r the normalized Reynolds-number coordinate, Eq. (13) gives the coefficient-optimized system that satisfied all verifier checks: z˙1 z_1 =1.143z3+0.0422r+0.201z2z3, =143z_3+0422r+201z_2z_3, (13) z˙2 z_2 =0.798z1, =798z_1, z˙3 z_3 =−0.995z1, =-995z_1, r˙ r =0. =0. The 1.143z31.143z_3 and −0.995z1-0.995z_1 terms define the leading linear oscillator in the (z1,z3)(z_1,z_3) coordinates, whose linearized angular frequency is 1.143×0.995=1.066 1.143× 0.995=1.066. The evolution of z2z_2 is coupled to the same cycle, while the quadratic contribution 0.201z2z30.201z_2z_3 feeds this coordinate back into z˙1 z_1, introducing a nonlinear deformation of the underlying oscillator. The additive 0.0422r0.0422r term shifts the latent vector field with Reynolds number, and r˙=0 r=0 preserves the selected parameter value throughout each rollout. Because z˙3=−(0.995/0.798)z˙2 z_3=-(0.995/0.798) z_2, the combination z3+1.247z2z_3+1.247z_2 is conserved to the precision of the reported coefficients. Consequently, for each fixed value of r, 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. ReRe Latent NMSE AE field (%) VG-to-AE field (%) End-to-end (%) 175 (I) 0.142 7.0 19.8 20.9 275 (I) 0.247 5.9 32.1 33.0 425 (I) 0.043 6.0 18.9 20.7 500 (E) 0.036 6.4 11.9 13.8 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 spatial L2L_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, generalization did not deteriorate monotonically from interpolation to extrapolation: the Re=500Re=500 extrapolation case produced the lowest latent NMSE (0.036) and end-to-end field error (13.8%), whereas the Re=275Re=275 interpolation case was the most difficult, with corresponding errors of 0.247 and 33.0%. Figure 4 provides important physical context for the latter value. Although displacement of the predicted vortices increases the spatial L2L_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. At Re=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. Figure 4: Simulation, autoencoder (AE) reconstruction, and VG-predicted vorticity at t∗=75t^*=75 for withheld Re=275Re=275 (interpolation, I) and Re=500Re=500 (extrapolation, E). All panels share a symmetric color scale. 6 Conclusions 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 (Marra et al. 2024; Cranmer et al. 2020), 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 T. Askham and J. N. Kutz (2018) Variable projection methods for an optimized dynamic mode decomposition. SIAM Journal on Applied Dynamical Systems 17 (1), p. 380–416. External Links: Document Cited by: §2. 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 Plasmas 27 (12), p. 123702. External Links: Document Cited by: §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 Cited by: §2. Z. Bastiani, R. M. Kirby, J. Hochhalter, and S. Zhe (2025) Diffusion-based symbolic regression. arXiv preprint arXiv:2505.24776. External Links: Link Cited by: §2. L. Biggio, T. Bendinelli, A. Neitz, A. Lucchi, and G. Parascandolo (2021) Neural symbolic regression that scales. In Proceedings of the 38th International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 139, p. 936–945. Cited by: §2. J. Brence, L. Todorovski, and S. Džeroski (2021) Probabilistic grammars for equation discovery. Knowledge-Based Systems 224, p. 107077. External Links: Document Cited by: §2. 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 Sciences 113 (15), p. 3932–3937. External Links: Document Cited by: §2. R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud (2018) Neural ordinary differential equations. In Advances in Neural Information Processing Systems, Vol. 31, p. 6571–6583. Cited by: §2. 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. In Advances in Neural Information Processing Systems, Vol. 33, p. 17429–17442. Cited by: §6. M. Cranmer (2023) Interpretable machine learning for science with PySR and SymbolicRegression.jl. arXiv preprint arXiv:2305.01582. Cited by: §2. S. d’Ascoli, S. Becker, A. Mathis, P. Schwaller, and N. Kilbertus (2024) ODEFormer: symbolic regression of dynamical systems with transformers. In The Twelfth International Conference on Learning Representations, External Links: Link Cited by: §2, §2, §3, §5.1. J. D. Eldredge (2021) ViscousFlow.jl: a framework for simulating viscous incompressible flows. Note: Zenodo External Links: Document, Link Cited by: §5.2. J. D. Eldredge (2022) A method of immersed layers on cartesian grids, with application to incompressible flows. Journal of Computational Physics 448, p. 110716. External Links: Document Cited by: §5.2. D. Etling (1990) Mesoscale vortex shedding from large islands: a comparison with laboratory experiments of rotating stratified flows. Meteorology and Atmospheric Physics 43, p. 145–151. External Links: Document Cited by: §1. F. Faraji, M. Reza, A. Knoll, and J. N. Kutz (2024) Dynamic mode decomposition for data-driven analysis and reduced-order modelling of E×BE× B plasmas: i. dynamics forecasting. Journal of Physics D: Applied Physics 57 (6), p. 065202. External Links: Document Cited by: §2. 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 Physics 137 (12), p. 123301. External Links: Document Cited by: §2, §5.2. 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 Physics 58 (10), p. 102002. External Links: Document Cited by: §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 Sciences 478 (2260), p. 20210904. External Links: Document Cited by: §2. 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 Sciences 481 (2307), p. 20240200. External Links: Document Cited by: §2. S. Greydanus, M. Dzamba, and J. Yosinski (2019) Hamiltonian neural networks. In Advances in Neural Information Processing Systems, Vol. 32, p. 15379–15389. Cited by: §2. M. Gruszecki, V. M. Nakariakov, T. Van Doorsselaere, and T. D. Arber (2010) Phenomenon of alfvénic vortex shedding. Physical Review Letters 105 (5), p. 055004. External Links: Document Cited by: §1. Japan Meteorological Agency (2019) Collection of images captured by himawari-8/9: karman vortex. Note: https://w.jma.go.jp/jma/jma-eng/satellite/introduction/image.htmlAccessed: 2026-07-22 Cited by: §1. P. Kamienny, S. d’Ascoli, G. Lample, and F. Charton (2022) End-to-end symbolic regression with transformers. In Advances in Neural Information Processing Systems, Vol. 35, p. 10269–10281. Cited by: §2. Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar (2021) Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, External Links: Link Cited by: §2. J. Loiseau and S. L. Brunton (2018) Constrained sparse galerkin regression. Journal of Fluid Mechanics 838, p. 42–67. External Links: Document Cited by: §2. 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 Computing 43 (6), p. B1105–B1132. External Links: Document Cited by: §2. G. Marra, S. Dumančić, R. Manhaeve, and L. De Raedt (2024) From statistical relational to neurosymbolic artificial intelligence: a survey. Artificial Intelligence 328, p. 104062. External Links: Document Cited by: §6. D. A. Messenger and D. M. Bortz (2021) Weak SINDy for partial differential equations. Journal of Computational Physics 443, p. 110525. External Links: Document Cited by: §2. 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 Learning 113, p. 7689–7721. External Links: Document Cited by: §2. 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. In International Conference on Learning Representations, External Links: Link Cited by: §2. 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 Physics 378, p. 686–707. External Links: Document Cited by: §2. Y. E. Şahin, N. Kilbertus, and S. Becker (2025) Predicting symbolic ODEs from multiple trajectories. In NeurIPS 2025 Workshop on Machine Learning and the Physical Sciences, External Links: Document, Link Cited by: §2. P. Shojaee, K. Meidani, A. Barati Farimani, and C. K. Reddy (2023) Transformer-based planning for symbolic regression. In Advances in Neural Information Processing Systems, Vol. 36, p. 45907–45919. External Links: Link Cited by: §2. L. Sirovich (1987) Turbulence and the dynamics of coherent structures. part i: coherent structures. Quarterly of Applied Mathematics 45 (3), p. 561–571. External Links: Document Cited by: §5.2. B. Taskin, W. Xie, and T. Lazebnik (2026) Knowledge integration for physics-informed symbolic regression using pre-trained large language models. Scientific Reports 16 (1), p. 1614. External Links: Document, Link Cited by: §2. S. Udrescu and M. Tegmark (2020) AI Feynman: a physics-inspired method for symbolic regression. Science Advances 6 (16), p. eaay2631. External Links: Document Cited by: §2. 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 Cited by: §2. Balth. van der Pol (1926) On “relaxation-oscillations”. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 2 (11), p. 978–992. External Links: Document Cited by: §5.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 Cited by: §2. C. H. K. Williamson and R. Govardhan (2004) Vortex-induced vibrations. Annual Review of Fluid Mechanics 36, p. 413–455. External Links: Document Cited by: §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, Link Cited by: §2. K. Yu, E. Chatzi, and G. Kissas (2026) Neuro-symbolic ODE discovery with latent grammar flow. arXiv preprint arXiv:2604.16232. External Links: Document, Link Cited by: §2. 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), detailed verifier definitions and their identifiability interpretation (Section S2), Van der Pol test-case details (Section S3), the fixed-Reynolds-number representation study (Section S4), the autoencoder sensitivity analysis (Section S5), and extended results from the multi-Reynolds-number test case (Section S6). Appendix S1 Sensitivity 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 reports the complete Van der Pol trajectory sets, selected equations, and final-test behavior. S1.1 Comparison 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μ=0.5 and 1.51.5, the same 12 model-construction trajectories were sampled over 0≤t≤200≤ t≤ 20. Each trajectory contributed an early window over 0≤t≤100≤ t≤ 10 and a late window over 10≤t≤2010≤ t≤ 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. 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 defines the verifier metrics and acceptance decisions, and Table S2 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.2 Validation performance and admissibility Table S1 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 μ, 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μ=1.5. The larger beam did not provide a consistent benefit: at temperature 0.1 no pre-optimization candidate was admissible for μ=0.5μ=0.5, while the post-optimization system for μ=1.5μ=1.5 no 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. Beam Temperature μ Pre-opt. candidate Pre-opt. NMSE Post-opt. NMSE Post-opt. admissible 20 0.1 0.5 Yes 0.288 7.97×−7.97× 10^-4 Yes 20 0.1 1.5 Yes 0.0145 2.74×−2.74× 10^-3 Yes 20 0.2 0.5 Yes 0.224 1.60×10−31.60× 10^-3 Yes 20 0.2 1.5 Yes 0.756 0.0255 Yes 50 0.1 0.5 No∗ 0.155∗ 0.0320∗ Yes∗ 50 0.1 1.5 Yes 0.331 0.0249 No 50 0.2 0.5 Yes 0.136 8.81×10−48.81× 10^-4 Yes 50 0.2 1.5 Yes 0.366 0.0205 Yes ∗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 and was carried forward to the cylinder-flow experiments. S1.3 Early- 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μ=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 to 7.39×10−47.39× 10^-4 and 3.19×10−43.19× 10^-4, respectively. For μ=1.5μ=1.5, the corresponding medians were 0.00537 and 0.00266 before optimization and 1.26×10−31.26× 10^-3 and 8.07×10−48.07× 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 S2 Verifier Definitions, Identifiability, and Sensitivity Considerations S2.1 Constrained model selection The verifier-guided workflow described in Section 4 of the paper separates dynamical admissibility from rollout-based ranking. Let KC_K denote the shortlist obtained after pooling, deduplicating, and ranking candidates across the discovery trajectories. For a candidate ∈Kh _K, verifier ℓ evaluates a metric mℓ()m_ (h) against a tolerance τℓ _ . An inequality-based decision can be written compactly as ℓ()=[mℓ()≤τℓ], pass_ (h)=I\! [m_ (h)≤ _ ], (S1) where [⋅]I[·] is the indicator function and τℓ _ 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:ℓ()=1for every applicable verifier ℓ,⋆=argmin∈ℒroll().A= \h _K: pass_ (h)=1\ for every applicable verifier \, = h \,L_roll(h). (S2) Equations (S1) and (S2) therefore define admissibility separately from ℒrollL_roll. Rollout error orders only the candidates that satisfy every applicable requirement. When A is empty, no admissible symbolic model is selected. S2.2 Verifier definitions and interpretation Local vector-field agreement. Let i(m)∈ℝdξ ξ_i^(m) ^d_ξ be the full state at sample i of discovery trajectory m. Its numerically estimated rate in evaluated coordinate j=1,…,dj=1,…,d is the scalar ξ˙ij(m) ξ_ij^(m), while hj(i(m))h_j( ξ_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() V_loc(h) =1d∑j=1dMSEm,i(hj(i(m)),ξ˙ij(m))Varm,i(ξ˙ij(m))+ϵ, = 1d _j=1^d MSE_m,i (h_j( ξ_i^(m)), ξ_ij^(m) )Var_m,i ( ξ_ij^(m) )+ε, (S3) loc() pass_loc(h) =[Vloc()≤2min~∈KVloc(~)], =I\! [V_loc(h)≤ 2 _ h _KV_loc( h) ], where ϵ=10−12ε=10^-12 prevents division by zero. For the directly observed and fixed-parameter systems, dξ=d_ξ=d; for the cross-ReRe system, dξ=d+1d_ξ=d+1 because the input also contains the normalized Reynolds coordinate r=(Re−300)/150r=(Re-300)/150, while the score is evaluated over the d latent-state equations. The variance normalization in Eq. (S3) 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) ξ_h^(s)(t) denote the trajectory obtained by integrating candidate h from observed starting state s, and let Bobs=maxm,i‖i(m)‖2B_obs= _m,i\| ξ_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≤cboundBobs]. pass_bound(h)=I\! [ ξ_h^(s)(t) is finite for all s,t\ \ _s,t \| ξ_h^(s)(t) \|_2≤ c_boundB_obs ]. (S4) In Eq. (S4), the norm and its reference scale are evaluated only over the latent coordinates and separately at each Reynolds number for the cross-ReRe case. The factor is cbound=1.5c_bound=1.5 for Van der Pol and cbound=3c_bound=3 for 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 coordinate j is estimated by the robust percentile half-range Aj()=12[Q0.99(ξj)−Q0.01(ξj)]A_j( ξ)= 12[Q_0.99( _j)-Q_0.01( _j)]. Let Aj(m)A_j^(m) be the amplitude of reference trajectory m, and let A^s,j(m) A_s,j^(m) be the corresponding amplitude of the candidate rollout from starting state s. If ℳvM_v indexes the reference cases used for verification and AJ_A the coordinates whose amplitudes are evaluated, then mamp() m_amp(h) =1|ℳv||A|∑m∈ℳv∑j∈A|mediansA^s,j(m)−Aj(m)||Aj(m)|+ϵ, = 1|M_v|\,|J_A| _m _v _j _A |median_s A_s,j^(m)-A_j^(m) ||A_j^(m)|+ε, (S5) amp() pass_amp(h) =[mamp()≤0.15]. =I\! [m_amp(h)≤ 0.15 ]. For Van der Pol, AJ_A contains the observed oscillatory coordinate, and ℳvM_v contains one pooled reference case whose target is formed from the discovery trajectories. For the fixed-ReRe flow case, ℳvM_v contains the single Re=300Re=300 reference trajectory; for the cross-ReRe case, it contains the three reference trajectories at Re=150,300,Re=150,300, and 450450. All retained coordinates belong to AJ_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). Period or dominant frequency. Let qj(m)q_j^(m) denote the characteristic temporal scale of coordinate j in reference case m, and let q^s,j(m) q_s,j^(m) denote the same quantity from the candidate rollout. Its relative discrepancy is mtime() m_time(h) =1|ℳv||q|∑m∈ℳv∑j∈q|mediansq^s,j(m)−qj(m)||qj(m)|+ϵ, = 1|M_v|\,|J_q| _m _v _j _q |median_s q_s,j^(m)-q_j^(m) ||q_j^(m)|+ε, (S6) time() pass_time(h) =[mtime()≤0.10]. =I\! [m_time(h)≤ 0.10 ]. For Van der Pol, q is the median interval between successive peaks after the transient and qJ_q contains the observed oscillatory coordinate. For each reduced flow coordinate, qjq_j is instead the dominant nonzero Fourier frequency, qj=argmaxν>0|ℱ[ξj−ξ¯j](ν)|2q_j= *arg\,max_ν>0|F[ _j- ξ_j](ν)|^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) tests preservation of the dominant oscillation rate. Starting-state consistency. Each candidate is integrated from three observed starting states. For a rollout property u∈A,qu∈\A,q\, where A is amplitude and q is period or dominant frequency, define Cu() C_u(h) =1|ℳv||u|∑m∈ℳv∑j∈ustds(u^s,j(m))means|u^s,j(m)|+ϵ, = 1|M_v|\,|J_u| _m _v _j _u std_s ( u_s,j^(m) )mean_s | u_s,j^(m) |+ε, (S7) cons() pass_cons(h) =[CA()≤0.10∧Cq()≤0.05]. =I\! [C_A(h)≤ 0.10\ \ C_q(h)≤ 0.05 ]. 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) excludes candidates for which those long-time properties depend strongly on the selected starting state. Parameter conservation and participation. For the cross-ReRe model, let hrh_r denote the candidate equation for the normalized Reynolds coordinate defined above and let hjh_j, j=1,…,dzj=1,…,d_z, denote the equations for the dzd_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]. pass_par(h)=I\! [h_r≡ 0\ \ _s,t|\, r_h^(s)(t)-r^(s)(0)\,|≤ 10^-8\ \ ∃\,j≤ d_z: ∂ h_j∂ r ≡ 0 ]. (S8) The first condition preserves the autonomous augmentation r˙=0 r=0 symbolically, and the second confirms the same conservation in numerical rollouts. The third requires the latent dynamics to depend explicitly on r; carrying an unused constant coordinate is not sufficient for a shared cross-parameter model. Table S2 summarizes the settings associated with Eqs. (S3)–(S8). The cross-ReRe 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. For the cross-ReRe integration, one model-time unit corresponds to 0.10.1 physical-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 were held constant within each experimental setting. Their sensitivity is examined for the representative fixed-ReRe case in Section S2.7. Table S2: Verifier settings used in the three experimental settings. Setting Van der Pol Fixed Re=300Re=300 Cross-ReRe Shortlist size 10 10 17∗17^* Local vector-field factor 2.0 2.0 2.0 Bound factor 1.5 3.0 3.0 Verification horizon 60 30 30∗30^** Discard interval 30 5 5∗5^** Amplitude tolerance 0.15 0.15 0.15 Period/frequency tolerance 0.10 0.10 0.10 Amplitude CV tolerance 0.10 0.10 0.10 Period/frequency CV tolerance 0.05 0.05 0.05 Parameter-drift tolerance NA NA 10−810^-8 ∗The two cross-ReRe shortlist branches yielded 17 distinct candidates. ∗Cross-ReRe horizon and discard values are reported in physical-time units. S2.3 Identifiability 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 ℋGH_G be the generated, deduplicated hypothesis class, R()R(f) its multi-trajectory rollout discrepancy, and ε a chosen rollout tolerance. The set of candidates that are indistinguishable at this rollout resolution is ℰR(ε)=∈ℋG:R()≤ε.E_R( )= \f _G:R(f)≤ \. (S9) Equation (S9) describes equivalence only with respect to the sampled rollout discrepancy. Let ℋA⊆ℋGH_A _G denote the candidates satisfying all applicable dynamical and physical-admissibility conditions. Verification then restricts the rollout-consistent set to ℰVG(ε)=ℰR(ε)∩ℋA.E_VG( )=E_R( ) _A. (S10) The intersection in Eq. (S10) 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.4 Local 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θ θ ^n_θ contain its nθn_θ free coefficients, ()y( θ) collect the sampled candidate rollouts, and ()m( θ) 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 VG()=[∂/∂/∂].J_VG( θ)= bmatrix /∂ θ\\ /∂ θ bmatrix. (S11) Each column of VGJ_VG in Eq. (S11) describes how the measured quantities change under a perturbation of one coefficient. Appending verifier sensitivities cannot reduce the rank available from the rollout samples: rankVG≥rank(∂/∂).rankJ_VG ( /∂ θ ). (S12) A strict inequality in Eq. (S12) 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.5 Implications 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 ℰϕE_ φ be an encoder with parameters ϕ φ, ϕ=ℰϕ()z_ φ=E_ φ( ω) its latent state, and g_ θ a symbolic latent vector field with coefficients θ. The sensitivity of reconstruction, latent rollout, and verifier quantities to both parameter sets could be represented by joint=∂(ϕ,)[ℛϕ(ϕ)^ϕ,(t)(ϕ,)],J_joint= ∂( φ, θ) bmatrixR_ φ(z_ φ)\\ z_ φ, θ(t)\\ m( φ, θ) bmatrix, (S13) where ℛϕR_ φ is the reconstruction map, ^ϕ,(t) z_ φ, θ(t) is the integrated latent trajectory, and m contains differentiable surrogates of the verifier properties. Equation (S13) 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.6 Parameter dependence in the cross-ReRe model In the cross-ReRe experiment, each decoding is performed for one trajectory whose appended normalized Reynolds coordinate r is constant. A single trajectory therefore cannot identify how the latent dynamics vary with Reynolds number. For example, on a trajectory with r=r0≠0r=r_0≠ 0, the terms αzj=αr0rzjwhen r=r0α z_j= αr_0\,rz_j r=r_0 (S14) are observationally indistinguishable. The occurrence of r in 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) is resolved only by evaluating a shared candidate across different values of r. Such terms can nevertheless arise because r is 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) additionally requires conservation of r and its participation in at least one latent equation. Parameter dependence in the selected model consequently emerges from candidate generation followed by cross-ReRe pooling, evaluation, and verification, rather than from simultaneous observation of all trajectories during an individual decode. S2.7 Sensitivity to verifier tolerances The numerical tolerances in Table S2 determine the boundary of the admissible set and therefore warrant a sensitivity check. We use the fixed-Re=300Re=300 cylinder 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-ReRe model. The same perturbation procedure can be applied to the Van der Pol and cross-ReRe 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 provides the full representation and selection comparison. Let γ∈0.75,1.0,1.5γ∈\0.75,1.0,1.5\ jointly scale the local vector-field agreement factor in Eq. (S3) and the amplitude, frequency, amplitude-CV, and frequency-CV tolerances. Thus γ<1γ<1 imposes stricter quantitative agreement, γ=1γ=1 recovers the nominal settings in Table S2, and γ>1γ>1 relaxes 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)γ=(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, 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-ReRe 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 S3 Van der Pol Test-Case Details and Extended Results Section S1 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.1 Trajectory Sets and Evaluation Roles For both μ=0.5μ=0.5 and 1.51.5, trajectories were integrated over 0≤t≤200≤ t≤ 20 at 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 j x0(0)x_0(0) x1(0)x_1(0) Candidate source 0 0.023643 0.900927 – 1 -0.711681 0.897299 late 2 -0.376337 -0.153347 early, late 3 0.655405 -0.181602 – 4 0.099187 -0.944882 – 5 0.507026 0.076287 late 6 -0.340537 0.576857 late 7 -0.393610 -0.093004 – 8 -0.731917 -0.193774 early, late 9 -0.593090 -0.475373 early 10 0.500729 -0.439182 – 11 -0.029618 0.961474 – Early window: 0≤t≤100≤ t≤ 10 (151 samples). Late window: 10≤t≤2010≤ t≤ 20 (151 samples). Beam–temperature validation j x0(0)x_0(0) x1(0)x_1(0) 0 -0.476776 -0.403018 1 0.628451 -0.816168 2 0.200201 0.457121 3 -0.624198 -0.889707 4 -0.450061 0.314866 5 0.124531 -0.699875 6 -0.134738 0.338595 7 -0.154431 0.266369 Final test j x0(0)x_0(0) x1(0)x_1(0) 0 -0.828702 -0.526379 1 0.602549 0.164324 2 -0.811743 -0.133746 3 -0.041897 -0.680522 4 0.469154 -0.772656 5 -0.217544 0.033480 6 -0.138744 0.173597 7 0.475676 0.912535 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. 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.2 Selected equations and verifier outcomes The canonical system has x˙0=x1 x_0=x_1 and x˙1=μx1−x0−μx02x1 x_1=μ x_1-x_0-μ x_0^2x_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 μ=0.5. Pre-optimization:x˙0 -optimization: x_0 =1.0239x1−0.0548x0, =0239x_1-0548x_0, (S15) x˙1 x_1 =0.2763x1−1.0463x0−0.0769x1(0.0736x0+3.6009x02). =2763x_1-0463x_0-0769x_1\! (0.0736x_0+3.6009x_0^2 ). Post-optimization:x˙0 -optimization: x_0 =0.99726x1−0.012400x0, =99726x_1-012400x_0, (S16) x˙1 x_1 =0.54166x1−0.99708x0−0.13251x1(0.029451x0+4.0969x02). =54166x_1-99708x_0-13251x_1\! (0.029451x_0+4.0969x_0^2 ). =1.5 μ=1.5. Pre-optimization:x˙0 -optimization: x_0 =0.0007(1+1.2715x0)2+0.9889x1, =0007 (1+1.2715x_0 )^2+9889x_1, (S17) x˙1 x_1 =1.4616x1−0.9972x0−0.1104x1(0.0690x0+12.5443x02). =4616x_1-9972x_0-1104x_1\! (0.0690x_0+12.5443x_0^2 ). Post-optimization:x˙0 -optimization: x_0 =0.00067573(0.99818+1.2726x0)2+0.99637x1, =00067573 (0.99818+1.2726x_0 )^2+99637x_1, (S18) x˙1 x_1 =1.43998x1−1.00151x0−0.11351x1(0.070367x0+12.8741x02). =43998x_1-00151x_0-11351x_1\! (0.070367x_0+12.8741x_0^2 ). 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 μ, 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 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. Table S3: Post-optimization verifier outcomes for the selected Van der Pol equations. Percentages are reported relative to the corresponding ground-truth oscillation statistic. Verifier Acceptance requirement μ=0.5μ=0.5 μ=1.5μ=1.5 Finite integration All three long-horizon rollouts complete 3/3 3/3 Boundedness Maximum state norm no greater than 1.51.5 times the observed maximum Pass Pass Local vector field Score no greater than the case-specific limit 8.72×10−48.72×10^-4 (limit 4.02×10−24.02×10^-2) 2.11×10−42.11×10^-4 (limit 4.60×10−34.60×10^-3) Amplitude error No greater than 15% 0.741% 0.678% Period error No greater than 10% 0.263% 0.236% Amplitude CV No greater than 10% 0.057% 0.016% Period CV No greater than 5% 0.370% Below numerical resolution Overall outcome Every applicable requirement satisfied Pass Pass The local vector-field scores are 2.17% and 4.59% of their respective limits for μ=0.5μ=0.5 and 1.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 separately assesses predictive generalization from initial conditions excluded from model construction. S3.3 Representative held-out rollouts Across the eight final-test initial conditions, coefficient optimization reduced the median VG rollout NMSE from 0.287 to 6.99×10−46.99× 10^-4 for μ=0.5μ=0.5, and from 0.0139 to 3.27×10−33.27× 10^-3 for μ=1.5μ=1.5. These values differ from the validation errors in Table S1: 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 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)x(0)=(-0.811743,-0.133746), for μ=0.5μ=0.5, and test trajectory 3, with (0)=(−0.041897,−0.680522)x(0)=(-0.041897,-0.680522), for μ=1.5μ=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. 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μ=0.5, both original ODEFormer equations begin to separate markedly from the ground truth after approximately t=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μ=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, 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 S4 Fixed-Reynolds-Number Representation and Selection Details This section expands the fixed-Re=300Re=300 experiment reported in Section 5.2. 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, S4.2, and S4.3 respectively examine the representation, configuration selection, and final-test behavior. S4.1 POD spectrum and mode-pair diagnostics POD was applied to fluctuations about the development-interval mean, rather than to the total vorticity fields. If (ti) ω(t_i) is the vectorized field at development time tit_i and ¯D ω_\!D is their temporal mean, the fluctuation snapshot matrix and its singular-value decomposition are ′ =[(t1)−¯D⋯(tND)−¯D]=, = bmatrix ω(t_1)- ω_\!D&·s& ω(t_N_D)- ω_\!D bmatrix=U V T, (S19) Er E_r =∑j=1rσj2∑j=1Nsσj2. = _j=1^r _j^2 _j=1^N_s _j^2. Here ND=250N_D=250, σj _j is the jjth singular value, NsN_s is the number of nonzero singular values, and ErE_r is the fraction of fluctuation energy represented by the first r modes. The mean field is restored only when the reduced coordinates are reconstructed in physical space. The near-equal singular values in Figure S2 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. Figure S2: POD singular-value spectrum and cumulative fluctuation energy at fixed Re=300Re=300. For a consecutive pair (2k−1,2k)(2k-1,2k), we quantify its relative energy imbalance by δk=|σ2k−12−σ2k2|σ2k−12+σ2k2. _k= | _2k-1^2- _2k^2| _2k-1^2+ _2k^2. (S20) The pair-energy column in Table S4 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 _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 therefore complements the cumulative spectrum by separating the energy contributed by each oscillatory pair. Table S4: POD mode-pair energy and imbalance diagnostics. Modes Pair energy (%) Cumulative energy (%) δk _k 1–2 73.44 73.44 0.1309 3–4 12.17 85.61 0.0379 5–6 8.91 94.52 0.0097 S4.2 POD-rank and decoding-window selection matrix The representation and decoding choices motivated in Section S4.1 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 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-ReRe 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. Rank Windows Adm. Window Continuous POD VG–POD End-to-end 4 8 10 1.30 1.54 20.87 4.57 21.39 4 24 10 1.11 2.03 20.87 5.88 21.75 6 8 6 4.36 3.95 12.70 7.57 14.93 6 24 10 4.61 7.11 12.70 12.84 18.28 Table S5 reports the admissible counts obtained at the nominal verifier settings. Section S2.7 examines how these counts change when the quantitative agreement tolerances are tightened or relaxed. For the selected rank-six, eight-window configuration, Eq. (S21) below is the pre-optimization equation selected at the nominal settings, and Eq. (S22) is the next-lowest-rollout admissible alternative. These are the two equations referenced in Section S2.7. Although the number of admissible rank-six candidates changes across the tested tolerance factors, Eq. (S21) remains selected in every case. Let =(z1,…,z6)z=(z_1,…,z_6) T denote the standardized POD coordinates, consistent with the notation used for the post-optimization system in Eq. (12). The pre-optimization system with the lowest development- window rollout error, whose structure was advanced to coefficient optimization, was z˙1 z_1 =1.0044z2−0.0700(12.3400−0.9594z2)−1, =0044z_2-0700(23400-9594z_2)^-1, z˙2 z_2 =−1.1105z1, =-1105z_1, (S21) z˙3 z_3 =−2.1602z4, =-1602z_4, z˙4 z_4 =2.1867z3, =1867z_3, z˙5 z_5 =3.4432z6−0.1382sin(0.1052+9.1998z2), =4432z_6-1382 (1052+1998z_2), z˙6 z_6 =−0.1399−2.8694z5. =-1399-8694z_5. The next-lowest-rollout admissible system provides a representative alternative: z˙1 z_1 =1.1208z2, =1208z_2, z˙2 z_2 =−1.0604z1, =-0604z_1, (S22) z˙3 z_3 =−2.2995z4, =-2995z_4, z˙4 z_4 =2.0672z3, =0672z_3, z˙5 z_5 =3.1395z6−0.2202sin(11.8700+148.5203z6), =1395z_6-2202 (18700+485203z_6), z˙6 z_6 =−0.1463−3.1228z5. =-1463-1228z_5. Both systems preserve three oscillator pairs and satisfy all applicable verifiers, but their nonlinear structures differ. In Eq. (S21), −0.0700(12.3400−0.9594z2)−1-0.0700(12.3400-0.9594z_2)^-1 adds a reciprocal correction to z˙1 z_1 in the first pair, while −0.1382sin(0.1052+9.1998z2)-0.1382 (0.1052+9.1998z_2) couples the first pair to z˙5 z_5 in the third pair. Equation (S22) leaves the first pair linear and instead introduces −0.2202sin(11.8700+148.5203z6)-0.2202 (11.8700+148.5203z_6) as a nonlinear correction within the third pair. Both systems also contain a small constant offset in z˙6 z_6. At the nominal settings, six of the ten shortlisted systems for this configuration satisfy every verifier. The coexistence of Eqs. (S21) and (S22) therefore illustrates the practical non-uniqueness discussed in Section S2.3: 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). 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). S4.3 Extended coordinate and field comparisons The final-test interval was evaluated only after the rank-six, eight-window configuration had been chosen. In Figures S3 and S4, aja_j denotes the plotted standardized POD coefficient and corresponds to zjz_j in Eqs. (S21) and (S22). 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(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. Figure S3: Final-test fixed-ReRe rollouts in the six standardized POD coordinates. Field error was decomposed into three contributions. For the simulated vorticity field (t) ω(t), its rank-six POD projection POD(t) ω_POD(t), and the field reconstructed from the symbolic rollout ^VG(t) ω_VG(t), the instantaneous relative errors are ePOD(t) e_POD(t) =‖POD(t)−(t)‖2‖(t)‖2, = \| ω_POD(t)- ω(t)\|_2\| ω(t)\|_2, (S23) eVG-POD(t) e_VG -POD(t) =‖^VG(t)−POD(t)‖2‖POD(t)‖2, = \| ω_VG(t)- ω_POD(t)\|_2\| ω_POD(t)\|_2, eend(t) e_end(t) =‖^VG(t)−(t)‖2‖(t)‖2. = \| ω_VG(t)- ω(t)\|_2\| ω(t)\|_2. Their time-averaged values are 12.53%, 8.77%, and 15.56%, respectively. Figure S4(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. Figure S4: Final-test fixed-ReRe 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. 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. Figure S5: Simulation, rank-six POD reconstruction, and VG prediction at three times in the final-test fixed-ReRe interval. Appendix S5 Autoencoder Representation–Discoverability Sensitivity S5.1 Controlled 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. The independently prepared autoencoder used for the cross-ReRe result in Section 5.3 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. 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 use these post-transient fields. One predeclared seed per architecture was carried forward to the symbolic discovery comparison. S5.2 Definition of the latent-representation diagnostics Three diagnostics characterize the temporal properties of each embedding before symbolic model construction. Let (m)∈ℝT×dzZ^(m) ^T× d_z contain the latent trajectory at the mmth Reynolds number, where T is the number of snapshots and dzd_z the latent dimension. Each coordinate is standardized using its mean μj _j and standard deviation σj _j on the concatenated trajectories: Z~t,j(m)=(Zt,j(m)−μj)/(σj+ϵ) Z_t,j^(m)=(Z_t,j^(m)- _j)/( _j+ε). With Δ and Δ2 ^2 denoting first and second temporal differences, respectively, the latent roughness is ℛz=1MR∑m=1MRmeant,j[(Δ2Z~t,j(m))2]meant,j[(ΔZ~t,j(m))2]+ϵ,R_z= 1M_R _m=1^M_R mean_t,j [( ^2 Z_t,j^(m))^2 ]mean_t,j [( Z_t,j^(m))^2 ]+ε, (S24) where MR=7M_R=7 is the number of Reynolds numbers. The denominator normalizes by the overall first-difference scale, so Eq. (S24) measures rapid changes in temporal slope rather than latent amplitude. Smaller values indicate smoother trajectories. For spectral concentration, let Pk,j(m)P_k,j^(m) be the discrete Fourier power of the centered standardized coordinate at positive, nonzero frequency bin k, and let 3(m,j)K_3^(m,j) contain its three largest-power bins. The reported diagnostic is z=1MRdz∑m=1MR∑j=1dz∑k∈3(m,j)Pk,j(m)∑k>0Pk,j(m)+ϵ.S_z= 1M_Rd_z _m=1^M_R _j=1^d_z _k _3^(m,j)P_k,j^(m) _k>0P_k,j^(m)+ε. (S25) Values of zS_z in Eq. (S25) 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 allZ_all and let λ1,…,λdz≥0 _1,…, _d_z≥ 0 be the eigenvalues of its sample covariance matrix. The covariance participation ratio is dPR=(∑j=1dzλj)2∑j=1dzλj2+ϵ.d_PR= ( _j=1^d_z _j )^2 _j=1^d_z _j^2+ε. (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)–(S26) describe complementary properties: temporal regularity, spectral organization, and covariance use of the available latent dimensions. S5.3 Latent-representation diagnostics evaluation Table S6 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 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 spatial L2L_2 error; the latent diagnostics are defined in Eqs. (S24)–(S26). Configuration Seed Field error (%) Roughness Spectral conc. Participation AE3S 17 5.636 0.0472 0.8830 1.975 AE3S 29 6.349 0.0543 0.8841 1.910 AE3S 43∗43^* 5.637 0.0332 0.8970 2.447 AE4S 17∗17^* 5.379 0.1145 0.7842 1.733 AE4S 29 5.825 0.1082 0.8065 2.664 AE4S 43 5.021 0.1143 0.7877 2.854 AE3D 17∗17^* 4.838 0.0347 0.9232 2.291 AE3D 29 5.644 0.0383 0.9087 2.370 AE3D 43 4.494 0.0465 0.9098 2.562 AE4D 17 5.091 0.0993 0.8646 3.677 AE4D 29∗29^* 4.625 0.0935 0.8753 3.204 AE4D 43 4.650 0.0733 0.8759 3.650 Independent shallow reference∗ – 5.826 0.0224 0.9264 2.149 ∗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. 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.4 Corresponding 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, 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 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. Table S7: Symbolic-discovery outcomes for the representations advanced from the controlled study and for the independent shallow reference. Representation Pre-opt. admissible Pre-opt. loss Post-opt. result AE3S, seed 43 No – Inadmissible AE4S, seed 17 No – Inadmissible AE3D, seed 17 Yes 2.169 Inadmissible AE4D, seed 29 No – Inadmissible Independent shallow reference Yes 1.042 0.200, admissible 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, in which dynamical verification can help distinguish latent embeddings with similar reconstruction fidelity but different compatibility with symbolic discovery. Appendix S6 Extended Results for the Cross-Parameter Cylinder-Flow Case S6.1 Representation and evaluation protocol The cross-ReRe study used the same shallow autoencoder as the result reported in Section 5.3; Section S5 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.Re∈\150,200,250,300,350,400,450\. The Reynolds number was appended as r=(Re−300)/150r=(Re-300)/150. The symbolic backbone was given the uniformly sampled trajectories at unit time increments, so its model time is τ=t∗/0.1τ=t^*/0.1, where consecutive physical snapshots are separated by Δt∗=0.1 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) ξ=(z_1,z_2,z_3,r) 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 approximately 50≤t∗≤10050≤ t^*≤ 100. The withheld cases comprise interpolation at Re=175,275,Re=175,275, and 425425 and extrapolation at Re=500Re=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 at Re=500Re=500, but the corresponding latent trajectory was excluded from symbolic model construction. Thus, “extrapolation” at Re=500Re=500 refers specifically to evaluating the Reynolds-conditioned equation at r=4/3r=4/3, beyond its construction interval −1≤r≤1-1≤ r≤ 1; it does not refer to extrapolation of the field representation. Candidates generated from these seven trajectories were verified using the cross-ReRe settings in Table S2. 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) and (S28) denote differentiation with respect to the model time τ. The pre-optimization and post-optimization systems are z˙1 z_1 =0.1070z3−0.0013r+0.0089z2z3, =0.1070z_3-0.0013r+0.0089z_2z_3, z˙2 z_2 =0.1003z1, =0.1003z_1, z˙3 z_3 =−0.1098z1, =-0.1098z_1, r˙ r =0, =0, (S27) z˙1 z_1 =0.114317z3+0.004217r+0.020148z2z3, =0.114317z_3+0.004217r+0.020148z_2z_3, z˙2 z_2 =0.079842z1, =0.079842z_1, z˙3 z_3 =−0.099544z1, =-0.099544z_1, r˙ r =0. =0. (S28) Since one unit of τ corresponds to 0.10.1 units of t∗t^*, multiplying the right-hand sides of Eq. (S28) by ten gives the physical-time equation reported in Eq. (13). S6.2 Latent-coordinate rollouts at withheld Reynolds numbers Figure S7 compares the latent trajectories obtained by encoding the simulated fields with the pre-optimization and post-optimization VG rollouts. We retain zjz_j for these nonlinear autoencoder coordinates, consistent with the notation in Section 5.3 and distinct from the POD coefficients aja_j in Figures S3 and S4. For an encoded trajectory (tn)z(t_n) and its prediction ^(tn) 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)]+ϵ,E_z= 13 _j=1^3 N^-1 _n=1^N ( z_j(t_n)-z_j(t_n) )^2Var_n\! [z_j(t_n) ]+ε, (S29) where N=500N=500 and ϵε is a small numerical constant. This metric weights each coordinate relative to its own variation; the conserved parameter coordinate is not included. Coefficient optimization reduces ℰzE_z at every withheld Reynolds number: • Re=175Re=175: 2.312→0.1422.312→ 0.142, a 93.9% reduction; • Re=275Re=275: 0.797→0.2470.797→ 0.247, a 69.0% reduction; • Re=425Re=425: 0.188→0.0430.188→ 0.043, a 77.2% reduction; • Re=500Re=500: 0.714→0.0360.714→ 0.036, a 94.9% reduction. All post-optimization rollouts remain finite and oscillatory over the complete interval. The largest remaining error occurs at Re=275Re=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), r enters explicitly through an additive term in z˙1 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 with ReRe. The strong result at Re=500Re=500 indicates that the post-optimization equation continues the same periodic shedding dynamics beyond the parameter range used to construct it. 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.3 Spatially 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.W= \(x/D,y/D):0.5≤ x/D≤ 5,\;|y/D|≤ 1.5 \. For a field ω∗(,t∗)=ωD/U∞ω^*(x,t^*)=ω D/U_∞, its spatial average is ⟨ω∗⟩(t∗)=1||∫ω∗(,t∗)d. ω^* _W(t^*)= 1|W| _Wω^*(x,t^*)\,dx. (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. 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 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 at Re=425Re=425 and 500500. At Re=175Re=175 the predicted phase gradually separates from the reconstructed field, while Re=275Re=275 exhibits 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.4 Field reconstruction and error decomposition The field comparison uses the same three-part error decomposition as Eq. (S23). Here the POD reconstruction POD ω_POD is replaced by the autoencoder reconstruction AE ω_AE. Accordingly, eAEe_AE and eVG-AEe_VG -AE replace ePODe_POD and eVG-PODe_VG -POD, respectively, while the end-to-end definition is unchanged. Table S8 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-ReRe decomposition, the end-to-end error is not the arithmetic sum of its two components. Table S8: Extended cross-ReRe latent and field errors. Field entries give the temporal mean with the median in parentheses, in percent. ReRe Regime Latent NMSE AE error VG–AE error End-to-end error Pre-opt. Post-opt. Mean (median) Mean (median) Mean (median) 175 Interpolation 2.312 0.142 6.96 (6.81) 19.85 (18.30) 20.86 (19.13) 275 Interpolation 0.797 0.247 5.90 (5.83) 32.07 (31.44) 33.01 (32.50) 425 Interpolation 0.188 0.043 6.00 (5.98) 18.88 (19.97) 20.74 (21.72) 500 Extrapolation 0.714 0.036 6.41 (6.46) 11.89 (13.10) 13.81 (14.26) Figure S9 consolidates all four withheld Reynolds numbers at t∗=75t^*=75, including the Re=275Re=275 and 500500 cases 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 at Re=275Re=275 is consistent with the accumulated phase difference in Figures S7 and S8; the closer latent rollouts at Re=425Re=425 and 500500 correspond to more closely aligned large-scale fields. As discussed in Section S6.2, 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. Figure S9: Simulation, autoencoder reconstruction, and VG prediction at t∗=75t^*=75 for all Reynolds numbers excluded from symbolic model construction. The first three rows are interpolation cases; the final row is extrapolation of the symbolic dynamics.