MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models

MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models Boxiao Wang

arXiv:2607.29561v1 [cs.LG] 31 Jul 2026

National Key Laboratory of Cognition and Decision Intelligence for Complex Systems Institute of Automation, Chinese Academy of Sciences Beijing, China wangboxiao2026@ia.ac.cn

Runxiang Wang

Chongming Li School of Astronomy and Space Science, University of the Chinese Academy of Sciences Beijing, China lichongming24@mails.ucas.ac.cn

Zhiwei Chen

Yifan Zhang

Jian Cheng

Relativistic Astrophysics, Institute for Theoretical Physics Goethe University Frankfurt Frankfurt am Main, Germany zchen@itp.unifrankfurt.de

National Key Laboratory of Cognition and Decision Intelligence for Complex Systems Institute of Automation, Chinese Academy of Sciences Beijing, China yfzhang@nlpr.ia.ac.cn

National Lab of Pattern Recognition Institute of Automation, Chinese Academy of Sciences Beijing, China jcheng@nlpr.ia.ac.cn

Abstract Discovering compact and interpretable equations from observational data is essential for understanding complex scientific systems. Symbolic regression (SR) provides a general computational framework for this task. While recent Large Language Model (LLM) based SR approaches show promise, they face two key limitations. First, they lack dedicated data analysis mechanisms for uncovering variable dependencies, which reduces the efficiency of equation discovery. Second, most methods rely on single-objective evaluation focused solely on fitting error. This neglect of structural complexity and generalization often causes models to converge prematurely to local optima, limiting their ability to explore the broader equation space. We propose Multi-Objective Tool-augmented Symbolic Regression (MOT-SR1 ), a unified framework that integrates external analytical tools to extract structural priors and guide equation generation, while jointly optimizing for accuracy, complexity, and generalization via a multi-objective evaluation module that maintains a dynamic Pareto front. MOT-SR employs two collaborative LLM modules: a Meta Strategy Generator, which selects tools and synthesizes structural optimization strategies based on Pareto-optimal equations, and an Equation Generator, which produces new candidate equations accordingly. The system operates in a closed-loop manner, continuously refining both strategies and equation structures. Across 40 standard tasks, MOT-SR outperforms existing SR methods in accuracy, generalization, and efficiency. We further validate MOT-SR on extreme mass-ratio inspiral (EMRI)

1

Kai Li

University of the Chinese National Key Laboratory of Academy of Sciences Cognition and Decision School of Advanced Intelligence for Complex Interdisciplinary Sciences Systems Beijing, China Institute of Automation, wangrunxiang2023@ia.ac.cn Chinese Academy of Sciences Beijing, China kai.li@ia.ac.cn

Code is available at https://github.com/wswbx/MOT-SR.

orbital modeling, an important problem in space-based gravitationalwave astronomy where small local errors can accumulate substantially over long-term evolution. The discovered interpretable correction achieves the lowest trajectory-level integration error on held-out configurations. These results demonstrate the potential of MOT-SR to enable reliable modeling of long-horizon scientific dynamics.

Keywords Symbolic Regression, Extreme Mass-Ratio Inspirals, Gravitational Wave Modeling, Large Language Models

1 Introduction Symbolic Regression (SR) [34] aims to discover underlying mathematical equations from data and has long been recognized as a key methodology in scientific discovery. It has been widely applied across disciplines, from identifying physical laws [33, 45] and modeling chemical systems [7, 12], to analyzing dynamic processes in biological or economic systems [51, 60]. By generating compact and interpretable equations, SR enables both accurate prediction and deep insight into system behavior. SR has long been recognized as an NP-hard problem [59], motivating diverse algorithmic developments. Early approaches based on genetic programming [9, 49] evolve equations via mutation and crossover. Reinforcement learning [44] models SR as a sequential decision-making process. Recently, Transformer-based models have enabled end-to-end learning from data to equations [4, 20, 63]. With

Wang et al.

Iterative Loop (Generation 𝑡 → 𝑡 + 1)

Observed Data 𝑋, 𝑦

Meta Strategy Generator

Equation Generator & Evaluation

Residual Calculation

Pareto Frontier Buffer ( 𝒫𝓉)

Scientific Toolbox

𝑒1 , 𝑒2 , … , 𝑒𝑁

Diversity-Aware Parent Selection

LLM Generation Correlation FFT/Wavelet Causal

Data Analysis Prompting

LLM Reasoning

Population of Equations

Search Strategy

Update Frontier ( 𝒫𝓉+1)

Multi-Objective Evaluation NMSE

(train-ID)

NMSE

(train-OOD)

AST Length (Complexity)

Non-Dominated Sorting & Pruning

Figure 1: Overview of MOT-SR. Given observed data (𝑋 , 𝑦), MOT-SR maintains a Pareto-frontier buffer 𝒫𝑡 of candidate equations and multi-objective scores. At iteration 𝑡 , a Meta Strategy Generator analyzes residuals of Pareto-optimal candidates, applies scientific tools, and outputs search strategy. Following this search strategy, an Equation Generator & Evaluation module performs diversity-aware parent selection, generates new equations, and evaluates them by NMSE ID and 𝐷 OOD (training-derived in-domain (ID)/out-ofon 𝐷train train domain (OOD) split) and by AST length for complexity. Nondominated sorting and pruning update the frontier to 𝒫𝑡+1 , forming a closed-loop scientific equation discovery process. the rise of large language models (LLMs), methods such as LLMSR [52] and LaSR [18] leverage LLM’s in-context learning capabilities and scientific priors to perform symbolic reasoning and equation generation. Despite encouraging progress, existing LLM-based SR methods face two key limitations. First, they typically lack systematic analysis of variable dependencies and data distributions, relying instead on problem descriptions as context. This often results in poorly constrained search spaces, which limits both the efficiency and directionality of equation exploration. Second, most approaches adopt a single-objective evaluation, typically minimizing fitting error, while overlooking other critical factors such as equation complexity, generalization, and diversity. This can lead to overfitting and premature convergence to locally optimal solutions. To address these challenges, we propose Multi-Objective Toolaugmented Symbolic Regression (MOT-SR), a unified framework (Figure 1) inspired by the human scientific modeling process. Scientists typically begin by analyzing data using a variety of tools, evaluate candidate hypotheses from multiple perspectives, and iteratively refine their modeling direction accordingly. MOT-SR operationalizes this human-like modeling paradigm through two core mechanisms and a cooperative, LLM-driven evolution system. At its core, MOT-SR first employs a tool-augmented analysis mechanism that invokes a suite of analytical tools to extract variable relationships from multiple complementary dimensions, such as linearity, periodicity, and causality. These insights are converted into interpretable priors and structural constraints that guide equation generation. In parallel, MOT-SR adopts a multi-objective evaluation mechanism that jointly assesses candidate equations across three dimensions: accuracy on the in-domain (ID) training subset

ID , generalization on the out-of-domain (OOD) training subset 𝐷train OOD 𝐷train , and structural complexity measured via abstract syntax tree (AST) length. A Pareto front is maintained using non-dominated sorting to preserve high-quality equations that represent optimal trade-offs. To realize these mechanisms in an adaptive and iterative manner, MOT-SR incorporates two cooperating LLMs. Specifically, the Meta Strategy Generator analyzes residual patterns of the current Pareto-optimal candidates and leverages external scientific tools (e.g., correlation, FFT/wavelet, and causal discovery) to produces a data-driven search strategy. Guided by this strategy, the Equation Generator & Evaluation module performs diversity-aware parent selection and generates new candidate equations, which are evaluated and integrated into the Pareto-frontier buffer via nondominated sorting and pruning. We evaluate MOT-SR on five benchmarks spanning physics, chemistry, biology, and materials science, using both the open-source LLaMA-3.1 [21] and the commercial GPT-4o mini [41]. Across all tasks, MOT-SR consistently outperforms traditional SR methods and recent LLM-based baselines in terms of accuracy, generalization, and equation compactness. Beyond benchmark evaluation, we apply MOT-SR to extreme mass-ratio inspiral (EMRI) orbital evolution, a challenging problem in space-based gravitational-wave astronomy where small local errors accumulate over many orbital cycles. Using 58 configurations for equation discovery and 30 held-out configurations for a posteriori evaluation, MOT-SR identifies a compact symbolic correction that achieves the lowest trajectory-level integration error— approximately three orders of magnitude lower than a neural residual baseline and 26.8 times lower than LLM-SR. These results demonstrate that MOT-SR can recover transferable correction structures that remain reliable across unseen physical configurations and long horizons, highlighting the potential of interpretable equation discovery for demanding scientific dynamical systems.

2 Preliminaries In SR, the learning task typically starts with a dataset consisting of input-output pairs:

𝐷 = {(x𝑖 , 𝑦𝑖 )}𝑛𝑖=1 ,

x𝑖 ∈ ℝ𝑑 ,

𝑦𝑖 ∈ ℝ,

where x𝑖 denotes a 𝑑 -dimensional input vector and 𝑦𝑖 is the corresponding scalar output. The goal is to discover an explicit analytical equation 𝑓 (⋅) such that the predicted outputs 𝑦𝑖̂ = 𝑓 (x𝑖 ) closely approximate the ground truth 𝑦𝑖 . To assess the quality of a candidate equation, the normalized mean squared error (NMSE) is defined as 𝑛

NMSE(𝑓 , 𝐷) =

2

𝑓 (x𝑖 ) − 𝑦𝑖 1 ∑( ) . 𝑛 𝑖=1 𝜎𝑦

(1)

where 𝜎𝑦 is the standard deviation of the target values across the dataset 𝐷 . This metric reflects the equation’s predictive accuracy, normalized by the variance of the outputs. Beyond fitting accuracy, SR also values simplicity and generalization, seeking equations that are not only accurate but also compact and transferable to unseen domains.

MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models

Our work builds on LLM-SR, a framework that leverages LLMs to generate symbolic equations through iterative optimization. Its core pipeline includes:

each variable conditioned on others. Together, these tools guide whether nonlinear or higher-order terms should be introduced. Time-frequency analysis tools help detect periodicity and transient dynamics, which frequently occur in oscillatory and mul• Equation skeleton generation: Structured prompts contain tiscale systems. MOT-SR employs two complementary methods: task-specific information (e.g., variable names, optimizaFast Fourier Transform [8] identifies dominant global frequency tion goals, example equations), guiding the LLM to procomponents, while the wavelet transform energy spectrum [11] capduce physically plausible equation skeletons. tures localized, non-stationary fluctuations. These insights support • Parameter optimization: The skeletons’ free parameters are the inclusion of periodic terms (e.g., sin) in candidate equations. optimized (e.g., via BFGS [14] ) and scored using NMSE. Causal inference tools are used to identify whether one vari• Feedback: High-quality equations are retained and reused able may influence another in a predictive or explanatory sense. as in-context examples, enabling iterative refinement through MOT-SR adopts Granger causality [16] for detecting temporal causal feedback-driven generation. influence in linear time-series, and Convergent Cross Mapping [55] While promising, LLM-based SR methods face key limitations: for identifying latent causality in nonlinear systems with potential they lack systematic modeling of variable dependencies, leading to delays. These methods provide structural signals that enhance the structurally under-informed equations; and they rely on a singleinterpretability and explanatory power of generated equations. objective evaluation focused solely on fitting error, neglecting facDynamic complexity tools help assess intrinsic system richtors such as complexity and generalization. These shortcomings ness and redundancy, guiding the pruning of over-specified comdiminish their effectiveness in solving more challenging tasks. ponents in equations. MOT-SR employs three methods: the Lyapunov exponent [61], which measures sensitivity to initial condi3 Method tions and indicates chaotic behavior; the correlation dimension [17], To address the above limitations, we propose MOT-SR (Multi-Objective which estimates the system’s effective degrees of freedom; and DyTool-augmented Symbolic Regression), a unified framework that namic Time Warping (DTW) [48], which evaluates time-shifted simintegrates external data analysis tools, multi-objective evaluation, ilarity between variable trajectories. Together, they offer structural and cooperative LLMs to enhance equation quality and search efcues for constructing compact and robust equations. ficiency. MOT-SR first extracts structural priors by analyzing variDistribution consistency tools evaluate whether input variable relationships with diverse analytical tools. It then introduces a ables behave uniformly across different input regions. MOT-SR multi-objective evaluation mechanism that jointly considers fitting uses the Kolmogorov–Smirnov (KS) test [35] to detect distributional error, equation complexity, and generalization, dynamically mainshifts between subdomains, informing the use of region-dependent taining a Pareto front through non-dominated sorting. Finally, two structures to reflect local variations in the data distribution. complementary LLMs work in tandem: one generates structural Rather than introducing new tools, the key innovation of MOTrefinement strategies, while the other synthesizes candidate equaSR lies in enabling LLMs to autonomously invoke and coordinate tions accordingly, forming a closed-loop process that continuously these tools to extract dependency patterns among variables. These guides and improves equation discovery. insights are distilled into concise guidance that informs variable selection and function composition, thereby improving responsive3.1 Tool-Augmented Variable Analysis ness to data characteristics and enhancing the scientific plausibilThis component constructs structural priors by quantifying diverse ity of generated equations. By adaptively combining outputs from variable relationships using a suite of carefully designed analytical heterogeneous analyses, for instance by linking correlation meatools. Details are provided in Appendix I. sures with periodicity detection, MOT-SR supports more targeted Linear correlation tools are essential for identifying domiequation discovery. The impact of this tool-augmented analysis on nant variables and constructing interpretable regression structures. equation generation is substantiated by the case studies in AppenMOT-SR integrates several complementary methods to assess lindix J.3. Looking ahead, MOT-SR opens the possibility for LLMs to ear dependencies from multiple statistical angles. The Pearson corsynthesize new tools, further expanding the scope of SR research. relation coefficient [42] quantifies pairwise linear associations, especially effective for Gaussian-like data. Simple linear regression 3.2 Multi-Objective Evaluation and residual variance analysis [38] evaluate predictive capacity and To enhance search efficiency and model quality, MOT-SR adopts error stability. PCA-based explained variance [19] identifies the key a multi-objective evaluation scheme that jointly considers predicdirections of structural variance. These tools jointly establish a tive accuracy, generalization, and structural simplicity. A Paretosolid basis for linear trend detection. based selection strategy maintains a diverse set of non-dominated Nonlinear dependency tools are employed to capture comcandidate equations, improving robustness and exploration of the plex interactions essential for modeling nonlinear systems. MOTsolution space. We define the ID/OOD regions within the training SR integrates three complementary methods: the Spearman rank set 𝐷train as a spatial split (Figure 2). The detailed configuration of correlation [54] measures monotonic associations based on rank, the ID-OOD partition for 𝐷train is provided in Appendix C. offering robustness to noise and non-Gaussian data; mutual information [50] measures the overall statistical dependency between variables without assuming any parametric form; and mutual inEvaluation Metrics. Unlike prior LLM-based SR methods that optiformation regression [43] quantifies the marginal contribution of mize only fitting error, MOT-SR employs a three-fold metric:

Wang et al. 22'Q  ,'Q 

22'Q  ,'Q 

