arXiv CSopen access

Agentic Symbolic Search: Characterizing PDEs Beyond Hand-crafted Expressions, Meshes, and Neural Networks

Agentic Symbolic Search: Characterizing PDEs Beyond Hand-crafted Expressions, Meshes, and Neural Networks Liu Yang∗ National University of Singapore [email protected]

arXiv:2606.20467v1 [cs.LG] 18 Jun 2026

Zongmin Yu National University of Singapore [email protected]

Abstract Mathematicians understand a PDE solution through mathematical structures rather than tables of computed values. Historically, this has been the product of mathematical analysis, carried out by hand for each problem individually. Neither numerical simulation nor neural networks produce those structures directly. We propose Agentic Symbolic Search (ASYS), a prior-guided framework in which an agent translates PDE theory, public problem constraints, and accumulated search experience into testable differentiable symbolic programs. The mathematical forms are refined under evolutionary search, while their continuous parameters are fit by gradient-based optimization. This makes the search an automated form of inductive-bias injection rather than blind symbolic regression. For problems with known analytical forms, ASYS recovers these forms naturally; for other problems, ASYS constructs analytical approximations which can guide mathematicians toward further analysis. In our experiments, across five problems spanning bounded dynamics, finite-time blow-up, and free-boundary focusing, ASYS produces interpretable representations, including a geometric interface formula for Allen–Cahn 2D dynamics and a nine-parameter contraction law for Keller–Segel chemotactic blow-up, in settings where no closed-form description was previously available. ASYS shows the possibility of a new paradigm for characterizing PDE solutions, beyond handcrafted analytical solutions, mesh-based numerical solutions, and neural network approximations.

1

Introduction

Partial differential equations (PDEs) constitute the core mathematical language for describing the dynamic evolution of the natural world. In the long history of science, the highest pursuit of understanding physical systems has always been the search for analytical solutions. From d’Alembert’s traveling-wave solution to Fourier’s expansion of the heat equation (Fourier, 1822), the analytical solution remains the “gold standard” of mathematical understanding precisely because of its readability. It directly and elegantly exposes how a system depends on parameters, initial data, and geometry, making asymptotic behavior, scaling laws, and singularities immediately visible. However, structural obstructions prevent the vast majority of nonlinear systems from admitting closed-form analytical solutions, severely limiting the universality of this paradigm. The advent of the electronic computer catalyzed the birth of the numerical solution paradigm. Beginning with the first numerical weather forecasts on ENIAC (Charney et al., 1950), techniques such as finite-difference, finite-element, and spectral methods made it possible to compute equations beyond the reach of analytic solution, thereby inaugurating an era of empirical discovery in computational science. Yet, the output of a numerical simulation is intrinsically a table of values, not an explanation. It merely records “what happened,” leaving the core mathematical structures that explain “why,” such as scaling laws and similarity exponents (Barenblatt, 1996; Choptuik, 1993), buried within massive amounts of discrete data, failing to automate the leap from “computation” to “understanding”. ∗ Corresponding author.

1

A Preprint

Agentic Symbolic Search

In recent years, the rise of deep learning has spawned a new generation of methods, such as physics-informed neural networks (PINNs) (Raissi et al., 2019) and neural operators (Lu et al., 2021; Li et al., 2021), which attempt to bridge this gap by replacing the computational grid with a trainable continuous function. Although these methods dispense with mesh generation and demonstrate significant potential in high-dimensional problems, they invariably lock the physical system’s solution within the opaque weights of a neural network. Crucial asymptotic forms, scaling exponents, and similarity structures remain implicit in the weights rather than explicit in the representation. In essence, the neural era removed the grid but reintroduced the very gap it set out to close: a network, like a grid, yields continuous values rather than explicit structure. Beyond opacity, neural solvers face a fundamental limitation regarding universal approximation: the theoretical guarantees that motivate them do not apply to the finite-width architectures actually deployed. The universal approximation theorems establish density for the family of feedforward networks taken over arbitrary width (Cybenko, 1989; Hornik, 1991; Pinkus, 1999); once an architecture is fixed, the realizable functions form a restricted class whose expressive power is provably bounded (Telgarsky, 2016; Barron, 1994; Lu et al., 2017), and training only adjusts parameters within that class. Finite-time blow-up makes the gap concrete: a network assembled from bounded activations is itself bounded, and cannot represent a solution that diverges in finite time, however its parameters are tuned. The obstruction is one of representation class, not of optimization. Classical approximation theory anticipates the remedy: adapting the form of the approximant in a problem-dependent way, rather than fitting parameters within a fixed basis, attains what no fixed linear space can (DeVore, 1998). A critical challenge remains in identifying a mechanism capable of performing this structural adaptation. Such a mechanism must search over mathematical forms rather than parameter values, and it must know where to look: for PDEs, that knowledge lives in the literature of profiles, scalings, and similarity structures, not in any fixed basis. Concurrently, rapid progress in Large Language Models (LLMs) is opening new avenues for scientific computation. Recent works using LLMs as optimizers and coding agents (e.g., FunSearch and AlphaEvolve (Romera-Paredes et al., 2024; Novikov et al., 2025)) have shown that language models, when coupled with evaluators, can propose algorithms and mathematical constructs that go beyond fixed expression grammars (Shojaee et al., 2025; Xia et al., 2025). Because these models operate over language and code, they can also act on mathematical knowledge expressed in natural language, from known singular profiles to scaling arguments, and turn it into executable hypotheses. This should not be confused with blind discovery from an empty prior. In ASYS, access to literature-level mathematical knowledge is the mechanism by which prior structure enters the search, just as a human applied mathematician begins from profiles, scalings, and asymptotic regimes rather than an unconstrained expression grammar. Along this line, agentic systems have recently been applied to PDE solving itself: ATHENA and its successor GRAFT-ATHENA (Toscano et al., 2025; 2026) deploy multi-agent teams that autonomously construct and debug numerical solvers and neural surrogates, while a self-evolving agent has been shown to discover interpretable fluid-control policies directly from physical simulation (Sun et al., 2026). These systems produce executable artifacts; ASYS differs in that it searches for the symbolic representation of the PDE solution itself, whose mathematical structure is the deliverable. The code generation and structural exploration capabilities of large language models provide an opportunity to address the interpretability and expressivity limitations described above. Rather than fitting parameters within a fixed network representation, one can use a coding agent to propose and revise the solution’s representation within a broader space of mathematical structures. We pursue this idea and propose Agentic Symbolic Search (ASYS). This framework employs a coding agent to translate prior knowledge from PDE theory and physical constraints, supplied in natural language, into differentiable symbolic programs. Fixed architectures admit no such channel: a network absorbs prior structure only insofar as an expert hand-encodes it into the architecture or the loss. By doing so, ASYS recovers at the level of the search the expressive freedom that universal approximation promises only at the level of the family, and automates the passage from numerical fitting to interpretable structure that previously required a human expert to supply the ansatz (Wang et al., 2023). The goal is not amortized zero-shot inference over many initial conditions, where neural operators are the natural tool, but structural characterization of a difficult evolution. Table 1 places ASYS among the preceding paradigms: explicit mathematical structure, previously obtainable only by a mathematician working case by case, is here produced by an automated search.

2

A Preprint

Agentic Symbolic Search

Table 1: Four paradigms for representing a PDE solution. The comparison concerns what object is returned, whether its mathematical structure is explicit, and how that representation is obtained, not an accuracy ordering. Explicit structure has historically required a human expert working case by case; ASYS is the only paradigm that produces it through an automated search.

Returned object

Analytical

Numerical solver

Neural network

ASYS

closed-form expression

values on a mesh

network weights

symbolic program

Mathematical structure

explicit

absent

implicit

explicit

Obtained by

a mathematician, case by case

a fixed numerical scheme

fitting a fixed architecture

an automated structural search

This comparison is therefore about representation, not about replacing every PDE solver. Neural operators and related surrogates are designed to learn maps across families of inputs. ASYS instead spends computation on one challenging trajectory in order to return an interpretable mechanism, such as a scaling law, interface geometry, or self-similar coordinate structure, that a mathematician can inspect.

2

Method

ASYS takes as input a known PDE together with its initial data, boundary conditions, and problem constants. The output is a differentiable symbolic program whose mathematical structure can be read directly. The method operates through two nested loops with strictly decoupled optimization objectives. The outer loop is scientifically a hypothesis generator over an ansatz space; the fact that each hypothesis is implemented as code is an execution detail. This loop proposes and revises the mathematical form of the ansatz: the coordinate transformations, the decomposition into branches, and the choice of trainable parameters; its progression is guided by a fixed, multi-dimensional evaluator score that measures mathematical correctness. An inner loop fits the continuous parameters of each proposed ansatz by quasi-Newton optimization under a fixed time budget, driven by an agent-defined scalar training loss which the agent is free to engineer. As summarized in Table 1, a classical solver (Courant et al., 1928; Charney et al., 1950) returns values on a mesh, a PINN (Raissi et al., 2019) or neural operator (Lu et al., 2021; Li et al., 2021) moves the solution into network weights, and ASYS makes the representation itself the search target. 2.1

Differentiable Symbolic Programs

Each representation is a differentiable program: a composition of coordinate transformations, analytic building blocks, piecewise branches, and trainable parameters 𝜃, written so that automatic differentiation can compute 𝜕𝑢 𝜃 /𝜕𝜃 for gradient-based fitting. The search space is therefore discrete at the level of program structure and continuous at the level of parameter values. This is broader than a symbolic-regression expression tree (which cannot naturally represent multi-branch decompositions or learned coordinate rescalings) but narrower than arbitrary executable code: programs must remain closed-form differentiable expressions that do not embed a numerical solver or access the reference solution. 2.2

Iterative Ansatz Refinement

Given this search space, ASYS iteratively refines the representation through an outer loop over program structure. At each iteration 𝑘, the agent receives three kinds of context: problem guidance, the highest-scoring previous representations, and a diagnostic summary of where the current best ansatz succeeds or fails. The guidance states the PDE, the admissible ansatz class, and relevant mathematical background, such as known profiles or scaling structures, in natural language. From this context the agent proposes a revised program: a new mathematical form, a custom scalar training loss for fitting its parameters, and optimization settings. Because of the decoupled objectives introduced above, the agent controls the structure and the formulation of the training objective, but not the gradient steps or the evaluation criteria. The evaluator then fits the trainable parameters 𝜃 by L-BFGS (Nocedal, 1980) under a fixed wall-clock budget and scores the result using its own fixed, independent criteria (Section 2.3). Because these symbolic programs may contain smooth gates, branch blends, singular scalings, and high-order nonlinear terms, the fitted loss landscape is generally non-convex and can become ill-conditioned near singularities. ASYS does not assume that the inner loop finds a global optimum. Instead, the outer loop perturbs the hypothesis family itself, giving the parameter optimizer new basins to explore and exposing when the limiting factor is the representation rather than the continuous fit. 3

A Preprint

Agentic Symbolic Search

