Paper deep dive
Bi-FORK: Generative Modeling of High-Dimensional Bifurcating Systems
Anna Zimmel, Fleur Hendriks, Markus Holzleitner, Florian Sestak, Martin Weichselbaumer, Vlado Menkovski, Johannes Brandstetter
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 10/11/2026, 2:26:08 AM
Summary
The paper introduces Bi-FORK, a generative deep learning framework designed to model high-dimensional bifurcating physical systems where a single input condition leads to multiple valid solution branches. Unlike traditional surrogates that average solutions, Bi-FORK uses latent flow matching to generate complete, coherent trajectories and employs repulsion-guided sampling to recover distinct solution modes in a single pass. The method is evaluated on buckling beams, mechanical metamaterials, and Allen-Cahn phase separation, demonstrating scalability to systems with up to 260,000 degrees of freedom.
Entities (8)
Relation Signals (7)
Bi-FORK → uses → Latent Flow Matching
confidence 95% · Bi-FORK generates complete trajectories through latent flow matching
Bi-FORK → uses → Repulsion-Guided Sampling
confidence 95% · uses repulsion-guided sampling to recover distinct solution branches
Bi-FORK → evaluatedon → Beam3D
confidence 90% · We evaluate Bi-FORK on buckling beams... Beam3D
Bi-FORK → evaluatedon → Mechanical Metamaterials
confidence 90% · We evaluate Bi-FORK on... mechanical metamaterials
Bi-FORK → evaluatedon → Allen-Cahn
confidence 90% · Allen-Cahn describes how a mixture separates... We evaluate Bi-FORK on... Allen-Cahn
Bi-FORK → usesmetric → Jensen-Shannon Divergence
confidence 90% · measuring orbit coverage... via the Jensen-Shannon divergence
Bi-FORK → usesmetric → MCon-MAE
confidence 90% · Mode Conditioned MAE (MCon-MAE) quantifies trajectory-level precision... We adopt two complementary metrics
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Bifurcations are ubiquitous in physical systems, from structural buckling to fluid and climate dynamics, yet they remain largely unexplored in deep learning. At a symmetry-breaking bifurcation, a single input admits multiple equally valid solutions, violating the one-to-one assumption underlying most learned physical surrogates. We introduce Bi-FORK, a generative framework for learning these one-to-many solution maps in high-dimensional systems. Bi-FORK generates complete trajectories through latent flow matching, preserving space and time coherence, and uses repulsion-guided sampling to recover distinct solution branches in a single amortized pass. We evaluate Bi-FORK on buckling beams, mechanical metamaterials, and Allen-Cahn phase separation, spanning continuous, discrete, and field-valued bifurcations with discretizations up to 260,000 points. Bi-FORK recovers the multimodal solution structure while scaling several orders of magnitude beyond prior approaches, opening generative modeling to high-dimensional bifurcating physical systems.
Tags
Links
- Source: https://arxiv.org/abs/2610.12449v1
- Canonical: https://arxiv.org/abs/2610.12449v1
Trouble viewing inline? Open PDF directly →
Full Text
101,956 characters extracted from source content.
Expand or collapse full text
Bi-FORK: Generative Modeling of High-Dimensional Bifurcating Systems Anna Zimmel Fleur Hendriks Markus Holzleitner Florian Sestak Affiliation: ELLIS Unit, LIT AI Lab, Institute for Machine Learning, JKU Linz Affiliation: Department of Mathematics and Computer Science, Eindhoven University of Technology Email: zimmel@ml.jku.at Martin Weichselbaumer Vlado Menkovski Johannes Brandstetter* Equal contribution as first authors, † Equal contribution as last authors Affiliation: ELLIS Unit, LIT AI Lab, Institute for Machine Learning, JKU Linz Affiliation: Department of Mathematics and Computer Science, Eindhoven University of Technology Affiliation: Mistral, Paris Email: f.hendriks@differ.nl Affiliation: Department of Mechanical Engineering, Eindhoven University of Technology Affiliation: DIFFER – Dutch Institute for Fundamental Energy Research Abstract Bifurcations are ubiquitous in physical systems, from structural buckling to fluid and climate dynamics, yet they remain largely unexplored in deep learning. At a symmetry-breaking bifurcation, a single input admits multiple equally valid solutions, violating the one-to-one assumption underlying most learned physical surrogates. We introduce Bi-FORK, a generative framework for learning these one-to-many solution maps in high-dimensional systems. Bi-FORK generates complete trajectories through latent flow matching, preserving space and time coherence, and uses repulsion-guided sampling to recover distinct solution branches in a single amortized pass. We evaluate Bi-FORK on buckling beams, mechanical metamaterials, and Allen-Cahn phase separation, spanning continuous, discrete, and field-valued bifurcations with discretizations up to 260,000 points. Bi-FORK recovers the multimodal solution structure while scaling several orders of magnitude beyond prior approaches, opening generative modeling to high-dimensional bifurcating physical systems. 1 Introduction Bifurcations are fundamental to nonlinear physical systems and occur across structural mechanics, fluid dynamics, climate dynamics, and phase-transition phenomena [49, 54, 47]. A bifurcation is a qualitative change in the number or stability of solutions as a control parameter varies. When a bifurcation breaks a symmetry, it produces a one-to-many mapping from control parameter to system response: multiple, equally valid solutions rather than one. For example, a compressed column is stable up to a critical load, then bifurcates into symmetric buckled branches [7]. Classical numerical methods, including finite-element and finite-volume methods, can compute these solutions accurately. However, for high-dimensional systems, recovering the full bifurcation structure requires stability and mode analyses together with repeated continuation or deflation solves [10, 18], which is computationally expensive. Therefore, machine learning (ML) methods are an attractive alternative as they amortize the cost of building the training set over the many fast evaluations of unseen system conditions [1, 4, 3, 20, 55]. However, because the system is one-to-many, training discriminative ML models averages over the solution space and, as a consequence, produces nonphysical predictions [50, 36, 46]. This has motivated the development of probabilistic methods based on GANs, diffusion models, and more recently flow matching, typically framed around uncertainty or the loss of fine-scale structure [19, 60, 12]. Most probabilistic physical surrogates model epistemic or unresolved uncertainty; symmetry-breaking bifurcations instead produce genuinely multiple, equally valid solutions structured by the system’s symmetry. Only a few recent approaches treat the bifurcation structure directly as the cause of solution multiplicity [17, 44, 61, 30, 25, 48]. Nevertheless, these approaches recover systems of at most a few hundred degrees of freedom with small, known, discrete symmetry groups. Figure 1: Learning bifurcating solution distributions. Left: Symmetry breaking results in multiple valid trajectories for the same control parameter and initial state. Center: Latent flow matching jointly samples trajectories and resolves distinct branches in one forward pass. Right: Unlike mode-averaging baselines, Bi-FORK recovers the target distribution. The main challenge Scaling to larger systems is challenging because recovering a bifurcating distribution requires every part of the system to agree on which branch it is realizing in both space and time. A beam buckling to the left at one end and right at the other is an invalid solution, as is a beam that first buckles to the left and then suddenly switches to buckling to the right halfway through a trajectory. Moreover, bifurcations are typically rare events in trajectories. Training a frame-by-frame autoregressive model to produce trajectories would face significant challenge of capturing this key aspect of the system’s behavior while being underrepresented in the training data. Predicting the full trajectory allows the model to instead represent the branch choice globally in latent space, for all steps and elements at the same time. Mesh-based GNN surrogates communicate through local message passing, making such global coordination increasingly difficult as systems grow due to oversmoothing and oversquashing [43, 13, 11, 5]. Evaluation is similarly nontrivial. Pointwise metrics such as rank histograms and CRPS [23, 19, 28, 42, 39] assess pointwise calibration but cannot determine whether samples are coherent trajectories, while best-of-K scores ignore most of the predicted distribution. Generic two-sample distances can detect mode collapse but become difficult to estimate in high dimensions and do not exploit the symmetry structure of bifurcating solutions [26, 9]. We address these challenges with Bi-FORK: Bifurcation Flow-matching Over Repulsion-guided Kinematics, a latent flow-matching framework for symmetry-breaking bifurcations (see Figure 1). Bi-FORK generates full trajectories jointly, keeping branch choices consistent across space and time rather than learning them through local transitions. Because full-resolution execution is intractable, flow matching operates in a compact latent space, decoupling global coordination from discretization. Although compression can average modes, branch identity dominates trajectory variance and is therefore preserved by the reconstruction objective without supervision. At inference, particle guidance [58, 14] repels concurrent samples, acting as a parallelized analog of classical deflation [10, 18] to recover distinct solutions. We assess the resulting distribution by matching predictions to their nearest ground truth and, per input, measuring orbit coverage under the relevant symmetry group via the Jensen-Shannon divergence, thereby detecting regime collapse and quantifying coverage of the underlying distribution. We evaluate the approach on three symmetry-breaking bifurcation problems: buckling beams, buckling mechanical metamaterials, and phase separation described by Allen-Cahn. These systems cover two- and three-dimensional grid- and mesh-based problems, various symmetry groups (both discrete and continuous), trajectories of 32 to 200 timesteps, and discretizations from a handful of points to 260,000. To our knowledge, this is the first demonstration of solution distribution recovery for one-to-many problems at that scale, and we show that our method outperforms existing approaches in both accuracy and efficiency. Our contributions are: 1. A latent flow-matching framework that models whole trajectories, making branch coherence structural, and scales several orders of magnitude beyond prior work. 2. A particle-guidance sampling mechanism recovering multiple distinct solutions in one amortized pass, as a parallel analog of deflation. 3. Diagnostics based on symmetry-orbit matching and per-input Jensen-Shannon divergence that detect regime collapse and extend distributional verification to continuous and infinite symmetry groups. 2 A Blueprint for Bifurcating Regimes 2.1 Bifurcation Systems Figure 2: Exemplary pitchfork bifurcation for the buckling metamaterials use case with increasing control protocol ℓ Recall that a bifurcation is the loss of uniqueness of a stable equilibrium as a control parameter crosses a critical value, with buckling as the canonical mechanical example [37, 6, 30]. Setting We represent trajectories as either point clouds or scalar fields on a regular grid. In both cases, t∈ℝN×CM^t ^N× C denotes the modeled response at N spatial locations and timestep t∈1,…,Tt∈\1,…,T\, and 1:T∈ℝT×N×CM^1:T ^T× N× C forms the complete response trajectory. For point clouds, the N locations are material points with reference coordinates 1∈ℝN×DX^1 ^N× D, where D∈2,3D∈\2,3\. The response tM^t is their displacement or microfluctuation field, so C=DC=D; deformed coordinates tX^t are recovered from 1X^1 and tM^t using the prescribed loading and dataset-specific kinematics. For regular-grid data, the N locations are fixed grid sites and tM^t is the scalar field itself, so C=1C=1. Datasets We consider three physical systems: (i) Beam3D models a slender beam that is compressed from one end until it buckles in arbitrary directions. (i) Mechanical Metamaterials captures the buckling behavior of 2D porous materials with different periodic patterns, leading to distinct stable shapes. (i) Allen-Cahn describes how a mixture separates into different phases over time, with a shared homogeneous initial state that evolves into several possible stable configurations. Modes as a symmetry orbit A problem instance is a condition =(1,ℓ)c=(X^1, ): a reference geometry together with a control protocol ℓ=(ℓ1,…,ℓT) =( ^1,…, ^T) that is fixed in advance and not influenced by the system’s response. The set of all feasible conditions we denote by C. For the mechanical metamaterials datasets, for example, ℓ is a quasi-static loading schedule generated by a single terminal load (tF^t interpolating toward TF^T), and the index t enumerates states and load increments simultaneously: tM^t is the equilibrium microfluctuation response to the input ℓt ^t (see Figure 2). For Allen-Cahn, the protocol is constant, ℓt≡(ϵ,μ) ^t≡(ε,μ), where ϵε is a length scale and μ determines local phase structure. The systems are multistable: under a fixed condition, several physically distinct outcomes—modes—are admissible. We denote their mode set by ℳM_c11 1 Here ‘Mode’ means a distinct admissible outcome, not an eigenmode: two symmetry-related metamaterial solutions stem from one buckling eigenmode but count as two modes. The usage agrees with statistics—outcomes follow p(⋅∣)p(· ), and for finite mode sets each branch is an isolated local maximum of this distribution; for the beam the distribution is uniform along the orbit, which is then its maximizer set., and write ℳtM_c^t for the set of their responses at timestep t. Symmetry organizes this set: with G=g∣g⋅= G_c=\\,g g·c=c\ \,\ (1) the symmetry group of the condition, any response M yields further solutions g⋅g·M, so the mode set contains full orbits, G⋅t⊆ℳt,Ht=g∈G∣g⋅t=t, G_c·M_c^t _c^t, H_c^t=\\,g∈ G_c g·M_c^t=M_c^t\,\, (2) where tM_c^t is a representative response and HtH_c^t its isotropy subgroup; each orbit is isomorphic to G/HtG_c/H_c^t with size [G:Ht][G_c:H_c^t]. For Beam3D and the mechanical metamaterials, the inclusion is an equality (a single orbit); for Allen-Cahn, the mode set is a union of orbits, one per morphology. Table 1 instantiates G_c, H_c, and ℳM_c per dataset; note e.g. the beam, where G=O(2)G_c=O(2) and the buckled state retains one reflection, giving the circular orbit O(2)/ℤ2≅SO(2)O(2)/Z_2 (2). Table 1: Symmetry structure per dataset: group G_c of the condition, isotropy H_c of an output state, resulting mode set ℳ≅G/HM_c G_c/H_c (union of orbits for Allen-Cahn). Dataset G_c H_c ℳM_c |ℳ||M_c| Beam3D O(2)O(2) ℤ2Z_2 circle of directions ϕφ ∞ Metamaterials wallpaper subgroup fixing ℓ pattern stabilizer translated patterns 1,2,41,2,4 Allen-Cahn domain syms ×u↦±u×\u ± u\ morphology symmetries orbits of morphologies varies Bifurcation time We define the bifurcation time as the first step at which the solution loses symmetry, tb=mint∈2,…,T∣Ht⊊Ht−1. t_b= \t∈\2,…,T\ H_c^t H_c^t-1\. (3) In practice, we estimate tbt_b as the first step at which trajectories under the same condition separate beyond a small threshold. This also accommodates Allen-Cahn, where all trajectories share the same homogeneous initial state but macroscopic modes emerge only later. The choice of representation, the practical estimator, and secondary bifurcations are discussed in Appendix B. Probabilistically, tbt_b marks the onset of distinct branches in the distribution of outcomes conditioned on the context c at timestep t, connecting physical symmetry breaking to the need for distribution-valued prediction. Three regimes of multimodality Our datasets span discrete (|ℳ|∈1,2,4|M_c|∈\1,2,4\, mechanical metamaterials of 29), continuous (|ℳ|=∞|M_c|=∞, Beam3D), and field-valued (Allen-Cahn) mode sets; see Appendix C. 2.2 Learning task Given a condition c, the task is to learn the conditional distribution over all feasible trajectories p(1:T∣),p\! (M^1:T ), (4) where the support of the target distribution is the mode set ℳM_c. The training datasets differ in their coverage of this support: the mechanical metamaterials enumerate all modes per condition, whereas Beam3D and Allen-Cahn provide samples from the mode set. We assume unknown perturbations that push the system toward one solution within an orbit have no preferred direction. Under this assumption, the condition’s symmetry implies that modes within an orbit are equally likely, so for single-orbit datasets the target is uniform over ℳM_c. Unlike standard trajectory prediction, the mapping from c to the system evolution is one-to-many: multiple physically valid trajectories arise under the same condition. During inference, the model receives c and generates one or more trajectory samples from a learned distribution pθflow(1:T∣)p_ _flow(M^1:T ). 2.3 Modeling Figure 3: Encoder-decoder structure for point cloud datasets. Input tokens (reference positions 1X^1 and displacements) are encoded into latent representation by cross-attending to L learned latent query tokens. The decoder reconstructs the displacements conditioned on the reference positions. Autoencoder We adopt a Perceiver-style cross-attention architecture [33, 4, 51] for all point cloud datasets. For each timestep t, the encoder maps the reference configuration 1X^1 and response field tM^t to L latent tokens of dimension DzD_z, i.e., t∈ℝL×DzZ^t ^L× D_z, using cross-attention with learned latent queries, followed by self-attention on the latent set. Importantly, the decoder reconstructs tM^t by querying these tokens at the fixed reference positions 1X^1. Thus, the latent representation tZ^t encodes all timestep-dependent states, while spatial reconstruction is always anchored to the same reference geometry 1X^1. This enables the independent analysis and prediction of latent trajectories tt=1T\Z^t\_t=1^T, which can then be decoded at the reference positions to recover tM^t. For regular grid datasets, we use a 3D Vision Transformer (ViT) [16, 45], detailed in the respective section (Section 3.3). Latent flow matching For the approximator, we draw inspiration from the flow-based model framework in [51]. Given TobsT_obs known states 1:TobsZ^1:T_obs and the condition c, it generates the remaining states Tobs+1:TZ^T_obs+1:T. The flow-based model then learns a path between noise and data. Let 1:=1:T∈ℝT×L×DzZ_1:=Z^1:T ^T× L× D_z be the true latent trajectory and ∼(0,I) ξ (0,I) be noise. For every τ∈[0,1]τ∈[0,1], a probability path is specified by τ=τ1+(1−τ), _τ=τ\,Z_1+(1-τ)\, ξ, (5) such that pure noise is obtained with τ=0τ=0, and the true trajectory at τ=1τ=1. We train a network with parameters θflow _flow to reverse the applied corruption: given the noised trajectory, the flow time τ, observed states and the control protocol, the clean trajectory is predicted by fθflow:(τ,τ,1:Tobs,)↦^1∈ℝT×L×Dz. f_ _flow:(Z_τ,τ,Z^1:T_obs,c) Z_1 ^T× L× D_z. (6) We adapt the standard flow matching approach through optimal transport-based coupling of in-batch pairs to minimize the path length between the noise and data distributions [57]. Additionally, we require the model to predict the bifurcation point, represented by a binary mask ^∈0,1T b∈\0,1\^T. At inference, ^1 Z_1 is recovered by solving the probability-flow ODE associated with Eq. 5, starting from pure noise. All timestep-dependent information for the prediction is expressed by τZ_τ and its conditioning signals. The formulation is model agnostic to the specific network architecture and training procedure used for fθflowf_ _flow (see Appendix D.1 for details). Inference using particle guidance Due to the multimodal nature of the problem, the control protocol and initial state do not determine which mode the system settles into, e.g., the orbit coordinate (e.g. the buckling direction ϕφ) is left free by c. While the generative model predicts a single valid mode rather than averaging across them, sampling alone does not control the distribution over modes. Figure 4: Particle guidance during sampling. Trajectories evolve identically up to bifurcation point tbt_b. From tbt_b, the repulsion term acts on the masked frames t≥tbt≥ t_b, pushing the K trajectories apart, yielding a multimodal distribution. We address this at inference with an approach inspired by particle guidance [14]. Rather than drawing K trajectories independently, we sample K candidate trajectories τ(1),…,τ(K)Z_τ^(1),…,Z_τ^(K) in a shared batch and integrate the probability-flow ODE for all of them jointly, leaving fθflowf_ _flow unmodified from training. To maintain diversity, we add an inter-sample repulsion term to the velocity field at each step (Figure 4). Each ^1(k) Z_1^(k) and its ordinary velocity field are computed independently via Eq. 6, as at training time. We additionally compute a repulsive term restricted to bifurcation-relevant frames, elected by the shared predicted mask ^∈0,1T b∈\0,1\^T, where b^t=1 b^t=1 indicates a predicted post-bifurcation frame. For each trajectory, the predicted values at these frames are masked and flattened into a single point; pairwise distances between trajectories’ points are weighted by a Gaussian RBF kernel with bandwidth set to the median pairwise distance. Each trajectory’s displacement force(k)force^(k) is the resulting kernel-weighted average of unit vectors pointing away from its neighbors, so trajectories that are similar at the bifurcation frames are pushed apart while already-distinct trajectories barely interact. The repulsive displacement modifies the velocity field v^θflow=(^1−τ)/(1−τ) v_ _flow=( Z_1-Z_τ)/(1-τ) of the probability-flow ODE for sample k. Schematically, v^(k)=wvelv^θflow(τ(k),τ)+wrep^⊙force(k), v^(k)=w_vel\, v_ _flow (Z_τ^(k),τ )+w_rep\, b ^(k), (7) where wvelw_vel weights the model velocity, wrepw_rep controls repulsion strength, and b selects the predicted post-bifurcation frames. The mask is broadcast over latent tokens and channels. This simplified expression omits the time-dependent scaling and normalization which we describe in Appendix D.5. As a result, all K trajectories are integrated jointly rather than independently: at likely branch points, samples that would otherwise collapse onto the same mode are pushed apart, so the sample set spreads across distinct branches instead of concentrating on one. 2.4 Metrics Evaluating generative models of bifurcating dynamical systems requires metrics that jointly capture trajectory-level precision and the distributional fidelity of predicted mode probabilities. Standard metrics such as MSE are inadequate in this setting, as they neither reflect the multimodal structure of the output space nor penalize mode collapse. We adopt two complementary metrics, adapted to systems with |ℳ|∈1,2,4|M_c|∈\1,2,4\ or |ℳ|=∞|M_c|=∞ responses depending on the condition c. Mode Conditioned MAE (MCon-MAE) quantifies trajectory-level precision. For each condition c, we draw K samples from the model and assign each to its nearest ground-truth mode by Euclidean distance. MCon-MAE is the resulting reconstruction error, averaged over samples and conditions: MCon-MAE=1||∑∈1K∑k=1Kmean(|^(k)−,m⋆(k)|)MCon-MAE= 1|C| _c 1K _k=1^Kmean\! ( | X_c^(k)-X_c,m_ ^(k) | ) (8) with mode assignment m⋆(k)=argminm∈ℳ‖^(k)−,m‖F2m_ ^(k)= _m _c \| X_c^(k)-X_c,m \|_F^2, where ,mX_c,m highlights the dependence on the condition c and mode m. We further report the Jensen-Shannon divergence (JSD) between the predicted and ground-truth mode distributions, estimated from K samples drawn at inference for each condition. This measures the extent to which a model represents each outcome with the correct relative frequency, rather than favoring a subset of modes. 3 Experiments 3.1 Beam3D Figure 5: Beam3D result comparison. Each panel shows an evaluated method, with an illustrative sample (colored) plotted over the actual trajectory (black). Their endpoints are shown against the full population and the ideal uniform dashed circle. Metrics are aggregated over all Beam3D test conditions using five random seeds. Bifurcating System and Learning Task The 3D buckling beam dataset represents multiple slender beams with a geometry 1X^1, which varies between data points, under axial compression. Each beam is fixed at one end with an incrementally increasing downward displacement d at the other. The input condition is therefore =(1,ℓ)c=(X^1, ). Below a critical threshold, the beam stays in axial compression. Beyond the threshold, the rotational symmetry breaks: the axisymmetric state remains an equilibrium but is no longer stable, and the beam buckles in an azimuthal direction ϕ∈[0,2π]φ∈[0,2π], yielding a continuous SO(2)SO(2) orbit of solutions. Given the condition c, the model is required to predict a representative set of K trajectories ^(k)k=1K\ X^(k)\_k=1^K, sampled from this continuous orbit. Modeling To model the bifurcating system, we use the architecture described in Section 2.3. In addition to tip displacement, we input per-node rotational stiffness and element-level attributes, including stiffness and length. At inference, we use particle guidance to distribute the samples more evenly over the admissible post-buckling manifold ℳM_c. We additionally evaluate Bi-FORK without particle guidance; we denote this variant Bi-FORK w/o PG. We sample 180180 times to estimate the distribution accuracy. Details on the chosen architecture parameters, such as latent dimension or hidden size, are provided in Appendix D. Metrics In addition to the metrics described in Section 2.4, we report the rejection rate (%\%): the fraction of generated trials whose mean Euclidean displacement from their matched mode exceeds a threshold. MCon-MAE and JSD are computed over all generated trials, including those that would have been rejected by the threshold. Results Bifurcating dynamical systems are naturally formulated as interacting particle systems, motivating models that operate on sets of arbitrary size without retraining. We compare against two such baselines: GeoTDM [27] and STFlow [8]. Figure 5 reports the MCon-MAE, JSD and rejection rate on the held-out test dataset. Among the models, GeoTDM achieves the lowest JSD due to its built-in equivariance, but rejects 99%99\% of its samples, with the worst MCon-MAE. STFlow reduces the trajectory rejection rate while maintaining low error and JSD values. Notably, Bi-FORK substantially improves MCon-MAE with a 0%0\% trajectory rejection rate and achieves a highly competitive JSD, actually closing the gap with the equivariant GeoTDM model. Visual inspection of the end-point distributions (Figure 5) confirms the quantitative findings: GeoTDM’s samples cluster off the target circle manifold, and STFlow samples are scattered noisily around it. Both Bi-FORK w/o PG and Bi-FORK trace the target manifold accurately. Table 11 summarizes performance across five random seeds, reporting the mean and standard deviation. We also compare inference runtime, shown in Table 15. Bi-FORK offers significant runtime advantages over graph-based baselines. It generates 180 samples for a single condition in 3.95 seconds, making it approximately 18× faster than STFlow and 1,966× faster than GeoTDM. 3.2 Mechanical Metamaterials Figure 6: Mechanical metamaterials result comparison. Ground truth modes and generated samples are shown for representative conditions with |ℳ|=2|M_c|=2 and |ℳ|=4|M_c|=4 buckling modes. Colors indicate the mapped ground-truth mode; missing modes are shown in gray. Affine baseline results in mean prediction. Reported MCon-MAE is averaged over the full test set, while JSD and full coverage are reported for |ℳ|∈2,4|M_c|∈\2,4\. Metrics reported over five random seeds on the test dataset. Bifurcating System and Learning Task The dataset of Hendriks et al. [29] contains buckling trajectories of periodic porous metamaterials spanning all 17 wallpaper groups. Depending on geometry and loading, symmetry breaking produces 1, 2, or 4 translation-related solutions, i.e., |ℳ|∈1,2,4|M_c|∈\1,2,4\: one mode when no observable symmetry breaking occurs, and two or four modes after isotropy reduction, with four modes potentially arising after a secondary bifurcation. The condition is therefore =(1,ℓ)c=(X^1, ) with ℓt=t ^t=F^t, with tF^t the applied macroscopic deformation gradient at time t. Since tF^t is known, the predicted ^tt=1T\ M^t\_t=1^T are the microfluctuation fields instead of the full deformation field. The deformed coordinates are then reconstructed as ^t=1(t)⊤+^t X^t=X^1(F^t) + M^t. Unlike the beam, the mode distribution here is discrete: the model must predict the correct number of modes for a given condition, because it should generate each physically plausible branch in ℳM_c. Modeling Following the architecture in Section 2.3, we use the first frame as the conditioning frame for each geometry. Additionally, we add node type information and the crystallographic group. Inference is executed similarly to the Beam3D dataset, with one difference for Bi-FORK: when no bifurcation time point is detected (|ℳ|=1|M_c|=1), the sampling algorithm discards the repulsive term in the diffusion process. In this setting, we sample 1616 times to estimate the potential distribution for each condition. Metrics We employ the metrics from Section 2.4. Additionally, we report coverage (%\%), defined as the percentage of conditions for which every ground-truth mode is represented by a valid generated sample. We report JSD and coverage only for cases with |ℳ|∈2,4|M_c|∈\2,4\. MCon-MAE and JSD are computed over all generated trials, including those that would have been rejected by the threshold. Table 2: Mode distributions of the two illustrated geometries in Figure 6 with |ℳ|∈2,4|M_c|∈\2,4\. Each row aggregates five seeds, with K=16K=16 samples per seed. Sample Method K Mode 1 Mode 2 Mode 3 Mode 4 2 modes Bi-FORK w/o PG 16 80.0% 20.0% – – Bi-FORK 16 56.2% 43.8% – – 4 modes Bi-FORK w/o PG 16 18.2% 0.0% 72.7% 9.1% Bi-FORK 16 23.2% 15.4% 36.3% 25.1% Results GeoTDM and STFlow were not viable baselines for this problem setting because by default, both rely on a fully connected graph, which is infeasible with thousands of nodes. They can also use sparse connectivity based on a distance threshold, but then they would need a number of message-passing steps equal to the graph diameter to achieve all-to-all communication, which is infeasible too. On Beam3D (⩽11 11 nodes for T=200T=200), both models showed no improvement and, being GNN-based, were already slow to train and sample from, requiring up to three orders of magnitude more inference time than Bi-FORK (Table 15). With roughly 160×160× more nodes per timestep and all timesteps processed jointly (Section 1), both models become computationally infeasible. We therefore compare against the affine baseline, which propagates 1X^1 using only the deformation gradient. Figure 6 shows that Bi-FORK w/o PG already improves over the affine baseline in accuracy and in distributional coverage: the per-mode JSD shows that it recovers multiple modes rather than collapsing to a mean prediction. With the full potential of our method, we further improve coverage, retrieving all four modes in roughly 50%50\% of the sampled test cases. Table 12 provides the quantitative summary, including standard deviation. Table 2 shows the diversity improvement for the plotted samples with two and four modes, producing a nearly uniform distribution with Bi-FORK. Finally, we ablate the bifurcation head on Beam3D and mechanical metamaterials. The results in Table 14 show that the learned bifurcation masks are predicted accurately on both datasets. 3.3 Allen-Cahn Figure 7: Allen-Cahn result comparison. Ground truth (black open circles) and predictions (teal and pink circles) are shown across a parameter sweep μ∈[−0.1,1.0]μ∈[-0.1,1.0] at fixed ϵ=0.1ε=0.1. The vertical axis shows the final-frame mean value. We evaluate metrics on the test dataset and average over five random seeds; lower values indicate better performance. Bifurcating System and Learning Task The 3D Allen-Cahn dataset describes the evolution of a phase-field variable governed by the Allen-Cahn equation, which models phase separation in binary mixtures on a periodic domain. Because all trajectories share the same initial homogeneous state, the condition is only =(ϵ,μ)c=(ε,μ), which is fixed over time. This initial state is invariant to G_c, which contains phase exchange u↦−u -\,u, periodic translations and the cubic point group. The initial state can evolve into many patterns, which can be symmetry-related through G_c or morphologically distinct. The number of morphologically distinct modes rises as ϵ−3ε^-3: shrinking ϵε opens more unstable directions. Modeling Given Allen-Cahn’s grid structure, the first stage uses a ViT instead of the Perceiver-style encoder. The periodic scalar field is compressed into a compact spatial-latent-token grid. Implementation details are in Appendix D. The second stage follows Section 2.3. As input conditions, we only use the physical parameters ϵε and μ; since the full trajectory is generated starting at t=1t=1, no additional conditioning frame is needed. The output is the full scalar-field rollout. We apply symmetry-aware augmentations (periodic translations and sign flips) that match the domain’s periodicity and point-group structure. At inference, the sampling follows the previously described pipeline, with one exception: we apply particle guidance at every valid frame by setting b^t=1 b^t=1 throughout the sampling process. Therefore, we do not train a bifurcation classification head for this problem setting. Metrics We extend the metrics described in Section 2.4 to the Allen-Cahn solution space. A raw distance between a generated field and any single reference is uninformative when the reference could equally be shifted, reflected, rotated, or sign-inverted and still represent the same solution. We therefore first match each generated trial to the closest branch among a predefined set of ground-truth solutions using an FFT-based descriptor invariant to these symmetries, then align it to that branch via an exhaustive search over the point group. Full implementation details are given in Appendix D.8. Estimating the JSD on the joint distribution over (branch, rotoreflection, shift) is infeasible, as the number of possible combinations exceeds the number of available samples. We therefore split the JSD metrics into three separate components: mode-coverage JSD over the predefined solutions (including their sign-inverted counterparts); a translation JSD that similarly fixes a branch and tests whether the recovered periodic shift is uniform; and a rotoreflection JSD that fixes a single matched branch and tests whether the recovered rotations and reflections relative to it are uniform. MCon-MAE compares the full predicted trajectory with its relaxed counterpart, obtained by a Gauss–Newton projection onto the Allen–Cahn PDE. It thus measures the distance to a physically valid solution rather than to precomputed references. Figure 8: Visual comparison of ground truth vs predicted samples for μ=0.86μ=0.86 and ϵ=0.094ε=0.094. Results Figure 7 summarizes the results, and Table 13(b) reports the mean and standard deviation over five seeds. Bi-FORK reduces the MCon-MAE from 0.052 (without particle guidance) to 0.045, showing that the generated trajectories require smaller corrections to satisfy the PDE. It also lowers the mode JSD from 0.131 to 0.089, reflecting improved coverage of the solution branches. The translation and rotoreflection JSDs are identical for both variants. To illustrate mode coverage, we fix ϵ=0.1ε=0.1 and vary only μ. Without particle guidance, the model samples mostly from the positive branch and misses several solutions on the negative branch. Bi-FORK samples both branches with similar density over the full range of μ. Figure 8 shows samples from both methods next to the ground truth. Even when the samples differ, they showcase valid solutions. 4 Discussion Limitations Our experiments show that Bi-FORK applies to a diverse set of bifurcating systems; however, some limitations remain. First, particle guidance can trade accuracy for diversity: the repulsion acts on all samples, including those already in the same valid mode, and can push them off the solution manifold. Second, the Allen–Cahn evaluation is necessarily incomplete. The full distribution over morphologies is unavailable, so the mode-coverage JSD relies on four simulated realizations per condition and their sign flips. In addition, separate JSDs for modes, translations, and rotoreflections do not guarantee coverage of the joint distribution. Third, although symmetry is central to the problem, the model learns it from data and augmentation rather than enforcing it architecturally. Conclusion We introduce Bi-FORK, a framework for learning symmetry-breaking bifurcations. It maps solution multiplicity via symmetry orbits and bifurcation time. The task is framed as learning a conditional distribution over complete trajectories using latent flow matching and particle guidance. We evaluate models using symmetry-aware metrics that detect mode collapse. Across buckling beams, mechanical metamaterials, and Allen-Cahn phase separation, Bi-FORK recovers multimodal solutions at discretizations up to 260,000 points, several orders of magnitude beyond prior work. Because no existing stochastic method scales to this size, we evaluate the largest systems against a deterministic affine baseline and our own ablations. Ultimately, this opens generative modeling to large-scale bifurcating physical systems. Acknowledgments We wish to thank Andreas Fürst for helpful discussions and feedback. The ELLIS Unit Linz, the LIT AI Lab, the Institute for Machine Learning, are supported by the Federal State Upper Austria. We thank the European Union’s Horizon Europe research and innovation programme under grant agreement number 101214398 (ELLIOT) and the Austrian Science Fund (FWF) 10.55776/COE12. We thank Audi AG, Silicon Austria Labs (SAL), Merck Healthcare KGaA, GLS (Univ. Waterloo), TÜV Holding GmbH, Software Competence Center Hagenberg GmbH, dSPACE GmbH, TRUMPF SE + Co. KG. We thank the projects FWF AIRI FG 9-N (10.55776/FG9), AI4GreenHeatingGrids (FFG-899943), FWF Bilateral Artificial Intelligence (10.55776/COE12). AI use statement In this work, we used generative AI tools for three tasks requiring disclosure: implementing methods, cleaning and reformatting datasets, and formulating mathematical claims, where AI served only as a source of inspiration. We used Codex to implement the baseline methods and supporting evaluation code in Section 3; two authors reviewed, verified, and tested this code, and we implemented our proposed method ourselves. We used Codex and Claude to clean and reformat the datasets and checked the processed data against the original sources. Codex and Claude also offered ideas for the mathematical claims in Section 2, but the final claims were formulated and verified by the authors. We have not used generative AI tools to generate synthetic datasets; to develop the problem setting, problem framing, or proposed architecture; to propose or refine hypotheses; to design or give feedback on experiments; or to interpret results. Translation and qualitative or thematic data analysis are not applicable to this work. Additionally, we used generative AI tools to identify relevant literature, to polish the manuscript for readability, and to make cosmetic improvements to figures. We did not use them for other recommended-disclosure tasks such as drafting text, suggesting hyperparameters, identifying research gaps, formatting references, or proposing a title. We have reviewed all AI-assisted work: suggested papers were read by the authors before being cited, text edits were checked to ensure that the technical content and claims were unchanged, and figure edits were limited to styling, with the underlying numbers never modified by AI. We take responsibility for the final content of this work, including text, claims, and artifacts produced with the aid of generative AI. Reproducibility statement We provide a detailed description of the datasets, including the evaluation splits used, in Appendix C. Appendix D also elaborates on the model architecture, training pipeline, loss functions, and hyperparameters for each dataset. Additionally we provide a detailed overview on the evaluation pipeline for each dataset in Appendix D.7 for Beam3D and mechanical metamaterials, and in Appendix D.8 for the Allen-Cahn dataset. The code is available at: https://github.com/ml-jku/Bi-FORK. References [1] J. Abramson, J. Adler, J. Dunger, R. Evans, T. Green, A. Pritzel, O. Ronneberger, L. Willmore, A. J. Ballard, J. Bambrick, et al. (2024) Accurate structure prediction of biomolecular interactions with alphafold 3. Nature 630 (8016), p. 493–500. Cited by: §1. [2] M. Albergo, N. M. Boffi, and E. Vanden-Eijnden (2025) Stochastic interpolants: a unifying framework for flows and diffusions. Journal of Machine Learning Research 26 (209), p. 1–80. Cited by: §D.1. [3] B. Alkin, M. Bleeker, R. Kurle, T. Kronlachner, R. Sonnleitner, M. Dorfer, and J. Brandstetter (2025) AB-upt: scaling neural cfd surrogates for high-fidelity automotive aerodynamics simulations via anchored-branched universal physics transformers. Transactions on Machine Learning Research. Cited by: §1. [4] B. Alkin, A. Fürst, S. Schmid, L. Gruber, M. Holzleitner, and J. Brandstetter (2024) Universal physics transformers: a framework for efficiently scaling neural operators. Advances in Neural Information Processing Systems 37, p. 25152–25194. Cited by: §D.1, §1, §2.3. [5] U. Alon and E. Yahav (2021) On the bottleneck of graph neural networks and its practical implications. International Conference on Learning Representations. Cited by: §1. [6] S. S. Antman (2005) Nonlinear problems of elasticity. 2 edition, Applied Mathematical Sciences, Vol. 107, Springer, New York. Cited by: Appendix B, §2.1. [7] Z. P. Bažant and L. Cedolin (2003) Stability of structures: elastic, inelastic, fracture, and damage theories. Courier Corporation. Cited by: §1. [8] K. Bennema ten Brinke, K. Minartz, and V. Menkovski (2026) STFlow: data-coupled flow matching for geometric trajectory simulation. In International conference on machine learning, Cited by: §3.1. [9] S. Bischoff, A. Darcher, M. Deistler, R. Gao, F. Gerken, M. Gloeckler, L. Haxel, J. Kapoor, J. K. Lappalainen, J. H. Macke, et al. (2024) A practical guide to sample-based statistical distances for evaluating generative models in science. Transactions on Machine Learning Research. Cited by: §1. [10] K. M. Brow and W. B. Gearhart (1971) Deflation techniques for the calculation of further solutions of a nonlinear system. Numerische Mathematik 16 (4), p. 334–342. Cited by: §1, §1. [11] C. Cai and Y. Wang (2020) A note on over-smoothing for graph neural networks. arXiv preprint arXiv:2006.13318. Cited by: §1. [12] G. Calzolari and W. Liu (2025) Accelerating large eddy simulations of urban airflow with generative adversarial networks. Building and Environment 286, p. 113622. Cited by: §1. [13] D. Chen, Y. Lin, W. Li, P. Li, J. Zhou, and X. Sun (2020) Measuring and relieving the over-smoothing problem for graph neural networks from the topological view. In Proceedings of the AAAI conference on artificial intelligence, Vol. 34, p. 3438–3445. Cited by: §1. [14] G. Corso, Y. Xu, V. De Bortoli, R. Barzilay, and T. Jaakkola (2024) Particle guidance: non-iid diverse sampling with diffusion models. In International Conference on Learning Representations, Vol. 2024, p. 22480–22507. Cited by: §D.5, §D.5, §1, §2.3. [15] M. G. Crandall and P. H. Rabinowitz (1971) Bifurcation from simple eigenvalues. Journal of Functional Analysis 8 (2), p. 321–340. Cited by: Appendix B. [16] A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly, et al. (2021) An image is worth 16x16 words: transformers for image recognition at scale. International Conference on Learning Representations. Cited by: §2.3. [17] G. Fabiani, F. Calabrò, L. Russo, and C. Siettos (2021) Numerical solution and bifurcation analysis of nonlinear partial differential equations with extreme learning machines. Journal of Scientific Computing 89 (2), p. 44. Cited by: §1. [18] P. E. Farrell, A. Birkisson, and S. W. Funke (2015) Deflation techniques for finding distinct solutions of nonlinear partial differential equations. SIAM Journal on Scientific Computing 37 (4), p. A2026–A2045. Cited by: §1, §1. [19] Z. Gao, X. Shi, B. Han, H. Wang, X. Jin, D. Maddix, Y. Zhu, M. Li, and Y. B. Wang (2023) Prediff: precipitation nowcasting with latent diffusion models. Advances in Neural Information Processing Systems 36, p. 78621–78656. Cited by: §1, §1. [20] M. Giacomini and P. Díez (2026) Surrogates for physics-based and data-driven modelling of parametric systems: review and new perspectives. Archives of Computational Methods in Engineering 33, p. 8519–8551. Cited by: §1. [21] M. Golubitsky, I. Stewart, and D. G. Schaeffer (1988) Singularities and groups in bifurcation theory, volume i. Applied Mathematical Sciences, Vol. 69, Springer, New York. Cited by: Appendix B. [22] M. Golubitsky and I. Stewart (2002) The symmetry perspective: from equilibrium to chaos in phase space and physical space. Progress in Mathematics, Vol. 200, Birkhäuser, Basel. Cited by: Appendix B. [23] P. Grönquist, C. Yao, T. Ben-Nun, N. Dryden, P. Dueben, S. Li, and T. Hoefler (2021) Deep learning for post-processing ensemble weather forecasts. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 379 (2194), p. 20200092. Cited by: §1. [24] J. Guckenheimer and P. Holmes (1983) Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Applied Mathematical Sciences, Vol. 42, Springer, New York. Cited by: Appendix B, Appendix B. [25] L. F. Guilhoto, A. Kaushal, and P. Perdikaris (2026) Multimodal scientific learning beyond diffusions and flows. arXiv preprint arXiv:2602.00960. Cited by: §1. [26] T. M. Hamill (2001) Interpretation of rank histograms for verifying ensemble forecasts. Monthly Weather Review 129 (3), p. 550–560. Cited by: §1. [27] J. Han, M. Xu, A. Lou, H. Ye, and S. Ermon (2024) Geometric trajectory diffusion models. Advances in Neural Information Processing Systems 37, p. 25628–25662. Cited by: §3.1. [28] S. He, H. Liang, Y. Zhang, and X. Yuan (2025) Skillful high-resolution ensemble precipitation forecasting with an integrated deep learning framework. arXiv preprint arXiv:2501.02905. Cited by: §1. [29] F. Hendriks, V. Menkovski, M. Doškář, M. G. Geers, K. Verbeek, and O. Rokoš (2025) Wallpaper group-based mechanical metamaterials: dataset including mechanical responses. Scientific Data 12, p. 1880. Cited by: §C.2, §2.1, §3.2. [30] F. Hendriks, O. Rokoš, M. Doškář, M. G. Geers, and V. Menkovski (2025) Equivariant flow matching for symmetry-breaking bifurcation problems. arXiv preprint arXiv:2509.03340. Cited by: §1, §2.1. [31] D. Hendrycks and K. Gimpel (2016) Gaussian error linear units (gelus). arXiv preprint arXiv:1606.08415. Cited by: §D.1. [32] J. Ho, A. Jain, and P. Abbeel (2020) Denoising diffusion probabilistic models. Advances in neural information processing systems 33, p. 6840–6851. Cited by: §D.7. [33] A. Jaegle, F. Gimeno, A. Brock, A. Zisserman, O. Vinyals, and J. Carreira (2021) Perceiver: general perception with iterative attention. In International conference on machine learning, p. 4651–4664. Cited by: §D.1, §2.3. [34] H. Kielhöfer (2012) Bifurcation theory: an introduction with applications to partial differential equations. 2 edition, Applied Mathematical Sciences, Vol. 156, Springer, New York. Cited by: Appendix B. [35] L. Klein, A. Krämer, and F. Noé (2023) Equivariant flow matching. Advances in Neural Information Processing Systems 36, p. 59886–59910. Cited by: §D.8. [36] D. Kochkov, J. Yuval, I. Langmore, P. Norgaard, J. Smith, G. Mooers, M. Klöwer, J. Lottes, S. Rasp, P. Düben, et al. (2024) Neural general circulation models for weather and climate. Nature 632 (8027), p. 1060–1066. Cited by: §1. [37] W. T. Koiter (1945) On the stability of elastic equilibrium. Ph.D. Thesis, Delft University of Technology. Note: English translation NASA T F-10833 (1967) / AFFDL-TR-70-25 (1970) Cited by: Appendix B, §2.1. [38] Y. A. Kuznetsov (2004) Elements of applied bifurcation theory. 3 edition, Applied Mathematical Sciences, Vol. 112, Springer, New York. Cited by: Appendix B. [39] E. Larsson, J. Oskarsson, T. Landelius, and F. Lindsten (2026) CRPS-lam: probabilistic regional weather forecasting with continuous ranked probability score. arXiv preprint arXiv:2510.09484. Cited by: §1. [40] I. Loshchilov and F. Hutter (2019) Decoupled weight decay regularization. International Conference on Learning Representations. Cited by: §D.3. [41] N. Ma, M. Goldstein, M. S. Albergo, N. M. Boffi, E. Vanden-Eijnden, and S. Xie (2024) Sit: exploring flow and diffusion-based generative models with scalable interpolant transformers. In European Conference on Computer Vision, p. 23–40. Cited by: §D.1, §D.1. [42] A. Mahesh, W. Collins, B. Bonev, N. Brenowitz, Y. Cohen, J. Elms, P. Harrington, K. Kashinath, T. Kurth, J. North, et al. (2025) Huge ensembles part i: design of ensemble weather forecasts using spherical fourier neural operators. Geoscientific Model Development 18, p. 5575–5603. Cited by: §1. [43] K. Oono and T. Suzuki (2020) Graph neural networks exponentially lose expressive power for node classification. International Conference on Learning Representations. Cited by: §1. [44] F. Pichi, F. Ballarin, G. Rozza, and J. S. Hesthaven (2023) An artificial neural network approach to bifurcating phenomena in computational fluid dynamics. Computers & Fluids 254, p. 105813. Cited by: §1. [45] C. Prabhakar, H. Li, J. Yang, S. Shit, B. Wiestler, and B. Menze (2024) Vit-ae++: improving vision transformer autoencoder for self-supervised medical image representations. In Medical Imaging with Deep Learning, p. 666–679. Cited by: §2.3. [46] I. Price, A. Sanchez-Gonzalez, F. Alet, T. R. Andersson, A. El-Kadi, D. Masters, T. Ewalds, J. Stott, S. Mohamed, P. Battaglia, et al. (2025) Probabilistic weather forecasting with machine learning. Nature 637 (8044), p. 84–90. Cited by: §1. [47] M. Rietkerk, R. Bastiaansen, S. Banerjee, J. van de Koppel, M. Baudena, and A. Doelman (2021) Evasion of tipping in complex systems through spatial pattern formation. Science 374 (6564), p. eabj0359. Cited by: §1. [48] I. Romero and M. Ortiz (2026) A note on data-driven methods for mechanical problems with non-unique solutions. Meccanica 61 (1), p. 26. Cited by: §1. [49] D. Ruelle and F. Takens (1971) On the nature of turbulence. Les rencontres physiciens-mathématiciens de Strasbourg-RCP25 12, p. 1–44. Cited by: Appendix B, §1. [50] C. Rupprecht, I. Laina, R. DiPietro, M. Baust, F. Tombari, N. Navab, and G. D. Hager (2017) Learning in an uncertain world: representing ambiguity through multiple hypotheses. In 2017 IEEE international conference on computer vision (ICCV), p. 3611–3620. Cited by: §1. [51] F. Sestak, A. P. Toshev, A. Fürst, G. Klambauer, A. Mayr, and J. Brandstetter (2025) LaM-SLide: latent space modeling of spatial dynamical systems via linked entities. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, Cited by: §D.1, §2.3, §2.3. [52] R. Seydel (2010) Practical bifurcation and stability analysis. 3 edition, Interdisciplinary Applied Mathematics, Vol. 5, Springer, New York. Cited by: Appendix B. [53] J. Song, C. Meng, and S. Ermon (2021) Denoising diffusion implicit models. International Conference on Learning Representations. Cited by: §D.7. [54] W. Steffen, J. Rockström, K. Richardson, T. M. Lenton, C. Folke, D. Liverman, C. P. Summerhayes, A. D. Barnosky, S. E. Cornell, M. Crucifix, et al. (2018) Trajectories of the earth system in the anthropocene. Proceedings of the national academy of sciences 115 (33), p. 8252–8259. Cited by: §1. [55] H. Tang, L. Shi, Y. Zhang, X. Liu, J. Tang, and J. Lu (2026) Learning structure, energy, and dynamics: a survey of artificial intelligence for protein dynamics. arXiv preprint arXiv:2604.25244. Cited by: §1. [56] J. M. T. Thompson and G. W. Hunt (1973) A general theory of elastic stability. Wiley, London. Cited by: Appendix B. [57] A. Tong, K. Fatras, N. Malkin, G. Huguet, Y. Zhang, J. Rector-Brooks, G. Wolf, and Y. Bengio (2024) Improving and generalizing flow-based generative models with minibatch optimal transport. Transactions on Machine Learning Research. Cited by: §2.3. [58] T. Unterthiner, B. Nessler, C. Seward, G. Klambauer, M. Heusel, H. Ramsauer, and S. Hochreiter (2018) Coulomb gans: provably optimal nash equilibria via potential fields. International Conference on Learning Representations. Cited by: §1. [59] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin (2017) Attention is all you need. Advances in neural information processing systems 30. Cited by: §D.1. [60] D. Yoon, M. Seo, D. Kim, Y. Choi, and D. Cho (2025) Probabilistic weather forecasting with deterministic guidance-based diffusion model. p. 108–124. Cited by: §1. [61] Z. Zou, Z. Wang, and G. Em Karniadakis (2025) Learning and discovering multiple solutions using physics-informed neural networks with random initialization and deep ensemble. Proceedings of the Royal Society A 481 (2325), p. 20250205. Cited by: §1. Appendix Contents Appendix A Notation The notation used throughout the paper is summarized in Table 3. Table 3: Overview of used symbols and notations. Superscripts t index trajectory frames, subscripts τ flow time, and parenthesized superscripts (k)(k) generated samples. Dataset-specific quantities are defined where they are used (Sections 3 and C). Symbol Meaning Space Data and trajectories N, n number of spatial locations (nodes); node index ℕN T, t number of timesteps; timestep index ℕN D, C spatial dimension; channels of the response 2,3\2,3\, 1,D\1,D\ 1X^1, tX^t reference geometry; deformed coordinates at timestep t ℝN×DR^N× D tM^t response (signal field) at timestep t ℝN×CR^N× C 1:TM^1:T response trajectory (1,…,T)(M^1,…,M^T) ℝT×N×CR^T× N× C Conditions, modes, and symmetry c, C condition (1,ℓ)(X^1, ), or (ϵ,μ)(ε,μ) for Allen-Cahn; set of conditions ℓ control protocol (loading schedule) (ℓ1,…,ℓT)( ^1,…, ^T) ℳM_c, ℳtM_c^t mode set (admissible trajectories); responses at timestep t ⊂ℝT×N×C ^T× N× C G_c symmetry group of the condition tM_c^t representative response at timestep t ∈ℳt _c^t HtH_c^t isotropy subgroup of tM_c^t ≤G≤ G_c tbt_b bifurcation time: first strict decrease of HtH_c^t (∞ if none) 2,…,T∪∞\2,…,T\∪\∞\ Distributions p(1:T∣)p(M^1:T ) ground-truth distribution of trajectories, supported on ℳM_c pθflow(1:T∣)p_ _flow(M^1:T ) learned distribution (latent flow model, decoded) Latent representation ℰϕencE_ _enc, ϕdecD_ _dec encoder (1,t)↦t(X^1,M^t) ^t; decoder (t,1)↦^t(Z^t,X^1) M^t tZ^t latent state at timestep t ℝL×DzR^L× D_z L, DzD_z number of latent tokens; token dimension ℕN Latent flow matching τ flow time (00: noise, 11: data) [0,1][0,1] 1Z_1 clean latent trajectory 1:TZ^1:T ℝT×L×DzR^T× L× D_z ξ Gaussian noise (0,I)N(0,I) τZ_τ noisy latent trajectory τ1+(1−τ) _1+(1-τ) ξ ℝT×L×DzR^T× L× D_z TobsT_obs number of observed (conditioning) frames 0,…,T\0,…,T\ fθflowf_ _flow, θflow _flow predictor (τ,τ,1:Tobs,)↦^1(Z_τ,τ,Z^1:T_obs,c) Z_1; its parameters ^1 Z_1 predicted clean latent trajectory ℝT×L×DzR^T× L× D_z v^θflow v_ _flow velocity field of the probability-flow ODE, induced by ^1 Z_1 ℝT×L×DzR^T× L× D_z Sampling and particle guidance K, k number of jointly generated samples; sample index ℕN τ(k)Z_τ^(k) latent state of sample k ℝT×L×DzR^T× L× D_z ptp^t, b predicted post-bifurcation probability; mask b^t=[pt≥1/2] b^t= 1[p^t≥ 1/2] [0,1][0,1], 0,1T\0,1\^T force(k)force^(k) repulsive force on sample k ℝT×L×DzR^T× L× D_z wvelw_vel, wrepw_rep weights of model velocity and repulsion ℝ≥0R_≥ 0 v^(k) v^(k) guided velocity of sample k ℝT×L×DzR^T× L× D_z Training and evaluation ℒrecL_rec, ℒflowL_flow, ℒbifL_bif reconstruction, latent-flow, and bifurcation losses btb^t, rtr^t bifurcation label [t≥tb] 1[t≥ t_b]; valid-frame indicator 0,1\0,1\ ^(k) X_c^(k), ,mX_c,m generated trajectory k; ground-truth trajectory of mode m m⋆(k)m_ ^(k) ground-truth mode matched to sample k MCon-MAE, JSD mode-conditioned MAE; Jensen–Shannon divergence ℝ≥0R_≥ 0 Appendix B Additional bifurcation context Estimating the bifurcation time In the datasets, tbt_b is obtained from simulated trajectories rather than from the isotropy subgroups directly. For two trajectories i,ji,j under the same condition, we compute the pairwise divergence dijt=∥it−jt∥d_ij^t= _i^t-M_j^t (the load-prescribed part cancels) and take the first timestep at which it exceeds a small threshold ϑ>0 >0. Before the bifurcation, the equilibrium is unique and G_c-invariant, so all trajectories coincide; afterwards, trajectories on distinct branches separate. For the quasi-static datasets, this estimate therefore recovers tbt_b up to the resolution ϑ . For Allen-Cahn, where small noise perturbs the shared homogeneous initial state, it marks the onset of macroscopic domain formation. Non-bifurcating conditions If no timestep satisfies Ht⊊Ht−1H_c^t H_c^t-1, the condition admits a single mode, and we set tb:=∞t_b:=∞. The bifurcation labels are then bt=[t≥tb]≡0b^t= 1[t≥ t_b]≡ 0, and the estimator above detects no separation. Accordingly, the bifurcation head is trained to predict ^= b=0 for such conditions, which deactivates the repulsion term in Eq. 7. Secondary bifurcations A loading path can exhibit several isotropy drops tb(1)<tb(2)<⋯t_b^(1)<t_b^(2)<·s, i.e., hierarchical bifurcations; the four-mode metamaterials arise from two such drops (Figure 10(b)). Throughout, tbt_b denotes the first one. Choice of representative The definition of tbt_b in Section 2.1 refers to a representative response tM_c^t at every timestep. Two responses on the same orbit have conjugate isotropy subgroups, Hg⋅=gHg−1H_g·M=g\,H_M\,g^-1 for g∈Gg∈ G_c, and conjugate subgroups are in general not nested. The inclusion Ht⊊Ht−1H_c^t H_c^t-1 is therefore to be read along a single trajectory: tM_c^t denotes frame t of one fixed admissible trajectory 1:T∈ℳM_c^1:T _c. With this convention, tbt_b does not depend on the trajectory chosen. Replacing 1:TM_c^1:T by g⋅1:T=(g⋅1,…,g⋅T)g·M_c^1:T=(g·M_c^1,…,g·M_c^T) conjugates every HtH_c^t by the same g; conjugation by a fixed element preserves inclusions and their strictness, so the first timestep at which Ht⊊Ht−1H_c^t H_c^t-1 holds is unchanged. For Allen-Cahn, where ℳM_c is a union of orbits, this defines tbt_b per orbit; the estimator above applies unchanged. Bifurcation, not chaos Unlike chaotic systems, where nearby trajectories diverge exponentially and long-run outcomes are unpredictable [49, 24], multistable systems branch into a small, structured set of symmetry-related outcomes. This makes the problem well posed for learning: the target is a multimodal distribution over a few valid branches, not an unpredictable attractor. Relation to bifurcation theory For the mechanical datasets, t counts load increments, so the “bifurcation time” is really a bifurcation load: the classical setting of equilibrium bifurcation theory, where a solution branch splits once the equilibrium stops being unique—mechanically, where the tangent stiffness matrix becomes singular [37, 56, 6, 15, 34, 52]. In symmetric problems, the new equilibria have strictly less symmetry and appear in symmetry-related copies, which is exactly the isotropy reduction defining tbt_b in Section 2.1 and the content of equivariant bifurcation theory [21, 22]. Allen-Cahn is instead genuinely time-dependent and covered by the dynamical version of the theory [24, 38]. We consider the generic supercritical scenario in which branches grow continuously out of the symmetric state; sudden jumps (subcritical or imperfect buckling) are out of scope. Appendix C Datasets C.1 Beam3D The dataset describes a 3D beam-buckling process, containing simulations with 33 to 1111 nodes. In the problem setting, the beam is fixed at the base and aligned along the z-axis in its undeformed state, with a specified downward displacement applied incrementally to the tip. This displacement sequence is the control protocol ℓ=(ℓ1,…,ℓT) =( ^1,…, ^T) in the condition =(1,ℓ)c=(X^1, ). Below a critical displacement threshold, the beam undergoes pure axial compression; beyond it, the system loses stability and deflects in any direction ϕφ drawn uniformly from [0,2π][0,2π] in the x-y-plane. As every deflection direction is equally probable, the system admits |ℳ|=∞|M_c|=∞ distinct post-buckling modes under identical loading. Dataset generation The beam contains NsegN_seg compressible elements, each characterized by a reference length LibeamL_i^beam, axial stiffness kiaxk_i^ax, and rotational spring constant kirotk_i^rot resisting bending between neighboring elements. The mechanical state is parametrized by axial strains ei=ln(Lidef/Libeam)e_i= (L_i^def/L_i^beam), where LidefL_i^def is the deformed element length, and qiq_i represents the bending angle from the vertical. The resulting node positions are expressed relative to the reference geometry as the displacement trajectory tt=1T\M^t\_t=1^T used by the model. Equilibrium configurations minimize the strain energy Ebeam(,)=12[k1rotq12+∑i=2Nsegkirot(qi−qi−1)2]+12∑i=1NsegLibeamkiaxei2E_beam(q,e)= 12 [k_1^rotq_1^2+ _i=2^N_segk_i^rot(q_i-q_i-1)^2 ]+ 12 _i=1^N_segL_i^beamk_i^axe_i^2 (9) subject to d=∑i=1NsegLibeam[1−exp(ei)cosqi]d= _i=1^N_segL_i^beam[1- (e_i) q_i], enforced via a Lagrange multiplier. The dataset contains 997997 trajectories of up to T=200T=200 timesteps, simulated to full deflection d, with system parameters NsegN_seg, LibeamL_i^beam, kirotk_i^rot, kiaxk_i^ax, and ϕφ sampled across instances. The sampling distributions are summarized in Table 4. Table 4: Sampling distributions across used parameters. Parameter Range Add. Information NsegN_seg 3–11 nodes (2–10 segments) per trajectory LibeamL_i^beam log-uniform[0.5,2.00.5,2.0] per segment kiaxk_i^ax log-uniform[0.5,2.00.5,2.0] per segment kirotk_i^rot log-uniform[0.5,2.00.5,2.0] per node ϕφ uniform on [0,2π0,2π) per trajectory The dataset is well-suited to our analysis since the buckling direction ϕφ is continuous, yielding infinitely many symmetry-related modes under the same loading condition. Examples Figure 9 visualizes representative 3D beam trajectories for increasing node counts, showing how the initially vertical beam progressively buckles into different directions over time. Figure 9: Example trajectories from the 3D beam dataset with different buckling directions ϕφ for diverse node counts over timestep t. Dataset Split Table 5 illustrates the distribution across different data splits. We consider a single sample per condition in the table, whereas during training the full range of buckling directions, ϕ∈[0,2π]φ∈[0,2π], is applied across the dataset. Table 5: Summary of 3D beam dataset splits. Split Trajectories Share (%) Training 697 70 Validation 150 15 Test 150 15 Total 997 100 C.2 Mechanical Metamaterials For the second problem setting, we use the publicly available dataset of Hendriks et al. [29]22 2 https://zenodo.org/records/18172947. It represents buckling trajectories of 2D porous elastic metamaterials, covering all 17 wallpaper groups. Each geometry is simulated on a 2×22× 2 representative volume element under 12 loading paths, giving 12,240 trajectories from 1,020 geometries, with mesh sizes between 1,000 - 36,000 per timestep. In our notation, the prescribed loading path is the control protocol ℓ ; for this dataset each increment is a known deformation gradient, ℓt=t ^t=F^t. The trajectories contain up to T=32T=32 timesteps, but can terminate earlier if the system runs into contact and reaches a steady state. For the underlying dataset the deformed configuration at timestep t is obtained by combining the macroscopic deformation with the displacement field tM^t: t=1(t)⊤+t. ^t=X^1(F^t) +M^t. (10) The deformation gradient is constructed by linearly interpolating from the identity matrix I to a prescribed terminal deformation gradient TF_T: t=+t−1T−1(T−),t∈1,…,T.F^t=I+ t-1T-1(F_T-I), t∈\1,…,T\. (11) The dataset is well-suited to our analysis because geometries can end in |ℳ|∈1,2,4|M_c|∈\1,2,4\ distinct post-buckling modes under the same condition. Before bifurcation, symmetry-related modes should remain close in latent space; after bifurcation, the model must represent distinct modes without losing their shared trajectory-level structure. Examples Figure 10 summarizes two examples of the metamaterials dataset with distinct modes. Figure 10(a) shows a trajectory with 2 modes. The starting configuration and the last timestep before buckling are displayed on the left in gray. Under progressive loading, the structure branches into two distinct post-buckling configurations: modes 1 and 2, which, in this case, are mirror reflections of one another. Both modes are shown at two successive timesteps, capturing how the transformed pattern becomes more pronounced with increasing load. (a) Example trajectory with one bifurcation point resulting in two modes. (b) Example trajectory with two bifurcation points resulting in four modes. Figure 10: Example mechanical metamaterials trajectories with distinct post-buckling mode structures. Figure 10(b) shows a trajectory with 4 modes. Starting from the undeformed configuration on the left, the structure first reaches a bifurcation point where two modes emerge. By moving toward the terminal deformation gradient TF_T, a second bifurcation point is reached, resulting in four distinct modes. This hierarchical branching behavior demonstrates how repeated symmetry-breaking events can generate increasingly varied deformed configurations from a single initial geometry. Dataset Split To prevent leakage between subsets, we split the data by geometry: each of the 1,020 geometries is assigned to exactly one of the training, validation, or test sets. All 12 loading paths of a geometry are kept together in the same subset, so no test geometry is seen during training. We provide a detailed overview of the training, validation, and test split in Table 6 and their mode distributions. Overall, the validation and test data share a similar number of modes across categories. Table 6: Summary of dataset splits. Split 1 mode 2 modes 4 modes Trajectories Share (%) Training 6,529 1,777 262 8,568 70 Validation 1,408 377 51 1,836 15 Test 1,391 354 91 1,836 15 Total 9,328 2,508 404 12,240 100 C.3 Allen-Cahn Dataset The Allen-Cahn dataset describes phase separation in a three-dimensional periodic domain Ω , see Figure 11. The state is a scalar phase field u(,t):Ω×[0,50]→ℝu(x,t) ×[0,50] governed by ∂u∂t=ϵ2∇2u−(u3−μu). ∂ u∂ t=ε^2∇^2u- (u^3-μ u ). (12) Here, ϵε sets the characteristic width and smoothness of the interfaces, and μ determines the local phase structure. For μ>0μ>0, small perturbations of the mixed state grow and form spatial domains approaching u=±μu=± μ. A small subset of the data uses μ<0μ<0, a non-separating regime in which the homogeneous state is stable. We discretize u on an Nx×Ny×NzN_x× N_y× N_z grid and denote the discrete state at timestep t by t∈ℝNx×Ny×Nzu^t ^N_x× N_y× N_z. Dataset generation We generated 1,0001,000 parameter conditions on a 64364^3 grid over the periodic domain [0,2π]3[0,2π]^3. For every condition, ϵε was sampled log-uniformly as ϵ=10zε=10^z with z∼(−1.5,−1.0)z (-1.5,-1.0), i.e. ϵ∈[10−1.5,10−1]≈[0.0316,0.1]ε∈[10^-1.5,10^-1]≈[0.0316,0.1]. Of the sampled μ values, 98%98\% were drawn uniformly from [0.05,1.0][0.05,1.0] and 2%2\% from [−0.1,0][-0.1,0]. Figure 11: Allen-Cahn example for two distinct ϵ,με,μ settings. All trajectories start from the same homogeneous state u0()≡0u^0(x)≡ 0. However, because this state is a valid solution but unstable, uniform random noise with an amplitude of 0.1 is added to the initial state in order to perturb it away from the unstable all-zero solution. The equation is integrated for 5,0005,000 steps with Δt=0.01 t=0.01, using a semi-implicit Fourier-spectral scheme: the nonlinear reaction term is evaluated explicitly, whereas the diffusion term is solved implicitly by pointwise division in Fourier space. Periodic boundary conditions are applied along all three axes. The initial state and every 100100th integration step are retained, giving T=51T=51 snapshots per trajectory over the time interval [0,50][0,50]. Consequently, the stored phase-field tensor before splitting has shape 4,000×51×64×64×644,000× 51× 64× 64× 64. Modes For Allen-Cahn the control protocol is constant, ℓt≡(ϵ,μ) ^t≡(ε,μ). For each sampled (ϵ,μ)(ε,μ) condition, four trajectories were simulated from the shared homogeneous initial state. These trajectories share the same physical parameters and provide samples of different phase-separation morphologies under the same condition c. In addition, the Allen-Cahn equation is invariant under the analytic transformation u↦−u -u. Applying this sign flip to each realization provides its phase-exchanged counterpart, yielding eight candidate modes per condition (four simulated realizations and their four sign-reversed counterparts). Dataset Split The data are split at the level of parameter conditions rather than individual trajectories. Thus, all four trajectories belonging to the same (ϵ,μ)(ε,μ) pair remain in the same subset, preventing information leakage between training and evaluation. The conditions were shuffled and divided into 70%70\% training, 15%15\% validation, and 15%15\% test data, as summarized in Table 7. Table 7: Summary of the Allen-Cahn dataset splits. Each condition contains four simulated trajectories and up to eight modes after analytic sign augmentation. Split Conditions Simulated trajectories Candidate modes Share (%) Training 700 2,800 5,600 70 Validation 150 600 1,200 15 Test 150 600 1,200 15 Total 1,000 4,000 8,000 100 Appendix D Experimental Details D.1 Architecture Autoencoder The encoder ℰϕencE_ _enc and decoder ϕdecD_ _dec form a Perceiver-style cross-attention autoencoder [33, 4, 51] for the 2D mechanical metamaterials and Beam3D datasets. The encoder maps N node points into a fixed, learned latent array of L tokens. The decoder reconstructs per-point deformations via cross-attention from positional queries. All nonlinear activations are GELU [31], and we additionally embed input and query positions via the sine-cosine position embedding from transformers [59]. A schematic of the architecture is provided in Figure 3; dataset-specific autoencoder hyperparameters are listed in Table 8. The Allen-Cahn Autoencoder uses a 3D Vision Transformer (ViT3D). The encoder patchifies the input voxel grid to 8×8×88× 8× 8 patches, which are run through self-attention/MLP transformer blocks to get a 4×4×44× 4× 4 latent representation. The decoder mirrors the encoder where the latent representation is upsampled to the original voxel grid. We ensure latent stability by normalizing the physical fields and conditioning parameters using dataset-specific statistics, while layer normalization is used throughout the attention-based encoders and latent flow network. This stabilizes the latent representations without introducing an additional variational regularization term. Flow Matching Model We employ a two-stage latent flow matching model. In the first stage, an encoder maps each physical state to a compact latent representation t=ℰϕenc(1,t),Z^t=E_ _enc(X^1,M^t), (13) while a decoder reconstructs the original state. The latent representation of a complete trajectory is 1:=(1,…,T)∈ℝT×L×Dz.Z_1:= (Z^1,…,Z^T ) ^T× L× D_z. (14) In the second stage, the autoencoder is frozen and a Scalable Interpolant Transformer (SiT) [41, 2] is trained to model the distribution of latent trajectories. A Gaussian sample ∼(0,I) ξ (0,I) is coupled to each data trajectory. At a randomly sampled flow time τ∈[0,1]τ∈[0,1], the interpolated latent state is τ=ατ1+στ,Z_τ= _τZ_1+ _τ ξ, (15) where Beam3D and the mechanical metamaterials use the generalized variance-preserving (GVP) [41] path with ατ=sin(12πτ) _τ= ( 12πτ) and στ=cos(12πτ) _τ= ( 12πτ), and Allen-Cahn uses the linear path with ατ=τ _τ=τ and στ=1−τ _τ=1-τ. Moreover, Eq. 5 in the main text shows the linear path. For a general path, the velocity field of Eq. 7 is v^θflow=α˙τ^1+σ˙τ(τ−ατ^1)/στ v_ _flow= α_τ Z_1+ σ_τ\,(Z_τ- _τ Z_1)/ _τ, which reduces to (^1−τ)/(1−τ)( Z_1-Z_τ)/(1-τ) for the linear path. The configured models directly estimate the clean latent trajectory, ^1=fθflow(τ,τ,1:Tobs,). Z_1=f_ _flow(Z_τ,τ,Z^1:T_obs,c). (16) The predicted latent trajectory is mapped back to the physical signal using the frozen decoder: ^t=ϕdec(^1t,1). M^t=D_ _dec ( Z_1^t,X^1 ). (17) D.2 Implementation Details Beam3D For Beam3D, each state contains the three-dimensional displacement field of the beam. The Perceiver encoder represents a state using six latent tokens of dimension 128128, producing a fixed-size representation that is independent of the number of mesh nodes. The flow model is conditioned on the prescribed scalar tip displacement and observes the initial trajectory frame in the standard configuration. Because the direction of beam buckling is rotationally symmetric, random azimuthal rotations are applied during training to expose the model to the complete SO(2) orbit of equivalent buckling directions. Mechanical Metamaterials Each mechanical metamaterial state is represented by a two-dimensional nodal displacement field together with the node positions, node types, and wallpaper-group information. A Perceiver autoencoder compresses each state into 128128 latent tokens of dimension 128128. The flow model is conditioned on the prescribed 2×22× 2 deformation gradient, represented as a four-dimensional vector. We group solution branches belonging to the same physical condition during training, and restrict optimal-transport coupling to samples within the same group. This encourages consistent matching between alternative bifurcation modes without coupling trajectories generated under different physical conditions. Allen–Cahn Each Allen–Cahn state is a scalar phase field defined on a periodic 64364^3 grid. A three-dimensional vision transformer divides the field into volumetric patches and encodes them into a structured latent grid. The flow model is conditioned on the physical parameters (ϵ,με,μ). Temporal and spatial rotary embeddings allow the transformer to distinguish both trajectory time and spatial latent-token position. The implementation supports periodic translations, sign inversion, axis reflections, and axis permutations as symmetry-preserving transformations of the Allen–Cahn system. D.3 Training Training is divided into two stages: (i) autoencoder training and (i) latent flow training. The first stage involves training the autoencoder to reconstruct individual physical states. The encoder ℰE and decoder D are jointly optimized to reconstruct individual physical states. Invalid or padded spatial elements are excluded from the reconstruction objective. The second stage involves training the latent flow model to learn the distribution of latent trajectories. The pretrained autoencoder is frozen, and only the SiT flow model and physical-parameter conditioner are optimized. Gaussian noise with the same shape as the latent trajectory is sampled and paired with the data using either independent coupling or a Hungarian optimal-transport assignment based on latent-space squared distance. When mode groups are available, the assignment is restricted to trajectories belonging to the same physical condition. For clean-trajectory prediction, the flow time is sampled according to τ∼(τmin,1−τmin),τmin=10−3τ ( _ ,1- _ ), _ =10^-3. The interpolated latent trajectory τZ_τ, flow time τ, control protocol, and observed-frame condition are then passed to the model. Invalid physical timesteps are replaced by zero before coupling and are excluded from both attention and loss computation. The standard conditional models observe the first trajectory frame. The Allen–Cahn transfer configuration uses unconditional generation by setting the number of observed frames to zero. Optimization The second-stage models are optimized with AdamW [40]. The standard learning rate is 10−410^-4, with cosine learning-rate decay for the mechanical metamaterial and standard Allen–Cahn models and a constant learning rate for Beam3D. Gradient norms are clipped at 11, and optimizer updates containing non-finite gradients are skipped. Training is step-based. Checkpoints are selected using the validation loss, with additional checkpoints maintained for generation metrics such as mode coverage. Data augmentation No generic spatial augmentation is applied to the mechanical metamaterial trajectories because their geometric, node-type, and wallpaper-group information is explicitly represented. Beam3D trajectories are randomly rotated around the longitudinal beam axis, reflecting the continuous rotational symmetry of the buckling problem. Allen–Cahn augmentation can combine periodic translations, sign inversion, axis reflection, and axis permutation. Each transformation is applied consistently to all frames of a trajectory to preserve its temporal dynamics. D.4 Loss Function The first-stage reconstruction model and the second-stage latent flow model are optimized separately. Reconstruction loss For a physical signal ∈ℝN×CM ^N× C and its reconstruction M, the first-stage reconstruction loss is the mean squared error over nodes: ℒrec=1N∑n=1N∥^n−n∥22.L_rec= 1N _n=1^N m_n-m_n _2^2. (18) Here, nm_n and ^n m_n denote the signal and its reconstruction at node n. Padded or invalid nodes, where present, are excluded from both the sum and the count N. For Metamaterials, n∈ℝ2m_n ^2 is the nodal displacement. For Beam3D, n∈ℝ3m_n ^3 is the three-dimensional nodal displacement. For Allen-Cahn, n∈ℝm_n is the scalar phase-field value. Latent flow loss The principal second-stage objective is the frame-masked mean-squared error between the predicted and target latent trajectories. With ^1=fθflow(τ,τ,1:Tobs,), Z_1=f_ _flow\! (Z_τ,τ,Z^1:T_obs,c ), (19) the loss is ℒflow=∑t=1Trt‖^1t−1t‖F2LDz∑t=1Trt.L_flow= _t=1^Tr^t \| Z_1^\,t-Z_1^\,t \|_F^2LD_z _t=1^Tr^t. (20) Here, τ denotes flow time and t indexes trajectory frames. The tensors ^1t,1t∈ℝL×Dz Z_1^\,t,Z_1^\,t ^L× D_z are the predicted and target clean latents at frame t, respectively. The indicator rt∈0,1r^t∈\0,1\ specifies whether frame t contributes to the loss. The denominator normalizes the squared error by the total number of latent entries across contributing frames. Bifurcation classification loss The Metamaterial and Beam3D models optionally include a temporal attention head that predicts whether each trajectory frame lies before or after bifurcation. The head is trained using masked binary cross-entropy: ℒbif=−∑t=1Trt[btlogpt+(1−bt)log(1−pt)]∑t=1Trt.L_bif=- _t=1^Tr^t [b^t p^t+(1-b^t) (1-p^t) ] _t=1^Tr^t. (21) Here, ptp^t is the predicted post-bifurcation probability, btb^t is the corresponding ground-truth target, and b^t=[pt≥1/2] b^t= 1[p^t≥ 1/2] is its binary prediction. The complete second-stage objective is ℒ=ℒflow+λbifℒbif.L=L_flow+ _bifL_bif. (22) The standard Allen–Cahn model does not use the auxiliary classification head and is optimized only using ℒflowL_flow. D.5 Adaptations to particle guidance In their notation, Corso et al. [14] add the gradient of a joint cross-particle potential to the reverse-diffusion drift. They evaluate this interaction on the noisy state xτx_τ. We make four adaptations for our latent data-prediction model. Denoised interaction We compute repulsion on the estimate ^1(k) Z_1^(k). It is restricted by the shared predicted mask b in Eq. 7, which is computed once before guidance begins and held fixed throughout guided sampling. The original method computes it directly on xτx_τ. We scale the resulting displacement by (στ2+δ)−1( _τ^2+δ)^-1 before adding it to the velocity update in Eq. 7; the constant δ>0δ>0 prevents divergence as στ→0 _τ→ 0. The interaction strength therefore follows the model’s denoising confidence. Normalized force The original force is an unnormalized gradient of a scalar potential. Our force is a kernel-weighted average of unit pairwise directions. We divide by the total kernel mass ∑jkij _jk_ij. This keeps its scale stable as the number and weight of neighboring particles change. Self-interaction For a translation-invariant kernel, the original self term has zero gradient and vanishes. Our self term is excluded from the numerator because its direction is undefined. It remains in the normalizing denominator and therefore damps the interaction. Frame and time gating The original potential acts on the full particle state with a generic coefficient. We apply repulsion only to the predicted bifurcation selected by ^=1 b=1. It is zero elsewhere. We use the bounded schedule γ(τ)=max(τκ,γmin)γ(τ)= (τ^κ, _ ). The schedule strengthens toward the data endpoint and never switches off. Here κ>0κ>0 is the schedule exponent and γmin∈(0,1] _ ∈(0,1] is its floor; wrepw_rep is the global repulsion weight. No additional per-frame weight is used beyond b. The bandwidth uses a median heuristic over squared, dimension-scaled pairwise distances. This choice follows the synthetic-experiment code released with Corso et al. [14]. Their paper leaves the kernel and bandwidth domain-specific and does not prescribe a general formula. We provide the latent-flow and guidance-schedule hyperparameters used for the reported models in Table 9. D.6 Hyperparameters Table 8: Hyperparameters of the autoencoder (first stage) used to produce the reported outputs, per dataset. Parameter Beam3D Metamaterials Allen-Cahn Architecture Perceiver AE Perceiver AE ViT3D + spatial token bottleneck Params 1.7×1061.7× 10^6 1.7×1061.7× 10^6 6.1×1066.1× 10^6 Input shape node states (N×CN× C) point set (17 groups, 8 node types) 64×64×64×164×64×64×1 voxel grid Patch size – – 8×8×88×8×8 (512 tokens) Latent shape (L×DzL× D_z) 6×1286× 128 128×128128× 128 4×4×44×4×4 (64×38464× 384) Encoder blocks (cross/self-attn) 1/ 11\,/\,1 1/ 11\,/\,1 ViT depth 66 (heads 6, embed 192) Decoder blocks (cross/self-attn) 1/ 11\,/\,1 1/ 11\,/\,1 ViT depth 44 (heads 6, embed 192) Head dims (enc cross/latent/dec cross/latent) 16/16/16/1616/16/16/16 16/16/16/1616/16/16/16 – Heads (enc cross/latent/dec cross/latent) 32/16/32/1632/16/32/16 32/16/32/1632/16/32/16 66 (enc) / 6/\,6 (dec) Decoder query width 160160 160160 – Activation GELU GELU GELU Dropout (latent / query) 0.0/ 0.050.0\,/\,0.05 0.0/ 0.050.0\,/\,0.05 0.00.0 (encoder) Loss masked MSE masked MSE masked MSE Optimizer AdamW AdamW AdamW Learning rate 10−410^-4 10−410^-4 2×10−42× 10^-4 Weight decay 0.00.0 0.010.01 0.010.01 Scheduler Cosine (eta_min 10−610^-6) Cosine (eta_min 1.5×10−51.5× 10^-5) none Gradient norm clipping 1.01.0 1.01.0 1.01.0 Batch size 128128 128128 128128 Precision bf16-mixed bf16-mixed bf16-mixed Table 9: Hyperparameters of the latent-space (second-stage) flow model used to produce the reported outputs, per dataset. Parameter Beam3D Metamaterials Allen-Cahn Params (non-AE) 7.1×1067.1× 10^6 17.2×10617.2× 10^6 31.8×10631.8× 10^6 in_dim / cond_dim 128/ 128128\,/\,128 128/ 128128\,/\,128 384/ 192384\,/\,192 Hidden size 256256 512512 512512 Depth 55 44 88 Heads 88 1616 1616 MLP ratio 2.02.0 2.02.0 2.02.0 Bifurcation head attention, λbif=0.1 _bif=0.1 attention, λbif=0.1 _bif=0.1 none Conditioner (in/dim/hidden) 1/ 128/ 5121\,/\,128\,/\,512 4/ 128/ 2564\,/\,128\,/\,256 2/ 192/ 5122\,/\,192\,/\,512 Transport path / target GVP / data GVP / data linear / data Sampler DDIM, η=0.0η=0.0 DDIM, η=0.0η=0.0 (eval η=1.0η=1.0) DDIM, η=0.0η=0.0 (eval η=1.0η=1.0) Loss flow MSE (frame-masked) flow MSE (frame-masked) flow MSE (frame-masked) Optimizer AdamW AdamW AdamW Learning rate 10−410^-4 10−410^-4 3×10−53× 10^-5 Weight decay 0.00.0 0.010.01 0.010.01 Scheduler none Cosine (eta_min 10−610^-6) Cosine (Tmax=500T_max=500k, eta_min 10−510^-5) Batch size 6464 1616 88 Precision bf16-mixed bf16-mixed fp32 PG settings (wvelw_vel/wrepw_rep/κ) 1.0/ 1.0/ 0.01.0\,/\,1.0\,/\,0.0 1.0/ 50/ 0.251.0\,/\,50\,/\,0.25 1.0/ 75/ 0.0251.0\,/\,75\,/\,0.025 D.7 Evaluation Details: Beam3D and Metamaterials We apply a mode-aware evaluation pipeline with task-specific adaptations to the underlying mode set. For discrete-mode settings, we match generated rollouts to their nearest ground-truth mode using a nearest-mode distance. For metamaterials, this uses an absolute distance threshold of 0.10.1 over the whole trajectory. We compute the mean absolute error over the whole trajectory and sample 1616 trajectories per condition, using 5050 integration steps. For both Beam3D and Allen-Cahn, the mode space is continuous. Rollouts for the Beam3D dataset are therefore evaluated by binning the buckling azimuth and matching against the nearest orbit member with a relative-distance rejection criterion. We reject a generated trajectory if its relative distance to the matched orbit member exceeds 50%50\% of its ground-truth displacement magnitude. Here we sample 180180 trials with fewer integration steps (1010 steps, DDIM [53]), and compare against GeoTDM (10001000 steps, DDPM [32]) and STFlow (1010 steps, Euler integration) using the same number of evaluation samples. D.8 Evaluation Details: Allen-Cahn The Allen-Cahn dataset dynamics are invariant under periodic translations, sign flips and point-group symmetries of the cubic domain, making pointwise comparison to a single reference rollout not meaningful. We evaluate in two steps: (i) each generated sample is matched to a ground-truth branch using canonical alignment, (i) and three different JSD metrics test distributional uniformity over the possible orbit operations. Matching For a generated field M, we compute the time-averaged radial power spectrum S()S(M) as the magnitude of the 3D FFT at each timestep, binned by integer radial frequency and averaged over time. This value is invariant to translation and the octahedral point group, making it suitable for identifying the best mode independent of orientation and shift. The matched branch is m⋆(k)=argmaxmcos_sim(S((k)),S((m)))m_ ^(k)= *arg\,max_mcos\_sim (S(M^(k)),\,S(R^(m)) ), where (m)R^(m) is the selected ground-truth reference for branch m. We count sign-inverted counterparts separately for mode-coverage JSD, but disregard sign for translation and rotoreflection JSD. Given m⋆(k)m_ ^(k), we recover translation and orientation by an exhaustive search over the 48-element octahedral group GoctG_oct (axis reflections and permutations). For each g∈Goctg∈ G_oct, we apply g to the sample’s mean field and phase-correlate it against the matched mode’s mean field using the normalized cross-power spectrum. The peak location gives the best periodic shift, and the height gives an alignment score. We only keep the g⋆g_ with the highest peak, resulting in one match result (m⋆(k),⋆(k),g⋆(k))(m_ ^(k), ψ_ ^(k),g_ ^(k)) per sample, where ⋆(k)∈ℤ643 ψ_ ^(k) _64^3 is the optimal periodic grid shift. The approach mirrors the alignment step in [35]. In contrast to that work, and because the symmetry group is small and finite, the exact minimum is found directly rather than approximated. MCon-MAE For each generated rollout ^(k) M^(k), let ^rel,(k) M^rel,(k) be its Gauss–Newton relaxation onto an Allen-Cahn PDE-consistent trajectory. We compute MCon-MAE=1K∑k=1Kmean(|^rel,(k)−^(k)|),MCon-MAE= 1K _k=1^Kmean\! ( | M^rel,(k)- M^(k) | ), (23) This dataset-specific metric compares each raw generation with its own relaxed projection; no saved reference rollout, matched branch, or branch label enters the calculation. It is therefore distinct from the nearest-reference MCon-MAE in Eq. 8. JSD Metrics Each JSD metric evaluates one component of the matching approach above over K trials. The results are then aggregated into a categorical distribution, and we compute the Jensen-Shannon divergence against the uniform distribution, as summarized in Table 10. Sampling. To estimate the distribution over the large number of potential solutions, we sample 100100 times for each condition c, using 3030 denoising steps. Table 10: Summary of the Allen-Cahn JSD metrics, including their categories and tests. Metric Categories Tests Mode Coverage JSD mapping to ℳM_c ground truth branches Are all solution branches reached equally often? Translation JSD grid shifts, coarsened to 8 bins per axis (83=5128^3=512) Given one ground truth sample, how often are periodic shifts recovered? Rotoreflection JSD 48 octahedral elements Given one ground truth sample, is the recovered orientation/reflection uniform? Appendix E Results E.1 Quantitative results We present a summary of our results for all three datasets in the tables below. We divide the results into two sections: an overview of the JSD and the MCon-MAE results. Section 3 explains additional metrics, which we include for each dataset. For the JSD tables, we also incorporate upper and lower bounds as references. The upper bound represents total mode collapse, while the lower bound reflects the expected JSD when sampling from a random uniform distribution. For all tables, we used five different evaluation seeds. Beam3D Table 11: Beam3D prediction performance. Lower MCon-MAE, JSD, and rejection rate indicate better reconstruction accuracy, better agreement with the target mode distribution, and fewer physically invalid rollouts. (mean ± std over five evaluation seeds) (a) JSD Method JSD GeoTDM 0.036±0.0010.036± 0.001 STFlow 0.074±0.0010.074± 0.001 Bi-FORK w/o PG 0.043±0.0030.043± 0.003 Bi-FORK 0.041±0.0010.041± 0.001 Upper bound∗ 0.644 Lower bound∗ 0.040 ∗ assuming total mode collapse, ∗ expected JSD when sampling from a random uniform distribution (b) MCon-MAE and rejection rate Method MCon-MAE Rejection GeoTDM 1.040±0.0011.040± 0.001 99.2±0.099.2± 0.0 STFlow 0.114±0.0000.114± 0.000 0.5±0.10.5± 0.1 Bi-FORK w/o PG 0.018±0.0010.018± 0.001 0.1±0.10.1± 0.1 Bi-FORK 0.018±0.0000.018± 0.000 <0.01<0.01 Metamaterials For the dataset Bi-FORK achieves JSD values below the lower-bound reference. This reference represents the expected JSD for a finite set of independent samples from the uniform target distribution, rather than a strict lower bound. Particle guidance can create dependence between samples that promotes low-discrepancy sampling, leading to more balanced mode frequencies and JSD values below the reference. Full coverage requires a valid sample for every mode; JSD includes all generated samples. Table 12: Result comparison of Bi-FORK (including w/o particle guidance) against the affine baseline, evaluated with MCon-MAE, JSD, and coverage for each mode count. (mean ± std over five evaluation seeds) (a) JSD Method JSD |ℳ|=1|M_c|=1 |ℳ|=2|M_c|=2 |ℳ|=4|M_c|=4 Affine Baseline 0.0000.000 0.2160.216 0.3800.380 Bi-FORK w/o PG 0.000±0.0000.000± 0.000 0.050±0.0020.050± 0.002 0.087±0.0080.087± 0.008 Bi-FORK 0.000±0.0000.000± 0.000 0.013±0.0010.013± 0.001 0.023±0.0020.023± 0.002 Upper bound∗ 0 0.216 0.380 Lower bound∗ 0 0.008 0.027 ∗ assuming total mode collapse, ∗ expected JSD when sampling from a random uniform distribution (b) MCon-MAE and full coverage Method MCon-MAE Full Coverage % |ℳ|=2|M_c|=2 |ℳ|=4|M_c|=4 Affine Baseline 0.0210.021 0 0 Bi-FORK w/o PG 0.016±0.0000.016± 0.000 38±0.938± 0.9 13±3.413± 3.4 Bi-FORK 0.018±0.0000.018± 0.000 74±1.074± 1.0 48±3.748± 3.7 Allen-Cahn Table 13: Result comparison of Bi-FORK (including w/o particle guidance), evaluated with relaxed MCon-MAE and three JSD metrics. (mean ± std over five evaluation seeds) (a) JSD Method Mean JSD JSD Mode-Coverage JSD Translation JSD Rotoreflection JSD Bi-FORK w/o PG 0.217±0.0040.217± 0.004 0.131±0.0100.131± 0.010 0.446±0.0000.446± 0.000 0.076±0.0010.076± 0.001 Bi-FORK 0.204±0.0010.204± 0.001 0.089±0.0010.089± 0.001 0.446±0.0000.446± 0.000 0.076±0.0020.076± 0.002 Upper bound∗ - 0.4970.497 0.6860.686 0.6420.642 Lower bound∗ - 0.0090.009 0.4460.446 0.0760.076 ∗ assuming total mode collapse, ∗ expected JSD when sampling from a random uniform distribution (b) MCon-MAE Method MCon-MAE Bi-FORK w/o PG 0.0515±0.0010.0515± 0.001 Bi-FORK 0.0453±0.0000.0453± 0.000 E.2 Ablations Bifurcation Head Prediction Accuracy To effectively use particle guidance, we predict a bifurcation mask corresponding to the bifurcation time point. We assess the accuracy of the predicted masks on the Beam3D and mechanical metamaterials datasets. The Allen-Cahn method does not use a bifurcation head. Instead, we set the guidance mask to b^t=1 b^t=1 for every valid frame, allowing guidance throughout the trajectory. We compare the results based on accuracy, precision, recall, F1 score, and balanced accuracy, as shown in Table 14. Table 14: Bifurcation-head classification on the complete test splits. Reporting the mean and standard deviation over five seeds. Dataset Accuracy Precision Recall F1 Balanced acc. Beam3D 0.989±0.0000.989± 0.000 0.988±0.0000.988± 0.000 1.000±0.0001.000± 0.000 0.994±0.0000.994± 0.000 0.950±0.0010.950± 0.001 Metamaterials 0.970±0.0000.970± 0.000 0.940±0.0020.940± 0.002 0.871±0.0020.871± 0.002 0.904±0.0020.904± 0.002 0.930±0.0010.930± 0.001 Wall Time We analyze the wall time of each tested method on the Beam3D dataset. For the underlying experiments, we used a B300 SXM6 with 288 GB evaluated on one input condition over four random draws. We summarize the resulting times in Table 15. Table 15: Wall-clock time taken to generate 180 samples of the Beam3D dataset over four random seeds. Method Samples Conditions Seconds Taken GeoTDM 180 1 7766.367766.36 STFlow 180 1 72.3272.32 Bi-FORK w/o PG 180 1 2.942.94 Bi-FORK 180 1 3.953.95 E.3 Qualitative Results Beam3D Figure 12: Beam bending process: showing predicted (green) vs actual (black) trajectory evolution for 180 samples. Node count ranging from 33 to 1111 nodes Allen-Cahn Figure 13: Sample comparison for increasing μ and ϵε, showcasing the diversity of potential solutions compared to the ground truth. Mechanical Metamaterials Figure 14: Visualization of one-mode samples generated by Bi-FORK with K=16K=16. The ground truth is compared with the closest predicted mode. Figure 15: Visualization of two-mode samples generated by Bi-FORK with K=16K=16. The ground truth is compared with the closest predicted mode. Figure 16: Visualization of a four-mode sample generated by Bi-FORK with K=16K=16. The ground truth is compared with the closest predicted mode.