Figure 2: Illustration of the training-derived ID/OOD split ID defined used in MOT-SR. The shaded region denotes 𝐷train by per-dimension percentile bounds, and the white region OOD ID , shown on two representative denotes 𝐷train = 𝐷train ∖ 𝐷train tasks. (1) In-Domain (ID) Accuracy. Measures interpolation performance within the central region of the input space. The in-domain ID is selected from the middle intervals along each input subset 𝐷train ID , dimension (Figure 2). Accuracy is quantified using NMSE on 𝐷train denoted as NMSEID . (2) Out-of-Domain (OOD) Generalization. Measures extrapOOD olation performance in the peripheral regions 𝐷train = 𝐷train ∖ ID 𝐷train , which correspond to the white regions in Figure 2. NMSE on this subset is denoted as NMSEOOD and reflects robustness to distributional shifts, a key requirement for scientific applications. It is worth noting that NMSEOOD is calculated on the training-derived OOD region 𝐷train and does not involve any held-out test data. (3) Equation Complexity. To promote interpretability, equation complexity is measured via the size of its abstract syntax tree (AST) [40]. In MOT-SR, each equation is represented as a Python function and parsed into an AST, where each node corresponds to a variable, operator, or function call. The total node count provides a low-overhead estimate of structural complexity, guiding the search toward compact equations. Although Python AST length may disproportionately penalize composite function calls, we adopt it as a computationally inexpensive proxy for expression-tree complexity; its empirical ranking consistency and associated trade-offs are discussed in Appendix L. Pareto-Guided Multi-Objective Optimization. Based on the three metrics, MOT-SR formulates equation discovery as a multi-objective optimization problem over a candidate set 𝒫 = {𝑓1 , 𝑓2 , … , 𝑓𝑀 }, aiming to balance fitting accuracy, extrapolation, and structural simplicity: ID ), NMSE(𝑓 , 𝐷 OOD ), ASTLen(𝑓 )) 𝑓 ∗ = arg min𝑓 ∈𝒫 (NMSE(𝑓 , 𝐷train train

A function 𝑓𝑎 dominates 𝑓𝑏 if it performs no worse across all objectives and better on at least one. The Pareto front 𝒫∗ consists of all non-dominated candidates in 𝒫. To improve the quality of the front, MOT-SR filters out overly simple candidates with low equation complexity that may mislead the search by applying ID ) and NMSE(𝑓 , 𝐷 OOD ). upper-bound thresholds on NMSE(𝑓 , 𝐷train train Remaining equations are ranked via non-dominated sorting to construct the current Pareto-optimal set. Details of the dynamic threshold construction are provided in Appendix C.4.

3.3

Cooperative LLMs for Equation Evolution

To integrate variable-aware analysis with multi-objective optimization, MOT-SR employs a cooperative framework consisting of two LLMs: the Meta Strategy Generator 𝜋stg and the Equation Generator 𝜋eq . The former extracts structural search strategies, while the latter generates candidate equations accordingly, jointly driving the population toward Pareto-optimality (Algorithm 1). Meta Strategy Generator. 𝜋stg formulates equation search strategies by synthesizing variable-level insights and structural abstraction. At each iteration, given the current Pareto front 𝒫∗𝑡−1 , 𝜋stg autonomously selects a subset of tools 𝒜 ⊆ Toolbox, applies them to the dataset 𝐷train , and summarizes key relationships as natural language descriptions ℛvarrel to guide downstream generation. In parallel, 𝜋stg performs structural abstraction over 𝒫∗𝑡−1 , producing complementary guidance ℛstruct through: 1) Commonality Extraction, which identifies frequent substructures that define prevailing symbolic motifs; 2) Disparity Analysis, which diagnoses regional residual patterns to uncover structural weaknesses; and 3) Blind Spot Discovery, which detects unexplored symbolic components to encourage diversity. The final strategy 𝒮𝑡 = (ℛvarrel , ℛstruct ) is passed to equation generator 𝜋eq to steer the next round of generation. Equation Generator. Upon receiving the strategy 𝒮𝑡 from 𝜋stg , 𝜋eq is responsible for synthesizing a new batch of candidate equations 𝒫𝑡 . To promote structural diversity and prevent premature convergence, MOT-SR incorporates a structure-diversity–guided parent sampling module, which explicitly prioritizes structurally distinct candidates during equation generation. Specifically, each equation 𝑓 ∈ 𝒫∗𝑡 on the Pareto front is parsed into an abstract syntax tree (AST), from which a set of symbolic subtrees 𝑆(𝑓 ) is extracted. The pairwise structural dissimilarity between two equations is computed via a subtree-overlap metric defined as: |𝑆(𝑓𝑖 ) ∩ 𝑆(𝑓𝑗 )| . SyntaxDiv(𝑓𝑖 , 𝑓𝑗 ) = − |𝑆(𝑓𝑖 )| This score reflects the relative uniqueness of 𝑓𝑖 ’s structure compared to 𝑓𝑗 . For each equation, we compute the average structural diversity score: Scorediv (𝑓𝑖 ) = 𝑁 1−1 ∑𝑗≠𝑖 SyntaxDiv(𝑓𝑖 , 𝑓𝑗 ).

A softmax sampling procedure is then applied to construct a parent set 𝒫parent , biased toward structurally diverse equations. These parents serve as in-context examples, which, together with the strategy 𝒮𝑡 , are fed into 𝜋eq to generate new equation candidates: 𝑓 ′ ∼ 𝜋eq (𝒮𝑡 , 𝒫parent ). This procedure facilitates the generation of structurally diverse and generalizable expressions, driving the continuous evolution and expansion of the equation population 𝒫.

4 Experiments 4.1 Datasets To assess MOT-SR’s performance, we adopt two sets of challenging datasets. The first includes four standard benchmarks from LLMSR, spanning nonlinear oscillatory systems (Oscillation 1 & 2), where Oscillation 1 focuses on periodic signal composition and Oscillation 2 introduces cross-variable interactions and non-periodic

MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models

Oscillation 1 Oscillation 2 E. coli growth Stress-Strain Accavg-0.001 (%)↑ NMSE↓ Accavg-0.001 (%)↑ NMSE↓ Accavg-0.1 (%)↑ NMSE↓ Accavg-0.1 (%)↑ NMSE↓ GPlern 0.11 0.0972 0.05 0.2000 0.76 1.0023 28.43 0.3496 PySR 3.80 0.0003 7.02 0.0002 2.80 0.4068 70.60 0.0347 1.49e-6 0.43 0.0282 2.04 0.2754 76.28 0.0282 RAG-SR 39.47 uDSR 1.78 0.0002 0.36 0.0856 1.12 0.5059 59.15 0.0639 LaSR (Llama-3.1) 2.79 0.7485 1.09 0.0310 3.44 0.1349 71.84 0.0320 LLM-SR (Llama-3.1) 12.67 2.55e-5 8.20 4.70e-5 1.36 0.5815 76.21 0.0333 LLM-SR (4o-mini) 11.12 2.07e-5 8.66 4.51e-5 3.24 0.0863 71.28 0.0491 MOT-SR (Llama-3.1) 100.00 1.27e-15 99.45 1.70e-10 6.60 0.0208 85.02 0.0150 MOT-SR (4o-mini) 99.99 1.42e-13 99.57 4.25e-10 6.32 0.0178 86.33 0.0144 Table 1: Overall performance of MOT-SR and baseline methods on four benchmarks. Model

disturbances to increase modeling difficulty, an E. coli Growth task modeling multivariate biological dynamics with nonlinear couplings [37, 46], and a Stress-Strain task from materials science featuring piecewise nonlinear deformation behavior [1]. The second is LSR-Synth–Chemistry from LLM-SRBENCH [53], which consists of 36 tasks derived from a chemical kinetics base equation with progressively increasing symbolic complexity. It is specifically designed to evaluate a model’s ability to generalize across nested, unseen, and semantically rich expressions. In this work, we focus on LSR-Synth–Chemistry rather than the complete LLM-SRBENCH suite because the other three datasets are constructed variants of the four benchmark problems already included in our evaluation. By concentrating on LSR-Synth–Chemistry, which is both novel and complementary, we ensure a comprehensive yet non-redundant assessment of model performance while also taking into account the practical constraints of available hardware resources. Full dataset descriptions are provided in Appendix F for completeness.

4.2

Baselines

We compare MOT-SR against a range of representative baselines from both classical and LLM-based SR methods. For the four standard tasks in the LLM-SR benchmark, we include GPleaRn, a classical genetic programming-based SR method; PySR [18], which combines evolutionary search with symbolic compression; uDSR [27], which replaces DSR’s RNN policy with a pretrained Transformer and neural-guided decoding; RAG-SR [63], which augments equation generation with structure retrieval; and LLM-SR [52]. On the more challenging LSR-Synth–Chemistry, we compare MOT-SR with leading LLM-enhanced methods, including SGA [32], which combines LLM-based hypothesis generation with physics-informed parameter optimization via bilevel search, and LaSR [18], which extracts abstract symbolic concepts from prior equations to guide hybrid LLM-evolutionary equation generation.

4.3

Evaluation Metrics

We evaluate different methods using three metrics: (1) Accuracy to tolerance 𝜏 , denoted as Accall (𝜏 ) and Accavg (𝜏 ), and (2) Normalized Mean Squared Error (NMSE). Accall (𝜏 ) measures task-level correctness by requiring all test points to satisfy the relative error bound 𝑦̂ −𝑦

𝜏 , i.e., Accall (𝜏 ) = 1 (max1≤𝑖≤𝑁test | 𝑖 𝑦 𝑖 | ≤ 𝜏 ). Accavg (𝜏 ) computes 𝑖

Model Accall-0.1 (%)↑ NMSE↓ SGA 8.33 0.0458 LaSR 27.77 2.77e-04 LLM-SR 66.66 8.01e-06 MOT-SR 86.11 3.85e-07 Table 2: Comparison on LSR-Synth–Chemistry. the proportion of test points that meet the same criterion, defined 𝑁

𝑦̂ −𝑦

test as Accavg (𝜏 ) = 𝑁1 ∑𝑖=1 1 (| 𝑖 𝑦 𝑖 | ≤ 𝜏 ). test

𝑖

We further evaluate the quality of the final Pareto front using Hypervolume (HV) [65] and Inverted Generational Distance (IGD) [64], two standard indicators in multi-objective optimization (see Appendix H for details). HV measures the volume dominated by the obtained solution set with respect to a fixed reference point, capturing both convergence and diversity. It reflects the overall coverage of the objective space and favors solution sets that are both wellconverged and diverse. IGD computes the average distance from the ground truth equations to its nearest counterpart in the generated front, emphasizing approximation accuracy with respect to the true Pareto-optimal equations. For LLM-SR, we construct a pseudo Pareto front by extracting the nondominated set at each generation and aggregating them over 100 generations, yielding a hindsight-aggregated (best-of-run) front for a stringent comparison.

4.4

MOT-SR Configuration

For fair comparison, we adopt the same LLMs across all methods: LLaMA-3.1-8B and GPT-4o-mini, covering both lightweight and high-performance scenarios. In each iteration, MOT-SR generates four candidate equations for evaluation. The total number of iterations is set to 2000 for standard benchmarks and 1000 for LSRSynth–Chemistry, following the LLM-SRBench protocol. Traditional baselines are allowed more iterations to ensure convergence. Additional details, including prompt design and sampling configurations, are provided in Appendix M and D.

5 Findings 5.1 MOT-SR Achieves the Best Overall Performance As shown in Table 1, MOT-SR consistently outperforms both classical and LLM-based SR baselines, achieving significantly lower NMSE and higher accuracy across all benchmarks. For example, on

MOT-SR (Llama-3.1)

10 2 10 4 10 6 10 8 10 10 10 12 10 14

100

10 4 10 7 10 10 10 13 0

MOT-SR

LLM-SR ID OOD

10 1 10 2

Bactgrow

Stress-Strain

(a) Generalization across domains under ID and OOD.

1000

Iteration

1500

2000

10 3 10 5 10 7 10 9 10 11

0

500

E. coli Growth

100

0

500

1000

Iteration

1000

Iteration

1500

2000

1500

2000

Stress-Strain

10 1

10 2

NMSE Comparison on LSR Synth Chemistry

500

LLM-SR (GPT-4o-mini)

Oscillation 2

NMSE (log scale)

PySR RAG-SR

NMSE (log scale)

Oscillator2

LLM-SR (Llama-3.1) 10 1

NMSE (log scale)

NMSE (log scale)

Oscillator1

MOT-SR (GPT-4o-mini)

Oscillation 1

10 1

NMSE (log scale)

NMSE (log scale)

Wang et al.

1500

2000

10 1

10 2

0

500

1000

Iteration

Figure 4: Training convergence comparison.

LLM-SR MOT-SR ID OOD 10 14

10 12

10 10

10 8

10 6 NMSE (log scale)

10 4

10 2

100

102

(b) Boxplot of generalization on 36 chemistry tasks.

Figure 3: Comparison of generalization performance.

Oscillator 1 with LLaMA, MOT-SR reaches an NMSE of 1.27×10−15 , far surpassing LLM-SR’s 2.55 × 10−5 . With both LLaMA and GPT4o-mini, MOT-SR attains over 90% accuracy on several datasets. On the more challenging LSR-Synth–ChemistRy (Table 2), MOTSR achieves the best overall performance with an average NMSE of 3.85 × 10−7 and 86.1% accuracy, significantly outperforming other LLM-based models. These results validate the effectiveness of MOT-SR’s tool-augmented multi-objective framework in enabling the discovery of more accurate equations.

5.2

MOT-SR Shows Stronger Generalization

We evaluate generalization on the held-out test set under in-domain and out-of-domain conditions. In this section (Figure 3), “ID/OOD” ID and 𝐷 OOD , obtained by applying denotes the test partitions 𝐷test test percentile-based bounds from 𝐷train (Sec. 3.2) to the test inputs. These partitions are used only for final reporting and are not involved in parameter fitting or any discovery-loop objectives. As shown in Figure 3a, MOT-SR consistently achieves the lowest ID and 𝐷 OOD , NMSE across all four standard tasks under both 𝐷test test significantly outperforming LLM-SR and other baselines. For inOOD stance, on the test OOD split of Oscillator 1 (𝐷test ), MOT-SR reaches −14 an NMSE of 6.20 × 10 , nearly eleven orders of magnitude lower than LLM-SR (1.4 × 10−3 ). Figure 3b further shows its superior median NMSEs across 36 chemistry tasks in LSR-Synth, under both test ID and test OOD conditions. MOT-SR’s generalization advantage stems from two key mechanisms. First, our multi-objective optimization explicitly incorporates performance on a training-derived OOD region, guiding the

search toward expressions that maintain low error under distributional shifts. Second, our meta-strategy module integrates variablelevel analysis and structural diagnostics: the former identifies stable input dependencies via analytical tools, while the latter detects effective substructures and potential failure modes across the population. Together, these components help MOT-SR uncover symbolic relationships that generalize beyond distribution-specific patterns.

5.3

MOT-SR Improves Discovery Efficiency

We compare the convergence behavior of MOT-SR and LLM-SR across four benchmark tasks. As shown in Figure 4, MOT-SR not only reduces error more rapidly but also converges to lower final NMSE values. In most cases, it outperforms LLM-SR’s best results (at 2000 iterations) within the first 1000 iterations. This efficiency gain stems from MOT-SR’s meta-strategy design. By extracting key variable dependencies through analytical tools and summarizing structural patterns from the Pareto front, MOT-SR narrows the search space and avoids redundant exploration, steering the model toward high-quality equations with fewer iterations. These results highlight its efficiency in navigating the symbolic search space while maintaining both accuracy and generalization. We further provide the computational cost analysis in Appendix J.2.

5.4

MOT-SR Exhibits Superior Multi-objective Optimization Quality

We compare MOT-SR and LLM-SR on Oscillator 1 using HV and IGD. As shown in Figure 5a, MOT-SR achieves faster growth in HV compared to LLM-SR, indicating earlier and more effective discovery of diverse, high-quality non-dominated solutions. Similarly, Figure 5b shows that MOT-SR consistently reduces IGD at a faster rate, with better convergence toward the true Pareto front. LLM-SR, in contrast, plateaus earlier with higher IGD, reflecting a tendency to remain in suboptimal regions. These results demonstrate that MOT-SR significantly outperforms baselines in terms of convergence speed, solution diversity, and optimization quality.

5.5

Ablation Study

To assess the contributions of key components in MOT-SR, we conduct ablation experiments on the Oscillator 1 benchmark with

MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models

HV

1.0

MOT-SR LLMSR

0.5

0.0

0

500

1000

Iterations

IGD

1.0

1500

2000

MOT-SR LLMSR

0.5

0.0

0

500

(a) HV

1000

Iterations

1500

2000

(b) IGD

5.6

Figure 5: Comparison of HV and IGD of different models.

NMSE (log scale)

10 6

1.53e-05

ID OOD

3.99e-07

9.76e-07

1.18e-10

10 10 10 12 10 14

1.70e-07

1.55e-08

10 8

6.20e-14 1.27e-15

Full MOT-SR

w/o Data

w/o Struct

These findings demonstrate that the data-driven and structureguided feedback mechanisms offer complementary strengths in guiding MOT-SR’s equation generation. Evaluations of the ablation settings on other benchmark datasets, are provided in Appendix G.1, and further ablation results on the role of different tools in the toolset are reported in Appendix G.2.

w/o Strategy

Figure 6: Ablation results on the Oscillation 1 problem. the LLaMA backbone. Our analysis focuses on two mechanisms: the multi-objective optimization mechanism and the meta strategy generator. 5.5.1 Impact of Multi-objective Optimization. We assess the impact of multi-objective optimization by replacing it with a singleobjective variant (w/o MultiObj) that optimizes only fitting error, excluding complexity and generalization. In contrast, w/o MultiObj recovers only 1 of 9 symbolic terms, compared to 4 of 6 for MOT-SR: Ground Truth 𝑣 ̇ = 0.8 sin(𝑥) − 0.5𝑥 ⋅ 𝑣 − 0.5𝑣 3 − 0.2𝑥 3 − 𝑥 cos(𝑥).