Figure 1: ASYS alternates between an outer loop over mathematical structure and an inner loop over continuous parameters. The agent reads problem guidance, previous high-scoring representations, and diagnostic summaries, then revises the representation, the training objective, and the optimizer configuration. The evaluator fits the resulting ansatz by quasi-Newton optimization and scores it using PDE-residual and public-constraint diagnostics. Note the separation of objectives: the inner loop optimizes the agent’s training loss, while the outer loop selects based on the evaluator’s fixed score.

ASYS runs the outer loop inside Evolutionary Ensemble of Agents (EvE) (Yu and Liu, 2026), an evolutionary framework that maintains a population of scored candidates (Figure 2). At each iteration, new programs are generated from problem guidance, diagnostic feedback, and examples drawn from stronger prior representations. After evaluation, the new programs are added to the population and can influence later iterations. Selection is therefore gradual: the search is biased toward representations that perform well across the score dimensions, without committing to a single best program or a fixed scalar objective. The same feedback loop also adapts the high-level search strategies used to propose subsequent representations. Appendix A.1 gives a self-contained overview of this outer-loop framework. 2.3

Scoring

The agent controls both the mathematical form of the representation and the scalar training loss used to fit its parameters, which may combine PDE residuals, initial and boundary terms, and regularization in any proportion. The score used to compare representations across iterations, however, is structurally separate and fixed by the evaluator. This separation of objectives allows the agent to freely engineer different training landscapes while the evaluation standard remains absolute. Evaluation in ASYS is driven entirely by the governing equation and its publicly specified constraints, strictly excluding any observed solution data. To rigorously assess each candidate, the evaluator constructs a four-dimensional score vector:
𝑠(𝑢 𝜃 ) = 𝑠phys , 𝑠ic , 𝑠bc , 𝑠comp . The first three dimensions correspond to the standard constraints of a well-posed PDE problem: 4

A Preprint

Agentic Symbolic Search

Figure 2: ASYS uses EvE as its outer-loop search mechanism. EvE maintains the pool of representations, samples stronger ones more often as later context, and drives selection; ASYS specializes this framework to differentiable symbolic programs for PDE solution characterization.

Physics residual The physics dimension measures equation satisfaction through the PDE residual evaluated on held-out collocation points. The specific residual metric is adapted to the scale structure of each PDE: a raw mean-squared residual for bounded problems, and a scale-free geometric-mean relative residual for problems with singular dynamics. Definitions are given in Appendix A.4. Initial value This dimension evaluates the spatial matching at 𝑡 = 0, defined as the mean squared error between the representation and the analytic initial data: ℓic =

2 1 ∑︁ 𝑢 𝜃 (𝑥𝑖 , 0) − 𝑢 0 (𝑥𝑖 ) . 𝑁ic 𝑖

Boundary consistency The boundary dimension is the squared residual of the prescribed boundary condition (Dirichlet, Neumann, or periodic, depending on the problem). Collocation counts and per-case boundary definitions are collected in Appendix A.4. While standard physics-informed frameworks typically rely solely on the triplet above, matching the initial spatial profile 𝑢 𝜃 (𝑥, 0) = 𝑢 0 (𝑥) and satisfying the boundary conditions does not uniquely determine the direction of temporal evolution. An ansatz that is constant in 𝑡 or trapped in a trivial branch such as 𝑢 ≡ 0 can satisfy the initial data exactly while departing along an incorrect solution trajectory. To address this critical structural obstruction, ASYS introduces a fourth dimension: Compatibility condition This score evaluates the kinematic consistency of the representation at the initial time. In our experiments, we implement this constraint by comparing a one-sided finite-difference time derivative against the PDE right-hand side applied to the initial data:  2 1 ∑︁ 𝑢 𝜃 (𝑥 𝑖 , Δ𝑡) − 𝑢 𝜃 (𝑥𝑖 , 0) ℓcomp = − 𝐹 [𝑢 0 ] (𝑥𝑖 ) . 𝑁ic 𝑖 Δ𝑡

5

A Preprint

Agentic Symbolic Search

This dimension acts as a vital phase-space constraint. By explicitly penalizing discrepancies in the initial derivative, it enforces that the solution leaves the initial manifold in a direction strictly compatible with the PDE vector field. 2.4

Selection and Offline Validation

The four-dimensional score drives the selection of examples shown to the agent during search. Final solution accuracy is evaluated separately after the run against independent numerical references. To select which previous representations appear as examples in later iterations, the population is ranked independently in each score dimension. Each representation receives an exponentially decaying weight per dimension based on its rank (best rank receives the highest weight), and the weights are summed across the four dimensions to produce a single sampling probability. Representations are then drawn without replacement, so stronger ones across any combination of dimensions are more likely to appear as context. No representation is ever discarded from the population. Appendix A.5 gives the precise weighting formula. Relative 𝐿 2 against independent numerical references (Appendix A.8) enters only as an offline validation check of whether the resulting representation tracks the solution, not as part of the score that guides the search. 2.5

Computational Cost

All experiments were run on a single Apple M4 laptop CPU with no GPU acceleration. Each representation is fitted by L-BFGS with a wall-clock budget of 30 seconds for bounded cases and 300 seconds for blow-up and free-boundary cases, with early stopping. Bounded problems run for five outer-loop iterations; blow-up and free-boundary cases run for ten. The dominant cost is the outer-loop structural search, not large-scale numerical training. Agent-side sessions were executed through a subscription-covered interactive coding-agent interface rather than a per-token API workflow; the corresponding billing convention is documented with the run configuration in Appendix A.2.

3

Results and Discussion

We evaluate ASYS on five PDE problems arranged by structural difficulty (Table 2). Each run is initialized from a minimal candidate and iteratively replaces it with increasingly structured representations. Evaluation protocol. During structural evolution and inner-loop parameter fitting, candidates are scored only by public mathematical constraints: the governing PDE residual, prescribed initial and boundary conditions, and the compatibility condition (Appendix A.4). No numerical reference solution, observed trajectory, or reference-derived metric is exposed to the search. After the search, we compute validation relative 𝐿 2 errors against independently generated numerical references that remain hidden from the coding agent throughout evolution (Appendix A.8). These offline validation errors determine the ASYS entries in Table 2 and the orange trajectories in the dual-axis figures. Wherever applicable, a dotted horizontal line denotes the validation error of a specialized self-similar baseline on the same right-hand axis. Specifically, for singular problems exhibiting a known similarity structure, we implement a self-similar PINN (SS-PINN) baseline adapted from Wang et al. (2025); full architecture, optimization, and profile-equation details are given in Appendix A.9. Table 2 summarizes the trainable parameter counts for both ASYS and SS-PINN. Crucially, for the SS-PINN baseline, this count reflects solely the weights within the profile-approximating MLP; the prescribed coordinate transformation is treated as a rigid mathematical prior rather than a learnable parameter.

6

A Preprint

Agentic Symbolic Search

Table 2: Reported validation relative 𝐿 2 errors and trainable parameter counts across the five PDE cases. All errors are computed offline against hidden independent numerical references (Appendix A.8). SS-PINN denotes a specialized self-similar coordinate transformation plus MLP, reported only where the known similarity variables define the baseline. A dash indicates that no SS-PINN baseline is reported for that case. SS-PINN parameter counts include trainable profile-MLP weights only; the hand-specified coordinate transform is not counted as trainable. ASYS is this method. Case

Type

NLS Allen–Cahn 2D Keller–Segel Graveleau gCLM

bounded bounded blow-up free boundary blow-up

3.1

SS-PINN

SS-PINN 𝑛params

ASYS

ASYS 𝑛params

– – 0.258 – 0.196

– – 49,921 – 49,921

0.0059 0.0107 0.188 0.00132 0.465

5 23 9 13,133 13,060

Bounded Demonstrations: NLS and Allen–Cahn 2D

The two bounded cases establish the basic mechanism before the harder singularity problems by testing whether the search can recover an explicit mathematical form when the solution remains 𝑂 (1) and admits compact symbolic structure. The nonlinear Schrödinger equation (benchmark setup following Raissi et al. (2019)),

NLS

1 𝑖𝜓𝑡 + 𝜓 𝑥 𝑥 + |𝜓| 2 𝜓 = 0, 2

𝑥 ∈ [−5, 5],

𝑡 ∈ [0, 𝜋/2],

has periodic boundary conditions and initial condition 𝜓(𝑥, 0) = 2 sech(𝑥). This initial condition is the classical 𝑁=2 higher-order soliton (Satsuma and Yajima, 1974), whose exact breather form is known. In this run, neural-network layers are excluded and the search targets purely symbolic programs. ASYS recovers the exact breather core and adds a five-parameter split endpoint correction for the periodic domain. The best representation has the form h i
2 𝜓(𝑥, 𝑡) = 𝜓core (𝑥, 𝑡) + 𝑎 𝜕𝑥 𝜓core (5, 𝑡) 𝑇3 (𝑥; 𝑚 3 ) − 𝑐 𝑇3 (𝑥; 𝑚 3 ) − 𝑇5 (𝑥; 𝑚 5 ) 𝑒 −𝑟𝑡 , where 𝑎, 𝑐, 𝑟, 𝑚 3 , and 𝑚 5 are trainable scalars. The closed-form core is
4𝑒 𝑖𝑡/2 cosh(3𝑥) + 3 cos(4𝑡) cosh 𝑥 + 3𝑖 sin(4𝑡) cosh 𝑥 𝜓core = . cosh(4𝑥) + 4 cosh(2𝑥) + 3 cos(4𝑡) The core is the Satsuma–Yajima breather. The endpoint correction uses a broader primary tail 𝑇3 and a sharper initial tail 𝑇5 , with exponential relaxation in time so the initial profile is preserved while the periodic derivative mismatch is flattened. The physics loss and the validation 𝐿 2 decrease together across iterations (Figure 3b), reaching relative 𝐿 2 = 0.00585 within five iterations. t = 0.00 4

t = 1.57

t = 0.79

candidate reference

|ψ(x, t)|

3

2

1

0 −4

−2

0

2

4

−4

−2

0

x

2

4

−4

−2

x

0

x

(a) Breather snapshots for | 𝜓 ( 𝑥, 𝑡 ) | at three times. Solid curves are the 𝐿 2 -best candidate and dashed curves are the independent reference.

Figure 3: NLS: breather profile snapshots for the 𝐿 2 -best symbolic candidate.

7

2

4

A Preprint

Agentic Symbolic Search

NLS physics loss (raw MSE) best-so-far validation L 2

100

physics loss (raw MSE)

L 2-best candidate

10−2

10−1

10−3

10−4 10−2 10

best-so-far validation L 2

10−1

−5

0

1

2

3

4

5

outer-loop iteration (b) Dual-axis trajectory: physics loss and best-so-far validation 𝐿 2 decrease together.

Figure 3: NLS (continued): residual scoring reaches a purely symbolic soliton-form candidate, and offline validation confirms the resulting profile.

While NLS is a case of analytic recall, Allen–Cahn 2D requires geometric construction: there is no classical formula for the evolving interface. Allen–Cahn 2D

The two-dimensional Allen–Cahn equation (Allen and Cahn, 1979), 𝑢 𝑡 = 𝜀 2 Δ𝑢 + 𝑢 − 𝑢 3 ,

