What Should a Large Language Model See? Physical Invariants as a Data Representation for PDE Discovery

# What Should a Large Language Model See? Physical Invariants as a Data Representation for PDE Discovery
Source: [https://arxiv.org/html/2608.25189](https://arxiv.org/html/2608.25189)
Matt ThomsonAffiliation:Division of Biology and Biological EngineeringAffiliation:California Institute of TechnologyAffiliation:Pasadena, CA 91125

###### Abstract

Understanding how molecular interactions govern macroscopic behaviour is a central challenge in molecular sciences\. However, conventional theory building cannot keep pace with the vast datasets modern experimentation routinely produces\. Large language models offer a promising route to automating theory construction, but a spatiotemporal field cannot be directly placed in a prompt\. Existing models generally learn about the data only through a score measuring how well each proposal fits it\. Here we introduce data interpretation, a stage that measures the field into the quantities a theorist would consult and supplies them to the model as a direct input\. On a benchmark of simulated fields, interpretation nearly triples the accuracy of recovered equations relative to showing the raw data, at negligible computational cost and without any training\. By allowing a language model to read field data as a theorist does, data interpretation offers a practical route to automated field theory construction that can coevolve with experimentation\.

## 1Introduction

Today, deriving a theoretical model for a complex chemical system typically requires months to years of specialized expertise, while experimentalists can routinely generate thousands of distinct systems by varying and recombining molecular components\. This gap arises from the bottleneck of human bandwidth\. Conventional modeling approaches fall into two categories: bottom\-up methods that coarse\-grain microscopic interactions through kinetic theory into continuum field equations, and top\-down methods that enumerate all symmetry\-permitted terms without reference to molecular details\. Both approaches are slow, expertise\-dependent, and unable to match the volume and complexity of data that modern experimentation now routinely produces\. AI, with its capacity for pattern recognition, feature extraction, and learning complex mappings across scientific domains, offers a compelling path toward automating this theory\-building process — at the scale and speed that molecular sciences now critically demand\.

In this paper, we focus on automating the top\-down approach: recovering partial differential equation \(PDE\)\-based field theories from data\. This is an important task in scientific machine learning, and the dominant formulation has been sparse regression\. SINDy[Brunton et al\. \(2016\)](https://arxiv.org/html/2608.25189#bib.bib9)and PDE\-FIND[Rudy et al\. \(2017\)](https://arxiv.org/html/2608.25189#bib.bib8)assemble a library of candidate terms into a design matrix, apply numerical regression to find the best\-fit candidates, and use a sparsity penalty to limit the length of PDEs\. However, this approach is limited by the library of derivative candidates, as the PDE terms must lie within the library written down in advance\.

Symbolic regression removes the library constraint at the cost of a combinatorial search\. Recent progress shows that pretrained large language models \(LLMs\), equipped with strong scientific reasoning capabilities, have great potential as the proposal mechanism for symbolic regression\. FunSearch[Romera\-Paredes et al\. \(2024\)](https://arxiv.org/html/2608.25189#bib.bib7)established the scheme of an LLM proposing candidates, a numerical evaluator scoring them, and an evolutionary algorithm driving iterative refinement\. LLM\-SR[Shojaee et al\. \(2025\)](https://arxiv.org/html/2608.25189#bib.bib6)applied this LLM\-guided evolutionary approach to scientific equation discovery\. Since then, many methods have been proposed to refine either the search or the feedback signal: in addition to evolving the equations, LASR[Grayeli et al\. \(2024\)](https://arxiv.org/html/2608.25189#bib.bib5)also evolves a library of hypotheses that can guide the equation discovery; DrSR[Wang et al\. \(2025\)](https://arxiv.org/html/2608.25189#bib.bib4)adds a data analysis LLM to learn data insights, and reflective feedback to store reusable strategies; IGSR[Saveliev et al\. \(2026\)](https://arxiv.org/html/2608.25189#bib.bib3)expands the evaluation metrics by scoring per\-term influence inside a tree search; LLM\-ACES[Abhyankar et al\. \(\)](https://arxiv.org/html/2608.25189#bib.bib2)incorporated active data acquisition to better evaluate and evolve candidate hypotheses; MARICL[Rezaei and Krishnan \(2026\)](https://arxiv.org/html/2608.25189#bib.bib1)applies agents to read high\-error examples from a base model and propose improvements\. These methods share a core architecture, shown in black in Fig\.[1](https://arxiv.org/html/2608.25189#S1.F1)\. A pre\-trained LLM proposes an equation, a numerical evaluator evaluates it, and the resulting score returns to the LLM to generate new candidates in a loop\. Within this architecture, the LLM proposes the symbolic equation structures, while numerical regression fits only the coefficients of the proposed PDE structure\. Such a division of labor eliminates the need for a predefined library for PDE term candidates\. Effort in the field has mainly concentrated on the evaluator and the search, through improved differentiation schemes, residual formulations, credit assignment, and exploration strategies\.

We focus on a distinct aspect that has received less attention: how the data is represented to the LLM\. We are motivated by how human theorists build theories from field data — by identifying salient structures in the data\. For example, an experienced theorist knows that a shock wave indicates a strong convection effect, and that a Turing pattern signals a diffusion\-reaction system\. Neither inference requires numerical values of the field, only an interpretation of its pattern\. We therefore propose supplying such an interpretation stage to the model as an input in its own right, rather than solely relying on the evaluator\.

![Refer to caption](https://arxiv.org/html/2608.25189v1/fig1.png)Figure 1:Architecture of symbolic regression to discover PDE\-based field theory from data\.Established methods pair an LLM with an evaluator: the LLM proposes a candidate equation, and the evaluator fits its coefficients and returns a score\. Because the field data typically cannot be placed in a prompt, that score is what informs the LLM about the data\. We contribute two interpretation stages, shown in green\. Data interpretation measures the observed field into physical quantities a theorist would compute\. In our ongoing work, we apply residual interpretation to analyze what a proposal failed to explain\. No training is required in this pipeline\.We introduce a stage called data interpretation, whose output is supplied directly to the LLM, as shown in green in Fig\.[1](https://arxiv.org/html/2608.25189#S1.F1)\. Feeding the raw field data to the LLM is prohibitively expensive, so most existing research uses it only within the evaluator\. Data interpretation instead lets the LLM read the data as a theorist does\. A set of measurements is computed and supplied as evidence, chosen to be the quantities a theorist would consult when encountering an unfamiliar field\. These diagnostics are cheap and compact: they take a small fraction of the time for a single model call, and occupy fewer tokens than the raw data they replace\. They are also interpretable, so what the LLM was told can be audited by both human and AI\. The present work implements a deliberately small set of diagnostics, which are sufficient to test the efficacy of the proposal rather than to exhaust it\. In our future work, we will extend the interpretation module with more comprehensive diagnostics\.

## 2Methodology

The pipeline of PDE discovery consists of four stages in a loop: “Data interpretation→\\toLLM proposing an equation→\\toParsing and canonicalisation→\\toEvaluation", as shown in Fig\.[1](https://arxiv.org/html/2608.25189#S1.F1)\. All experiments in this paper use QwQ\-32B\. The model is frozen throughout\.

### 2\.1Data interpretation

We apply spectral analysis as the main component of data interpretation\. We consider a scalar fieldu⁡\(𝒙,t\)u\(\\bm\{x\},t\)on a periodic domain, governed by a constant\-coefficient PDE of first or second order in time\. The Fourier transform

u^𝒌\(t\)=∫u\(𝒙,t\)e−i𝒌⋅𝒙d𝒙\\hat\{u\}\_\{\\bm\{k\}\}\(t\)=\\int u\(\\bm\{x\},t\)e^\{\-i\\bm\{k\}\\cdot\\bm\{x\}\}\\mathrm\{d\}\\bm\{x\}\(1\)can decouple the field datau⁡\(𝒙,t\)u\(\\bm\{x\},t\)into a series of spatial modesu^𝒌​\(t\)\\hat\{u\}\_\{\\bm\{k\}\}\(t\)at wavevector𝒌\\bm\{k\}\. For a linear PDE, each modeu^𝒌​\(t\)\\hat\{u\}\_\{\\bm\{k\}\}\(t\)evolves independently, obeying an ordinary differential equation whose parameters are the coefficients we seek\. The LLM can therefore infer the coefficients directly from how each mode grows, decays, or oscillates\. On the contrary, if energy transfers between modes, the decoupling will break down, which indicates a nonlinear PDE\.

The interpretation addresses four questions, each with corresponding measurements, as listed below\.

Q1: First or second order in time?For a PDE of first order in time, e\.g\.,∂u∂t=ν​∇2u\+σ​u\\frac\{\\partial u\}\{\\partial t\}=\\nu\\nabla^\{2\}u\+\\sigma u, whereν\\nuandσ\\sigmaare coefficients, each Fourier mode satisfies

d​u^𝒌d​t=\(σ−ν​\|𝒌\|2\)​u^𝒌,sou^𝒌​\(t\)=u^𝒌​\(0\)​e\(σ−ν​\|𝒌\|2\)​t,\\frac\{\\mathrm\{d\}\\hat\{u\}\_\{\\bm\{k\}\}\}\{\\mathrm\{d\}t\}=\(\\sigma\-\\nu\|\\bm\{k\}\|^\{2\}\)\\hat\{u\}\_\{\\bm\{k\}\},\\qquad\\mbox\{so\}\\quad\\hat\{u\}\_\{\\bm\{k\}\}\(t\)=\\hat\{u\}\_\{\\bm\{k\}\}\(0\)e^\{\(\\sigma\-\\nu\|\\bm\{k\}\|^\{2\}\)t\},\(2\)Therefore, the amplitude\|u^𝒌​\(t\)\|\|\\hat\{u\}\_\{\\bm\{k\}\}\(t\)\|is monotonic in time\.

For a PDE of second order in time, e\.g\.,∂2u∂t2=c2​∇2u\+σ​u\+γ​∂u∂t\\frac\{\\partial^\{2\}u\}\{\\partial t^\{2\}\}=c^\{2\}\\nabla^\{2\}u\+\\sigma u\+\\gamma\\frac\{\\partial u\}\{\\partial t\}, each mode is instead a damped oscillator

d2​u^𝒌d​t2−γ​d​u^𝒌d​t\+\(c2​\|𝒌\|2−σ\)​u^𝒌=0,\\frac\{\\mathrm\{d\}^\{2\}\\hat\{u\}\_\{\\bm\{k\}\}\}\{\\mathrm\{d\}t^\{2\}\}\-\\gamma\\frac\{\\mathrm\{d\}\\hat\{u\}\_\{\\bm\{k\}\}\}\{\\mathrm\{d\}t\}\+\\left\(c^\{2\}\|\\bm\{k\}\|^\{2\}\-\\sigma\\right\)\\hat\{u\}\_\{\\bm\{k\}\}=0,\(3\)and the amplitude\|u^𝒌​\(t\)\|\|\\hat\{u\}\_\{\\bm\{k\}\}\(t\)\|is oscillatory in time\.

We measure the time evolution of\|u^𝒌​\(t\)\|\|\\hat\{u\}\_\{\\bm\{k\}\}\(t\)\|to infer the temporal order of the PDE\. However, this test is not exhaustive\. For example, an overdamped second\-order equation, withγ2≥4​\(c2​\|𝒌\|2−σ\)\\gamma^\{2\}\\geq 4\(c^\{2\}\|\\bm\{k\}\|^\{2\}\-\\sigma\)decays monotonically and is indistinguishable from a first\-order one\. This shows that data interpretation constrains the space of plausible structures rather than determining one\.

Q2: What are the linear coefficients?We measure the decay rater⁡\(𝒌\)=−d​ln⁡\|u^𝒌\|d​tr\(\\bm\{k\}\)=\-\\frac\{\\mathrm\{d\}\\ln\|\\hat\{u\}\_\{\\bm\{k\}\}\|\}\{\\mathrm\{d\}t\}to probe the PDE coefficients\. For \([2](https://arxiv.org/html/2608.25189#S2.E2)\),r⁡\(𝒌\)=ν​\|𝒌\|2−σr\(\\bm\{k\}\)=\\nu\|\\bm\{k\}\|^\{2\}\-\\sigma\. One regression across modes𝒌\\bm\{k\}therefore returns both linear coefficientsν\\nuandσ\\sigma\. For oscillatory modes, we additionally estimate the frequencyω⁡\(𝒌\)\\omega\(\\bm\{k\}\)from the spacing of successive minima inu^𝒌​\(t\)\\hat\{u\}\_\{\\bm\{k\}\}\(t\)\. For \([3](https://arxiv.org/html/2608.25189#S2.E3)\),ω2​\(𝒌\)≈c2​\|𝒌\|2−σ\\omega^\{2\}\(\\bm\{k\}\)\\approx c^\{2\}\|\\bm\{k\}\|^\{2\}\-\\sigma, and the regression can recover the coefficientc2c^\{2\}\. Sinceω⁡\(𝒌\)\\omega\(\\bm\{k\}\)is only an estimate, we do not use the intercept to determineσ\\sigma\.

Q3: Is the PDE linear or nonlinear?We measure two quantities\. The first is the fraction of spectral energy above a wavenumber threshold, compared between the start and end of a time window\. The second is the coefficient of determinationR2R^\{2\}of the rate regression in Q2\. Both quantities can probe the nonlinearity of the PDE\.

A nonlinear term generally couples different modes\. For example, the Fourier transform of an advection termu​∂u∂xu\\frac\{\\partial u\}\{\\partial x\}is

ℱ​\[u​∂u∂x\]𝒌=∑𝒌1\+𝒌2=𝒌i​𝒌2,x​u^𝒌1​u^𝒌2,\\mathcal\{F\}\\left\[u\\frac\{\\partial u\}\{\\partial x\}\\right\]\_\{\\bm\{k\}\}=\\sum\_\{\\bm\{k\}\_\{1\}\+\\bm\{k\}\_\{2\}=\\bm\{k\}\}i\\bm\{k\}\_\{2,x\}\\hat\{u\}\_\{\\bm\{k\}\_\{1\}\}\\hat\{u\}\_\{\\bm\{k\}\_\{2\}\},\(4\)where𝒌2,x\\bm\{k\}\_\{2,x\}is thexx\-component of𝒌2\\bm\{k\}\_\{2\}\. Therefore mode𝒌\\bm\{k\}is driven by products of other modes\. Energy consequently migrates toward high wavenumbers, which the first measurement detects\. And because the regression in Q2 presumes each mode decays independently, nonlinearity degrades its fit, which the coefficient of determinationR2R^\{2\}detects\. A lowR2R^\{2\}therefore indicates the PDE is likely nonlinear\.

Q4: Is there advection?We measure the gain in spatial correlation under translation\. For two time framest0<t1t\_\{0\}<t\_\{1\}, we compute

g=max𝒔⁡\[corr​\(u⁡\(𝒙,t0\),u⁡\(𝒙\+𝒔,t1\)\)−corr​\(u⁡\(𝒙,t0\),u⁡\(𝒙,t1\)\)\]g=\\max\_\{\\bm\{s\}\}\\left\[\\text\{corr\}\\left\(u\(\\bm\{x\},t\_\{0\}\),u\(\\bm\{x\}\+\\bm\{s\},t\_\{1\}\)\\right\)\-\\text\{corr\}\\left\(u\(\\bm\{x\},t\_\{0\}\),u\(\\bm\{x\},t\_\{1\}\)\\right\)\\right\]\(5\)together with the maximising shift𝒔∗\\bm\{s\}^\{\*\}\. Advection translates structure without deforming it, so a field under transport correlates poorly with itself at zero shift but strongly once shifted by the displacement, giving a largeggand a nonzero𝒔∗\\bm\{s\}^\{\*\}\. Diffusion and growth act without preferred direction, leaving the maximum at𝒔∗=𝟎\\bm\{s\}^\{\*\}=\\bm\{0\}andggnear zero\. The measurement thus separates directional transport from isotropic evolution, and𝒔∗\\bm\{s\}^\{\*\}also estimates the advection velocity given the time intervalt1−t0t\_\{1\}\-t\_\{0\}\.

### 2\.2LLM proposing an equation

The diagnostics of §[2\.1](https://arxiv.org/html/2608.25189#S2.SS1)are rendered as text, approximately 300 tokens, and appended to a prompt with two further specifications\. First, the symbols an equation may be built from: the fielduu, and its temporal and spatial derivatives\. Any algebraic combination is admissible\. We deliberately do not include a library of candidate terms\. Second, the output form: the LLM reasons in prose, then states its PDE proposal in a delimited block with explicit numerical coefficients\. This makes the answer separable from a reasoning thread that typically consists of thousands of tokens\.

### 2\.3Parsing and canonicalisation

The proposed PDE is extracted and rewritten into a single canonical form, so that expressions differing only in notation score identically\. A language model does not adhere to one convention: the same Laplacian may arrive asd​2​u/d​x​2\+d​2​u/d​y​2d2u/dx2\+d2u/dy2,∇2u\\nabla^\{2\}u, orΔ​u\\Delta u\. Canonicalisation therefore ensures that a correct proposal is scored on its content rather than on its notation\.

### 2\.4Evaluation

Given a proposed PDE structure, the coefficients are fitted by regularized least squares and the equation is scored by both a weak form residual and a sparsity penalty\. This regression only fits coefficients for a proposed PDE structure, and never selects the structure itself\. We note that the residual alone would be insufficient, because least squares can assign a near zero coefficient to a redundant term\. The sparsity penalty can select the shorter PDE when both explain the data equally well, as a theorist would\.

## 3Experiments

### 3\.1Data and benchmark

Fields are generated by numerically solving PDEs assembled from a library of eight terms\. For a scalar fieldu⁡\(𝒙,t\)u\(\\bm\{x\},t\)on a periodic domain, a first\-order sample takes the form∂u∂t=ν​∇2u\+𝒂⋅∇u\+σ​u\+β​u​∂u∂x\+r⁡\(u−u2\)\\frac\{\\partial u\}\{\\partial t\}=\\nu\\nabla^\{2\}u\+\\bm\{a\}\\cdot\\bm\{\\nabla\}u\+\\sigma u\+\\beta u\\frac\{\\partial u\}\{\\partial x\}\+r\(u\-u^\{2\}\), and a second\-order example∂2u∂t2=c2​∇2u\+σ​u\+γ​∂u∂t\\frac\{\\partial^\{2\}u\}\{\\partial t^\{2\}\}=c^\{2\}\\nabla^\{2\}u\+\\sigma u\+\\gamma\\frac\{\\partial u\}\{\\partial t\}\. During data generation, each spatial term is included or omitted at random and its coefficient drawn from a fixed range\. These terms are the standard building blocks of transport and kinetics in the molecular sciences\. Their combinations recover the reaction\-diffusion systems underlying pattern formation, the convection\-diffusion equations of flow dynamics, and the damped wave equations of acoustic and elastic responses\. Each field is integrated on a 64 by 64 grid over 50 frames\.

Table 1:Benchmark strata and accuracy results\.Strata vary along two dimensions, the temporal order and the nonlinear terms within\.*Self advection, suppressed*and*wave, weakly identifiable*are deliberately unfavourable controls where the data does not determine the answer\.Accuracy is quantified by the overlap between the proposed and true term sets\. DenotingTTas the terms of the generating equation andPPfor those of the proposal, and\|⋅\|\|\\cdot\|for the number of elements in a set, the accuracy is measured byF1=2​\|T∩P\|\|T\|\+\|P\|F\_\{1\}=2\\frac\{\|T\\cap P\|\}\{\|T\|\+\|P\|\}\. We report two metrics\. MeanF1F\_\{1\}measures how much of each equation was identified\. Exact recovery is the fraction of exact\-matching samples withP=TP=T\. We report both because at this sample size exact recovery is coarse, whileF1F\_\{1\}can be inflated by a guess containing whichever terms are common\. The fixed\-term set floor quantifies that inflation\.

Samples are drawn in fixed proportions across six strata\. The strata vary along two dimensions, the temporal order of the PDE and which nonlinear terms it contains, and are listed with the accuracy in Table[1](https://arxiv.org/html/2608.25189#S3.T1)\. Two are deliberately unfavourable\. In “self advection, suppressed", the field decays exponentially, so a term scaling asu2u^\{2\}in \([4](https://arxiv.org/html/2608.25189#S2.E4)\) is negligible for most of the observation window whether or not it appears in the equation\. In “wave, weakly identifiable", the PDE contains two terms that are difficult to distinguish in data\.

### 3\.2Results and discussion

We compare three inputs to the LLM, with results given in Table[2](https://arxiv.org/html/2608.25189#S3.T2)\. The first is the data interpretation of section[2\.1](https://arxiv.org/html/2608.25189#S2.SS1), approximately 310 tokens for each sample\. The second is a permuted interpretation, computed from a different sample, which controls against the model inferring a plausible equation without consulting its input\. The third replaces the interpretation with ten 1D slices of the raw field, approximately 1,230 tokens\. As a floor we also report the best fixed term set: the single equation achieving the highest mean accuracy across all samples, which uses no information from any individual field and requires no model call\. Table[2](https://arxiv.org/html/2608.25189#S3.T2)shows that data interpretation exceeds the floor by 0\.329 at paired Wilcoxonp=2\.6×10−5p=2\.6\\times 10^\{\-5\}and doubles exact recovery atp=6\.6×10−3p=6\.6\\times 10^\{\-3\}\.

Table 2:Accuracy by inputs\.All results use the same benchmark, model, token budget and parser, and differ only in what the model is shown about the field\. The floor is the single fixed term set achieving the highest meanF1F\_\{1\}across all samples\.Interpreting the field is substantially more effective than showing the raw data\. Data interpretation reachesF1=0\.720F\_\{1\}=0\.720against 0\.225 for raw field slices, with exact recovery on 14 of 44 samples against 2\. The gain in accuracy comes from the interpretations\. Permuting the assignment between interpretations and fields, so that each sample receives another sample’s measurements, reduces accuracy from 0\.720 to 0\.244 \(p<3×10−7p<3\\times 10^\{\-7\}\)\. This collapse shows the LLM reads the interpretations and follows them, even to a wrong answer\. In Table[2](https://arxiv.org/html/2608.25189#S3.T2), both controls fall below the optimized fixed equation, which is consistent with the LLM not relying on common terms when its input is uninformative\.

Accuracy is highest where the interpretation is most complete\. In Table[1](https://arxiv.org/html/2608.25189#S3.T1), the “linear only" stratum is recovered perfectly, consistent with the dispersion fit in Q2\. The shortfall concentrates on the nonlinear strata, where the interpretation establishes that a nonlinear term is present but cannot determine its form\.

## 4Concluding remarks

Limits\. The limits of this work are mainly on the small scale of datasets and on the insufficient characterization tools in data interpretation\. On scale, this is a workshop paper: 44 samples, one frozen model, and noise free simulated fields drawn from a library of eight terms\. That being said, the differences we report are large, and the paired tests are exact rather than asymptotic, so the sample size is adequate to establish them\. On characterization, the data interpretation defines the space of PDE structures that can be discovered\. Our characterizations establish that a nonlinear term is present without indicating its form, and the model consequently identifies those terms less accurately\. A richer set of characterization tools would expand the space of discoverable PDEs; the ones reported here were chosen to test the idea rather than to exhaust it\. The characterizations are also fragile to noise, so robust differentiation and noise filtering are required before being applied to experimental data\.

Conclusions\. Our results show that an LLM given a small set of physical measurements of a field identifies its governing equation substantially better than one given the field itself\. The measurements cost a fraction of a second per field, use fewer tokens than the raw data they replace, and require no training\. Therefore, data interpretation is a compact and useful input for LLM\-driven symbolic regression\. More broadly, this pipeline is promising to automate the top\-down route to field theory construction\. Identifying which mechanisms underlie field data, and quantifying whether a candidate equation accounts for it, are the two key steps that conventionally take most time and effort for a theorist\. In our LLM\-driven symbolic regression, the data interpretation performs the first and the evaluator the second\. Because neither is trained or tailored to a particular system, a new field can be analyzed in minutes rather than months, and every proposal is tested against the data as it is made\. Overall, this pipeline opens a path to theory construction that scales with experiment: as a laboratory varies molecular components across thousands of systems, each new dataset could automatically yield a candidate theory, already validated against the data that produced it\.

Future work\. A direct extension of our work is to apply data interpretation to the residual, as shown in dotted lines in Fig\.[1](https://arxiv.org/html/2608.25189#S1.F1)\. Given a proposed equation, the part of the dynamics it fails to account for is itself a field, and interpreting that field would indicate what the proposal is missing\. The LLM could then revise its answer, and the process could iterate until the equation residual is minimal\. Three further directions follow\. Robust differentiation would test whether the interpretation survives experimental noise, a precondition for any use on experimental data\. Prepending the data interpretation to an existing evolutionary search would test whether it improves such frameworks without modification to the search itself\. And a richer set of characterization tools would widen the space of PDEs that can be reached\.

## References

- \[1\]S\. L\. Brunton, J\. L\. Proctor, and J\. N\. Kutz\(2016\)Discovering governing equations from data by sparse identification of nonlinear dynamical systems\.Proceedings of the National Academy of Sciences113,pp\. 3932–3937\.Cited by:[§1](https://arxiv.org/html/2608.25189#S1.p2.1)\.
- \[2\]S\. H\. Rudy, S\. L\. Brunton, and J\. N\. Kutz\(2017\)Data\-driven discovery of partial differential equations\.Science Advances3,pp\. e1602614\.Cited by:[§1](https://arxiv.org/html/2608.25189#S1.p2.1)\.
- \[3\]B\. Romera\-Paredes, M\. Barekatain, A\. Novikov, M\. Balog, M\. P\. Kumar, E\. Dupont, F\. J\. R\. Ruiz, J\. S\. Ellenberg, P\. Wang, O\. Fawzi, P\. Kohli, and A\. Fawzi\(2024\)Mathematical discoveries from program search with large language models\.Nature625,pp\. 468–475\.Cited by:[§1](https://arxiv.org/html/2608.25189#S1.p3.1)\.
- \[4\]P\. Shojaee, K\. Meidani, S\. Gupta, A\. B\. Farimani, and C\. K\. Reddy\(2025\)LLM\-SR: scientific equation discovery via programming with large language models\.arXiv preprint arXiv:2404\.18400\.Cited by:[§1](https://arxiv.org/html/2608.25189#S1.p3.1)\.
- \[5\]A\. Grayeli, A\. Sehgal, O\. Costilla\-Reyes, M\. Cranmer, and S\. Chaudhuri\(2024\)Symbolic regression with a learned concept library\.arXiv preprint arXiv:2409\.09359\.Cited by:[§1](https://arxiv.org/html/2608.25189#S1.p3.1)\.
- \[6\]R\. Wang, B\. Wang, K\. Li, Y\. Zhang, and J\. Cheng\(2025\)DrSR: LLM based scientific equation discovery with dual reasoning from data and experience\.arXiv preprint arXiv:2506\.04282\.Cited by:[§1](https://arxiv.org/html/2608.25189#S1.p3.1)\.
- \[7\]E\. S\. Saveliev, S\. Holt, N\. Seedat, D\. L\. Bentley, J\. Weatherall, and M\. van der Schaar\(2026\)Influence\-guided symbolic regression: scientific discovery via LLM\-driven equation search with granular feedback\.arXiv preprint arXiv:2605\.29184\.Cited by:[§1](https://arxiv.org/html/2608.25189#S1.p3.1)\.
- \[8\]N\. Abhyankar, S\. Li, S\. Kabra, N\. Ramakrishnan, Y\. Gel, and C\. K\. ReddyLLM\-aces: closed\-loop discovery of dynamical systems with llm\-guided adaptive search\.arXiv preprint arXiv:2606\.25039\.Cited by:[§1](https://arxiv.org/html/2608.25189#S1.p3.1)\.
- \[9\]M\. R\. Rezaei and R\. G\. Krishnan\(2026\)From residuals to reasons: llm\-guided mechanism inference from tabular data\.arXiv preprint arXiv:2605\.22897\.Cited by:[§1](https://arxiv.org/html/2608.25189#S1.p3.1)\.