MOT-SR w/o MultiObj 𝑣 ̇ = 𝑎0 tanh(𝑎1 𝑥) − 𝑎2 tanh(𝑎3 𝑣) − 𝑎4 𝑥 − 𝑎5 𝑣 − 𝑎6 𝑥𝑣 − 𝑎7 (𝑥 − 𝑣) − 𝑎8 (𝑥 + 𝑣) − 𝑎9 𝑣 − (1 − 𝑎9 )𝑥.

MOT-SR 𝑣 ̇ = −𝑎0 𝑥 + 𝑎1 𝑣 − 𝑎2 𝑥 3 + 𝑎3 𝑣 3 − 𝑎4 𝑥𝑣 + 𝑎5 sin(𝑎6 𝑥) .

Case Study: EMRI Orbital Evolution Correction

Problem and scientific objective. Space-based gravitational-wave detection is a major frontier in modern astronomy, opening the millihertz band to compact systems inaccessible to ground-based detectors [47]. Extreme mass-ratio inspirals (EMRIs), in which a compact object inspirals around a much more massive black hole, are key sources in this band and probe the spacetime and physical processes near the central black hole [2]. Their many orbital cycles require accurate evolution models, since small local errors can cause substantial long-term trajectory and phase deviations. We present MOT-SR as the first application of SR to EMRI orbitalevolution modeling and use this task to examine whether symbolic equation discovery can address a practical modeling problem in gravitational-wave astrophysics. Pn5AAK is widely used for its computational efficiency but may lose accuracy over long evolutions, while self-force calculations are more accurate but costly and currently limited in scope. We study eccentric equatorial orbits in Schwarzschild spacetime using Fast Self-forced Inspirals (FSI) [57] as the reference and Pn5AAK [22], denoted PN5 below, as the approximation. This comparison quantifies accumulated orbital deviations and identifies where PN5 remains reliable. We focus on correcting the evolution of the dimensionless semilatus rectum 𝑝 . Given the eccentricity 𝑒 and mass ratio 𝜂, the learning target is the discrepancy between the local evolution rates prė 𝑒, 𝜂) = 𝑝FSI ̇ (𝑝, 𝑒, 𝜂)−𝑝PN5 ̇ (𝑝, 𝑒, 𝜂). dicted by FSI and PN5, defined as Δ𝑝(𝑝, The discovered residual is added to the PN5 evolution rate and integrated over time to obtain the corrected trajectory 𝑝corr (𝑡). This setting evaluates whether an equation fitted to evolution-rate data remains reliable when repeatedly applied within a dynamical solver. Full details of the residual construction, unit conversion, equation search, and numerical integration are provided in Appendix E.

This indicates multi-objective optimization is essential for discovDataset construction. We generate Schwarzschild eccentric-EMRI ering compact and interpretable equations. In contrast, single-objective trajectories with central black hole masses between 104 and 107 𝑀⊙ , optimization often leads to redundancy and overfitting. secondary masses between 5 and 100𝑀⊙ , initial eccentricities 0 ≤ 5.5.2 Effectiveness of the Meta Strategy Generator. To further in𝑒0 ≤ 0.2, and initial semi-latus recta satisfying 6 + 2𝑒0 + 0.05 ≤ vestigate the contribution of the meta strategy generator, we con𝑝0 ≤ 12. Each case is simulated with both FSI and PN5, and the duct a series of ablation studies by progressively removing its core resulting evolution rates are aligned at the same orbital states becomponents: the data analysis submodule (w/o Data), the structure fore the residual labels are computed. Among the 100 generated analysis submodule (w/o Struct), and the entire module (w/o Stratcases, Cases 89–100 contain only one or two valid states and are egy). removed before evaluation. The remaining 88 cases are divided at As shown in Figure 6, removing the data analysis module imthe trajectory level: 58 complete trajectories are used for equation pairs strategy refinement based on variable relationships, leading discovery and parameter fitting, and 30 complete trajectories are to a notable performance drop. Excluding the structure module reheld out for testing. No trajectory contributes samples to both subsults in a noticeable decline in generalization and accuracy, as the sets. The equation structure, fitted constants, and candidate selecmodel can no longer extract structural patterns from prior equation are fixed before evaluation on the 30 held-out cases. tions to guide generation. Removing the entire meta strategy generator results in prompt-only generation without feedback, proBaselines and evaluation. We compare against NN, an adapted ducing the lowest accuracy and stability among all variants. neural residual baseline motivated by the neural UDE formulation

Wang et al. FSI NIT MOT-SR

case42

LLM-SR

8

LLM-SR NN

1.62 1.20 1.48

3.48e-03 2.65e-01 1.17e-03 3.13e-02 2.69e-03 1.86e-01

(b) Held-out cases by trajectory length Medium (500–2,000) 4.47e-01 1.32e-04 1.09e-02 8.18e-05 4.19e-03 Long (2,000–10,000) 1.47 Extra-long (> 10,000) 2.81 1.40e-03 8.56e-02

Table 3: Mean case-level NMSE of the integrated 𝑝(𝑡) trajectories. Panel (a) reports results on the discovery, held-out, and complete sets. Panel (b) groups the held-out cases by trajectory length.

introduced by Keith et al. [23]. NN uses a feed-forward ResidualMLP with hidden dimensions [128, 128, 64] and SiLU activations to predict Δ𝑝 ̇ from (𝑝, 𝑒, 𝜂). The predicted residual is added to the PN5 evolution rate before integration. We also report LLM-SR as an alternative LLM-based symbolic correction. Following the standard a posteriori evaluation protocol for scientific equation discovery, we assess the complete integrated trajectory rather than point-wise residual fitting. For each case, we apply the NMSE defined in Eq. (1) to the predicted trajectory 𝑝corr (𝑡) and the FSI reference 𝑝FSI (𝑡). Results. Table 3 reports the mean case-level NMSE. Lower values indicate closer agreement with the FSI reference trajectory. On the 30 held-out trajectories, MOT-SR achieves a mean NMSE of 1.17 × 10−3 , approximately three orders of magnitude lower than NN and 26.8 times lower than LLM-SR. MOT-SR also consistently achieves the lowest mean NMSE among the three methods on the 58 discovery cases and across all 88 valid cases, demonstrating its accuracy under complete-trajectory integration. More importantly, the evaluation split is constructed at the case level: each held-out trajectory represents a distinct EMRI configuration that is entirely excluded from equation discovery. The shared symbolic correction discovered from 58 cases remains effective when integrated over 30 unseen configurations. This result provides evidence that MOT-SR generalizes across EMRI cases, instead of relying on samples drawn from trajectories encountered during discovery. To examine error accumulation over longer integration horizons, Table 3 groups the held-out trajectories with more than 500 valid integration points by trajectory length. The length-stratified results reveal clear differences in long-term integration behavior. The mean NMSE of NN increases monotonically from 4.47 × 10−1 on medium trajectories to 1.47 on long trajectories and 2.81 on extra-long trajectories. This trend indicates that small errors in the neural residual accumulate when the fitted correction is repeatedly applied during integration. In comparison, MOT-SR maintains a mean NMSE between 8.18 × 10−5 and 1.40 × 10−3 and achieves the lowest error in every trajectory-length group.

FSI NIT MOT-SR

9 case53

LLM-SR NN

p

8

7

7

10 1 10 5

100 10 5

0.000

0.001

0.002

t (yr)

(a) Case 42

0.003

|p pFSI |/|pFSI |

MOT-SR

p

(a) Overall evaluation Discovery cases (58) Held-out cases (30) All valid cases (88)

NN

|p pFSI |/|pFSI |

Evaluation group

0.000

0.005

0.010

t (yr)

0.015

0.020

(b) Case 53

Figure 7: A posteriori integration for two held-out EMRI configurations, comparing FSI with NN, LLM-SR, and MOT-SR. Quantitative results are based on all 30 held-out cases.

This behavior is consistent with the training design of MOT-SR. MOT-SR divides 𝒟train into ID and OOD subsets and evaluates candidate equations for both ID accuracy and OOD generalization during discovery. This internal generalization criterion discourages equations whose accuracy is restricted to a narrow subset of the discovery configurations, yielding a correction that remains reliable under repeated long-term integration. Compared with LLMSR, MOT-SR reduces the mean NMSE by factors of 82.6, 51.2, and 61.1 on the medium, long, and extra-long trajectories, respectively. The consistent advantage over LLM-SR shows that producing a symbolic expression alone does not fully realize the potential of SR for this task; the ID/OOD-aware discovery process of MOT-SR further exploits the ability of symbolic equations to capture correction structures shared across EMRI configurations. Collectively, these results show that MOT-SR preserves low integration error across different trajectory lengths and unseen EMRI configurations. Figure 7 provides two illustrative examples. In both cases, the trajectory produced by MOT-SR remains closer to the FSI reference over the plotted integration horizon, while the neural residual leads to a visibly larger deviation. Together with the aggregate results, these examples provide evidence of trajectory-level transfer across unseen EMRI configurations within the considered Schwarzschild equatorial eccentric regime.

6 Conclusion We introduced MOT-SR, a unified symbolic regression framework combining tool-augmented variable analysis, multi-objective evaluation, and cooperative LLMs for scientific equation discovery. Across standard benchmarks, MOT-SR consistently improves predictive accuracy, generalization, equation compactness, and search efficiency over classical and LLM-based baselines. Beyond benchmark tasks, MOT-SR discovers a compact symbolic correction for EMRI orbital evolution that substantially reduces long-horizon trajectorylevel integration error and remains effective across held-out configurations. These results demonstrate the potential of tool-augmented, multi-objective equation discovery for building interpretable and transferable models of complex scientific dynamics.

Limitations and Ethical Considerations MOT-SR currently relies on a fixed set of evaluation objectives and manually specified analytical tools. Future work will introduce domain-specific objectives, automate tool discovery and refinement, and extend MOT-SR to high-dimensional scientific data

MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models

and more complex dynamical systems. This work involves no human participants or personal data. The benchmark and astrophysical data are used solely for scientific equation discovery and model evaluation.

References [1] B. S. Aakash, JohnPatrick Connors, and Michael D Shields. 2019. Stress-strain data for aluminum 6061-T651 from 9 lots at 6 temperatures under uniaxial and plane strain tension. Data in Brief 25 (Aug 2019), 104085. doi:10.1016/j.dib.2019. 104085 [2] Pau Amaro-Seoane, Jeff Andrews, Manuel Arca Sedda, Abbas Askar, Quentin Baghi, Razvan Balasov, Imre Bartos, Simone S Bavera, Jillian Bellovary, Christopher PL Berry, et al. 2023. Astrophysics with the laser interferometer space antenna. Living Reviews in Relativity 26, 1 (2023), 2. [3] Johan Andersson. 2000. A survey of multiobjective optimization in engineering design. Department of Mechanical Engineering, Linktjping University. Sweden (2000), 38. [4] Luca Biggio, Tommaso Bendinelli, Alexander Neitz, Aurelien Lucchi, and Giambattista Parascandolo. 2021. Neural Symbolic Regression that Scales. arXiv:2106.06427 [cs.LG] https://arxiv.org/abs/2106.06427 [5] L. Biggio*, T. Bendinelli*, A. Neitz, A. Lucchi, and G. Parascandolo. 2021. Neural Symbolic Regression that Scales. In Proceedings of 38th International Conference on Machine Learning (ICML 2021) (Proceedings of Machine Learning Research, Vol. 139). PMLR, 936–945. *equal contribution. https://proceedings.mlr.press/ v139/biggio21a.html [6] Sandra C Cerda-Flores, Arturo A Rojas-Punzo, and Fabricio Nápoles-Rivera. 2022. Applications of multi-objective optimization to industrial processes: a literature review. Processes 10, 1 (2022), 133. [7] Jindou Chen, Jidong Tian, Liang Wu, ChenXinWei, Xiaokang Yang, Yaohui Jin, and Yanyan Xu. 2025. KinFormer: Generalizable Dynamical Symbolic Regression for Catalytic Organic Reaction Kinetics. In International Conference on Representation Learning, Y. Yue, A. Garg, N. Peng, F. Sha, and R. Yu (Eds.), Vol. 2025. 67058–67080. https://proceedings.iclr.cc/paper_files/paper/2025/file/ a76b693f36916a5ed84d6e5b39a0dc03-Paper-Conference.pdf [8] James W. Cooley and John W. Tukey. 1965. An Algorithm for the Machine Calculation of Complex Fourier Series. Math. Comput. 19 (1965), 297–301. doi:10. 1090/S0025-5718-1965-0178586-1 [9] Miles Cranmer. 2023. Interpretable Machine Learning for Science with PySR and SymbolicRegression.jl. arXiv:2305.01582 [astro-ph.IM] https://arxiv.org/ abs/2305.01582 [10] Laure Crochepierre, Lydia Boudjeloud-Assala, and Vincent Barbesant. 2022. Interactive Reinforcement Learning for Symbolic Regression from Multi-Format Human-Preference Feedbacks. In IJCAI 2022- 31st International Joint Conference on Artificial Intelligence. Vienne, Austria. https://hal.science/hal-03695471 [11] Ingrid Daubechies. 1992. Ten Lectures on Wavelets. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, USA. [12] Song Deng, Junjie Wang, Li Tao, Su Zhang, and Hongwei Sun. 2023. EV charging load forecasting model mining algorithm based on hybrid intelligence. Computers and Electrical Engineering 112 (2023), 109010. doi:10.1016/j.compeleceng. 2023.109010 [13] Mengge Du, Yuntian Chen, and Dongxiao Zhang. 2023. DISCOVER: Deep identification of symbolically concise open-form PDEs via enhanced reinforcementlearning. arXiv:2210.02181 [cs.LG] https://arxiv.org/abs/2210.02181 [14] Roger Fletcher. 1987. Practical Methods of Optimization (2nd ed.). John Wiley & Sons, Chichester, New York. [15] Jenna C Fromer and Connor W Coley. 2023. Computer-aided multi-objective optimization in small molecule discovery. Patterns 4, 2 (2023). [16] C. W. J. Granger. 1969. Investigating Causal Relations by Econometric Models and Cross-spectral Methods. Econometrica 37, 3 (1969), 424–438. http://www. jstor.org/stable/1912791 [17] Peter Grassberger and Itamar Procaccia. 1983. Measuring the strangeness of strange attractors. Physica D: Nonlinear Phenomena 9, 1 (1983), 189–208. doi:10. 1016/0167-2789(83)90298-1 [18] Arya Grayeli, Atharva Sehgal, Omar Costilla-Reyes, Miles Cranmer, and Swarat Chaudhuri. 2024. Symbolic Regression with a Learned Concept Library. arXiv:2409.09359 [cs.LG] https://arxiv.org/abs/2409.09359 [19] I. T. Jolliffe. 2002. Principal Component Analysis (2nd ed.). Springer‑Verlag, New York, NY, USA. doi:10.1007/b98835 [20] Pierre-Alexandre Kamienny, Stéphane d’Ascoli, Guillaume Lample, and François Charton. 2022. End-to-end symbolic regression with transformers. arXiv:2204.10532 [cs.LG] https://arxiv.org/abs/2204.10532 [21] Paul Kassianik, Baturay Saglam, Alexander Chen, Blaine Nelson, Anu Vellore, Massimo Aufiero, Fraser Burch, Dhruv Kedia, Avi Zohary, Sajana Weerawardhena, Aman Priyanshu, Adam Swanda, Amy Chang, Hyrum Anderson, Kojin Oshiba, Omar Santos, Yaron Singer, and Amin Karbasi. 2025. Llama-3.1-FoundationAI-SecurityLLM-Base-8B Technical Report. arXiv:2504.21039 [cs.CR] https://arxiv.org/abs/2504.21039 [22] Michael L Katz, Alvin JK Chua, Lorenzo Speri, Niels Warburton, and Scott A Hughes. 2021. Fast extreme-mass-ratio-inspiral waveforms: New tools for millihertz gravitational-wave data analysis. Physical Review D 104, 6 (2021), 064047. [23] Brendan Keith, Akshay Khadse, and Scott E. Field. 2021. Learning orbital dynamics of binary black hole systems from gravitational wave measurements.

Wang et al.

Phys. Rev. Research 3 (Nov 2021), 043101. Issue 4. doi:10.1103/PhysRevResearch. 3.043101 [24] Michael Kommenda, Andreas Beham, Michael Affenzeller, and Gabriel Kronberger. 2015. Complexity measures for multi-objective symbolic regression. In International Conference on Computer Aided Systems Theory. Springer, 409–416. [25] J.R. Koza. 1990. Genetically breeding populations of computer programs to solve problems in artificial intelligence. In [1990] Proceedings of the 2nd International IEEE Conference on Tools for Artificial Intelligence. 819–827. doi:10.1109/TAI.1990. 130444 [26] Jiří Kubalík, Erik Derner, and Robert Babuška. 2021. Multi-objective symbolic regression for physics-aware dynamic modeling. Expert Systems with Applications 182 (2021), 115210. [27] Mikel Landajuela, Chak Shing Lee, Jiachen Yang, Ruben Glatt, Claudio P Santiago, Ignacio Aravena, Terrell Mundhenk, Garrett Mulcahy, and Brenden K Petersen. 2022. A Unified Framework for Deep Symbolic Regression. In Advances in Neural Information Processing Systems, S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (Eds.), Vol. 35. Curran Associates, Inc., 33985–33998. https://proceedings.neurips.cc/paper_files/paper/2022/file/ dbca58f35bddc6e4003b2dd80e42f838-Paper-Conference.pdf [28] Mikel Landajuela, Brenden K. Petersen, Soo K. Kim, Claudio P. Santiago, Ruben Glatt, T. Nathan Mundhenk, Jacob F. Pettit, and Daniel M. Faissol. 2021. Improving exploration in policy gradient search: Application to symbolic optimization. arXiv:2107.09158 [cs.LG] https://arxiv.org/abs/2107.09158 [29] Wenqiang Li, Weijun Li, Linjun Sun, Min Wu, Lina Yu, Jingyi Liu, Yanjie Li, and Song Tian. 2023. Transformer-based model for symbolic regression via joint supervised learning. In International Conference on Learning Representations. https://api.semanticscholar.org/CorpusID:259298765 [30] Xi Lin, Zhiyuan Yang, and Qingfu Zhang. 2022. Pareto Set Learning for Neural Multi-objective Combinatorial Optimization. arXiv:2203.15386 [cs.LG] https: //arxiv.org/abs/2203.15386 [31] Thibaut Lust and Jacques Teghem. 2012. The multiobjective multidimensional knapsack problem: a survey and a new approach. International Transactions in Operational Research 19, 4 (2012), 495–520. arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1475-3995.2011.00840.x doi:10.1111/j.1475-3995.2011.00840.x [32] Pingchuan Ma, Tsun-Hsuan Wang, Minghao Guo, Zhiqing Sun, Joshua B. Tenenbaum, Daniela Rus, Chuang Gan, and Wojciech Matusik. 2024. LLM and Simulation as Bilevel Optimizers: A New Paradigm to Advance Physical Scientific Discovery. In Proceedings of the 41st International Conference on Machine Learning (Proceedings of Machine Learning Research, Vol. 235), Ruslan Salakhutdinov, Zico Kolter, Katherine Heller, Adrian Weller, Nuria Oliver, Jonathan Scarlett, and Felix Berkenkamp (Eds.). PMLR, 33940–33962. https://proceedings.mlr.press/ v235/ma24m.html [33] Nour Makke and Sanjay Chawla. 2024. Data-driven discovery of Tsallis-like distribution using symbolic regression in high-energy physics. PNAS Nexus 3, 11 (10 2024), pgae467. arXiv:https://academic.oup.com/pnasnexus/articlepdf/3/11/pgae467/60816181/pgae467.pdf doi:10.1093/pnasnexus/pgae467 [34] Nour Makke and Sanjay Chawla. 2024. Interpretable scientific discovery with symbolic regression: a review. Artificial Intelligence Review 57 (01 2024). doi:10. 1007/s10462-023-10622-0 [35] Frank J. Massey. 1951. The Kolmogorov-Smirnov Test for Goodness of Fit. J. Amer. Statist. Assoc. 46, 253 (1951), 68–78. http://www.jstor.org/stable/2280095 [36] Matteo Merler, Katsiaryna Haitsiukevich, Nicola Dainese, and Pekka Marttinen. 2024. In-Context Symbolic Regression: Leveraging Large Language Models for Function Discovery. In Proceedings of the 62nd Annual Meeting of the Association for Computational Linguistics (Volume 4: Student Research Workshop). Association for Computational Linguistics, 589–606. doi:10.18653/v1/2024.acl-srw.49 [37] Jacques Monod. 1949. THE GROWTH OF BACTERIAL CULTURES. Annual Review of Microbiology 3, Volume 3, 1949 (1949), 371–394. doi:10.1146/annurev. mi.03.100149.002103 [38] Douglas C. Montgomery, Elizabeth A. Peck, and G. Geoffrey Vining. 2013. Introduction to Linear Regression Analysis (5 ed.). Wiley, Hoboken, NJ. [39] T. Nathan Mundhenk, Mikel Landajuela, Ruben Glatt, Claudio P. Santiago, Daniel M. Faissol, and Brenden K. Petersen. 2021. Symbolic Regression via Neural-Guided Genetic Programming Population Seeding. arXiv:2111.00053 [cs.NE] https://arxiv.org/abs/2111.00053 [40] Iulian Neamtiu, Jeffrey S. Foster, and Michael Hicks. 2005. Understanding source code evolution using abstract syntax tree matching. In Proceedings of the 2005 International Workshop on Mining Software Repositories (St. Louis, Missouri) (MSR ’05). Association for Computing Machinery, New York, NY, USA, 1–5. doi:10.1145/1083142.1083143 [41] OpenAI. 2024. GPT‑4o mini: advancing cost‑efficient intelligence. https: //openai.com/index/gpt-4o-mini-advancing-cost-efficient-intelligence/. Accessed: 2025‑07‑29. [42] Karl Pearson. 1895. Note on Regression and Inheritance in the Case of Two Parents. Proceedings of the Royal Society of London 58 (1895), 240–242. http: //www.jstor.org/stable/115794

[43] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, Jake Vanderplas, Alexandre Passos, David Cournapeau, Matthieu Brucher, Matthieu Perrot, and Édouard Duchesnay. 2011. Scikit‑learn: Machine Learning in Python. Journal of Machine Learning Research 12 (2011), 2825–2830. [44] Brenden K. Petersen, Mikel Landajuela, T. Nathan Mundhenk, Claudio P. Santiago, Soo K. Kim, and Joanne T. Kim. 2021. Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. arXiv:1912.04871 [cs.LG] https://arxiv.org/abs/1912.04871 [45] Julia Reuter, Hani Elmestikawy, Fabien Evrard, Sanaz Mostaghim, and Berend van Wachem. 2023. Graph Networks as Inductive Bias for Genetic Programming: Symbolic Models for Particle-Laden Flows. In Genetic Programming, Gisele Pappa, Mario Giacobini, and Zdenek Vasicek (Eds.). Springer Nature Switzerland, Cham, 36–51. [46] L Rosso, J. R. Lobry, S Bajard, and J. P. Flandrois. 1995. Convenient Model To Describe the Combined Effects of Temperature and pH on Microbial Growth. Applied and Environmental Microbiology 61, 2 (Feb 1995), 610–6. doi:10.1128/ aem.61.2.610-616.1995 [47] Wen-Hong Ruan, Zong-Kuan Guo, Rong-Gen Cai, and Yuan-Zhong Zhang. 2020. Taiji program: Gravitational-wave sources. International Journal of Modern Physics A 35, 17 (2020), 2050075. [48] H. Sakoe and S. Chiba. 1978. Dynamic programming algorithm optimization for spoken word recognition. IEEE Transactions on Acoustics, Speech, and Signal Processing 26, 1 (1978), 43–49. doi:10.1109/TASSP.1978.1163055 [49] Michael Schmidt and Hod Lipson. 2009. Distilling Free-Form Natural Laws from Experimental Data. Science 324, 5923 (2009), 81– 85. arXiv:https://www.science.org/doi/pdf/10.1126/science.1165893 doi:10.1126/ science.1165893 [50] C. E. Shannon. 1948. A mathematical theory of communication. The Bell System Technical Journal 27, 3 (1948), 379–423. doi:10.1002/j.1538-7305.1948.tb01338.x [51] Hao Shi, Weili Song, Xinting Zhang, Jiahe Shi, Cuicui Luo, Xiang Ao, Hamid Arian, and Luis Seco. 2024. AlphaForge: A Framework to Mine and Dynamically Combine Formulaic Alpha Factors. arXiv:2406.18394 [q-fin.CP] https://arxiv. org/abs/2406.18394 [52] Parshin Shojaee, Kazem Meidani, Shashank Gupta, Amir Barati Farimani, and Chandan K Reddy. 2025. LLM-SR: Scientific Equation Discovery via Programming with Large Language Models. arXiv:2404.18400 [cs.LG] https://arxiv.org/ abs/2404.18400 [53] Parshin Shojaee, Ngoc-Hieu Nguyen, Kazem Meidani, Amir Barati Farimani, Khoa D Doan, and Chandan K Reddy. 2025. LLM-SRBench: A New Benchmark for Scientific Equation Discovery with Large Language Models. arXiv:2504.10415 [cs.CL] https://arxiv.org/abs/2504.10415 [54] C. Spearman. 1904. The Proof and Measurement of Association between Two Things. The American Journal of Psychology 15, 1 (1904), 72–101. http://www. jstor.org/stable/1412159 [55] George Sugihara, Robert May, Hao Ye, Chih hao Hsieh, Ethan Deyle, Michael Fogarty, and Stephan Munch. 2012. Detecting Causality in Complex Ecosystems. Science 338, 6106 (2012), 496– 500. arXiv:https://www.science.org/doi/pdf/10.1126/science.1227079 doi:10.1126/science.1227079 [56] Mojtaba Valipour, Bowen You, Maysum Panju, and Ali Ghodsi. 2021. SymbolicGPT: A Generative Transformer Model for Symbolic Regression. arXiv:2106.14131 [cs.LG] https://arxiv.org/abs/2106.14131 [57] Maarten van de Meent and Niels Warburton. 2018. Fast self-forced inspirals. Classical and Quantum Gravity 35, 14, Article 144003 (July 2018), 144003 pages. arXiv:1802.05281 [gr-qc] doi:10.1088/1361-6382/aac8ce [58] Martin Vastl, Jonáš Kulhánek, Jiří Kubalík, Erik Derner, and Robert Babuška. 2024. SymFormer: End-to-End Symbolic Regression Using Transformer-Based Architecture. IEEE Access 12 (2024), 37840–37849. doi:10.1109/ACCESS.2024. 3374649 [59] Marco Virgolin and Solon P. Pissis. 2022. Symbolic Regression is NP-hard. arXiv:2207.01018 [cs.NE] https://arxiv.org/abs/2207.01018 [60] Ylva Wahlquist, Jesper Sundell, and Kristian Soltesz. 2024. Learning pharmacometric covariate model structures with symbolic regression networks. Journal of Pharmacokinetics and Pharmacodynamics 51, 2 (2024), 155–167. doi:10.1007/ s10928-023-09887-3 [61] Alan Wolf, Jack B. Swift, Harry L. Swinney, and John A. Vastano. 1985. Determining Lyapunov exponents from a time series. Physica D: Nonlinear Phenomena 16, 3 (1985), 285–317. doi:10.1016/0167-2789(85)90011-9 [62] Shunyu Yao, Fei Liu, Xi Lin, Zhichao Lu, Zhenkun Wang, and Qingfu Zhang. 2025. Multi-objective Evolution of Heuristic Using Large Language Model. arXiv:2409.16867 [cs.AI] https://arxiv.org/abs/2409.16867 [63] Hengzhe Zhang, Qi Chen, Bing XUE, Wolfgang Banzhaf, and Mengjie Zhang. 2025. RAG-SR: Retrieval-Augmented Generation for Neural Symbolic Regression. In The Thirteenth International Conference on Learning Representations. https://openreview.net/forum?id=NdHka08uWn

MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models

[64] Eckart Zitzler, Kalyanmoy Deb, and Lothar Thiele. 2000. Comparison of Multiobjective Evolutionary Algorithms: Empirical Results. Evolutionary Computation 8, 2 (2000), 173–195. doi:10.1162/106365600568202 [65] E. Zitzler and L. Thiele. 1999. Multiobjective evolutionary algorithms: a comparative case study and the strength Pareto approach. IEEE Transactions on Evolutionary Computation 3, 4 (1999), 257–271. doi:10.1109/4235.797969

Wang et al.

Appendix A Related Work A.1 Symbolic Regression Traditional SR methods mainly rely on evolutionary algorithms, reinforcement learning [44], and Transformers [5]. For instance, genetic programming [25] formulates equation discovery as an evolutionary search over tree-based representations, refining structures via mutation and crossover. Reinforcement learning-based symbolic regression was first introduced by Petersen et al. [44] and has since developed into various policy-optimization frameworks [10, 13, 28, 39]. More recently, Transformer-based models [20, 29, 56, 58, 63] have been applied to SR, leveraging large-scale pretraining to improve equation generation. However, these models typically lack mechanisms to incorporate physical priors. With advances in natural language processing, LLM-based SR methods such as LLM-SR [52], LaSR [18], and ICSR [36] have emerged. LLM-SR utilizes scientific priors encoded in LLMs to generate plausible equation forms, followed by data-driven parameter fitting. LaSR introduces abstract concept generation to guide hypothesis construction, while ICSR formats training data as in-context prompts to induce function generation. Yet these approaches still rely heavily on pre-trained knowledge and lack explicit modeling of variable relationships or structured reasoning over data, which limits both their search efficiency and generalization ability.

A.2

Multi-objective Optimization

Multi-objective optimization is a fundamental research direction in the optimization community and has been widely applied in areas such as combinatorial optimization [30, 31, 62] and industrial design [3, 6, 15]. In the context of SR, Kommenda et al. systematically compared various complexity measures for multi-objective SR and proposed a new metric that preserves semantic information while improving search efficiency [24]. Kubalík et al. introduced a physics-aware multi-objective approach [26] that jointly optimizes model accuracy and physical consistency, enhancing interpretability and reliability. MOT-SR is the first framework to integrate LLMs into multi-objective SR, enabling data-driven equation discovery with improved accuracy, simplicity, and generalization.

B

Algorithmic Pseudocode for MOT-SR

Algorithm 1 MOT-SR 1: Input: Training set 𝐷 ; Maximum iterations 𝑇 ; 𝑛 samples per iteration; Initial population 𝒫0 (optional); Meta Strategy Generator 𝜋stg ;

Equation Generator 𝜋eq ; Toolbox: Data analysis toolset ∗ 2: Output: Approximate Pareto front 𝒫 3: Initialize 𝒫0

4: for 𝑡 = 1, … , 𝑇 do

𝑅𝑒𝑠 ← EvaluateEquationResiduals(𝒫∗𝑡−1 , 𝐷) 6: ℛvarrel ← 𝜋stg .DataAnalysis(Toolbox, 𝑅𝑒𝑠, 𝒫∗𝑡−1 ) 7: ℛstruct ← 𝜋stg .GenerateStructuralPrompts(𝒫∗𝑡−1 ) 8: 𝒮𝑡 ← (ℛvarrel , ℛstruct ) 9: 𝒫parent ← ParentSelection(𝒫∗𝑡−1 ) 10: 𝒫𝑡 ← 𝒫∗𝑡−1 11: for 𝑖 = 1, … , 𝑛 do 12: 𝑓 ∼ 𝜋 eq (𝒮𝑡 , 𝒫parent ) 13: 𝑠𝑐𝑜𝑟𝑒 ← MultiObjectiveEvaluation(𝑓 , 𝐷) 14: 𝒫𝑡 ← 𝒫𝑡 ∪ {𝑓 , 𝑠𝑐𝑜𝑟𝑒} 15: end for 16: 𝒫∗𝑡 ← PopulationManagement(𝒫𝑡 ) 17: end for ∗ 18: 𝒫 ← 𝒫𝑇 5:

C Data Partition and Thresholding Protocol for ID/OOD Evaluation C.1 Motivation and Protocol Definition In scientific discovery, a valid symbolic model must not only fit the observed data (interpolation) but also capture the underlying physical laws to predict system behaviors in unobserved regions (extrapolation). To rigorously evaluate this capability, we adopt a percentile-based masking strategy to partition the global data distribution 𝒟total into In-Domain (ID) and Out-of-Domain (OOD) subsets. Unlike random splitting, which results in independent and identically distributed (i.i.d.) subsets, our protocol enforces a strictly defined distributional shift between training (ID) and testing (OOD), thereby assessing the model’s robustness against extrapolation risks.

MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models

C.2

Partition Methodology

For each task, let x ∈ ℝ𝑑 denote the input variables and 𝑦 the target. The protocol is defined based on per-dimension marginal percentiles of the training inputs in 𝐷train . Percentile candidates. We iterate over percentile values

𝑝 ∈ {10, 11, … , 49}.

(2)

Per-dimension central intervals. For each input dimension 𝑗 ∈ {1, … , 𝑑}, we compute lower and upper bounds from the training inputs:

ℓ𝑗 (𝑝) = percentile(𝑋𝑗 , 𝑝),

𝑢𝑗 (𝑝) = percentile(𝑋𝑗 , 100 − 𝑝),

(3)

where 𝑋𝑗 denotes the set of training values of the 𝑗 -th coordinate. Training-derived ID region as a Cartesian product. We define the training-derived ID subset as the Cartesian product of these perdimension central intervals: ID (𝑝) = x, 𝑦) ∈ 𝐷 𝐷train (4) {( train ∶ ℓ𝑗 (𝑝) ≤ 𝑥𝑗 ≤ 𝑢𝑗 (𝑝), ∀𝑗}. ID (𝑝) forms a central hyper-rectangle in the input space (Figure 2). Geometrically, 𝐷train

Training-derived OOD region as the complement within 𝐷train . The training-derived OOD subset is defined as the spatial complement of the ID region within the training set: OOD ID (𝑝). 𝐷train (𝑝) = 𝐷train ∖ 𝐷train (5) ID (𝑝). This subset corresponds to peripheral regions along at least one dimension, inducing an explicit distributional shift relative to 𝐷train

Selecting the split percentile. Among candidate percentiles, we select 𝑝 ⋆ to balance the sample counts in the two regions:

𝑝 ⋆ = arg

max

𝑝∈{10,…,49}

OOD ID (𝑝)|, min (|𝐷train (𝑝)|), |𝐷train

which encourages a non-trivial OOD subset while retaining sufficient ID samples for stable fitting. Finally, we set OOD OOD ⋆ ID ∶= 𝐷 ID (𝑝 ⋆ ), 𝐷train 𝐷train ∶= 𝐷train (𝑝 ). train

C.3

(6)

(7)

Usage in MOT-SR and No-Leakage Guarantee

MOT-SR uses the above training-derived split as follows: ID . • Parameter fitting (discovery loop). Equation constants are fitted only on 𝐷train OOD • Generalization feedback (discovery loop). Extrapolation feedback is computed on 𝐷train to guide multi-objective selection, but this subset is never used for parameter fitting. ID and 𝐷 OOD ) are held out and used only for final reporting. • Final evaluation (locked test). Benchmark test sets (e.g., 𝐷test test OOD Therefore, the “OOD” objective used inside the discovery loop refers strictly to the training-derived region 𝐷train (Eq. 5), and does not involve any held-out test data.

C.4

Dynamic NMSE Thresholding for Pre-Filtering

To avoid flooding the Pareto selection with clearly underfitting candidates, we apply a loose upper-bound NMSE threshold in each generation. For every candidate 𝑖, we compute its worst-case training error across the two training-derived regions: (𝑖)

(𝑖)

𝑒𝑖 = max (NMSEID , NMSEOOD ),

(8)

and define the current best worst-case error in the population as

𝑒 ⋆ = min 𝑒𝑖 .

(9)

𝑖

We then construct a dynamic upper bound ⋆

𝑇 = 100.5 log10 (𝑒 ) = √𝑒 ⋆ ,

(10)

so that 𝑇 is on the order of √𝑒 ⋆ in log-scale. Any candidate satisfying (𝑖)

(𝑖)

max (NMSEID , NMSEOOD ) > 𝑇

(11)

is treated as clearly underfitting; in implementation we assign it a very large objective value so it cannot enter the Pareto front. This step does not penalize low-NMSE solutions, and only removes candidates whose errors are orders of magnitude worse than the current best.

Wang et al.

D MOT-SR Configuration and Language Model Details D.1 MOT-SR Configuration We implement MOT-SR with both open-source and commercial LLM backbones: LLaMA-3.1-8B-Instruct and GPT-4o-mini. The LLaMA3.1 model is quantized and deployed locally on NVIDIA H100 80GB GPUs for efficient inference. And GPT-4o-mini is accessed via the OpenAI API, providing high-quality reasoning without requiring local resources. Following the LLM-SRBench protocol, the maximum number of optimization iterations is set to 2000 for the four main benchmark datasets, and 1000 for the LSR-Synth–Chemistry dataset. In each iteration, the Meta Strategy Generator selects 3 tools from the tool set to analyze variable-level relationships and generates a natural language strategy prompt based on both data characteristics and symbolic structural patterns. Both the Meta Strategy Generator and Equation Generator in MOT-SR utilize LLM decoding with top-𝑘 = 30, top-𝑝 = 0.3, and temperature 0.6. In each iteration, the Equation Generator produces 4 candidate expressions. During Pareto frontier construction, we apply an NMSE-based filtering mechanism to eliminate trivial expressions: for each candidate, we compute the base-10 logarithm of its minimal NMSE (across ID and OOD), and discard it if it exceeds 10⌈0.5⋅min log10 NMSE⌉ , where the minimum is taken over the current candidate pool.

D.2 Large Language Models MOT-SR employs two backbone LLMs: LLaMA-3.1-8B-Instruct and GPT-4o-mini, covering both open-source and API-accessible commercial models. The LLaMA-3.1 model is locally quantized to 4-bit precision and deployed on NVIDIA H100 80GB GPUs, requiring approximately 8GB of VRAM during inference, thus allowing execution on consumer-grade hardware with appropriate quantization. The GPT-4o-mini model is accessed via OpenAI’s API, offering high-quality reasoning with minimal deployment overhead.

E Symbolic Equation Discovery for EMRI Orbital Evolution E.1 EMRI Orbital-Evolution Dataset Here, we provide the preprocessing details used to construct the symbolic regression dataset. Each case is simulated using both FSI and PN5. The resulting evolution rates are aligned at the same orbital states (𝑝, 𝑒, 𝜂) before the residual labels are computed. If the output points of the two models do not coincide, they are aligned through interpolation or recomputation. The original data retain time and evolution rates in seconds. To express the rates with respect to the dimensionless time 𝑡/𝑀 , we define

𝐺𝑀 . 𝑐3 For 𝑋 ∈ {FSI, PN5}, the dimensionless semi-latus-rectum evolution rate is 𝑀sec =

(12)

𝑝𝑋̇ = 𝑀sec 𝑝𝑋̇ sec ,

(13)

where 𝑝𝑋̇ sec = 𝑑𝑝/𝑑𝑡sec . The symbolic regression target is then constructed as

̇ sec − 𝑝PN5 ̇ sec ) . Δ𝑝 ̇ = 𝑀sec (𝑝FSI

(14)

The resulting regression samples have the form

𝜇 . (15) 𝑀 Among the 100 generated cases, Cases 89–100 contain only one or two valid states and are removed before evaluation. The remaining 88 cases are split at the trajectory level: 58 complete trajectories are used for equation discovery and parameter fitting, while 30 complete trajectories are held out for testing. No trajectory contributes samples to both subsets. The equation structure, fitted constants, and candidate selection are fixed before evaluation on the held-out cases. Only 𝑝 , 𝑒 , and 𝜂 are used as inputs to the symbolic regression model. The masses 𝑀 and 𝜇 determine the mass ratio 𝜂, while 𝑀 additionally defines the time-conversion factor 𝑀sec . Extrinsic parameters, including redshift, sky location, source orientation, distance, and initial phases, are excluded because they do not determine the orbit-averaged evolution rate considered here. (𝑝, 𝑒, 𝜂) ⟼ Δ𝑝,̇

E.2

𝜂=

Formulating EMRI Orbital Evolution as Symbolic Regression

Predicting a local correction instead of a trajectory. Directly regressing the full trajectory 𝑝(𝑡) would entangle the local modeling error with numerical integration and the selected initial conditions. We instead search for a local and configuration-independent symbolic residual

𝑓𝜃 ∶ (𝑝, 𝑒, 𝜂) ⟼ Δ𝑝,̇

(16)

where 𝑓𝜃 is an explicit symbolic expression and 𝜃 denotes its fitted constants. The symbolic hypothesis is trained to approximate the rate discrepancy in Equation (14): ̇ (𝑝, 𝑒, 𝜂) − 𝑝PN5 ̇ (𝑝, 𝑒, 𝜂). 𝑓𝜃 (𝑝, 𝑒, 𝜂) ≃ 𝑝FSI (17) This formulation asks the symbolic model to identify the systematic part of the approximation error shared across EMRI configurations, rather than memorizing a separate time-dependent correction for each trajectory.

MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models

Corrected orbital-evolution rate. Given a candidate symbolic expression, the corrected PN5 evolution rate is reconstructed as

̇ (𝑝, 𝑒, 𝜂) = 𝑝PN5 ̇ (𝑝, 𝑒, 𝜂) + 𝑓𝜃 (𝑝, 𝑒, 𝜂). 𝑝corr

(18)

When the orbital evolution is integrated in physical time, the corresponding rate is

̇ sec = 𝑝corr

̇ 𝑝PN5 + 𝑓𝜃 (𝑝, 𝑒, 𝜂) . 𝑀sec

(19)

The symbolic correction is therefore repeatedly evaluated at the current orbital state during numerical integration, rather than being applied once as a post-processing adjustment. Symbolic fitting and multi-objective evaluation. For a set of residual samples 𝒟, the point-wise prediction error of a candidate expression is measured by 2 ∑𝑖 [𝑓𝜃 (𝑝𝑖 , 𝑒𝑖 , 𝜂𝑖 ) − Δ𝑝𝑖̇ ] . (20) NMSEΔ𝑝 ̇ (𝑓𝜃 , 𝒟) = 2 ̇ ] ∑𝑖 [Δ𝑝𝑖̇ − Δ𝑝 𝒟 Within the 58-trajectory discovery set, MOT-SR applies the same training-derived ID/OOD partition used in the general framework. The free constants are optimized using the ID subset, while the OOD subset supplies extrapolation feedback. Candidate expressions are then jointly evaluated according to (21) (NMSEID , NMSEOOD , ASTLen) , so that the search favors expressions that balance local fitting accuracy, cross-region generalization, and structural compactness. A posteriori trajectory evaluation. A low point-wise residual error does not necessarily guarantee an accurate long-term orbit, because small systematic errors can accumulate when the learned expression is repeatedly evaluated inside the dynamical solver. We therefore perform an a posteriori evaluation by inserting Equation (18) into the orbital-evolution procedure and integrating from the initial conditions of each complete trajectory. This produces the corrected trajectory 𝑝corr (𝑡). For each EMRI case, trajectory-level accuracy is measured against the FSI reference using 2

NMSE𝑝(𝑡) =

∑𝑘 [𝑝corr (𝑡𝑘 ) − 𝑝FSI (𝑡𝑘 )] ∑𝑘 [𝑝FSI (𝑡𝑘 ) − 𝑝 FSI ]

2

.

(22)

The equation structure, fitted constants, and final candidate selection are fixed before integration on the 30 held-out configurations. Hence, Equation (22) evaluates whether a single symbolic correction discovered from the 58 discovery cases remains reliable under repeated integration across unseen EMRI systems and different trajectory lengths.

F Datasets F.1 Nonlinear Oscillator Oscillatory systems with nonlinear damping are foundational in physics and engineering for modeling the motion of objects subjected to restoring and dissipative forces. These systems are governed by second-order differential equations of the form 𝑥 ̈ + 𝑓 (𝑡, 𝑥, 𝑥)̇ = 0, where the nonlinear function 𝑓 encapsulates dynamic interactions among position, velocity, and possibly time. To evaluate a model’s ability to recover such complex dynamics, two synthetic oscillator tasks are used. The first system is defined by:

𝑣 ̇ = 𝐹 sin(𝜔𝑥) − 𝛼𝑣 3 − 𝛽𝑥 3 − 𝛾 𝑥𝑣 − 𝑥 cos(𝑥) (𝐹 = 0.8, 𝛼 = 0.5, 𝛽 = 0.2, 𝛾 = 0.5, 𝜔 = 1.0) The second system follows:

𝑣 ̇ = 𝐹 sin(𝜔𝑡) − 𝛼𝑣 3 − 𝛽𝑥𝑣 − 𝛿𝑥 exp(𝛾 𝑥) (𝐹 = 0.3, 𝛼 = 0.5, 𝛽 = 1.0, 𝛿 = 5.0, 𝛾 = 0.5, 𝜔 = 1.0) with initial conditions 𝑥 = 0.5, 𝑣 = 𝑥 ̇ = 0.5, and simulation time 𝑡 ∈ [0, 50]. These equations exhibit rich nonlinear structures and variable couplings, making them ideal benchmarks for testing symbolic reasoning and generalization beyond simple oscillatory behavior.

F.2 Bacterial Growth Accurately capturing the dynamics of E. coli proliferation under varying environmental conditions is of critical importance in areas such as biotechnology and microbiological risk assessment. This benchmark reflects realistic yet challenging biological modeling, where growth depends on multiple interacting factors. The dataset models population change via: 𝑑𝐵 = 𝑓𝐵 (𝐵) ⋅ 𝑓𝑆 (𝑆) ⋅ 𝑓𝑇 (𝑇 ) ⋅ 𝑓pH (pH) 𝑑𝑡 where 𝐵 is bacterial density, 𝑆 is nutrient concentration, 𝑇 is temperature, and pH denotes acidity. Each term accounts for a distinct physiological influence on growth.

Wang et al.

The explicit expression used is:

(pH − pHmin )𝜋 tanh 𝑘(𝑇 − 𝑥0 ) 𝑑𝐵 𝑆 = 𝜇max 𝐵 ( )( ) ) exp (−|pH − pHopt |) sin ( 4 𝑑𝑡 𝐾𝑆 + 𝑆 pHmax − pHmin 1 + 𝑐(𝑇 − 𝑥decay )

2

This formulation introduces complex nonlinearities and multimodal interactions across environmental axes, challenging models to integrate biological structure while avoiding rote memorization.

F.3 Material Stress Behavior To evaluate symbolic models under realistic experimental conditions, this benchmark focuses on the stress-strain response of Aluminum 6061-T651 under thermal influence. The dataset records tensile strength measurements at six different temperatures, ranging from room temperature to 300℃, simulating diverse material states. In contrast to synthetic equations, this task lacks an explicit ground-truth formula, requiring data-driven inference of hidden physical laws. It serves as a test of a model’s capacity to identify empirical regularities in noisy, high-variance regimes. A widely used approximation of this behavior is given by:

𝜎 = (𝐴 + 𝐵𝜀 𝑛 ) (1 − (

𝑇 − 𝑇𝑟 𝑚 ) ) 𝑇𝑚 − 𝑇𝑟

where 𝜎 denotes stress, 𝜀 is strain, 𝑇 is the temperature, 𝑇𝑟 is a reference point, and 𝑇𝑚 is the melting point. The coefficients 𝐴, 𝐵, 𝑛, and 𝑚 are empirically determined for the alloy. This benchmark emphasizes the need for robust symbolic reasoning in the absence of prior symbolic templates, bridging theory-driven modeling with experimental data interpretation.

F.4 LSR-Synth–Chemistry The LSR-Synth–Chemistry dataset, introduced as part of the LLM-SRBench benchmark by Shojaee et al., is designed to evaluate symbolic regression models on chemically motivated yet synthetically constructed reaction kinetics. It comprises 36 differential equation discovery tasks, each modeling the time evolution of a reactant concentration 𝐴(𝑡) using novel, data-driven expressions. Each equation features a distinct combination of symbolic components—ranging from classic kinetic motifs (e.g., first- and second-order decay terms like −𝑘𝐴(𝑡), −𝑘𝐴(𝑡)2 ) to synthetic nonlinearities. These include exponential decays such as exp(−𝑘𝑠 𝑡), logarithmic forms like 𝐴(𝑡)2

log(𝐴(𝑡) + 1), square roots, and oscillatory terms such as sin(𝜔𝐴(𝑡)) and cos(⋅). In addition, rational expressions like 1+𝛽𝐴(𝑡)4 introduce

challenges in terms of singularity avoidance and numerical stability. The dataset emphasizes symbolic diversity and parametric variability, making it a rigorous testbed for evaluating model generalization, interpretability, and robustness across structurally distinct formulations. Each task was carefully validated for analytical solvability, numerical stability, and scientific plausibility via expert review.

G Additional Ablation Studies G.1 Ablation Study on Remaining Benchmarks We further validated the key components via ablations on the remaining three benchmarks, with the results shown in Figure 8 to Figure 10. The experimental results are highly consistent with the findings presented in the main text. First, for the single-objective optimization variant (w/o MultiObj), a significant degradation in out-of-domain (OOD) generalization ability was observed across all benchmarks when compared to the full model. Although the in-domain (ID) fitting error of this variant is comparable to that of the full MOT-SR on some tasks, its performance on OOD data systematically confirms the tendency of single-objective optimization to overfit the ID data. This tendency consequently hinders the discovery of universally applicable scientific laws. Second, the ablation experiments on the internal components of the Meta Strategy Generator show that removing either the data analysis module (w/o Data) or the structure analysis module (w/o Struct) leads to a distinct decline in model performance. More critically, the degree of performance degradation caused by removing both modules simultaneously (w/o Strategy) exceeds the impact of removing either single module. This phenomenon clearly reveals the functional complementarity and synergistic effect between the data-driven prior knowledge and the structure-driven feedback mechanism. Both are indispensable for achieving MOT-SR’s final performance.

G.2

Ablation Study on Data Analysis Tools

To evaluate the contributions of the data analysis tools in MOT-SR, we conducted a series of ablation experiments. Our framework integrates six categories of analysis tools, and we set up six control groups (w/o tool-1 to 6), each corresponding to the removal of one category of tools. As shown in the Figure 11, removing any single category of analysis tools in this test leads to a decline in the model’s final performance. It is noteworthy that although there is only one tool in the sixth category, its removal still resulted in a final error nearly an order of magnitude higher than that of the full MOT-SR. This result strongly demonstrates that external analysis tools can effectively enhance a

MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models

large language model’s insight into the underlying structure and variable relationships within data, enabling it to capture deeper trends behind the data points and thereby generate more accurate symbolic equations.

Oscillation 2

10 3

ID 7.72e-05

10 4

OOD 9.01e-05

NMSE (log scale)

10 5

2.32e-06

1.85e-06

1.63e-06

10 6

3.36e-07

2.16e-07 10 7

8.26e-09

8.10e-09

10 8 10 9

1.70e-10 10 10

Full MOT-SR

w/o Multi-Objective

w/o Struct

w/o Data

w/o Strategy

Figure 8: Ablation study on Oscillation 2.

E. coli Growth ID

1.04e+00

100

NMSE (log scale)

6.37e-01

10 1

8.48e-01

7.24e-01

OOD 9.23e-01

7.21e-02 3.06e-02 2.21e-02

2.08e-02 1.45e-02 10 2

Full MOT-SR

w/o Multi-Objective

w/o Struct

w/o Data

Figure 9: Ablation study on E. coli growth.

w/o Strategy

Wang et al.

Stress-Strain

10 1

ID

NMSE (log scale)

6.04e-02

1.98e-02 2.02e-02

1.54e-02

1.50e-02

2.01e-02

1.93e-02

OOD

2.33e-02

1.21e-02

10 2

2.04e-03

10 3

Full MOT-SR

w/o Multi-Objective

w/o Struct

w/o Data

w/o Strategy

Figure 10: Ablation study on Stress-strain.

3.13e-09 10 9

5.97e-11

NMSE (log scale)

10 10

5.55e-11

10 11

1.48e-12

1.38e-12

10 12

7.41e-14

10 13 10 14 10 15

1.27e-15

Full MOT-SR

w/o tool-1

w/o tool-2

w/o tool-3

w/o tool-4

w/o tool-5

w/o tool-6

Figure 11: Ablation study of the Data Analysis Tools on Oscillator 1.

H

HV and IGD Definition

F.1 HV. The hypervolume (HV) measures the volume in the objective space that is dominated by the approximate Pareto front and bounded by a reference point. It reflects the convergence and diversity of the solutions. The HV is formally defined as follows: HV(𝑃, 𝑟 ∗ ) = VOL ( ⋃ [𝑣1 , 𝑟1∗ ] × [𝑣2 , 𝑟2∗ ] × ⋯ × [𝑣𝑚 , 𝑟𝑚∗ ]) , v∈𝑃

where 𝑃 denotes the approximate Pareto front obtained by the symbolic regression algorithm, v = (𝑣1 , … , 𝑣𝑚 )⊤ represents an objective vector, VOL(⋅) indicates the Lebesgue measure, and 𝑟 ∗ = (𝑟1∗ , … , 𝑟𝑚∗ )⊤ is the reference point. To eliminate the impact of varying scales and units among objectives in HV calculation, we normalize each objective to the interval [0, 1]. For an objective value 𝑓𝑖 (𝑥), the normalized value is:

𝑓𝑖′ (𝑥) = where

𝑓𝑖 (𝑥) − 𝑧𝑖ideal 𝑧𝑖nadir − 𝑧𝑖ideal

𝑧𝑖ideal = min 𝑣𝑖 , v∈𝑃

,

𝑧𝑖nadir = max 𝑣𝑖 v∈𝑃

MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models

After normalization, each objective is scaled to the range [0, 1]. Based on this, the reference point is set as 𝑟 ∗ = (1.05, … , 1.05)⊤ . F.2 IGD. The Inverted Generational Distance (IGD) measures both the convergence and diversity of the predicted front by computing the distance from points on the true Pareto front to their nearest counterparts in the predicted set. The formula is as follows: IGD(𝑃, 𝑃 ∗ ) =

1 ∑ min 𝑑(p, p∗ ), |𝑃 ∗ | p∗ ∈𝑃 ∗ p∈𝑃

where 𝑃 is the predicted Pareto front, 𝑃 ∗ is a set of reference points sampled from the true Pareto front, and 𝑑(p, p∗ ) represents the Euclidean distance between vectors p and p∗ in objective space. The computation of the Inverted Generational Distance (IGD) metric conventionally depends on the availability of a known true Pareto front, which is not always accessible in practical scenarios. In the context of symbolic regression evaluation, however, the IGD metric benefits from an inherent advantage, as the true Pareto front can be explicitly defined by the ground truth equation 𝑓 true . This allows for an exact assessment of the quality of generated equations relative to the ground truth, without reliance on an approximate Pareto front. Accordingly, the IGD formulation in this study is simplified as follows: IGD(𝑃, {𝑓 true }) = min 𝑑(𝑓 , 𝑓 true ) 𝑓 ∈𝑃

Afterwards, all objectives are normalized to the [0, 1] range to ensure consistency with the true objective values.

I Data Analysis Tools in MOT-SR I.1 Linear Correlation Tools Linear correlation tools provide foundational insights into how variables are related through linear or approximately linear patterns. These tools are instrumental for identifying primary dependencies, detecting dominant axes of variation, and establishing structural priors for symbolic equation modeling. MOT-SR incorporates four linear-correlation-based modules:

• Pearson Correlation Coefficient: It quantifies the degree of linear association between two variables 𝑥 and 𝑦 using the classical Pearson correlation coefficient 𝑟 . The input consists of two real-valued vectors, the output includes the correlation score and the associated 𝑝 -value indicating statistical significance. A high absolute value of 𝑟 (close to 1) suggests a strong linear trend, while values near 0 indicate little to no linear dependency. This measure is especially suitable for Gaussian-like distributions and helps identify candidate variables for symbolic terms with additive or multiplicative linear effects. • Simple Linear Regression: This module fits a univariate linear model 𝑦 = 𝑎 + 𝑏𝑥 using least squares estimation and returns the regression coefficients, the 𝑅2 value (explained variance), and associated statistics such as 𝑝 -values and standard errors. The inputs are one-dimensional arrays 𝑥 and 𝑦 , and the output reflects how well a straight line fits the data. The 𝑅2 value is particularly informative, indicating how much of the variance in 𝑦 can be explained by 𝑥 . It provides a predictive perspective on the relationship beyond correlation, supporting model selection based on explanatory power. • Residual Variance: This tool computes the variance of residuals 𝑦 − 𝑦̂ after applying linear regression. The lower the residual variance, the better the linear model captures the relationship between 𝑥 and 𝑦 . It takes the same inputs as Simple Linear Regression and internally reuses the linear regression output. This metric emphasizes prediction error dispersion, offering a direct measure of the noise or unmodeled nonlinear structure in the data. It is particularly valuable for pruning noisy or unstable variables from candidate model terms. • PCA Explained Variance: Principal Component Analysis (PCA) is applied to the joint space of 𝑥 and 𝑦 , and the proportion of variance explained by the first principal component is reported. The input is a two-dimensional matrix formed by stacking 𝑥 and 𝑦 , the output is a scalar value indicating the strength of the shared linear structure. This approach is robust to scaling and rotation and can capture the dominant direction of variation when the relationship is more general than univariate regression. It aids in identifying latent linear couplings and guides the selection of structurally informative variables.

I.2

Nonlinear Dependency Tools

Nonlinear dependency tools are employed to capture complex, non-monotonic interactions between variables, which are essential for accurately modeling nonlinear systems. Unlike linear tools, these methods can reveal relationships that are not adequately described by simple linear transformations. MOT-SR incorporates three nonlinear-dependency-based modules:

• Spearman Rank Correlation: This tool measures the strength and direction of monotonic relationships between two variables using their rank values. It computes Spearman’s 𝜌 , a non-parametric counterpart to Pearson’s 𝑟 , by evaluating how well the relationship between 𝑥 and 𝑦 can be described by a monotonic function. The input consists of two numerical vectors, the output includes the correlation

Wang et al.

coefficient and its 𝑝 -value. This method is particularly effective when data exhibit nonlinear but monotonic trends, such as saturating growth or sigmoid-like behavior. • Mutual Information: This tool quantifies the total amount of information shared between two variables, regardless of the specific functional form of their relationship. It estimates mutual information by discretizing the input variables 𝑥 and 𝑦 into bins and calculating the joint entropy. The output is a scalar metric representing dependency strength. Unlike correlation coefficients, mutual information can capture both linear and highly nonlinear dependencies, making it suitable for detecting complex statistical associations in symbolic modeling. • Mutual Information Regression Score: This variant employs the scikit-learn implementation of mutual information to assess the relevance of 𝑥 in predicting 𝑦 using a regression-based formulation. The method internally estimates how much knowing 𝑥 reduces uncertainty about 𝑦 . It accepts one-dimensional arrays as input and outputs a scalar score. This score is robust to arbitrary nonlinearities and discontinuities, making it valuable for feature selection in nonlinear symbolic regression tasks.

I.3

Time-Frequency Analysis Tools

Time-frequency analysis tools are essential for detecting periodicities, oscillations, and transient dynamics that are characteristic of nonlinear and multi-scale systems. These tools enable symbolic regression models to incorporate periodic or time-varying terms where appropriate. MOT-SR employs two such modules:

• Fast Fourier Transform Frequency Difference: This module applies the discrete Fourier transform to the input signals 𝑥 and 𝑦 to extract their respective power spectra. It identifies the dominant frequency component for each variable and returns the absolute difference between them. The input consists of two one-dimensional arrays representing time-series data, and the output includes the dominant frequencies of both variables as well as their difference. This tool captures global periodic patterns and is useful for detecting whether the signals share synchronized or harmonically related structures. A small frequency difference suggests potential functional alignment via sinusoidal or oscillatory terms. • Wavelet Energy Correlation: This module leverages discrete wavelet decomposition to capture localized energy features of 𝑥 and 𝑦 at multiple temporal scales. It decomposes both signals into several levels using a specified wavelet basis (e.g., Daubechies 4), computes the energy at each level, and measures the Pearson correlation between the resulting energy vectors. The input includes two one-dimensional arrays and optional parameters for wavelet type and decomposition level. The output is a correlation coefficient representing how similarly the energy of the two signals is distributed across scales. Unlike Fourier-based methods, this tool excels at capturing transient and nonstationary dependencies, offering robust structural priors for equations involving local periodicity or bursts.

I.4

Causal Inference Tools

Causal inference tools aim to uncover directional relationships between variables, particularly whether the historical behavior of one variable contributes to or influences another. These tools help MOT-SR to build models that not only fit the data well but also respect temporal or structural causality, thereby improving interpretability and generalization. Two modules are implemented:

• Granger Causality: This method assesses whether past values of a variable 𝑥 improve the prediction of a variable 𝑦 in a multivariate time series setting. The inputs are two equal-length time series arrays, and the test evaluates multiple lagged regression models to determine if 𝑥 Granger-causes 𝑦 . The output includes 𝑝 -values for different lag orders, and the smallest 𝑝 -value is used as the primary metric. A statistically significant result implies that the past of 𝑥 contains information predictive of 𝑦 , suggesting a directional dependency. This method is especially valuable for systems where delayed effects are prominent, such as in control dynamics or feedback loops. • Convergent Cross Mapping: CCM is a nonlinear causal discovery technique grounded in dynamical systems theory. It tests whether the state of variable 𝑥 can be reconstructed from the historical trajectory of 𝑦 , indicating that 𝑥 leaves an imprint on 𝑦 . The input consists of two time series and parameters for embedding dimension and library sampling. CCM constructs a manifold from the delay embedding of 𝑦 , uses it to cross-predict 𝑥 , and evaluates the reconstruction accuracy (typically using a cross-map skill 𝜌 ). A high 𝜌 suggests that 𝑥 causally influences 𝑦 . Unlike Granger causality, CCM does not rely on linearity or temporal precedence, making it ideal for identifying nonlinear, feedback-driven causality in complex systems.

I.5

Dynamic Complexity Tools

Dynamic complexity tools capture nuanced structural and temporal irregularities in time-series data, particularly useful for distinguishing chaotic, nonlinear, or asynchronous behaviors. MOT-SR leverages three such tools:

• Lyapunov Exponent Difference: This tool estimates the maximal Lyapunov exponents (MLEs) of two time-series variables 𝑥 and 𝑦 , then computes their absolute difference. Lyapunov exponents characterize how sensitive a system is to initial conditions. A high MLE suggests chaotic behavior, whereas near-zero or negative values indicate stable dynamics. The input consists of two real-valued time series. The tool automatically determines the embedding dimension and time lag for phase space reconstruction. The resulting metric reflects how similarly (or differently) 𝑥 and 𝑦 behave in terms of dynamical predictability.

MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models

• Correlation Dimension Difference: This tool quantifies and compares the fractal (correlation) dimensions of 𝑥 and 𝑦 . The correlation dimension serves as a complexity measure, indicating how densely a system’s trajectory fills its phase space. The input includes two realvalued sequences and a fixed embedding dimension (usually 2). The metric is the absolute difference between their estimated correlation dimensions. This captures differences in dynamic complexity and is useful for identifying structural mismatches in multivariate nonlinear processes. • Dynamic Time Warping Distance: DTW measures the alignment cost between 𝑥 and 𝑦 by allowing local nonlinear stretching or compression in time. Unlike Euclidean distance, DTW is robust to phase shifts and unequal pacing in temporal evolution. The input is a pair of univariate time series, and the output is a scalar DTW distance. This tool is particularly effective in identifying temporal patterns that share similar shapes but occur at different rates or phases.

I.6

Distribution Consistency Tools

This category includes statistical tools that assess whether two variables share similar probability distributions. Such tools are especially useful when validating whether a generated or transformed signal preserves the underlying distributional structure of the original data.

• Kolmogorov–Smirnov Test: The Kolmogorov–Smirnov (KS) test is a non-parametric method that quantifies the maximum distance between the empirical cumulative distribution functions (ECDFs) of two datasets 𝑥 and 𝑦 . The input to this tool is a pair of real-valued vectors, and it returns the KS statistic (a scalar metric) along with the corresponding 𝑝 -value indicating statistical significance. The KS statistic captures how different the two distributions are—larger values indicate more significant deviations. This test is sensitive to both location and shape differences in the distributions and is thus useful for evaluating whether symbolic transformations (e.g., derived equations) maintain statistical fidelity to the data source.

J Results Presentation J.1 Equation Recovery and Pareto Front Evolution Fig 12a to Fig 12d are sampled from four test cases in which MOT-SR successfully identified the ground-truth equations on the LSR-Synth– Chemistry dataset. Each figure illustrates the evolution of NMSE of the Pareto front equation set across iterations. For visualization, we present the structural skeletons of the equations discovered by MOT-SR. Highlighted segments in red indicate terms that match exactly with the corresponding ground-truth equation. It is noteworthy that the size of the Pareto front typically expands rapidly during the early stages, forming a diverse equation population. As iterations proceed, this population undergoes continuous updates. Once sufficient structural experience has accumulated, the Pareto front enters a phase of rapid convergence and eventually stabilizes to the correct target equation. Below we present comparisons between MOT-SR and LLM-SR on the same tasks: CRK14: Ground truth: −𝑘𝐴(𝑡) + 𝑘𝑝 sin(𝜔𝐴(𝑡)) MOT-SR:

𝜃0 𝐴 + 𝜃1 sin(𝜃2 𝐴) LLM-SR:

𝜃0 𝐴 (1 + 1+ CRK19: Ground truth:

(𝜃1 𝐴)2 ) 𝜃2

1+ 𝐴𝜃 1+

2

5 (𝜃1 𝐴)2 𝜃2

+

𝜃4 𝐴 + 𝜃3 𝐴 𝜃5 + 𝐴

−𝑘𝐴(𝑡)2 + 𝑘𝑝 sin(𝜔𝐴(𝑡))

MOT-SR:

𝜃0 sin(𝜃1 𝐴) + 𝜃2 𝐴2 LLM-SR:

−𝜃0 𝐴 − 𝜃1 𝐴2 + 𝜃2 𝑡 + 𝜃3 𝑡 2 + 𝜃4 𝐴3 CRK21: Ground truth:

−𝑘𝐴(𝑡)𝑒 −𝑘𝑠 𝑡 + 𝑘𝑝 sin(𝜔𝐴(𝑡))

MOT-SR:

𝜃0 𝐴𝑒 −𝜃1 𝑡 + 𝜃2 sin(𝜃3 𝐴)

Wang et al.

LLM-SR:

𝜃0 𝑡 + 𝜃1 𝐴𝑒 −𝜃2 𝑡 +

𝜃5 𝐴𝜃6 𝜃3 𝐴 + 𝜃4 𝑡 + 1 1 + 𝜃7 𝐴𝜃6

CRK36: Ground truth:

−𝑘𝐴(𝑡) + 𝑘𝑞 𝐴(𝑡) log(𝛾 𝑡 + 1) MOT-SR:

𝜃0 𝐴 + 𝜃1 𝐴 log(𝜃2 𝑡 + 1) LLM-SR:

𝜃0 𝐴 + 𝜃1 𝑒 −𝜃2 𝑡 𝐴 +

𝜃5 𝐴 𝜃3 𝐴 + + 𝜃7 𝑡 + 𝜃8 𝐴 + 𝜃9 𝑡𝐴 𝜃4 + 𝐴 1 + 𝜃6 𝐴

These results clearly demonstrate that MOT-SR is capable of precisely capturing the exact components of the ground-truth equations, yielding structurally closer expressions than LLM-SR. We attribute this improvement to MOT-SR’s enhanced exploratory capability derived from multi-objective evaluation, as well as its strategy-guided generation process. NMSE Spread over Iterations

NMSE Spread over Iterations

10 4

Min-Max Range Min NMSE

10 2

eq1.

Nmse Test All (log scale)

Nmse Test All (log scale)

10 6

eq2.

10 8

Min-Max Range Min NMSE

100

eq3.

10 10

eq1.

10 4

eq2.

10 6 10 8 10 10

10 12 10 12 0

50

100

Iterations

150

200

250

0

(a) Pareto Front Performance in CRK14.

50

100

150

Iterations

200

(b) Pareto Front Performance in CRK19. NMSE Spread over Iterations

NMSE Spread over Iterations Min-Max Range Min NMSE

10 1

Min-Max Range Min NMSE

10 3

eq1.

10 5

Nmse Test All (log scale)

Nmse Test All (log scale)

eq2.

10 3 10 5 10 7 10 9

10 7 eq.1

10 9

eq.2

10 11

eq.3 eq.4

10 13

10 11 0

50

100

Iterations

150

(c) Pareto Front Performance in CRK21.

200

0

50

100

Iterations

150

200

250

(d) Pareto Front Performance in CRK36.

Figure 12: Pareto Front Performance results across four benchmarks.

J.2

Comparison of convergence dynamics.

Figure 13 presents the wall-clock time evolution of NMSE (log scale) for MOT-SR and the baseline LLM-SR across four representative benchmarks. The results demonstrate that MOT-SR consistently achieves substantially faster error reduction, while LLM-SR converges more slowly and plateaus at significantly higher error levels. This advantage is particularly pronounced in the Oscillation tasks, where

MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models

MOT-SR (Llama-3.1)

10 4 10 7 10 10 10 13 0

4

12

16

10 3 10 5 10 7 10 9 10 11

0

4

E. coli Growth

10 1

0

4

8

Time (hours)

8

Time (hours)

12

16

12

16

Stress-Strain NMSE (log scale)

100

NMSE (log scale)

8

Time (hours)

Oscillation 2

10 1

NMSE (log scale)

NMSE (log scale)

10 1

10 2

LLM-SR (Llama-3.1)

Oscillation 1

12

16

10 1

10 2

0

4

8

Time (hours)

Figure 13: Wall-clock time evolution of NMSE showing faster convergence of MOT-SR over LLM-SR.

MOT-SR rapidly identifies accurate dynamical structures, but is also evident in more challenging domains such as E. coli growth and stress– strain modeling, where the baseline stagnates prematurely. These findings highlight that the integration of tool-guided analysis and multiobjective optimization enables MOT-SR to explore the hypothesis space more efficiently, yielding superior equations under comparable computational budgets.

Wang et al.

Evolution of Discovered Equations

100

10 2

c1x + c0v + c2ve c3x

10 4

Normalized MSE

J.3

c1sin(c2x) + c4cos(c 5x) + c6x + c0v + c3ve c7x + c8xv

10 6

10 8

10 10

c0x + c1v c2x2 + c3v2 + c4x3 + c5v3 c6xv + c7sin(c8x) + c9v

10 12

c0x + c1v + c2x3 + c3v3 c4xv + c5sin(c6x)

10 14 0

250

500

750

1000

Iteration

1250

1500

1750

Figure 14: Iterative evolution of equation structures under MOT-SR (Oscillation1).

2000

MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models

100

10 2

c0t 2A + c1tA + c2A c3A 2 + c4t + c5t 2 + c6

Normalized MSE

10 4

10 6

c0t 2 + c1A + c2ec3t + c4sin(c5t)A + c6tA A c0 + c1e c2t + c3t c4 + c5t 2

10 8

(

)

10 10

10 12

Ac0 + Ac1log(c2t + 1)

10 14 0

250

500

Iteration

750

1000

Figure 15: Iterative evolution of equation structures under MOT-SR (CRK36). Evolution of Equation Structures for Oscillation 1. Figure 14 presents the iterative evolution of equation structures for the problem oscillation 1. In the early iterations, the generated expressions contained a wide range of nonlinear and polynomial terms, many of which represented exploratory hypotheses rather than meaningful components. Guided by the TOOL-AUGMENTED VARIABLE ANALYSIS module, the framework progressively refined these structures. When the analytical tools within this module (e.g., wavelet corr, FFT, and PCA mapping) revealed periodic features, sinusoidal terms such as sine and cosine functions were introduced into the candidate equations, thereby complementing the polynomial basis and improving alignment with the observed data. During the iteration, structurally inconsistent or statistically unsupported terms were systematically removed, while the retained nonlinear and oscillatory components captured the essential dynamics of the system. The process ultimately converged toward compact formulations that balance accuracy, generalization, and complexity. This refinement demonstrates how the iterative mechanism, augmented by TOOL-AUGMENTED VARIABLE ANALYSIS, effectively distills the equation space and produces models consistent with both empirical behavior and theoretical plausibility. Evolution of Equation Structures for CRK36. Figure 15 presents the iterative evolution of equation structures for the problem CRK36. At the beginning of the process, the candidate functions contained a variety of nonlinear combinations involving both the temporal variable and A. Under the guidance of two correlation coefficients, the framework progressively refined these structures by eliminating inconsistent nonlinearities and retaining only the dominant dependencies. As the iterations advanced, the generated functions converged toward a fully linear dependence on A.

K

Gaussian Noise Evaluation

To evaluate the robustness of MOT-SR in noisy environments, we conduct experiments on the Oscillator1 benchmark by injecting Gaussian noise with standard deviations 𝜎 = {0.001, 0.002} into the training data. Model performance is assessed under both in-domain (ID) and out-of-domain (OOD) test conditions. Unlike conventional symbolic regression approaches that rely solely on fitting accuracy, MOT-SR integrates structural priors—derived from analytical tools—and performs multi-objective evaluation that jointly considers accuracy, complexity, and generalization. This allows it to maintain high-fidelity symbolic recovery even when data are corrupted. Table 4 shows that MOT-SR consistently outperforms LLMSR across all metrics, achieving significantly lower NMSE and higher Accavg-0.01 scores under both noise settings. These gains stem from MOT-SR’s integration of tool-informed structural priors via the Meta Strategy Generator and its multi-objective Equation Generator, which jointly enable the model to identify robust, generalizable expressions even under noisy supervision.

Wang et al.

Table 4: Comparison of MOT-SR and LLMSR on Oscillator1 under Gaussian noise (𝜎 = {0.001, 0.002}). Metrics include NMSE and accuracy (Accavg-0.01 ) on ID and OOD test sets. All models use LLaMA-3.1 as backbone.

𝜎 = 0.001 Model

LLMSR MOT-SR

L

ID

𝜎 = 0.002 OOD

ID

OOD

NMSE↓

Accavg-0.01 ↑

NMSE↓

Accavg-0.01 ↑

NMSE↓

Accavg-0.01 ↑

NMSE↓

Accavg-0.01 ↑

1.59𝑒 -5 5.43e-11

89.20% 99.99%

0.0021 5.49e-7

10.46% 98.79%

2.68𝑒 -5 3.98e-8

82.08% 99.74%

0.0025 8.57e-7

10.50% 97.88%

Assessment of Python AST Length as a Complexity Proxy

Classical symbolic regression methods operate directly on mathematical expression trees, where structural complexity is commonly measured by counting expression-tree nodes. MOT-SR allows the Equation Generator to produce executable Python functions without a fixed operator dictionary or function whitelist. This program-level representation does not directly provide a canonical mathematical expression tree. Known syntactic bias. In our implementation, each candidate function is parsed with ast.parse, and its Python AST length is defined as the total number of nodes returned by ast.walk. This metric reflects the syntactic structure of Python code and introduces a known bias. For example, a composite mathematical call such as np.sin(x) can produce more Python AST nodes than a basic arithmetic operation, even when the two have comparable mathematical expression-tree complexity. Python AST length is thus treated as a computational proxy for mathematical expression complexity. Computational trade-off. A stricter alternative would convert each generated Python function into a normalized mathematical expression, construct the corresponding expression tree, and count its nodes. This procedure aligns more closely with conventional complexity measures in symbolic regression, but adds a conversion step for every candidate generated during the search. MOT-SR evaluates thousands of candidates, making repeated conversion a substantial computational burden. Python AST length can be computed deterministically and directly from the generated source code with negligible overhead. Empirical ranking consistency. Pareto-based selection depends primarily on the relative complexity ordering of candidate equations. We therefore assess whether Python AST length preserves the ordering induced by mathematical expression-tree complexity. Let

𝐶AST (𝑓 ) denote the Python AST node count of a candidate equation 𝑓 , and let

𝐶Expr (𝑓 ) denote the node count of its normalized mathematical expression tree. For a set 𝒬 of candidate pairs, we define the pairwise rankingdisagreement ratio as 1 ∑ 𝕀 [sign (𝐶AST (𝑓𝑖 ) − 𝐶AST (𝑓𝑗 )) ≠ sign (𝐶Expr (𝑓𝑖 ) − 𝐶Expr (𝑓𝑗 ))] , 𝑟dis = (23) |𝒬| (𝑓 ,𝑓 )∈𝒬 𝑖 𝑗

where comparisons tied under either metric are excluded. We retrospectively computed both complexity measures for valid candidate equations collected from the search logs and evaluated nearly 400,000 candidate pairs. The resulting pairwise disagreement ratio was approximately 5%, indicating that the two metrics produced the same relative ordering for approximately 95% of the evaluated pairs. Within the candidate distributions observed in our experiments, Python AST length therefore preserves the expression-tree ordering in most comparisons. Because the complexity objective is used to determine relative dominance and construct the Pareto front, this level of agreement indicates that the proxy has a limited effect on the practical selection order. Composite-function robustness. We further examine whether the additional syntactic cost of composite functions systematically suppresses candidates containing trigonometric or other non-polynomial terms. The analysis covers tasks whose target structures depend strongly on such functions, including the Oscillation 1 and Oscillation 2 benchmarks and trigonometry-heavy tasks in LSR-Synth–Chemistry. MOT-SR discovers compact equations containing terms such as sin(⋅) and cos(⋅), while achieving strong numerical and symbolic recovery performance on these tasks. These results indicate that the additional AST cost assigned to composite calls does not prevent the corresponding structures from entering or remaining on the Pareto front. Python AST length is a language-dependent and syntax-dependent measure, rather than an exact representation of mathematical complexity. We use it as a low-overhead engineering proxy that supports flexible Python-based equation generation. Within the candidate distributions examined in this work, the retrospective ranking analysis and the results on composite-function-heavy tasks show that this proxy provides a sufficiently consistent complexity ordering for Pareto-based search.

MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models

M

Prompt

J.1 Prompt Design. The prompts used to construct the meta strategy are shown in Fig. 16 and Fig. 17. Specifically, Fig. 16 presents the meta prompt used during initialization, while Fig. 17 corresponds to the prompts designed for iterative optimization. Together, these constitute the prompt inputs to the Meta Strategy Generator. J.2 Examples. The prompts used by the Equation Generator for the Oscillators 1 task are illustrated in Fig. 18 to Fig. 21. The prompts received by the Meta Strategy Generator for the same task are shown in Fig. 22 to Fig. 27. Fig. 28 provides examples of data analysis tool usage, drawn from tool invocations in the CRK36 task. Notably, through the integration of external tools, the language model successfully identified a strong linear relationship between the variable 𝑥0 (that is, 𝐴(𝑡)) and the target variable 𝑦 . As demonstrated in the final output, MOT-SR successfully recovered the ground-truth governing equation for the CRK36 problem.

Wang et al.

Meta Prompt SYSTEM ROLE: You are a **scientist-assistant LLM** supervising a symbolic-regression pipeline. You must analyze existing equations, select appropriate analysis tools, and produce guidance for a generator LLM. ## Allowed analysis tools {allowed_list} Choose **up to 3** tools—exclusively from the list above—based **solely on the diagnostic needs revealed by the residual data you receive**.If a desired tool is not listed, omit it. ## RESPONSE FORMAT (STRICT) Return exactly one JSON object (no markdown fences). It must contain: 1. "structure_insight_prompt_for_generator_LLM" – either a string, or an object with `guidance.prompt_for_generator_LLM` (string). 2. "analysis_tools" – a JSON array of 1–5 tool names from the allowed list. Escape newlines inside strings as `\n`. After the closing brace `}`, output nothing else.

Figure 16: Meta Prompt.

Refinement Prompt ## Context {context} ## Mathematical Pattern Analysis Task Analyze the following equations to identify **common mathematical patterns and shared structural elements**:

.join(blocks) ## Available Analysis Tools {tool_list} ## Pattern Analysis Instructions 1. **Examine Mathematical Structures**: Look across all equations for: - Recurring mathematical functions (trigonometric, exponential, polynomial) - Similar variable interaction patterns (x*y, x+y, x^n, etc.) - Consistent coefficient magnitudes or sign patterns - Common denominators or fraction structures - Similar nesting or grouping of terms 2. **Identify Shared Elements**: Find mathematical components that appear in multiple high-performing equations 3. **Interpret Mathematical Significance**: Explain what these patterns might indicate about: - Underlying physical relationships - Mathematical constraints or principles - Optimal solution characteristics

## Pattern Analysis Instructions 1. **Examine Mathematical Structures**: Look across all equations for: - Recurring mathematical functions (trigonometric, exponential, polynomial) - Similar variable interaction patterns (x*y, x+y, x^n, etc.) - Consistent coefficient magnitudes or sign patterns - Common denominators or fraction structures MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models - Similar nesting or grouping of terms 2. **Identify Shared Elements**: Find mathematical components that appear in multiple high-performing equations 3. **Interpret Mathematical Significance**: Explain what these patterns might indicate about: - Underlying physical relationships - Mathematical constraints or principles - Optimal solution characteristics 4. **Generate Guidance**: Create a detailed prompt that will guide future equation generation based on successful patterns 5. **Select Analysis Tools**: Choose up to 5 data analysis tools that could help validate the patterns you identified ## RESPONSE FORMAT Return ONLY JSON with keys: structure_insight_prompt_for_generator_LLM, analysis_tools. - structure_insight_prompt_for_generator_LLM: Detailed guidance based on your pattern analysis for generating better equations - analysis_tools: Array of up to 5 tool names from the available list Focus on mathematical commonalities and their significance in equation generation.

Figure 17: Prompts Designed for Iterative Optimizationt.

Oscillator1 You are a helpful assistant tasked with discovering mathematical function structures for scientific systems. Complete the 'equation' function below, considering the physical meaning and relationships of inputs.

# Scientific Analysis Based on Pareto Front Selection ## Core Generation Guidance ## Core Generation Guidance Generate equations that incorporate the identified mathematical patterns and shared structural elements. Emphasize leveraging variable interactions, consistent coefficient magnitudes, and common denominators. Prioritize equations that respect physical constraints and optimization characteristics. ## Generation Instructions Based on the above analysis, focus on: 1. Incorporating the identified recurring mathematical functions 2. Using the successful variable interaction patterns 3. Maintaining coefficient relationships that work well 4. Respecting the physical constraints and relationships 5. Building upon the shared mathematical elements

Figure 18: Example Prompt in Oscillator1 - For Equation Generator.

Wang et al.

## Data Analysis Results Applied 3 analysis tools to all variable pairs: **Tools used**: spearman_corr, residual_var, corr **Dataset overview**: - 2 independent variables analyzed - 10000 total data points ## Correlation Analysis: x0_y **Strong correlation detected** (r=-0.920) Consider incorporating a strong relationship between x0_y and output in your equation. The relationship appears to be negative - as one increases, the other decreases. ## Residual Variance Analysis: x0_y **Residual variance**: 0.000419 Low residual variance suggests the linear relationship captures most variation.

## Correlation Analysis: x0_y **Strong correlation detected** (r=-0.951) Consider incorporating a strong relationship between x0_y and output in your equation. The relationship appears to be negative - as one increases, the other decreases. ## Correlation Analysis: x1_y **Weak correlation detected** (r=0.103) The relationship between x1_y and output appears weak or non-linear. ## Residual Variance Analysis: x1_y **Residual variance**: 0.004359 Low residual variance suggests the linear relationship captures most variation. ## Correlation Analysis: x1_y **Weak correlation detected** (r=-0.032) The relationship between x1_y and output appears weak or non-linear.

Figure 19: Example Prompt in Oscillator1 - For Equation Generator.

## Generation Strategy Based on the analysis of Pareto front samples and data relationships: 1. **Build upon success**: Use the selected high-quality solution as foundation 2. **Targeted improvement**: Focus on reducing residuals in identified problem areas 4. **Data-driven insights**: Incorporate findings from variable relationship analysis 5. **Variable interactions**: Consider the significant correlations and patterns identified 6. **Maintain balance**: Preserve the multi-objective balance that led to Pareto selection 7. **Consider trade-offs**: Balance accuracy, generalization, and computational efficiency

## Generation Strategy Based on the analysis of Pareto front samples and data relationships: 1. **Build upon success**: Use the selected high-quality solution as foundation MOT-SR: Multi-Objective Tool-Augmented Scientific Equation DiscoveryFocus with Largeon Language Models residuals in identified 2. **Targeted improvement**: reducing problem areas 4. **Data-driven insights**: Incorporate findings from variable relationship analysis 5. **Variable interactions**: Consider the significant correlations and patterns identified 6. **Maintain balance**: Preserve the multi-objective balance that led to Pareto selection 7. **Consider trade-offs**: Balance accuracy, generalization, and computational efficiency

""" Find the mathematical function skeleton that represents acceleration in a damped nonlinear oscillator system with driving force, given data on position, and velocity. """

Figure 20: Example Prompt in Oscillator1 - For Equation Generator.

def equation(x: np.ndarray, v: np.ndarray, params: np.ndarray) -> np.ndarray: """ Mathematical function for acceleration in a damped nonlinear oscillator Args: x: A numpy array representing observations of current position. v: A numpy array representing observations of velocity. params: Array of numeric constants or parameters to be optimized Return: A numpy array representing acceleration as the result of applying the mathematical function to the inputs. """ # Driving force considering correlation with x0_y and output f = params[0] * np.exp(-params[1] * v**2) # Driving force with damping # Incorporating strong linear relationship between x0_y and output linear_term = params[2] * x # Linear coefficient for damping

# Damping effect considering correlation with x1_y and output damp = params[3] * (1 + 0.5 * np.sin(params[4] * x)) # Damping coefficient with trigonometric term # Nonlinear oscillator equation considering correlations a = -(damp + 0.5 * linear_term) * v # Additional non-linear terms for better model fit b = params[5] * np.sin(params[6] * x + params[7] * v) # Nonlinear term with hybrid sine function c = params[8] * v**2 + f * np.sin(params[9] * x) # Nonlinear term with exponenial trig function # Regularization technique to prevent overfitting return a + b + c def equation_v1(x: np.ndarray, v: np.ndarray, params: np.ndarray) -> np.ndarray: """Improved version of `equation_v0`."""

Figure 21: Example Prompt in Oscillator1 - For Equation Generator.

Wang et al.

Oscillator1 SYSTEM ROLE: You are a **scientist-assistant LLM** supervising a symbolic-regression pipeline. You must analyze existing equations, select appropriate analysis tools, and produce guidance for a generator LLM. ## Allowed analysis tools ccm_causality, corr, corr_dim_relation, dtw_distance, fft_cross_freq, granger_causality, ks_test_diff, lin_reg, lyapunov_relation, mutual_info, mutual_info_regression_score, pca_mapping, pearson_corr, residual_var, spearman_corr, wavelet_corr Choose **up to 3** tools—exclusively from the list above—based **solely on the diagnostic needs revealed by the residual data you receive**.If a desired tool is not listed, omit it. ## RESPONSE FORMAT (STRICT) Return exactly one JSON object (no markdown fences). It must contain: 1. "structure_insight_prompt_for_generator_LLM" – either a string, or an object with `guidance.prompt_for_generator_LLM` (string). 2. "analysis_tools" – a JSON array of 1–5 tool names from the allowed list. Escape newlines inside strings as `\n`. After the closing brace `}`, output nothing else. ## Context Symbolic regression problem: oscillator1. Dataset has 2 input variables and 1 target variable. Training data shape: (10000, 2) -> (10000,) Analyzing equations selected from Pareto front with residual data. ## Mathematical Pattern Analysis Task Analyze the following equations to identify **common mathematical patterns and shared structural elements**:

Figure 22: Example Prompt in Oscillator1 - For Meta Strategy Generator.

### Equation 1 Expression : """Mathematical function for acceleration in a damped nonlinear oscillator Args: x: A numpy array representing observations of current position. v: A numpy array representing observations of velocity. params: Array of numeric constants or parameters to be optimized params[0]: Spring constant, effectively the driving force coefficient params[1]: Damping coefficient params[2]: Amplitude, the maximum value of driving force params[3]: Angular frequency, c=driving force frequency params[4]: Mass of object, m params[5]: Friction or driving force damping coefficient params[6]: Non-linear damping exponent params[7]: Non-linear driving force frequency millers params[8]: Driving force input frequency peak frequency millers params[9]: Non-linear displacement spring constant magnetic field coefficients Returns: A numpy array representing acceleration as the result of applying the mathematical function to the inputs. """ # Driving force c = params[0] # Original driving force coefficient A = params[2] # Original amplitude

params[3]: Angular frequency, c=driving force frequency params[4]: Mass of object, m params[5]: Friction or driving force damping coefficient params[6]: Non-linear damping exponent params[7]: Non-linear driving force frequency millers params[8]: Driving force input frequency peak frequency millers params[9]: Non-linear displacement spring constant magnetic MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models field coefficients Returns: A numpy array representing acceleration as the result of applying the mathematical function to the inputs. """ # Driving force c = params[0] # Original driving force coefficient A = params[2] # Original amplitude # Frictionless surface parameters (assuming frictionless) m = params[4] # mass of object # Non-linear ODE (On-Off Oscillators) # x-dependent, non-linear spring constant k = params[9] + params[1] * np.exp(-np.abs(v)) constant

# Original spring

# Incorporating driving force non-linearity # exponential riding f = c * (A + np.tanh(params[6] * (v + params[7]))) * np.sin(params[8] * x) # Non-linear ODE application return -m * (k * v + f)

Figure 23: Example Prompt in Oscillator1 - For Meta Strategy Generator.

Wang et al.

ID Score : -0.00013475 (in-distribution performance) OOD Score : -5.38816e-05 (out-of-distribution performance) Eval Time : 10.6924s (computational efficiency) Avg Accuracy : -9.43156e-05 (overall accuracy) Gen Gap : 8.08681e-05 (generalization stability) Efficiency : 0.0855254 (time-normalized performance)

**In-Distribution (ID) Dataset Analysis:** Worst Region Statistics: - Sample Count: 520 samples - Threshold: 0.000958914 - Mean Residual: -0.0008611 - Mean Abs Residual: 0.00112713 - Max Abs Residual: 0.00148016 - Std Residual: 0.000750307 Input Variables in Worst Region: x0: [-0.481813, -0.419474] x1: [-0.124159, 0.113473] Variable (y) in Worst Region: y: [0.048175, 0.0951525] **Out-of-Distribution (OOD) Dataset Analysis:** Worst Region Statistics: - Sample Count: 480 samples - Threshold: 0.000740516 - Mean Residual: -0.000603334 - Mean Abs Residual: 0.000866861 - Max Abs Residual: 0.000997573 - Std Residual: 0.000629749 Input Variables in Worst Region: x0: [-0.520795, 0.489177] x1: [-0.0900758, 0.282136] Variable (y) in Worst Region: y: [-0.0614007, 0.060305]

Target

Target

Figure 24: Example Prompt in Oscillator1 - For Meta Strategy Generator.

MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models

### Equation 2 Expression : """ Mathematical function for acceleration in a damped nonlinear oscillator. Args: x: A numpy array representing observations of current position. v: A numpy array representing observations of velocity. params: Array of numeric constants or parameters to be optimized params[0]: Spring constant (k), driving force coefficient (omega) params[1]: Damping coefficient (gamma) params[2]: Amplitude (A) params[3]: Angular frequency (omega) Returns: A numpy array representing acceleration as the result of applying the mathematical function to the inputs. """ # Parameters for driving force omega = params[0] # Angular frequency omega_driving_force = omega # Parameters for friction or driving force damping gamma = params[1] # Damping coefficient gamma_friction_driving_damping = gamma # Mass and amplitude for non-linear oscillator m = params[4] # Mass A_driving_damping = params[2] # Amplitude # Calculate velocity induced force force = omega_driving_force * np.sin(params[3] * np.tanh(params[5] * v)) # Acceleration equation # Consider strong correlations between position and velocity acceleration = -m * (force + gamma_friction_driving_damping * v) return acceleration

Figure 25: Example Prompt in Oscillator1 - For Meta Strategy Generator.

ID Score : -1.16205 (in-distribution performance) OOD Score : -1.00194 (out-of-distribution performance) Eval Time : 0.698804s (computational efficiency) Avg Accuracy : -1.082 (overall accuracy) Gen Gap : 0.160103 (generalization stability) Efficiency : 0.58865 (time-normalized performance)

**In-Distribution (ID) Dataset Analysis:** Worst Region Statistics: - Sample Count: 520 samples - Threshold: 0.101862 - Mean Residual: 0.0907549 - Mean Abs Residual: 0.109957 - Max Abs Residual: 0.113557 - Std Residual: 0.0622102 Input Variables in Worst Region: x0: [-0.512046, 0.444764] x1: [0.11177, 0.212819] (y) in Worst Region: y: [-0.104672, 0.110646] **Out-of-Distribution (OOD) Dataset Analysis:** Worst Region Statistics: - Sample Count: 480 samples - Threshold: 0.0978116 - Mean Residual: 0.0446015 - Mean Abs Residual: 0.103459

Target Variable

**In-Distribution (ID) Dataset Analysis:** Worst Region Statistics: - Sample Count: 520 samples - Threshold: 0.101862 - Mean Residual: 0.0907549 - Mean Abs Residual: 0.109957 - Max Abs Residual: 0.113557 - Std Residual: 0.0622102 Input Variables in Worst Region: x0: [-0.512046, 0.444764] x1: [0.11177, 0.212819] Target Variable (y) in Worst Region: y: [-0.104672, 0.110646] **Out-of-Distribution (OOD) Dataset Analysis:** Worst Region Statistics: - Sample Count: 480 samples - Threshold: 0.0978116 - Mean Residual: 0.0446015 - Mean Abs Residual: 0.103459 - Max Abs Residual: 0.111418 - Std Residual: 0.0934521 Input Variables in Worst Region: x0: [-0.562765, 0.495206] x1: [0.0432888, 0.221801] Target Variable (y) in Worst Region: y: [-0.103036, 0.10915]

Figure 26: Example Prompt in Oscillator1 - For Meta Strategy Generator.

## Available Analysis Tools ccm_causality, corr, corr_dim_relation, dtw_distance, fft_cross_freq, granger_causality, ks_test_diff, lin_reg, lyapunov_relation, mutual_info, mutual_info_regression_score, pca_mapping, pearson_corr, residual_var, spearman_corr, wavelet_corr ## Pattern Analysis Instructions 1. **Examine Mathematical Structures**: Look across all equations for: - Recurring mathematical functions (trigonometric, exponential, polynomial) - Similar variable interaction patterns (x*y, x+y, x^n, etc.) - Consistent coefficient magnitudes or sign patterns - Common denominators or fraction structures - Similar nesting or grouping of terms 2. **Identify Shared Elements**: Find mathematical components that appear in multiple high-performing equations 3. **Interpret Mathematical Significance**: Explain what these patterns might indicate about: - Underlying physical relationships - Mathematical constraints or principles - Optimal solution characteristics 4. **Generate Guidance**: Create a detailed prompt that will guide future equation generation based on successful patterns 5. **Select Analysis Tools**: Choose up to 5 data analysis tools that could help validate the patterns you identified ## RESPONSE FORMAT Return ONLY JSON with keys: structure_insight_prompt_for_generator_LLM, analysis_tools. - structure_insight_prompt_for_generator_LLM: Detailed guidance based on your pattern analysis for generating better equations - analysis_tools: Array of up to 5 tool names from the available list Focus on mathematical commonalities and their significance in equation generation.

Figure 27: Example Prompt in Oscillator1 - For Meta Strategy Generator.

Wang et al.

MOT-SR: Multi-Objective Tool-Augmented Scientific Equation Discovery with Large Language Models

Partial Data Analysis Tools Usage in CRK36:

pearson_corr: ## Correlation Analysis: x0_y **Strong correlation detected** (r=1.000) Consider incorporating a strong relationship between x0_y and output in your equation. The relationship appears to be positive - as one increases, the other increases.

spearman_corr: ## Correlation Analysis: x0_y **Strong correlation detected** (r=1.000) Consider incorporating a strong relationship between x0_y and output in your equation. The relationship appears to be positive - as one increases, the other increases.

lin_reg: ## Linear Regression Analysis: x1_y **Strong linear relationship** (R²=0.996) Linear equation: output ≈ 0.000 × x1_y + 0.000 Consider incorporating linear terms in your equation.

residual_var: ## Residual Variance Analysis: x0_y **Residual variance**: 0.057119 Low residual variance suggests the linear relationship captures most variation.

mutual_info: ## Mutual Information Regression Analysis: x1_y **Mutual information score**: 8.294050 Mutual information regression measures the dependency between variables in a regression context.

Figure 28: Examples of Data Analysis Tools Usage.

Related documents

• Collaborative Principle Evolution via Evidence Transfer for Scientific Discovery#1129907
• SCIRIGOR:Evaluating Open-Ended Scientific Analysis Beyond Final Scores#973792
• Large language model-assisted discovery of cohorts from scientific literature#617033
• Large language model-assisted discovery of cohorts from scientific literature#616533
• Scalable scientific interest profiling using large language models.#242317