(𝑥, 𝑦) ∈ [−1, 1) 2 ,

𝑡 ∈ [0, 8],

𝜀 = 0.04,

starts from two overlapping diffuse disks whose union forms an asymmetric peanut-shaped region. As in the NLS case, neural-network layers are excluded. The phenomenon is geometric: Figure 4 shows a connected peanut-shaped positive phase whose neck smooths and whose outer boundary relaxes toward an oval. ASYS builds this peanut-to-oval dynamics from a geometric signed-distance formula: an early smooth union of two moving diffuse lobes, time-dependent curvature shrinkage, center relaxation toward the merged state, a sigmoid-gated blend to a rotating oval, a localized neck-fill distance shift, memory decay that preserves the exact initial phase field, and an initial right-hand-side jet. In compact form,  
(1 − 𝑏(𝑡))𝑆union + 𝑏(𝑡)𝑆oval + 𝑆neck 𝑢(𝑥, 𝑦, 𝑡) = tanh + 𝑀 (𝑡) 𝑢 0 − 𝑢 geom (0) + 𝐽 (𝑡)𝐹 [𝑢 0 ], √ 2𝜀 where 𝑆union is a softened maximum of the two lobe signed distances, 𝑆oval is an anisotropic oval distance in the lobe-axis frame, 𝑆neck fills the overlap neck, 𝑀 (𝑡) decays the exact initial-shape memory, and 𝐽 (𝑡)𝐹 [𝑢 0 ] supplies the initial PDE tangent. The term 𝑆neck is the explicit mechanism for the interface-fusion region: it fills the overlap neck instead of treating the two lobes as independent shrinking circles. The 23 trainable parameters have direct geometric meanings: lobe shrink rates, center-relaxation time and fraction, oval blend time, oval axes and shrink rates, union softness, neck-fill scales, memory decay, and tangent-jet strength. No MLP correction is needed. The best candidate reaches relative 𝐿 2 = 0.0107 with 23 trainable scalars.

8

A Preprint

Agentic Symbolic Search

y

t = 0.00

t = 4.21

t = 8.00

1

1

1

0

0

0

−1

−1 −1

0

1

−1 −1

0

x

1

−1

x

0

1

x

(a) Allen–Cahn 2D peanut-to-oval snapshots. Background color shows the reference solution (red: 𝑢 ≈ +1, blue: 𝑢 ≈ −1); the black contour at 𝑢 = 0 traces the candidate’s phase interface.

Figure 4: Allen–Cahn 2D: the 𝐿 2 -best symbolic candidate tracks the peanut-to-oval geometry.

Allen-Cahn 2D

10−2

physics loss (raw MSE)

physics loss (raw MSE)

L 2-best candidate

10−1 10−3

10−4

best-so-far validation L 2

best-so-far validation L 2

10−2

0

1

2

3

4

5

outer-loop iteration (b) Dual-axis trajectory: best-so-far validation 𝐿 2 drops early and the raw PDE residual continues to decrease through iteration 5.

Figure 4: Allen–Cahn 2D (continued): residual and best-so-far validation trajectories.

3.2

Keller–Segel Radial Blow-up

Beyond smooth solutions, we further evaluate ASYS on PDEs with finite-time blow-up and free-boundary singularities. The radial Keller–Segel system (Keller and Segel, 1970), 𝑢𝑡 =

∫ 𝑟

1 𝜕𝑟 (𝑟𝑢𝑟 + 𝑢𝑀 (𝑟)) , 𝑟

𝑀 (𝑟) =

𝑠 𝑢(𝑠) 𝑑𝑠, 0

2

models chemotactic aggregation on 𝑟 ∈ [0, 30] with 𝑢(𝑟, 0) = 9.5𝑒 −𝑟 and symmetry at the origin. The total mass 9.5𝜋 exceeds the critical threshold 8𝜋, so the density concentrates near the origin in finite time while retaining nontrivial outer mass. The optimal validation accuracy is achieved at iteration 8, where a parsimonious nine-parameter candidate yields a relative 𝐿 2 error of 0.188. For context, the specialized SS-PINN baseline (Appendix A.9) is explicitly endowed with the analytical similarity-coordinate transform (e.g., 𝜂 = 𝑟/𝜏 𝛽 ), which maps the severe finite-time blow-up in physical coordinates into a well-conditioned stationary profile-learning task for the MLP; it is trained directly on the resulting profile equation, the most favorable setting for the baseline. Leveraging this hand-crafted coordinate system, SS-PINN attains an 𝐿 2 error of 0.258. ASYS, by contrast, uncovers a lower-error representation (𝐿 2 = 0.188) without any a priori exposure to the global similarity reduction. The parameter landscape summarized in Table 2 highlights this architectural

9

A Preprint

Agentic Symbolic Search

Keller-Segel radial profile overlay 103 102

u(r, t)

101 100 10−1 10−2 10−3

candidate solid, reference dashed t/T * = 0.00

10−4 10−5

t/T * = 0.62

10−4

t/T * = 0.90

10−3

10−1

10−2

100

101

r

Figure 5: Radial profiles for the 𝐿 2 -best representation at iteration 8. Solid curves show the best representation found; dashed curves show the finite-volume reference.

contrast: while SS-PINN requires 49,921 learned weights to fit the transformed profile, ASYS condenses the underlying dynamic mechanism into nine interpretable trainable scalars within an explicit contraction law. The best candidate blends an early Taylor branch, anchored to the initial condition and its PDE tangent, with a late blow-up branch composed of a critical core and two Gaussian halos: 𝑢(𝑟, 𝑡) = (1 − 𝛼(𝑡))𝑢 T (𝑟, 𝑡) + 𝛼(𝑡)𝑢 B (𝑟, 𝑡), 𝑢B =

𝑀𝑐

+

2 𝜆(𝑡) 2 1 + (𝑟/𝜆(𝑡)) 2

𝑀in − (𝑟/𝜎ℎ ) 2 𝑀out − (𝑟/𝜎𝑜 ) 2 𝑒 + 2 𝑒 . 𝜎ℎ2 𝜎𝑜

The rational core 𝑀𝑐 /𝜆2 (1 + (𝑟/𝜆) 2 ) 2 reuses the classical critical Keller–Segel profile (Herrero and Velázquez, 1996). The contraction scale 𝜆(𝑡) = 𝑓𝑐 + (1 − 𝑓𝑐 ) ((𝑇 ∗ − 𝑡)/𝑇 ∗ ) 𝛾 parameterizes the pre-asymptotic path toward that profile with a learned power law, for which existing theory does not provide a closed form. The two Gaussian halos carry the supercritical excess mass 𝜋(9.5 − 𝑀𝑐 ) above the 8𝜋 threshold. The full representation compresses the initial 12,737-parameter neural candidate to nine interpretable scalars (core mass, halo widths, contraction exponent, blend time). Later iterations continue to reduce the physics loss but do not change the iteration 8 validation 𝐿 2 (Figure 6c). Rather than merely reflecting an optimization trace, this divergence isolates two distinct phenomena: the residual continues to reward lower equation residual, while the low-dimensional nine-parameter family reaches its useful profile-matching limit near the singular window. At the same time, the contraction-scale panels in Figure 6 show the physical mechanism extracted by the search: a finite-window power law for the inner scale 𝜆(𝑡). Thus the main scientific output is the explicit contraction mechanism, even though the validation profile error plateaus after iteration 8.

10

A Preprint

Agentic Symbolic Search

Keller-Segel contraction scale

Keller-Segel scale log-log fit 100

0.8

learned scale fit slope 0.96

learned core scale

learned core scale

1.0

0.6

0.4

0.2

variable exponent γ = 1.33

10−1 0.0

0.2

0.4

0.6

10−1

0.8

100

t/T *

τ = (T * − t)/T *

(a) Learned contraction scale over the scored window.

(b) Log–log fit exposes the power law encoded by the learned scale.

Keller-Segel

100

physics loss (ρgeo)

physics loss (ρgeo)

L 2-best candidate

100

10−1

iteration 8 best L 2

SS-PINN

10−2

0

1

2

3

4

5

6

7

8

9

best-so-far validation L 2

best-so-far validation L 2

10

outer-loop iteration (c) Dual-axis trajectory: iteration 8 gives the best validation 𝐿 2 value; the dotted line marks the specialized SS-PINN baseline.

Figure 6: Keller–Segel: the 𝐿 2 -best candidate appears at iteration 8 with relative 𝐿 2 = 0.188 and nine parameters; the trajectory makes the residual-versus-validation distinction explicit.

3.3

Graveleau PME Focusing

Keller–Segel shows progressive structural enrichment: the search builds its representation piece by piece across iterations. Graveleau tests a different mode: can ASYS instantiate a known mathematical framework as a differentiable program and fit its free parameters? The second-kind self-similar structure for the porous medium equation (𝑚 = 2) is classical (Aronson and Graveleau, 1993) and is provided as problem guidance, but the similarity exponent 𝛽 is a nonlinear eigenvalue that must be determined numerically. In the pressure variable 𝑣, the equation becomes 𝑣 𝑡 = 𝑣Δ𝑣 + |∇𝑣| 2 in radial geometry, with a centered hole in the initial pressure and zero-flux at the outer boundary. The best candidate appears at iteration 4 with relative 𝐿 2 = 0.00132 and a learned exponent 𝛽model = 0.928 (reference fit 𝛽ref = 0.877). It combines a Taylor anchor, a self-similar pressure profile, and an MLP correction:
𝑣(𝑟, 𝑡) = 𝑣 anchor + 𝑝 𝑚 𝑔(𝑡)𝜙(𝑟) 𝑣 sim − 𝑣 anchor + 𝑔(𝑡)𝜙(𝑟)𝜎𝑐 𝑁 𝜃 (𝑧), where 𝜙(𝑟) = (1 − (𝑟/6) 2 ) 2 tapers the correction at the outer boundary. The self-similar component is 𝑟 𝑣 sim = 𝐴 𝜏 2𝛽−1 𝐹 (𝜂), 𝜂= , 𝑅(𝑡) = 𝜏 𝛽 , 𝑅(𝑡) 11

A Preprint

Agentic Symbolic Search

𝑠 . 1+𝑠 The Taylor anchor matches 𝑣(𝑟, 0) and 𝜕𝑡 𝑣(𝑟, 0) exactly, ensuring the correct initial evolution. The self-similar branch encodes the second-kind scaling 𝑅(𝑡) ∼ 𝜏 𝛽 and the far-field power 𝐹 (𝜂) ∼ 𝜂2−1/𝛽 (Aronson and Graveleau, 1993); the search instantiates this structure as a differentiable program and fits 𝛽, the amplitude 𝐴, and three log-shape correction coefficients from the PDE-residual and public-constraint scores. A four-hidden-layer MLP (13121 parameters) corrects the transition between the anchor and self-similar regimes. 𝐹 (𝜂) = 𝛽𝑠(1 + 𝑠) 1−1/𝛽 exp(𝑐 1 𝑦 + 𝑐 2 𝑦 2 + 𝑐 3 𝑦 3 ),

𝑠 = (𝜂 − 1)+ ,

𝑦=

Figure 7 should be read as a scaling-law check, not only as a profile overlay. The log–log radius analysis is the standard way to expose a second-kind exponent, and here it gives 𝛽model = 0.928 against the reference fit 𝛽ref = 0.877. The 6% gap reflects the indirect determination: the evaluator scores PDE satisfaction and public constraints, not the free-boundary radius. Despite this gap, the low validation 𝐿 2 = 0.00132 confirms that the overall pressure profile is well tracked. The important point is that a residual-guided symbolic program recovered the right self-similar framework closely enough for the nonlinear similarity exponent to be inferred after the run. Figure 8 separates the agent-visible raw residual from the hidden offline validation metric. The initial high-capacity seed has the smallest raw residual on the public collocation score, but it lacks an explicit second-kind free-boundary law and therefore remains less accurate than later structured candidates in validation 𝐿 2 . By iteration 4, the search has injected the rigid self-similar framework above. This restriction reduces local fitting freedom and leaves a larger raw residual than the seed, but it aligns the macroscopic focusing geometry and gives the best validation error. Thus the figure is not a train–validation discrepancy; it shows why residual-only local scoring must be interpreted together with structural offline validation in free-boundary focusing. Graveleau pressure profile overlay 1.4

pressure v(r, t)

1.2 1.0 0.8 0.6 0.4 0.2

candidate solid, reference dashed t = 0.124

0.0 0.0

0.5

1.0

1.5

t = 0.589

t = 0.617

2.0

2.5

r

(a) Pressure profiles for the focusing free-boundary shape. Solid curves show the 𝐿 2 -best candidate; dashed curves show the finite-volume reference. Graveleau front radius

Graveleau beta analysis 100

reference β = 0.877

0.8

candidate β = 0.928

front radius R(t)

candidate β = 0.928

0.6

R(t)

reference β = 0.877

0.4

0.2

10−1

10−2 0.0 0.1

0.2

0.3

0.4

0.5

10−2

0.6

10−1

100

τ = (T − t)/T

t

(b) Free-boundary radius with 𝛽model = 0.928 and reference fit 𝛽ref = 0.877.

(c) Log–log analysis of the second-kind exponent scale.

Figure 7: Graveleau focusing: the 𝐿 2 -best candidate instantiates a second-kind self-similar free-boundary framework with 𝛽model = 0.928, about 6% above the reference fit 𝛽ref = 0.877.

12

A Preprint

Agentic Symbolic Search

Graveleau physics loss (raw residual)

physics loss (raw residual)

L 2-best candidate

4 × 10−3 3 × 10−3

2 × 10−3

10−3

best-so-far validation L 2

6 × 10−3

best-so-far validation L 2

10−3

0

1

2

3

4

5

6

7

8

9

10

outer-loop iteration Figure 8: Graveleau residual and offline validation trajectory. The lowest raw residual occurs at the initial candidate, but the best profile match appears later at iteration 4, after the second-kind self-similar free-boundary structure is introduced.

3.4

Stress Test: gCLM

To probe the expressivity limits of ASYS, we consider a final stress test where the underlying equation features strong nonlocal coupling. The generalized Constantin–Lax–Majda equation (Chen, 2020) with 𝑎 = 0.25, 𝜔𝑡 + 𝑎 𝑢 𝜔 𝑥 = 𝑢 𝑥 𝜔,

𝑢 𝑥 = 𝐻𝜔,

𝑥 ∈ [−8, 8],

has periodic boundary conditions and 𝜔(𝑥, 0) = 10 sin(𝜋𝑥/8). The Hilbert transform 𝐻 couples every spatial mode nonlocally, and the solution develops a cusp whose width contracts sharply near the focusing time. The search represents the nonlocal coupling explicitly inside its candidates and tries pole-type and Poisson-kernel forms across iterations, but no candidate captures the late-time concentration. This failure is structurally different from the bounded cases: the nonlocal Hilbert transform turns the search for a compact local ansatz into a global-mode problem, causing a dimensional explosion in the symbolic search space. The best candidate is a Taylor-anchored neural ansatz with a fixed blow-up clock. It tracks the early-time evolution but misses the cusp (Figure 9a). The physics loss decreases across iterations, yet the best validation 𝐿 2 remains at 0.465 (Figure 9b). The SS-PINN baseline (Appendix A.9) reaches 𝐿 2 = 0.196 once the known self-similar coordinate system is supplied. This diagnosis isolates the main missing ingredient in the ASYS run: the search did not construct the global self-similar reduction within ten iterations. Unlike Keller–Segel, where the critical profile provides a strong scaffold, this run lacked a comparable structural starting point. The parameter metrics in Table 2 should be interpreted in this precise context: although the SS-PINN profile MLP uses 49,921 trainable weights, its decisive conceptual advantage stems from the pre-supplied self-similar coordinate framework, an inductive bias that the current ASYS search space failed to construct autonomously within ten iterations.

13

A Preprint

Agentic Symbolic Search

ω(x, t)

t/T * = 0.00

t/T * = 0.69

10.0

20

7.5

15

5.0

10

2.5

5

0.0

0

−2.5

−5

−5.0

−10

−7.5

−15

−10.0

−20 −5

0

t/T * = 0.99

400

reference candidate

200 0 −200 −400

5

−5

0

x

5

−5

x

0

5

x

(a) gCLM 𝑎 = 0.25 profile for the 𝐿 2 -best candidate. The late-time concentration remains visibly mismatched.

Figure 9: gCLM stress test: the 𝐿 2 -best candidate misses the late-time concentration.

gCLM

100

physics loss (ρgeo)

physics loss (ρgeo)

L 2-best candidate

6 × 10−1

10−3

4 × 10−1 3 × 10−1 SS-PINN

2 × 10−1

best-so-far validation L 2

best-so-far validation L 2

10−4 0

1

2

3

4

5

6

7

8

9

10

outer-loop iteration (b) Dual-axis trajectory: the physics loss decreases, but best-so-far validation 𝐿 2 remains high; the dotted line marks the specialized SS-PINN baseline.

Figure 9: gCLM stress test (continued): the physics loss decreases while the offline validation error stays high.

4

Conclusion

Interpretable mathematical structure, historically the product of case-by-case mathematical analysis, can be obtained through automated search. The five case studies confirm this across a range of structural difficulty: ASYS recovers a known analytical form when one exists (NLS), constructs new mathematical representations where no closed form was previously available (Allen–Cahn, Keller–Segel, Graveleau), and exposes its current boundary on problems whose global coordinate structure remains out of reach (gCLM). Three limitations shape the current results. First, ASYS is a structural characterizer for difficult individual evolutions, not an amortized solver for families of initial conditions. When the goal is rapid inference across many inputs, neural operators and related surrogate models remain the appropriate tools. The value of ASYS is different: it spends search computation to expose an interpretable mechanism. Second, the quality of the search depends on the coding agent’s ability to propose useful ansatz families within a limited number of iterations, as well as on how the evolutionary framework organizes exploration across the population. The gCLM stress test illustrates this directly: a self-similar PINN performs well once the coordinate transform is supplied, while ASYS does not discover that global coordinate system within ten iterations. The Keller–Segel residual-versus-validation gap points to the same bottleneck. In both

14

A Preprint

Agentic Symbolic Search

cases the limitation lies in the agent’s structural reach, not in the scoring mechanism. High-dimensional, strongly coupled, or globally nonlocal PDEs are likely to require hybrid representations, for example a symbolic singular core coupled to a neural operator or spectral surrogate for the smoother background field. As the capabilities of the underlying language models advance, the structural search is likely to benefit directly. Third, the inner-loop optimizer (L-BFGS) operates on a non-convex landscape whose conditioning deteriorates near singularities; a more flexible optimization strategy, such as continuation, relaxed gradient flows, or Hessian-free updates, could stabilize parameter fitting in blow-up and free-boundary regimes. Several extensions follow naturally from the current framework. The scoring system introduced in Section 2.3 evaluates representations using only equation residuals and public constraints, with no observed solution data. Incorporating observed data as a third scoring signal would extend the method to inverse problems, where the governing equation is partially known and unknown terms or parameters must be recovered from measurements. Because the evaluator structure remains unchanged, this extension requires only an additional score dimension, not a redesign of the search. A second extension is methodological: ablate the specificity of the problem guidance, separating what comes from public mathematical prior knowledge from what is found by outer-loop exploration under weaker prompts. Beyond forward and inverse problems on canonical equations, applying the method to real-world systems introduces further challenges: noisy and partial observations, uncertain or approximate governing equations, and phenomena such as turbulence whose characteristic structures remain open mathematical questions. Finally, the representations produced here are starting points for rigorous analysis: an ansatz that captures the correct scale, exponent, or profile can inform conjectures, guide function-space design, or seed computer-assisted proofs (Wang et al., 2023; Chen and Hou, 2025). Connecting agent-generated representations to formal verification is a natural next step.

Acknowledgements Liu Yang acknowledges support from the National Research Foundation, Singapore, under the NRF fellowship (Project No. NRF-NRFF17-2025-0006).

References Samuel M. Allen and John W. Cahn. A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening. Acta Metallurgica, 27(6):1085–1095, 1979. doi: 10.1016/0001-6160(79)90196-2. D. G. Aronson and J. Graveleau. A self-similar solution to the focusing problem for the porous medium equation. European Journal of Applied Mathematics, 4(1):65–81, 1993. doi: 10.1017/S095679250000098X. Grigory Isaakovich Barenblatt. Scaling, Self-Similarity, and Intermediate Asymptotics. Cambridge Texts in Applied Mathematics. Cambridge University Press, 1996. doi: 10.1017/CBO9781107050242. Andrew R. Barron. Approximation and estimation bounds for artificial neural networks. Machine Learning, 14(1): 115–133, 1994. doi: 10.1007/BF00993164. Jule G. Charney, Ragnar Fjørtoft, and John von Neumann. Numerical integration of the barotropic vorticity equation. Tellus, 2(4):237–254, 1950. doi: 10.3402/tellusa.v2i4.8607. Jiajie Chen. Singularity formation and global well-posedness for the generalized Constantin–Lax–Majda equation with dissipation. Nonlinearity, 33(5):2502–2532, 2020. doi: 10.1088/1361-6544/ab74b0. arXiv:1908.09385. Jiajie Chen and Thomas Y. Hou. Singularity formation in 3D Euler equations with smooth initial data and boundary. Proceedings of the National Academy of Sciences, 122(27):e2500940122, 2025. doi: 10.1073/pnas.2500940122. Matthew W Choptuik. Universality and scaling in gravitational collapse of a massless scalar field. Physical Review Letters, 70(1):9–12, 1993. doi: 10.1103/PhysRevLett.70.9. Richard Courant, Kurt Friedrichs, and Hans Lewy. über die partiellen differenzengleichungen der mathematischen physik. Mathematische Annalen, 100(1):32–74, 1928. doi: 10.1007/BF01448839. George Cybenko. Approximation by superpositions of a sigmoidal function. Mathematics of Control, Signals and Systems, 2(4):303–314, 1989. doi: 10.1007/BF02551274. 15

A Preprint

Agentic Symbolic Search

Ronald A. DeVore. Nonlinear approximation. Acta Numerica, 7:51–150, 1998. doi: 10.1017/S0962492900002816. Joseph Fourier. Théorie analytique de la chaleur. Firmin Didot, Paris, 1822. Miguel A. Herrero and Juan J. L. Velázquez. Singularity patterns in a chemotaxis model. Mathematische Annalen, 306 (1):583–623, 1996. Kurt Hornik. Approximation capabilities of multilayer feedforward networks. Neural Networks, 4(2):251–257, 1991. doi: 10.1016/0893-6080(91)90009-T. Evelyn F. Keller and Lee A. Segel. Initiation of slime mold aggregation viewed as an instability. Journal of Theoretical Biology, 26(3):399–415, 1970. doi: 10.1016/0022-5193(70)90092-5. Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations (ICLR), 2015. Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, 2021. URL https://openreview.net/forum?id=c8P9NQVtmnO. Ziyang Liu, Xinyan Guo, Xuchen Wei, Han Hao, and Yang Liu. Escher-loop: Mutual evolution by closed-loop self-referential optimization. arXiv preprint arXiv:2604.23472, 2026. Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3:218–229, 2021. doi: 10.1038/s42256-021-00302-5. Zhou Lu, Hongming Pu, Feicheng Wang, Zhiqiang Hu, and Liwei Wang. The expressive power of neural networks: A view from the width. In Advances in Neural Information Processing Systems (NeurIPS), pages 6232–6240, 2017. arXiv:1709.02540. Jorge Nocedal. Updating quasi-Newton matrices with limited storage. Mathematics of Computation, 35(151):773–782, 1980. doi: 10.1090/S0025-5718-1980-0572855-7. Alexander Novikov, Ngân Vũ, Marvin Eisenberger, Emilien Dupont, Po-Sen Huang, Adam Zsolt Wagner, Sergey Shirobokov, Borislav Kozlovskii, Francisco J. R. Ruiz, Abbas Mehrabian, M. Pawan Kumar, Abigail See, Swarat Chaudhuri, George Holland, Alex Davies, Sebastian Nowozin, Pushmeet Kohli, and Matej Balog. AlphaEvolve: A coding agent for scientific and algorithmic discovery. arXiv:2506.13131, 2025. Allan Pinkus. Approximation theory of the MLP model in neural networks. Acta Numerica, 8:143–195, 1999. doi: 10.1017/S0962492900002919. Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686–707, 2019. doi: 10.1016/j.jcp.2018.10.045. Bernardino Romera-Paredes, Mohammadamin Barekatain, Alexander Novikov, Matej Balog, M Pawan Kumar, Emilien Dupont, Francisco JR Ruiz, Jordan S Ellenberg, Pengming Wang, Omar Fawzi, et al. Mathematical discoveries from program search with large language models. Nature, 625(7995):468–475, 2024. doi: 10.1038/s41586-023-06924-6. Junkichi Satsuma and Nobuo Yajima. Initial value problems of one-dimensional self-modulation of nonlinear waves in dispersive media. Progress of Theoretical Physics Supplement, 55:284–306, 1974. Parshin Shojaee, Kazem Meidani, Shashank Gupta, Amir Barati Farimani, and Chandan K. Reddy. LLM-SR: Scientific equation discovery via programming with large language models. In International Conference on Learning Representations, pages 1–22, 2025. doi: 10.48550/arXiv.2404.18400. Oral presentation; arXiv:2404.18400. Boai Sun, Wenjin Guo, Zongmin Yu, and Liu Yang. Self-evolving scientific agent discovers generalizable physicallyreasoned fluid control. arXiv preprint arXiv:2606.08405, 2026.

16

A Preprint

Agentic Symbolic Search

Matus Telgarsky. Benefits of depth in neural networks. In Conference on Learning Theory (COLT), pages 1517–1539, 2016. arXiv:1602.04485. Juan Diego Toscano, Daniel T. Chen, and George Em Karniadakis. ATHENA: Agentic team for hierarchical evolutionary numerical algorithms. arXiv preprint arXiv:2512.03476, 2025. doi: 10.48550/arXiv.2512.03476. Juan Diego Toscano, Zhaojie Chai, and George Em Karniadakis. GRAFT-ATHENA: Self-improving agentic teams for autonomous discovery and evolutionary numerical algorithms. arXiv preprint arXiv:2605.11117, 2026. doi: 10.48550/arXiv.2605.11117. Yongji Wang, Ching-Yao Lai, Javier Gómez-Serrano, and Tristan Buckmaster. Asymptotic self-similar blow-up profile for three-dimensional axisymmetric Euler equations using neural networks. Physical Review Letters, 130(24):244002, 2023. doi: 10.1103/PhysRevLett.130.244002. Yongji Wang, Mehdi Bennani, James Martens, Sébastien Racanière, Sam Blackwell, Alex Matthews, Stanislav Nikolov, Gonzalo Cao-Labora, Daniel S. Park, Martin Arjovsky, Daniel Worrall, Chongli Qin, Ferran Alet, Borislav Kozlovskii, Nenad Tomašev, Alex Davies, Pushmeet Kohli, Tristan Buckmaster, Bogdan Georgiev, Javier Gómez-Serrano, Ray Jiang, and Ching-Yao Lai. Discovery of unstable singularities. arXiv preprint arXiv:2509.14185, 2025. doi: 10.48550/arXiv.2509.14185. Shijie Xia, Yuhan Sun, and Pengfei Liu. SR-Scientist: Scientific equation discovery with agentic AI. arXiv preprint arXiv:2510.11661, 2025. doi: 10.48550/arXiv.2510.11661. Zongmin Yu and Yang Liu. Evolutionary ensemble of agents. arXiv preprint arXiv:2605.09018, 2026.

Appendix A.1

Outer-Loop Framework

ASYS uses EvE (Yu and Liu, 2026; Liu et al., 2026) as its outer-loop search framework. This section summarizes the aspects of EvE that are relevant to the reported experiments; full details are given in the EvE reference. EvE maintains two co-evolving populations (Figure 2). The candidate population stores every candidate program produced during the run, together with its four-dimensional score. No candidate is discarded; the population grows monotonically. The strategy population stores search strategies, each a guidance document describing how to propose new candidates (what mathematical directions to explore, what pitfalls to avoid). Strategies are scored by the quality of the candidates they produce, so effective strategies are sampled more often in later iterations. At each iteration, the framework runs the following steps: 1. Sample context. For each of the two parallel workers, sample up to four example candidates from the candidate population using rank-exponential-sum weights (Appendix A.5), and one strategy from the strategy population. 2. Propose. Each worker, a coding-agent session, reads the problem guidance, the sampled examples, and the sampled strategy, then writes a new candidate program (the representation, the training loss, and the optimizer settings). 3. Evaluate. The evaluator fits the candidate’s trainable parameters by L-BFGS under a fixed wall-clock budget, then scores the result on the four dimensions (Appendix A.4). 4. Update populations. The scored candidate is added to the candidate population. A new or revised strategy, reflecting what the worker learned during this iteration, is added to the strategy population. The iteration count, worker count, example count, wall-clock budget, and other run-level parameters are listed in Appendix A.2.

17

A Preprint

Agentic Symbolic Search

A.2

Run Configuration

Table 3 lists the outer-loop parameters shared across all five reported cases. The terms used are defined below. An iteration is one cycle of the outer loop. At each iteration, two workers run in parallel; each worker is an independent coding-agent session that proposes one new candidate program. Before proposing, each worker receives up to four examples: previous candidates sampled from the population by the rank-exponential-sum rule (Appendix A.5). After the worker submits a candidate, the evaluator fits its trainable parameters by L-BFGS (Nocedal, 1980) under a fixed wall-clock budget and scores it on the four dimensions (Appendix A.4). The scored candidate is then added to the population and becomes available as a future example. In parallel with the candidate population, a second population of search strategies is maintained. Each strategy is a guidance document that a worker follows when proposing candidates. Strategies are scored by the quality of the candidates they produce, so strategies that lead to lower-loss candidates are sampled more often in later iterations. In the reported experiments, two strategies are produced per iteration, matching the two workers. Table 3: Outer-loop configuration for all reported experiments. Parameter

Value

Workers per iteration Examples shown per worker Outer-loop iterations (bounded cases) Outer-loop iterations (blow-up / stress test) Inner-loop optimizer Wall-clock budget per candidate Score dimensions Selection weights (𝑤𝑑 , 𝑇𝑑 ) Population Hardware

2 4 5 10 L-BFGS 30 s (bounded) / 300 s (blow-up and free-boundary), with early stopping 4 (physics, ic, bc, comp) 1.0, 1.0 for all dimensions append-only (no discard) Apple M4 laptop CPU, no GPU

For cost accounting, the coding-agent sessions in these experiments were run through a subscription-covered interactive agent interface, not through a per-token API workflow. Thus the reported runs incurred no additional agent-side API billing beyond the subscription already used to access the interactive interface. Any API list-price or equivalent-token accounting used for cross-report comparison should be read only as a normalization convention, not as the invoice paid for these runs. Each case begins from guidance containing the PDE and public constraints, the editable representation files, the evaluation dimensions, and physical background or literature hints. The guidance permits rewriting the differentiable representation and the scalar training objective, while excluding external solvers, data loading, network access, and dense lookup tables. A.3

Agent Workspace Examples

Each search worker receives a structured workspace containing three components: a task description (README.md) that states the PDE, score dimensions, admissibility rules, and editable files; a strategy record (guidance/strategy.md) with search suggestions; and a reference record (guidance/references.md) with public mathematical facts. The worker may edit four files that define the candidate representation, the training loss, the optimizer configuration, and fitting hyperparameters. Numerical reference solutions used for offline validation are not present in the workspace. The two examples below reproduce sanitized workspace packets. They intentionally show the public mathematical priors exposed to the search, such as transition layers, curvature motion, critical mass, and self-similar clocks. These are the inductive biases ASYS is designed to operationalize, not hidden solution tables: absolute paths, run-management material, previous-candidate diagnostics, and all numerical validation references are omitted. The task statement, scoring dimensions, admissibility rules, editable surface, strategy record, and reference record are preserved.

18

A Preprint

Agentic Symbolic Search

Bounded example: Allen–Cahn 2D peanut merger. README.md Find a compact, trainable scalar field 𝑢(𝑥, 𝑦, 𝑡) for 𝑢 𝑡 = 𝜀 2 (𝑢 𝑥 𝑥 + 𝑢 𝑦 𝑦 ) + 𝑢 − 𝑢 3 ,

𝜀 = 0.04.

The domain is (𝑥, 𝑦) ∈ [−1, 1) 2 , periodic in both directions, with 0 ≤ 𝑡 ≤ 8. The model input is physical (𝑥, 𝑦, 𝑡) and the output is one column 𝑢(𝑥, 𝑦, 𝑡). The initial geometry is a connected, asymmetric two-source peanut. It is the smooth union of two diffuse disks: center (−0.15, 0.02) and radius 0.22 for the first blob; center (0.18, −0.03) and radius 0.18 for the second. The center distance is about 0.337, the radius sum about 0.400, and the overlap width about 0.063, roughly 1.5𝜀. The initial positive phase is already connected. This is not two independent one-dimensional fronts. The useful target is an explicit time-dependent signed-distance geometry for a peanut that relaxes toward a smoother oval or rounder interface while preserving the periodic square geometry. Allen–Cahn diffuse interfaces have a natural tanh profile across an 𝑂 (𝜀) transition layer. In the sharp-interface limit, curved interfaces move approximately by mean curvature. The overlap neck is therefore a genuine two-dimensional interaction zone: a candidate that treats the two disks as uncoupled shrinking circles is missing the main geometry. Scoring. Training inside the submitted code fits the trainable parameters of a fixed structure by minimizing a scalar loss chosen by the worker. The search score, however, evaluates the structure on four independent dimensions: physics, the strong-form Allen–Cahn residual; ic, the match to the initial peanut field at 𝑡 = 0; ic_compatibility, agreement of the initial time derivative with 𝜀 2 Δ𝑢 0 + 𝑢 0 − 𝑢 30 ; and bc, periodic value and first-derivative matching on opposite sides of the square. These dimensions are ranked independently and are not reduced to one weighted scalar. No external reference table is supplied to the evaluator. Symbolic representation constraint. This bounded run is restricted to explicit differentiable symbolic programs. The submitted representation may use trainable scalar parameters, analytic functions, coordinate transforms, piecewise branches, smooth gates, and closed-form distance constructions, but it may not introduce neural network layers such as fully connected, convolutional, recurrent, or transformer modules. The purpose of the run is to test whether the peanut-to-oval dynamics can be captured by explicit geometry, not by a black-box correction term. Admissibility rules. The workspace rules exclude shortcuts that would replace representation search with a different task: abandoning the initial condition after 𝑡 = 0, embedding a numerical time-marching solver, using a lookup table keyed to grid points or external validation data, or returning a field independent of trainable parameters. Closed-form transforms, analytic series from a reduced equation, and coefficients with an analytic source remain admissible. Editable surface. The worker may edit only problems/pinn/allen_cahn_2d/model.py, loss.py, optim.py, and config.yaml.

guidance/strategy.md The strategy record says the main difficulty is the neck: the initial condition is a connected asymmetric peanut built from overlapping diffuse disks, so simply shrinking two disks misses the curvature-driven relaxation into a smoother oval. It also warns that matching the value at 𝑡 = 0 is not enough; the initial time derivative must match 𝜀 2 Δ𝑢 0 + 𝑢 0 − 𝑢 30 . Periodicity matters as well: opposite sides of the square must agree in both value and derivative. The suggested directions are: build a low-parameter signed-distance scaffold for the connected peanut; blend from the peanut distance field to an oval or rounded connected interface with controls for neck smoothing, axis lengths, and center drift; add a curvature-aware early-time correction for the PDE tangent; use periodic coordinate features or boundary factors; and keep any secondary neural or Fourier correction compact so that it does not dominate the explicit geometry.

guidance/references.md The reference record keeps the structural facts short: tanh transition layers are natural for Allen–Cahn interfaces, and mean-curvature motion is the relevant sharp-interface background. It supplies public mathematical context, not a trajectory table, grid table, or numerical validation artifact.

19

A Preprint

Agentic Symbolic Search

Blow-up example: Keller–Segel 2D radial concentration. README.md Find a compact, trainable radial density 𝑢(𝑟, 𝑡) for 𝑢𝑡 =

∫ 𝑟

1 𝜕𝑟 𝑟𝑢𝑟 + 𝑢𝑀 (𝑟) , 𝑟

𝑀 (𝑟) =

𝑠 𝑢(𝑠, 𝑡) 𝑑𝑠. 0 2

The domain is 0 ≤ 𝑟 ≤ 30, with 0 ≤ 𝑡 ≤ 0.9 · 0.9628. The initial data are 𝑢(𝑟, 0) = 9.5𝑒 −𝑟 , with total mass 9.5𝜋, above the 8𝜋 critical mass. The boundary conditions are origin symmetry and outer decay at 𝑟 = 30. The model input is physical (𝑟, 𝑡) and the output is one density column. This is a supercritical radial Keller–Segel concentration benchmark. The chemotactic drift is nonlocal: 𝑀 (𝑟) is a cumulative radial mass at each time slice, and the residual is computed by a finite-volume radial operator on a geometric grid. A generic pointwise neural representation is not a strong hypothesis; useful candidates should make positivity, mass concentration, and a shrinking inner scale visible while preserving a halo that carries the excess mass. Scoring. Training fits parameters of a fixed structure by minimizing a worker-defined scalar loss. The search score separately 2 ranks four dimensions: physics, the finite-volume radial Keller–Segel residual; ic, the match to 9.5𝑒 −𝑟 at 𝑡 = 0; ic_compatibility, agreement of the one-sided initial time derivative with the finite-volume right-hand side at the initial condition; and bc, origin symmetry plus outer value/derivative decay. Mass conservation and positivity diagnostics are also reported as guards against degenerate density shapes. Admissibility rules. The workspace rules exclude shortcuts that would break the interpretation of the candidate as a continuous representation: matching 𝑢(𝑟, 0) and then switching to an unrelated later field; embedding a finite-volume, finite-difference, spectral, Runge–Kutta, ETD, split-step, or other time-marching solver; using a dense lookup table keyed to grid points or external validation data; returning a parameter-independent field; or tuning a correction only to evaluator probe points rather than defining one continuous physical field. Analytic self-similar forms, asymptotic profiles, low-rank corrections, and spatial projections with analytic source remain admissible. Editable files and residual helpers. The worker may edit only problems/blowup/ks_radial_2d/model.py, loss.py, optim.py, and config.yaml. The provided residual helpers expose the Keller–Segel residual, boundary loss, positivity penalty, initial condition, radial right-hand side, cumulative mass, total mass, and constants including the critical mass, scored time horizon, outer radius, and radial finite-volume grid measures.

guidance/strategy.md The strategy record emphasizes that the equation is nonlocal: a pointwise ansatz can match the initial condition but miss the cumulative-mass balance that drives the finite-volume residual. Since 9.5𝜋 > 8𝜋, the profile must represent finite-time concentration: the central density grows and the inner length scale contracts far below the initial Gaussian width. A broad halo or tail is still needed to carry the excess mass and protect outer decay; the core and halo play different roles. The suggested directions are: start from the exact Gaussian initial condition and preserve it by construction; add a smooth time factor whose first derivative approximates the finite-volume right-hand side at 𝑡 = 0; prefer positive bubble-plus-halo profiles with trainable scale, amplitude, and tail parameters; explore self-similar clocks with a shrinking inner length scale and a broad mass reservoir; use the provided mass diagnostics in the training loss; and allow low-rank radial corrections only when they are anchored to analytic functions such as Gaussians, rational tails, or smooth bump factors. Dense per-radius coefficients and interpolation surfaces are ruled out.

guidance/references.md The reference record keeps the problem facts explicit: the PDE and domain above, the 8𝜋 critical mass threshold, the independent physics/ic/ic_compatibility/bc score dimensions, and the fact that numerical reference solutions are reserved for offline validation, not design data.

A.4

Score Dimensions

Section 2.3 introduces the four score components used to compare candidates. This appendix gives the case-specific definitions, including the physics-residual normalization, collocation counts, and finite-difference steps. The four evaluation dimensions are
𝑠(𝑢 𝜃 ) = −ℓphys , −ℓic , −ℓbc , −ℓcomp , 20

A Preprint

Agentic Symbolic Search

each defined below. The physics-loss trajectories shown in the Results figures correspond to ℓphys ; equivalently, the score components are (𝑠phys , 𝑠ic , 𝑠bc , 𝑠comp ). Physics residual For bounded cases (NLS and Allen–Cahn 2D), the physics loss is the raw mean-squared PDE residual: 𝑁 1 ∑︁ 2 ℓphys = 𝑅 . 𝑁 𝑗=1 𝑗 For blow-up and singular-structure cases, the physics loss is a scale-free geometric-mean relative residual. For pointwise residuals, 𝑁 |𝑅 𝑗 | © 1 ∑︁ ª 𝜌𝑗 = , ℓphys = 𝜌geo = exp ­ log max(𝜌 𝑗 , 𝜀log ) ® , 𝐷 𝑗 + 𝜀 denom 𝑁 𝑗=1 « ¬ with 𝜀 denom = 10−10 and 𝜀log = 10−30 . Keller–Segel uses this idea per time slice, with radial finite-volume weights: ∥𝑅 𝑘 ∥ 𝑘 , ∥𝑢 𝑡 ,𝑘 ∥ 𝑘 + ∥ diff 𝑘 ∥ 𝑘 + ∥ chem 𝑘 ∥ 𝑘 ! 1 ∑︁ log max(𝜌 𝑘 , 10−12 ) , ℓphys = exp 𝑁𝑡 𝑘 𝜌𝑘 =

Í where ∥ 𝑓 ∥ 2𝑘 = 𝑗 Δ𝑚 𝑗 𝑓 𝑘,2 𝑗 and Δ𝑚 𝑗 = 12 (𝑟 2𝑗+1 − 𝑟 2𝑗 ). The reported Graveleau free-boundary run instead uses the raw residual mean-square loss, 𝑁
2 1 ∑︁ Grav 𝑣 𝑡 − 𝑣Δ𝑟 𝑣 − 𝑣 𝑟2 𝑗 . ℓphys = 𝑁 𝑗=1 Table 4: Physics residual definitions and denominators used for reported scores or diagnostics. Case

Physics score definition

NLS

Raw MSE over the real and imaginary residual components of the split NLS system. The relative denominator |𝑢𝑡 | + 21 |𝑣𝑥 𝑥 | + | | 𝜓 | 2 𝑣 | and |𝑣𝑡 | + 12 |𝑢 𝑥 𝑥 | + | | 𝜓 | 2 𝑢| is retained only as a diagnostic. Raw MSE of 𝑢𝑡 − 𝜀 2 Δ𝑢 − 𝑢 + 𝑢3 . The diagnostic denominator is |𝑢𝑡 | + | 𝜀 2 Δ𝑢| + |𝑢| + |𝑢3 | . Per time slice, ∥𝑢𝑡 ∥ + ∥ diff ∥ + ∥ chem ∥ with the finite-volume measure above. Diagnostic denominator |𝑣𝑡 | + |𝑣Δ𝑟 𝑣 | + |𝑣𝑟2 | ; the reported run uses the raw residual mean-square score above. | 𝜔𝑡 | + | 𝑎𝑢𝜔 𝑥 | + |𝑢 𝑥 𝜔 | , with 𝑢 𝑥 = 𝐻 [ 𝜔] and zero-mean velocity.

Allen–Cahn 2D Keller–Segel Graveleau gCLM

Initial value ℓic =

2 1 ∑︁ 𝑢 𝜃 (𝑥𝑖 , 0) − 𝑢 0 (𝑥 𝑖 ) , 𝑁ic 𝑖

with componentwise real and imaginary parts for complex NLS. The initial-value term uses 𝑁ic = 50 for NLS, 1024 for Allen–Cahn, and 256 for Keller–Segel, Graveleau, and gCLM. Compatibility condition The compatibility dimension compares a one-sided finite-difference time derivative against the PDE right-hand side applied to the analytic initial condition:  2 1 ∑︁ 𝑢 𝜃 (𝑥 𝑖 , Δ𝑡) − 𝑢 𝜃 (𝑥𝑖 , 0) ℓcomp = − 𝐹 [𝑢 0 ] (𝑥𝑖 ) . 𝑁ic 𝑖 Δ𝑡 The finite-difference steps are 1.571 × 10−3 for NLS, 8.0 × 10−3 for Allen–Cahn, 8.6652 × 10−4 for Keller–Segel, 6.1699 × 10−4 for Graveleau, and 2.38491 × 10−4 for gCLM. The right-hand sides are respectively the NLS split system applied to 2 sech 𝑥, the Allen–Cahn operator applied to the two-disk phase field, the finite-volume Keller–Segel right-hand side on the Gaussian initial condition, 𝑣 0 Δ𝑟 𝑣 0 + (𝑣 0′ ) 2 on the Graveleau shell, and 𝑢 𝑥 [𝜔0 ]𝜔0 − 𝑎𝑢[𝜔0 ]𝜔0,𝑥 for gCLM. 21

A Preprint

Agentic Symbolic Search

Boundary consistency The boundary loss is the squared residual of the public boundary condition: periodic value and derivative matching for NLS, periodic value and derivative matching in both coordinate directions for Allen–Cahn, radial symmetry at the origin plus outer decay for Keller–Segel, 𝜕𝑟 𝑣(6, 𝑡) = 0 for Graveleau, and periodic value matching for gCLM. A.5

Candidate Selection

When the outer loop selects which previous candidates to show the agent as examples, it ranks the full population independently in each of the four score dimensions. Let 𝑟 𝑑 (𝑐) denote the rank of candidate 𝑐 in dimension 𝑑 (rank 0 is best). Each candidate receives a sampling weight ∑︁
𝑤(𝑐) = 𝑤 𝑑 exp −𝑟 𝑑 (𝑐)/𝑇𝑑 , 𝑑 ∈ {phys, ic, bc, comp}

with 𝑤 𝑑 = 1 and 𝑇𝑑 = 1 for all four dimensions in the reported experiments, so no dimension is preferred. Candidates are then drawn without replacement using these weights, so candidates that rank well in any combination of dimensions are more likely to appear as context. The population is append-only: no candidate is discarded after evaluation. A.6

Problem Setups

Table 5 collects the PDE, domain, and public constraints for each of the five reported cases. These are the inputs given to the search; no additional problem data is supplied. Table 5: PDEs, domains, and public constraints used by the five reported cases. Case

PDE

Domain

Initial and boundary data

NLS

𝑖 𝜓𝑡 + 12 𝜓𝑥 𝑥 + | 𝜓 | 2 𝜓 = 0

𝑥 ∈ [ −5, 5] , 𝑡 ∈ [0, 𝜋/2]

Allen–Cahn 2D

𝑢𝑡 = 𝜀 2 Δ𝑢 + 𝑢 − 𝑢3 , 𝜀 = 0.04

( 𝑥, 𝑦) ∈ [ −1, 1) 2 , 𝑡 ∈ [0, 8]

𝜓 ( 𝑥, 0) = 2 sech 𝑥 ; periodic value and derivative matching at 𝑥 = ±5. Smooth union of two tanh disks centered at (−0.15, 0.02) and (0.18, −0.03) with radii 0.22 and 0.18; periodic value and derivative matching.

Keller–Segel

∫𝑟 𝑢𝑡 = 𝑟 −1 𝜕𝑟 (𝑟𝑢𝑟 + 𝑢𝑀 (𝑟 ) ) , 𝑀 (𝑟 ) = 0 𝑠𝑢(𝑠) 𝑑𝑠

Graveleau

𝑣𝑡 = 𝑣Δ𝑟 𝑣 + 𝑣𝑟2 , Δ𝑟 𝑣 = 𝑣𝑟𝑟 + 𝑟 −1 𝑣𝑟

gCLM

𝜔𝑡 + 𝑎𝑢𝜔 𝑥 = 𝑢 𝑥 𝜔 , 𝑢 𝑥 = 𝐻 [ 𝜔 ] , 𝑎 = 0.25

𝑟 ∈ [0, 30] , 𝑡 ∈ [0, 0.9𝑇 ∗ ] , 𝑇 ∗ = 0.9628 𝑟 ∈ [0, 6] , 𝑡 ∈ [0, 0.995𝑇 𝑓 ] , 𝑇 𝑓 = 0.6200901157 𝑥 ∈ [ −8, 8] , 𝑡 ∈ [0, 0.99𝑇 ∗ ] , 𝑇 ∗ = 0.2409

𝑢(𝑟 , 0) = 9.5𝑒 −𝑟 ; radial symmetry at 𝑟 = 0 and outer decay at 𝑟 = 30. 2 𝑣 (𝑟 , 0) = 1.6(1 − 𝑒 − ( (𝑟 −1) /0.6) )1𝑟 ≥1 ; outer Neumann condition at 𝑟 = 6. 𝜔 ( 𝑥, 0) = 10 sin( 𝜋 𝑥/8) ; periodic boundary conditions.

A.7

2

Candidate Formulas

The Results section presents the best candidate for each case in compact form. This appendix gives the complete mathematical expressions, including auxiliary definitions, MLP architectures, and the breakdown between analytic and fitted components. To keep the formulas readable, each expression is decomposed into its functional roles: initial anchor, structural branch, transition gate, boundary factor, and residual correction. Each formula is grounded line-by-line in the corresponding model.py source file. NLS

Let 𝐷 = cosh(4𝑥) + 4 cosh(2𝑥) + 3 cos(4𝑡), 𝑁 𝑅 = cosh(3𝑥) + 3 cos(4𝑡) cosh 𝑥, 𝑁 𝐼 = 3 sin(4𝑡) cosh 𝑥.

The analytic core is 𝜓core (𝑥, 𝑡) =

4𝑒 𝑖𝑡/2 (𝑁 𝑅 + 𝑖𝑁 𝐼 ) . 𝐷

With 𝐿 = 5, 𝐸 𝑝 (𝑥) =

 𝑥 2𝑝 , 𝐿

𝑉 𝑝 (𝑥) = −

𝐿 𝐸 𝑝 (𝑥), 2𝑝

22

𝑍 𝑝 (𝑥) =

  𝐿 𝑥2 1 − 2 𝐸 𝑝 (𝑥), 2 𝐿

A Preprint

Agentic Symbolic Search

𝑇 𝑝 (𝑥; 𝑚) = 𝑚𝑉 𝑝 (𝑥) + (1 − 𝑚)𝑍 𝑝 (𝑥). The candidate uses a split endpoint tail, h i
2 𝜓(𝑥, 𝑡) = 𝜓core (𝑥, 𝑡) + 𝑎 𝜕𝑥 𝜓core (𝐿, 𝑡) 𝑇3 (𝑥; 𝑚 3 ) − 𝑐 𝑇3 (𝑥; 𝑚 3 ) − 𝑇5 (𝑥; 𝑚 5 ) 𝑒 −𝑟𝑡 . The five trainable scalars are a bounded endpoint strength 𝑎, primary-tail mix 𝑚 3 , initial-tail mix 𝑚 5 , relaxation rate 𝑟, and compensation strength 𝑐. In the implementation these are smooth transforms of raw trainable parameters so that the endpoint strength, mixes, and relaxation rate remain in controlled ranges. Allen–Cahn 2D

For disk phases 

 𝑅 𝑘 − ∥ (𝑥, 𝑦) − 𝑐 𝑘 ∥ , 𝑝 𝑘 = tanh √ 2𝜀

𝑢0 =

1 ( 𝑝 1 + 𝑝 2 + 1 − 𝑝 1 𝑝 2 ), 2

and 𝐹 [𝑢 0 ] = 𝜀 2 Δ𝑢 0 + 𝑢 0 − 𝑢 30 . The evolving distance field is
𝑆 = (1 − 𝑏(𝑡)) 𝑆union + 𝑆neck + 𝑏(𝑡)𝑆oval . The union distance is a softened maximum of two moving lobe distances,   𝑆union = 𝑤 log 𝑒 𝑑1 /𝑤 + 𝑒 𝑑2 /𝑤 , where the lobe centers relax toward the merged center and the radii shrink by separate curvature rates. The oval distance is   √︁ 𝑆oval = 𝑎 min 1 − (𝑥𝑟 /𝐴) 2 + (𝑦 𝑟 /𝐵) 2 , with rotated coordinates in the two-lobe axis frame and time-dependent major and minor axes. A localized Gaussian distance shift 𝑆neck fills the overlap neck. The full candidate is  
2 2 𝑆 𝑢 = tanh √ + 𝑒 −𝜇 (𝑡/8) 𝑢 0 − 𝑢 geom (0) + 𝜂 𝑡𝑒 − (𝑡/𝜏 𝑗 ) 𝐹 [𝑢 0 ]. 2𝜀 The 23 trainable scalars control the two shrink rates, geometry start time, center relaxation, oval blend and axes, union softness, neck-fill amplitude and scales, initial-shape memory decay, and initial right-hand-side jet. There is no neural network component in this representation. Keller–Segel

The nine-parameter candidate is 𝑢(𝑟, 𝑡) = (1 − 𝛼(𝑡))𝑢 T (𝑟, 𝑡) + 𝛼(𝑡)𝑢 B (𝑟, 𝑡).

The Taylor branch is anchored to the Gaussian initial condition and its continuous right-hand side,    2 2 𝑢 T = 𝑢 0 exp 𝑞(𝑡) 4𝑟 2 − (4 + 𝐴IC ) + 2𝐴IC 𝑒 −𝑟 , 𝑢 0 = 9.5𝑒 −𝑟 . The blow-up branch is 𝑢B =

2 2 𝑀𝑐 𝑀in 𝑀out + 2 𝑒 − (𝑟/𝜎ℎ ) + 2 𝑒 − (𝑟/𝜎𝑜 ) , 𝜆2 (1 + (𝑟/𝜆) 2 ) 2 𝜎ℎ 𝜎𝑜

with 𝜆 = 𝑓𝑐 + (1 − 𝑓𝑐 ) ((𝑇 ∗ − 𝑡)/𝑇 ∗ ) 𝛾 . The coefficient 𝑀𝑐 corresponds to physical mass 𝜋𝑀𝑐 ; the halo coefficients split the excess over the near-critical core, and 𝛼(𝑡) is a learned smooth transition from the Taylor branch to the blow-up branch.

23

A Preprint

Agentic Symbolic Search

2

Graveleau Let 𝑣 0 = 1.6(1 − 𝑒 − ( (𝑟 −1)/0.6) )1𝑟 ≥1 , 𝐹1 = 𝑣 0 Δ𝑟 𝑣 0 + (𝑣 0′ ) 2 , and 𝐹2 = 𝜕𝑡2 𝑣| 𝑡=0 computed from the same pressure PDE. The anchor is 𝑣 anchor = 𝑣 0 + ℎ1 (𝑡; 𝑑1 )𝐹1 + 𝑏 2 ℎ2 (𝑡; 𝑑2 )𝐹2 , where 𝑏 2 scales the second-order Taylor term. The self-similar branch is 𝑣 sim = 𝐴𝜏 2𝛽−1 𝐹 (𝜂),

𝜏=

𝑇𝑓 − 𝑡 , 𝑇𝑓

𝐹 (𝜂) = 𝛽𝑠(1 + 𝑠) 1−1/𝛽 exp(𝑐 1 𝑦 + 𝑐 2 𝑦 2 + 𝑐 3 𝑦 3 ),

𝜂=

𝑟 , 𝜏𝛽

𝑠 = (𝜂 − 1)+ ,

𝑦=

𝑠 . 1+𝑠

The full candidate is 𝑣 = 𝑣 anchor + 𝑝 𝑚 𝑔(𝑡)𝜙(𝑟) (𝑣 sim − 𝑣 anchor ) + 𝑔(𝑡)𝜙(𝑟)𝜎𝑐 𝑁 𝜃 (𝑧),   𝑟 22 𝜙(𝑟) = 1 − . 6 The best representation has 𝛽 = 0.928055, while the reference free-boundary fit is 𝛽 = 0.877058. The candidate’s profile flexibility enters through the scalar amplitude, blend gates, three log-shape coefficients, and the additive MLP correction. gCLM The gCLM stress-test candidate is a Taylor-anchored neural tail. The first line below is the initial sinusoid and its first PDE tangent, the second line adds a short-time quadratic Taylor correction, and the last line appends a trainable antisymmetric tail plus higher-order Fourier Taylor modes. With 𝜃 = 𝜋𝑥/8, 𝜏 = 𝑡/𝑇max (where 𝑇max = 0.99𝑇 ∗ is the scored time horizon), and 𝐶 (𝜏) = (1 − 0.90𝜏) −1 , it uses 𝜔 = 10 sin 𝜃 − 37.5𝑡 sin(2𝜃)
1 + 𝑡 2 𝑞(𝑡) 70.3125 sin 𝜃 + 257.8125 sin(3𝜃) 2 + 𝜏 3 𝐶 (𝜏) 𝐴 𝜃 (𝑥, 𝑡) + 𝑃Tay (𝑥, 𝑡), where 𝑞(𝑡) is a short quadratic-onset ramp, 𝐴 𝜃 is an antisymmetrized MLP tail, and 𝑃Tay is a trainable combination of the stored 𝑡 3 , 𝑡 4 , 𝑡 5 Fourier Taylor modes. A.8

Numerical Reference Solutions

The validation relative 𝐿 2 errors reported in the Results section are computed against independently generated numerical reference solutions. These references are not supplied to the search; they are used only for offline validation after the run. Table 6 summarizes the numerical method, grid, and convergence checks for each case. Table 6: Numerical references used only for offline validation. Case

Method

Grid

Tolerances or step

Validation note The generator script is not archived; the stored metadata reports the method and tolerances above. Independent pseudospectral reference used for offline validation. 1024 versus 1600 grid check requires max weighted relative 𝐿 2 < 10 −2 . 800/1600/3200 grid comparison with frontwindow relative 𝐿 2 and Richardson estimate.

NLS

Periodic Fourier method of lines with DOP853

1024 internal modes, 512 stored points, 81 stored times

rtol 10 −10 , atol 10 −11

Allen–Cahn 2D

Fourier pseudospectral ETD-RK4

128 × 128, 20 stored snapshots

fixed Δ𝑡 = 0.004, 2000 steps

Keller–Segel

Conservative finite volume on a geometric radial grid, Radau time integration Uniform radial finite volume for density, BDF with sparse Jacobian, pressure stored as 𝑣 = 2𝑢 Fourier pseudospectral Hilbert transform and zero-mean velocity, Radau time integration

1600 radial cells, 40 stored times 1600 radial cells, 80 stored times

rtol 10 −9 , atol 10 −11

1024 periodic points, 400 stored times

rtol 10 −12 , atol 10 −14

Graveleau

gCLM

rtol 10 −8 , atol 10 −10

256/512/1024 spectral grid comparison by max absolute and relative-max differences.

The Graveleau exponent 𝛽ref = 0.877058 is not an analytic constant. It is obtained by fitting the late-time free-boundary radius over 20 time samples in [0.80, 0.995] of the event focus time to 𝑅(𝑡) = 𝐶 (𝑇 − 𝑡) 𝛽 , with 𝐶, 𝑇, and 𝛽 fitted by nonlinear least squares. 24

A Preprint

Agentic Symbolic Search

A.9

Self-Similar Baseline

For singular cases with a known self-similar structure, we implement a specialized self-similar PINN (SS-PINN) following the approach of Wang et al. (2025). Rather than learning the singular dynamics in physical spatiotemporal coordinates, SS-PINN is given the analytical self-similar coordinate transformation as a rigid mathematical prior, so that the MLP learns a well-conditioned stationary profile. The physical solution is then reconstructed from this profile and the prescribed scaling, and validation 𝐿 2 is computed on the full physical grid. This baseline is reported for Keller–Segel and gCLM, where a similarity reduction is known. It is not applied to bounded cases (NLS, Allen–Cahn), which do not exhibit singular self-similar structure, or to Graveleau, where the self-similar structure describes only the local focusing interface: even an oracle single-profile fit in similarity coordinates yields a full-grid 𝐿 2 worse than a standard MLP trained in physical coordinates. Architecture and optimization The profile 𝑈 (𝜂) is parameterized by a fully connected MLP with four hidden layers of 128 neurons each, using tanh activations. Inputs are affine-normalized to [−1, 1]; weights are initialized by Xavier normal and biases by zero. This yields exactly 49,921 trainable parameters in double precision (float64). Three architecture variants (4×64, 4×128, 6×128) are swept and the best validation 𝐿 2 is reported. Optimization proceeds in two phases under a wall-clock budget of 600 s. An Adam warm-up (Kingma and Ba, 2015) with learning rate 10−3 runs for 5,000 gradient steps, followed by L-BFGS (Nocedal, 1980) with strong-Wolfe line search and history size 100 until the time budget is exhausted. Gradients are clipped at ℓ 2 -norm 1,000. Profile equations Rather than routing the MLP output through the original-coordinate evaluator, SS-PINN trains the profile network directly on the analytically substituted stationary profile equation in similarity coordinates. This is the most favorable setting for the baseline: the profile equation is a well-conditioned ODE (or pseudo-ODE for the nonlocal gCLM case), free of the multi-scale blow-up dynamics that make the physical-coordinate PDE difficult. Keller–Segel. With 𝜏 = 𝑇 ∗ − 𝑡 (𝑇 ∗ = 0.9628), 𝛼 = 1, 𝛽 = 12 , 𝜂 = 𝑟/𝜏 1/2 , and 𝑢(𝑟, 𝑡) = 𝜏 −1𝑈 (𝜂), the stationary profile equation is ∫ 𝜂 1 + 𝑀 (𝜂) ′ 𝜂 ′ 𝑈 ′′ + 𝑈 − 𝑈 + 𝑈 2 − 𝑈 = 0, 𝑀 (𝜂) = 𝑠 𝑈 (𝑠) 𝑑𝑠. 𝜂 2 0 The cumulative mass 𝑀 (𝜂) is evaluated by trapezoid quadrature on 10,000 uniformly spaced 𝜂 points in [10−3 , 𝜂max ], with 𝜂max ≈ 96.7. The total loss is ℓpde + 1000 (ℓsym + ℓtail ) + 100 ℓnorm + ℓpos , where ℓsym enforces origin symmetry 𝑈 ′ (0) = 0, ℓtail penalizes nonzero far-tail values, ℓnorm anchors 𝑈 (0) to a reference-derived normalization estimate to prevent the trivial 𝑈 ≡ 0 solution, and ℓpos penalizes negative density. gCLM. With 𝜏 = 𝑇 ∗ − 𝑡 (𝑇 ∗ = 0.2409), 𝛼 = 1, 𝛽 = 34 , 𝑎 = 0.25, and the similarity coordinate centered at the periodic boundary (𝑥 center = 8) so that 𝜉 = (𝑥 − 𝑥 center )/𝜏 3/4 and 𝜔(𝑥, 𝑡) = 𝜏 −1 Ω(𝜉), the stationary profile equation is 𝛼 Ω + 𝛽 𝜉 Ω′ + 𝑎 U Ω′ − 𝐻 [Ω] Ω = 0,

U ′ = 𝐻 [Ω],

where 𝐻 denotes the Hilbert transform, computed by periodic FFT on a centered 𝜉 grid of 10,000 points spanning [−𝜉max , 𝜉max ]. The total loss is ℓpde + 1000 (ℓcenter + ℓtail ) + 1000 ℓnorm , where ℓcenter enforces Ω(0) = 0, ℓtail penalizes nonzero far-tail values, and ℓnorm anchors Ω(−1) to a reference-derived normalization target. Scope The prescribed coordinate transformation is the decisive inductive bias: it maps the finite-time blow-up into a bounded stationary-profile fitting task that a standard MLP can handle. In the paper comparisons, this transform is counted as a rigid mathematical prior and is not included in the trainable parameter count; the reported 49,921 parameters refer solely to the profile MLP weights. A.10

Detailed Result Metrics

Table 7 separates the two quantities that appear together in the main Results table. The 𝐿 2 -best representation is chosen only after the run, using the independent numerical reference. The physics-best point is the candidate with the smallest positive PDE-residual loss under the case-specific score definition. These points can differ, which is why the figures plot both the physics-loss trajectory and the validation 𝐿 2 trajectory. The iteration columns identify the outer-loop

25

A Preprint

Agentic Symbolic Search

iteration in which the selected candidate was evaluated; because an iteration may evaluate more than one candidate, equal iteration numbers do not imply the same candidate. Table 7: Detailed separation of offline 𝐿 2 validation and PDE-residual physics loss at the candidate level. Lower values indicate lower loss for both columns.

Case

NLS 1.342 Allen–Cahn 2D 0.182 Keller–Segel 0.991 Graveleau 0.00389 gCLM 0.567

3 4 8 4 10

Phys. loss Phys.-best at 𝐿 2 best

iter.

Min. Rel. 𝐿 2 phys. loss at phys. best

0.00585 1.71×10 −5 0.0107 8.41×10−5 0.188 4.49×10 −3 0.00132 2.13×10−3 0.465 3.06×10 −4

4 5 10 0 9

7.26×10−6 5.18×10−5 4.39×10−3 4.95×10−4 1.49×10−4

Initial 𝐿 2 -best Best rel. 𝐿 2 iter. rel. 𝐿 2

26

0.00630 0.0108 0.849 0.00389 0.482