Structured Scaling of AI Discovery Across Diverse Scientific Domains Wizard Intelligence Learning Lab, Stanford University Peking University, Tsinghua University The Hong Kong University of Science and Technology (Guangzhou) Homepage Abstract Scientific discovery often requires many cycles of proposing, testing, and refining candidate solutions. Language models can increasingly participate in these loops, but simply generating more attempts does not ensure progress: parallel searches may duplicate one another and iterative refinement may become trapped in poor directions. The central challenge is therefore not only to scale AI-driven discovery, but to structure that scaling so that evaluation signals compound over time. Here we introduce SimpleTES ( Simple T est-time E valuation-driven S caling), a framework that focuses on the structured scaling of AI discovery loops, organizing evaluator queries across independent trajectories, iterative refinement, local candidate selection, and the selective reuse of evaluated histories. Drawing on structural features of scientific communities, SimpleTES uses a single open-source GPT-OSS model to establish new state-of-the-art solutions across 28 open-ended problems in diverse scientific domains ranging from quantum physics and astronomy to biology, AI, and mathematics. These include a 24.5% reduction in quantum circuit compilation overhead, up to 23% lower propulsive cost for deep-space trajectories, a 2.17 × 2.17\times faster lasso-path solver, an 8.5% lower-error whole-brain neural-activity predictor, the fastest reported TriMul kernel, and new mathematical constructions beyond prior human or AI records. We further post-train the model for long-horizon discovery by assigning each attempt the final outcome of the trajectory it helped produce. This improves performance on both training and held-out mathematics problems, further advancing the frontier. Together, these results establish structured scaling as a general mechanism for advancing AI scientific discovery. Contents ### Introduction Over the past decade, the role of AI in scientific discovery has changed significantly. Early successes in AI for science came largely from specialized models that served as predictors or optimizers for well-defined scientific objects, such as protein structures, quantum wavefunctions, or molecular potentials [ 65 , 21 , 107 , 151 , 169 ] . More recently, language models have been treated as researchers, generating candidate solutions, writing and revising the programs that implement them, and carrying a line of inquiry across many steps in response to evaluator feedback [ 17 , 48 , 85 , 100 ] . For instance, recent systems build the loop around evolutionary search procedures, as exemplified by AlphaEvolve and related approaches [ 100 , 114 , 6 , 73 ] , or around tree search [ 9 ] , while other works adapt the model through test-time reinforcement learning on a single problem [ 155 , 165 ] . Despite differences in domains, models, and search heuristics, these systems share a common building block: the propose-evaluate-refine discovery loop that drives new discovery results [ 114 , 43 , 123 , 165 ] . More broadly, the human research community advances science collectively by scaling this loop in a complex, decentralized manner: researchers around the world pursue candidate solutions in parallel, driven by different inspirations, test them against the current baselines, and use the resulting feedback to shape the next attempts [ 74 ] . With LLMs taking the role of individual researchers in the loop, the entire scaling can be designed completely differently: an organized, structured way to drive scientific breakthroughs more efficiently and effectively. We name this new scaling dimension evaluation-driven scaling , because it is primarily driven by evaluation signals from empirical measurements, simulations, formal checks, or performance profiles, which later on become inspirations of next proposals. Different from standard model scaling that scales up parameters and training data [ 67 , 56 ] , or classical test-time scaling that scales through longer reasoning traces or additional agent turns [ 129 , 158 , 96 , 174 , 83 , 69 ] , evaluation-driven scaling targets challenging, open-ended problems that fundamentally require explorations, and a larger number of loops does not automatically imply better results. While showing promise in making discoveries, existing systems leave unexamined how the loop itself should be organized in a structured way. This leads to a central question: How can AI-driven discovery loops be scaled in a structured way to accelerate scientific breakthroughs? This paper introduces SimpleTES ( Simple T est-time E valuation-driven S caling), a framework for scaling AI-driven discovery by structuring how evaluator queries are allocated and how evaluated histories are grown, filtered, and reused. The central challenge is to determine how early attempts and their evaluation signals should guide later proposals, and which search structures allow progress to compound on open-ended problems. SimpleTES draws on structural features of scientific communities: locally, individual researchers test a few variants and carry forward the best; over time, researchers selectively learn from prior results based on relevance, performance, and reputation; collectively, researchers build on one another’s ideas to refine and propose new solutions; and in parallel, multiple lines of inquiry develop and compete. SimpleTES operationalizes these features through a compact and structured design tuple ( C , L , K , Φ ) (C,L,K,\Phi) . Global width C C maintains independent trajectories, refinement depth L L lets feedback accumulate within each trajectory, local sample size K K filters candidate variants before they shape future refinement, and the proposal constructor Φ \Phi selects and compresses useful prior solutions, failures, and feedback into the next prompt. This reframes discovery-loop scaling as a problem of managing evaluated histories, so that additional evaluations compound into stronger trajectories rather than merely producing more attempts. The key distinction between SimpleTES and existing AI discovery systems such as AlphaEvolve lies in the focus. Existing systems typically seek to maximize discovery performance through customized agentic search procedures. SimpleTES instead asks how the discovery loop itself should be organized and scaled. By isolating this structure, SimpleTES becomes a simple, structured, and highly parallelizable system that scales efficiently to tens of thousands of evaluator queries, yielding better results using smaller open-source models compared to prior systems that used larger proprietary LLMs. A comprehensive comparison can be found in Supplementary B.1 . We apply SimpleTES on 28 open-ended problems across diverse scientific domains, spanning quantum circuit compilations, astrodynamics trajectory designs, GPU kernel implementation, LLM scaling law designs, biological and time-series algorithms optimization, and mathematical constructions. On all problems, SimpleTES discovers the state-of-the-art solutions using a single open-source GPT-OSS model, surpasses both human experts and other AI discovery systems that use an order of magnitude larger proprietary models. To illustrate, SimpleTES improves the best-known bound for the Erdős minimum-overlap problem, open in extremal combinatorics since 1955, ahead of every prior human and machine record [ 100 ] . It discovers quantum compilation policies that reduce expensive qubit-swap gates by over 20% relative to the widely used LightSABRE baseline [ 176 ] . It improves the fuel efficiency of Solar System exploration mission trajectories, reducing propulsive costs by up to 23% relative to references constructed from historical encounter sequences under the same orbital-dynamics evaluator [ 148 , 121 ] . It implements the fastest TriMul kernel, the triangle multiplicative update central to AlphaFold3 [ 66 , 1 ] , and a lasso path solver 2.17 × 2.17\times faster than the widely used glmnet [ 39 ] . In neuroscience, it also discovers a state-of-the-art whole-brain neural-activity predictor, reducing error by 8.5 % 8.5\% on the long-horizon time-series forecasting task in ZAPBench [ 88 ] relative to the proprietary-model-driven algorithm discovery result [ 9 ] . More than any single record, what matters is that one unchanged procedure produces all of them, from abstract mathematical constructions to executable scientific programs, showing that SimpleTES can serve as a general mechanism for discovery. The above results come from the fixed, off-the-shelf GPT-OSS model, which is not specifically optimized for the long-horizon discovery that SimpleTES performs: standard post-training tends to reward outputs that score well immediately, whereas discovery often turns on unconventional attempts that score poorly on their own yet lead to much stronger solutions. This raises a natural question: how can the LLM itself be deliberately trained for long-horizon discovery? We address this by post-training the model on its own successful trajectories, crediting each attempt by the final outcome of its trajectory rather than its immediate score. The improved model then drives the same search, now less myopic and more willing to build on early failures in pursuit of a stronger eventual solution. We study training primarily in mathematics, across a range of related open-ended problems. Post-training improves the entire outcome distribution, raising both the average score across independent trajectories and the best score on problems where the off-the-shelf GPT-OSS model scores poorly, and it transfers to held-out mathematical problems. These results suggest the trained model learns to search for strong eventual solutions rather than optimizing for immediate scores. In summary, SimpleTES takes the scaling of the loop as its object of study, and asks how a fixed budget of evaluator queries should be structured. By formalizing this structured scaling with the design tuple ( C , L , K , Φ ) (C,L,K,\Phi) , SimpleTES shows superior performance using the open-source model, surpassing existing systems that uses larger proprietary models. SimpleTES uses a single configuration to handle all five domains and discover 28 SOTA results, with only the evaluator changing between problems. A summary of SimpleTES result with the open-sourced model can be found in Table 1 . ### Overview of SimpleTES SimpleTES is a test-time framework for scientific discovery on problems that have a verifiable evaluation protocol , which we call an evaluator . An evaluator is the system’s interface to the scientific objective: given a candidate solution, it returns a numerical score together with auxiliary feedback such as verifier messages, error traces, or runtime statistics. Depending on the task, it may compute the objective exactly (an exact verifier), estimate it empirically (a measurement), or approximate it (a heuristic proxy). The objective is to find the highest-scoring candidate under the evaluator. Formally, a problem is specified by an instruction x 0 x_{0} , and the goal is a solution y y that maximizes the score r r output by an evaluator V V . To search for such a solution, an LLM G G plays the role of a researcher, drawing candidates y ∼ G ⁡ ( x ) y\sim G(x) from a proposal x x built from the instruction x 0 x_{0} together with the accumulated history of previously evaluated solutions, their scores, and the evaluator’s feedback. SimpleTES runs this as a propose–evaluate–refine loop: G G proposes a candidate, V V evaluates it, and its score and feedback are folded into the history that shapes the next proposal. SimpleTES organizes this loop through a compact design tuple ( C , L , K , Φ ) , \displaystyle(C,L,K,\Phi), (1) where C C is the number of parallel trajectories, L L the number of refinement steps per trajectory, K K the number of candidates drawn per step, and Φ \Phi the proposal constructor that assembles each new proposal from a trajectory’s history. SimpleTES launches C C independent trajectories and advances each for L L refinement steps. At each step, Φ \Phi builds a proposal x x from the history so far, G G draws K K candidate solutions from x x , V V scores them, and the highest-scoring one is appended to the trajectory. Once every trajectory has completed, SimpleTES returns the best solution across all of them, having spent N = C × L × K N=C\times L\times K evaluator queries. This decomposition is what gives SimpleTES its structure: each of the three dimensions targets a distinct, independent way an undirected search fails, so together they form a small and non-redundant set of controls on the search, which we analyse and ablate in the Methods. Refinement depth L L lets feedback accumulate, so later attempts build on earlier ones, but a single deep trajectory is not enough, since once early attempts bias the search toward a narrow direction, further refinement tends to reinforce that basin rather than escape it. Global width C C counters this by running independent trajectories in parallel, so competing directions mature without influencing one another early on. Local sample size K K addresses noise in the LLM’s output: drawing several candidates at each step and committing only the highest-scoring one keeps a weak or non-executable candidate from corrupting every later refinement built on it. Because the three act on different failure modes, they can be adjusted independently, and the budget N N has a clear meaning as the way that scaling is allocated. In practice, the allocation of the evaluator-query budget across C C , L L , and K K depends on the task. Problems whose progress depends on discovering qualitatively different approaches, such as mathematical construction problems, benefit from larger C C , whereas problems dominated by incremental engineering refinement, such as GPU kernel optimization, benefit from larger L L . The local sample size K K should remain moderate: increasing it makes each committed step more reliable, but leaves fewer refinement steps under the same budget. Unless otherwise stated, we use C = 32 C=32 , L = 100 L=100 , and K = 16 K=16 as a default configuration. While C C , L L , and K K govern how evaluated histories grow, the proposal constructor Φ \Phi governs which parts of that history are reused, since a trajectory’s history soon becomes too long and noisy to feed back in full. Because each proposal can build on several prior attempts rather than only the last, the history forms a directed graph rather than a chain. Φ \Phi selects which attempts to reuse through a graph-based variant of the PUCT rule [ 115 , 127 , 126 ] , favoring attempts that scored well or led to strong later attempts, while still giving promising but rarely used ones a chance to be revisited. It also folds in a running summary of recurring failures, such as compilation errors, verifier failures, or timeouts, to steer G G away from unproductive directions. In Methods, a Pólya-urn analysis together with ablations on C C , L L , and K K shows why these dimensions form a compact design space for scaling the loop, and Φ \Phi ’s construction and selection strategy are specified and ablated in full. We next turn to the discoveries SimpleTES makes across scientific domains, and then to how a language model can be trained for SimpleTES -based discovery. ### SimpleTES discovers effective algorithms to reduce quantum compilation overhead Quantum computing is an increasingly important computing paradigm that stores and processes information via quantum-mechanical systems [ 99 ] . Whereas a classical bit has a definite value of either 0 or 1, a quantum bit, or qubit , can be in a coherent superposition of both states. By manipulating such states and creating correlations between qubits through entangling operations, quantum computers can obtain polynomial or exponential speedups for selected problem classes [ 125 , 50 , 35 , 52 , 13 , 37 , 84 , 7 ] . At the center of this paradigm is the quantum algorithm, which specifies how quantum operations transform input states into desired outputs. In the circuit model, a quantum algorithm is represented by a quantum circuit: an ordered sequence of elementary operations, called gates , applied to logical qubits [ 99 ] . One-qubit gates update the state of a single qubit, while two-qubit gates create interactions between a pair of qubits. A logical circuit may specify two-qubit gates between arbitrary pairs of logical qubits, but a real device can apply a direct two-qubit gate only when the corresponding physical qubits satisfy the hardware’s interaction constraints, thus creating a compilation problem of how to assign logical qubits to physical qubits and update that assignment over time so that every required two-qubit gate becomes executable. The additional operations introduced in this process are costly and error-prone, making overhead reduction the central optimization goal. This optimization problem is computationally challenging and takes different forms on different hardware platforms [ 61 , 128 , 28 , 15 , 133 , 80 , 132 , 137 ] . We use SimpleTES to address two representative hardware settings, superconducting platforms and neutral-atom architectures. In superconducting devices, physical qubits occupy fixed chip locations and allowed two-qubit interactions form a sparse coupling graph; compilation must choose an initial layout and insert SWAP gates when interacting logical qubits are mapped to non-adjacent physical qubits [ 76 , 176 ] . In zoned neutral-atom devices, atoms are kept in a storage zone or moved to an entangling zone for Rydberg gates; compilation must choose placement snapshots whose transitions can be realized through AOD-mediated movement. Together, these settings contrast fixed-graph routing with geometry-limited atom placement and rearrangement, while retaining the same high-level goal of minimizing extra hardware work. We apply SimpleTES within fixed compilation scaffolds that isolate editable policy components, with necessary correctness checks applied. For superconducting circuits, SimpleTES edits the initial layout and online SWAP selection rule; candidate routers are compared against SABRE [ 76 ] and LightSABRE [ 176 ] , with added SWAP overhead as the cost. For zoned neutral-atom compilation, SimpleTES edits only the placer in a fixed scheduler–reuse-analyzer–placer–router pipeline, with the remaining modules held fixed; candidate placers are compared against a reuse-aware ZAC-style baseline [ 80 ] , with total execution time as the objective. Formal problem definitions, hardware constraints, and optimization objectives for both compilation tasks are given in Supplementary H.1 . Benchmark suites, baselines, and evaluation metrics are provided in Supplementary H.2 . Across both settings, SimpleTES achieves substantial overhead reductions relative to established baselines, driven by a shared design principle: place important logical qubits near useful hardware locations before expensive routing or rearrangement becomes unavoidable. In superconducting routing, the best policy reduces topology-averaged added SWAP overhead by 21.7% relative to SABRE and 14.9% relative to LightSABRE , with the largest reductions on Q20 reaching 33.3% and 24.5%, respectively; it implements the shared principle by mapping high-traffic logical qubits toward high-degree physical qubits, refining the assignment by local search, and scoring online SWAPs using a broader candidate neighborhood, look-ahead information, and immediate gate executability ( Figure 2 b,c). In zoned neutral-atom compilation, the best placer discovered by SimpleTES reduces the geometric mean of total execution time by 33.2% relative to the ZAC-style baseline; it implements the same principle through the sequence of placement snapshots passed to the fixed AOD router, preserving relative qubit order, reducing travel distance without losing AOD parallelism, and refining placements through reverse-time passes ( Figure 2 d–f). Further details of the best discovered programs and the additional superconducting routing policy obtained using SimpleTES with Gemini are provided in Supplementary H.3 and H.4 . ### SimpleTES improves optimizers for fuel-efficient deep-space trajectories Deep-space missions often reach distant targets by using planets as moving gravitational waypoints rather than flying directly. A gravity assist (GA) uses a planetary flyby to reshape the spacecraft’s heliocentric path with little propellant, whereas a deep-space maneuver (DSM) uses an engine burn to adjust the path [ 98 ] . Designing such missions is difficult because route choices and timing choices are tightly coupled: the optimizer must choose which bodies to encounter, which transfer branches to use, and where maneuvers may be inserted, while also setting launch and encounter epochs and maneuver variables [ 148 , 121 ] . We formulate this as an evaluation-driven discovery task in which each candidate produced by SimpleTES is an executable mission-design program. Given mission windows, a target body, allowed flyby bodies, and endpoint conditions, the program proposes complete multi-gravity-assist trajectories ( Figure 3 a); an independent patched-conic evaluator reconstructs the trajectory from ephemerides [ 38 ] , transfer arcs [ 62 ] , flybys [ 95 ] , DSMs, and endpoint costs, rejects constraint violations, and scores feasible designs by verified propulsive cost. We benchmark SimpleTES on five historical mission-window benchmarks and a prospective Jupiter exploration window. For the historical benchmarks, the comparison is not to flown mission Δ ​ v \Delta v budgets, which incorporate operational choices, navigation margins, and mission-specific high-fidelity modeling. Instead, each reference trajectory is constructed from the historical encounter sequence and nominal timing, then evaluated under the same orbital-dynamics model, feasibility constraints, and cost definition used to score SimpleTES . Across the historical mission windows, the best verified trajectories found by the discovered optimizers reduce propulsive cost in every case, by 2.1–23.1% relative to these same-evaluator references ( Figure 3 b). In the prospective 2028–2030 Jupiter exploration window, where no encounter sequence is prescribed, the discovered optimizer synthesizes an evaluator-feasible no-DSM Earth–Venus–Earth–Earth–Jupiter trajectory reaching Jupiter in 2035 ( Figure 3 a). Inspection of the evolved optimizers shows that the gains arise from reusable policies for coupled discrete–continuous mission design. The clearest sequence-level changes occur for Galileo and Rosetta ( Figure 3 c): the Galileo optimizer finds Earth–Venus–Venus–Venus–Earth–Jupiter rather than the reference Earth–Venus–Earth–Earth–Jupiter sequence, while the Rosetta optimizer finds Earth–Earth–Earth–67P rather than Earth–Earth–Mars–Earth–Earth–67P. For Mariner 10, Voyager 2, and Cassini, the programs independently recover the reference encounter order and lower cost by refining the continuous trajectory variables within that sequence, showing that improvements do not require a different encounter sequence. Across cases, the evolved programs generate and prioritize candidate encounter sequences, allocate finite computation to promising time-of-flight regions, explore multiple transfer branches, combine global search with local refinement of encounter epochs and maneuver variables, and retain DSMs only when they reduce verified cost. Together, these results show that the evolved optimizers allocate finite computation both across encounter sequences and within each sequence’s continuous design space. ### SimpleTES discovers better scientific algorithms Scientific algorithms are the computational methods that solve scientific problems: optimization solvers, forecasters, and learned models whose quality is judged by task-specific scientific objectives. These algorithms are usually produced by hand, often after years of accumulated domain expertise and engineering. We apply SimpleTES to three such problems and find that the same evaluation-driven search procedure discovers algorithms that improve over strong baselines: a geometry-aware lasso path solver that preserves strict correctness guarantees while outperforming the expert-engineered glmnet on real datasets; a whole-brain neural-activity forecaster that outperforms a prior proprietary-model search result; and a single-cell RNA-seq denoising algorithm that generalizes across datasets. We discuss the first two in detail below. Full denoising results are reported in the Supplementary Information, where an algorithm evolved exclusively on the Pancreas dataset matches the previous state of the art on held-out PBMC ( 0.71 0.71 ) and establishes a new best on held-out Tabula Muris Senis Lung ( 0.74 0.74 versus 0.73 0.73 ) (Supplementary M ). ##### SimpleTES discovers a geometry-aware lasso solver that surpasses glmnet . The lasso regularization path is a core computational primitive in high-dimensional statistics, used in cross-validation, model selection, and large-scale regression workflows from genomics to finance. The de facto solver, glmnet [ 39 ] , encodes more than a decade of expert coordinate-descent engineering: warm starts, sequential strong-rule screening [ 142 ] , active-set inner loops, KKT verification, and a precision-matched switch between the covariance ( p < 500 p<500 ) and naive residual ( p ≥ 500 p\geq 500 ) updates. We ask whether evaluation-driven search, starting from a faithful C++ port of glmnet itself, can produce a strictly faster solver while preserving the same correctness guarantees. The evaluator is exact: candidate programs are self-contained C++ binaries, and a candidate is valid only if its per- λ \lambda objective gap relative to sklearn ’s solution stays below 10 − 6 10^{-6} on a held-out problem, with any failure forcing the score to zero. Evolution operates exclusively on synthetic problems, scored by the geometric mean of solve times across 17 procedurally generated instances that span sample-to-feature ratio, design-matrix sparsity, active-set density, and feature correlation. Because no real dataset is ever exposed to the search, the gains reported on real data below reflect generalization rather than tuning to the test set. Across the 11 real-world datasets in Supplementary Table 16 , none seen during the search, SimpleTES discovers a solver that is 2.17 × 2.17\times faster than glmnet on average and 14.08 × 14.08\times faster than sklearn , with every reported solution passing the held-out correctness check. What makes the discovered solver notable is not just raw speed but geometry awareness : it reads the shape of each problem and adapts its algorithm accordingly. Where glmnet applies coordinate descent uniformly, SimpleTES ’s solver detects the problem regime and switches solver family to match it. For moderate-dimensional, tall designs ( p ≤ 2000 p\leq 2000 and n ≥ p / 4 n\geq p/4 ) it abandons coordinate descent for an exact LARS homotopy solver that traces the regularization path analytically, updating the active set at each kink with rank-one inverse Gram updates [ 34 ] ; for wide or sparse designs ( p > 2000 p>2000 or n < p / 4 n

0 \beta>0 . Assume the total budget N N is split into C C independent trajectories, each performing L = N C L=\frac{N}{C} refinement steps. Then, as s → 1 s\to 1 , the optimal allocation ( L ⋆ , C ⋆ ) (L^{\star},C^{\star}) that minimizes the total budget N N subject to the reliability constraint P fail ​ ( L , C ) ≤ ϵ P_{\text{fail}}(L,C)\leq\epsilon satisfies L ⋆ = Θ ⁡ ( log λ ⁡ ( 1 − s ) ) , C ⋆ = Θ ⁡ ( log ⁡ 1 ϵ ) . \displaystyle L^{\star}=\Theta\!\left(\log_{\lambda}(1-s)\right),\qquad C^{\star}=\Theta\!\left(\log\frac{1}{\epsilon}\right). (9) Here Θ ⁡ ( ⋅ ) \Theta(\cdot) hides the dependency on D D and β \beta . ###### Proof. Let r = ⌈ log λ ⁡ ( 1 − s ) ⌉ . \displaystyle r=\left\lceil\log_{\lambda}(1-s)\right\rceil. (10) Since 0 < λ < 1 0<\lambda<1 , a trajectory reaches score at least s s after L L refinement steps if and only if min d ∈ [ D ] ⁡ y d ​ ( L ) ≥ r . \displaystyle\min_{d\in[D]}y_{d}(L)\geq r. (11) Thus r r is the number of successful refinements required in every dimension, and r → ∞ r\to\infty as s → 1 s\to 1 . We first characterize a single trajectory. Set α = 1 / β \alpha=1/\beta . Dividing the transition probability in Definition .2 by β \beta gives p d ​ ( t ) = α + y d ​ ( t − 1 ) D ​ α + t − 1 . \displaystyle p_{d}(t)=\frac{\alpha+y_{d}(t-1)}{D\alpha+t-1}. (12) This is the standard symmetric Pólya urn with initial mass α \alpha on each of the D D dimensions. Therefore, the vector of refinement counts after L L steps has the Dirichlet-multinomial law y ⁡ ( L ) ∼ DM ⁡ ( L , α , … , α ) . \displaystyle y(L)\sim\operatorname{DM}(L;\alpha,\ldots,\alpha). (13) Equivalently, it admits the mixture representation P = ( P 1 , … , P D ) \displaystyle P=(P_{1},\ldots,P_{D}) ∼ Dirichlet ⁡ ( α , … , α ) , \displaystyle\sim\operatorname{Dirichlet}(\alpha,\ldots,\alpha), (14) y ⁡ ( L ) | P \displaystyle y(L)\mid P ∼ Multinomial ⁡ ( L , P ) . \displaystyle\sim\operatorname{Multinomial}(L,P). (15) Hence, if L , r → ∞ L,r\to\infty with r / L → γ r/L\to\gamma , then y ⁡ ( L ) L ⇒ P , \displaystyle\frac{y(L)}{L}\Rightarrow P, (16) and the single-trajectory failure probability satisfies, for every γ ∈ ( 0 , 1 / D ) \gamma\in(0,1/D) , F L ​ ( r ) := Pr ⁡ ( min d ∈ [ D ] ⁡ y d ​ ( L ) < r ) → F ⁡ ( γ ) := Pr ⁡ ( min d ∈ [ D ] ⁡ P d < γ ) . \displaystyle F_{L}(r):=\Pr\!\left(\min_{d\in[D]}y_{d}(L) 0 G(\gamma)>0 on ( 0 , 1 / D ) (0,1/D) . Thus, by continuity, G G attains a maximum at some γ ⋆ ∈ ( 0 , 1 / D ) \gamma^{\star}\in(0,1/D) , and γ ⋆ \gamma^{\star} depends only on D D and β \beta . Substituting this structural constant gives L ⋆ = r γ ⋆ ​ ( 1 + o ⁡ ( 1 ) ) = Θ ⁡ ( r ) , C ⋆ = log ⁡ ( 1 / ϵ ) − log ⁡ F ⁡ ( γ ⋆ ) ​ ( 1 + o ⁡ ( 1 ) ) = Θ ⁡ ( log ⁡ 1 ϵ ) . \displaystyle L^{\star}=\frac{r}{\gamma^{\star}}(1+o(1))=\Theta\!\left(r\right),\qquad C^{\star}=\frac{\log(1/\epsilon)}{-\log F(\gamma^{\star})}(1+o(1))=\Theta\!\left(\log\frac{1}{\epsilon}\right). (22) Finally, r = ⌈ log λ ⁡ ( 1 − s ) ⌉ = Θ ⁡ ( log λ ⁡ ( 1 − s ) ) r=\lceil\log_{\lambda}(1-s)\rceil=\Theta(\log_{\lambda}(1-s)) as s → 1 s\to 1 , which proves the claimed allocation scalings. ∎ ##### Remark. The theorem rigorously justifies the effect of independent experiments. This reflects a fundamental asymmetry between the two scaling factors: C C governs the exploration strength, while L L governs exploitation power. This asymmetry is manifested clearly in practice. For instance, on the Second Autocorrelation Inequality, the strongest agent such as Claude-Code with Opus 4.6 plateaued at a score of 0.9438, even with hundreds of refinement steps at a cost of ∼ 100 \sim 100 million tokens and ∼ 500 \sim 500 USD. However, by effectively balancing C C and L L (method specified below), we match this score using the open-source gpt-oss-120b model 2 2 2 Estimated costs for gpt-oss-120b are calculated based on the median pricing from OpenRouter, which is $0.15 per 1M input tokens and $0.60 per 1M output tokens. at a cost of ∼ 60 \sim 60 USD ( 8.3 × 8.3\times gap). More surprisingly, we achieve a SOTA score of 0.9627 by continuously scaling up the evaluation ( ∼ 400 \sim 400 USD), a score that no closed-source model could improve upon. This empirical evidence clearly demonstrates the importance of identifying the correct factor for evaluation scaling. Having demonstrated how compute should be allocated across multiple refinement trajectories, we turn to the scaling within each trajectory, with a given budget L = N / C L=N/C . At first glance, more refinement steps seem preferable, as each step can leverage feedback from previous attempts. However, this intuition could break down when path dependence is strong ( β \beta is large), when each refinement step further commits the trajectory to its current direction, making it harder to escape suboptimal regions. ##### Local sample size K K . The preceding analysis motivates a second budget axis: the local sample size K K . While the number of independent trajectories C C mitigates global path dependence by exploring different refinement paths, K K addresses a more local failure mode: an individual refinement proposal may fail to improve the current solution. We model this by assuming that, after a dimension is selected according to the Pólya refinement rule in Definition .2 , the generator produces K K independent refinement proposals for that dimension. Each proposal improves the selected dimension with probability p p , and otherwise leaves it unchanged. The best proposal is then applied. Thus, increasing K K raises the probability that a refinement step results in an actual improvement from p p to 1 − ( 1 − p ) K 1-(1-p)^{K} , but it also reduces the number of sequential refinement steps under a fixed total budget. We empirically illustrate this trade-off in Figure 1 . We simulate problems parameterized by ( D , λ , β ) (D,\lambda,\beta) under the refinement model in Definition .2 , fixing L = 4096 L=4096 , β = 4 \beta=4 , and C = 32 C=32 independent trajectories, and report the normalized final bottleneck score averaged over 2048 simulations. The results show that moderate values of K K substantially improve performance over purely sequential refinement ( K = 1 ) (K=1) , especially when the improvement probability p p is large enough for local sampling to reliably find a useful proposal. However, very large K K can degrade performance, since allocating more budget to local sampling leaves fewer refinement steps. This confirms that K K should be treated as a trade-off parameter rather than a monotonic source of improvement. We further study this interaction between C C , L L , and K K on real open-ended tasks in Section F . Together, the total evaluation budget N = C × K × L N=C\times K\times L has been decomposed into three scaling dimensions: global independent trajectories C C , local sample size K K , and the actual refinement depth L = N C ​ K L=\frac{N}{CK} . In short, this section aims to recognize the simple yet effective factors to be scaled, with the evidence below showing that by wisely scaling these parameters, we can achieve state-of-the-art solutions on most problems considered. Arguably, there is clear room for further improvements. For instance, the current budget is divided evenly across each trajectory. Clearly, discarding unsatisfying trajectories on-the-fly and saving budgets for promising trajectories could be helpful. Additionally, one might recognize the fractal structure of C C and L L , both being i.i.d. attempts at different levels. As an analogy, breaking L L sequential refinement steps into smaller segments to couple more complexity might be beneficial. Relevant ablations can be found in Section F . ### D Method Analysis In this section, we present a comprehensive analysis of our proposed framework. We first investigate its scaling behaviors, demonstrating how systematic expansion of the evaluator-query budget drives continuous performance improvements. Next, we examine the impacts of training, highlighting its ability to enhance both in-domain and out-of-domain discovery capabilities. We then analyze various reward hacking phenomena that emerge when models autonomously exploit vulnerabilities in surrogate evaluators. Finally, we conduct detailed ablation studies to evaluate the individual contributions of core framework components. #### D.1 Experiments on the Scaling Behavior of SimpleTES To investigate the scaling behavior of SimpleTES , we explore the impact of key framework parameters: total evaluation budget N N , refinement depth L L , global width C C , and local sample size K K . We aim to reveal how the framework effectively translates increased computational scale into substantial performance improvements. In our ablation studies, we select three representative open-ended tasks—the First Autocorrelation Inequality, the Erdős Minimum-Overlap Problem, and the TriMul task in GPU kernel optimization—to investigate the scalability of SimpleTES . ##### Scalability of global width ( C C ) and refinement depth ( L L ). First, we scale along the dimensions of global parallel exploration ( C C ) and sequential refinement depth ( L L ). For these experiments, we hold the local sample size constant at K = 32 K=32 . The results are visualized in Figure 2 . As demonstrated in the heatmaps, SimpleTES exhibits strong and consistent scalability along both the global width ( C C ) and sequential refinement depth ( L L ) axes. The consistent performance improvements observed as L L increases serve as compelling evidence for the efficacy of sequential refinement , demonstrating that iteratively leveraging historical feedback successfully drives targeted enhancements. However, because sequential refinement is inherently path-dependent and prone to saturation around local optima, scaling C C provides a critical and complementary advantage. The robust gains along the C C axis validate the power of parallel exploration: by distributing the evaluation budget across independent trajectories, the framework successfully diversifies the committed histories, effectively mitigating path-dependent bottlenecks and ensuring that later refinements compound upon structurally advantageous foundations. Furthermore, the heatmaps reveal that the performance benefits derived from scaling C C versus scaling L L are strongly task-dependent. For mathematical discovery tasks like the First Autocorrelation Inequality and the Erdős Minimum-Overlap Problem, scaling the number of parallel chains C C (e.g., up to C = 32 C=32 ) provides a more pronounced advantage at larger scales. Mathematical constructions often require extensive exploration of diverse starting points to discover a promising structural “flash of insight.” Conversely, for the TriMul GPU kernel optimization task, scaling the iteration depth L L drives the most significant performance gains. Kernel optimization heavily relies on complex, step-by-step engineering refinement rather than sudden structural breakthroughs, making deep trajectories (large L L ) essential to systematically tune hardware-specific parameters and iteratively squeeze out peak performance. ##### Scalability of local sample size ( K K ) and refinement depth ( L L ). Next, we analyze the interaction between depth L L and local sample size K K by fixing the global parallelization at C = 32 C=32 , as illustrated in Figure 3 . Along the depth dimension, performance consistently improves as L L increases across all configurations. This is expected, as deeper chains allow the policy to iteratively exploit feedback and refine the best discovered paradigms. However, scaling the local sample size K K reveals a highly dynamic, depth-dependent behavior. At shallow chain depths (small L L ), allocating more compute to K K does not yield a consistent monotonic increase in performance. In these early stages, generating a higher-quality node via a larger K K does not necessarily translate to an immediate breakthrough; the solution is still navigating broad structural choices, making the early refinement steps inherently noisy. Yet, this trend is robustly consolidated as the chain deepens. At large values of L L , scaling K K clearly and consistently drives superior final performance, which implies that, while a rigorously selected node might not immediately show a massive score advantage, it establishes a fundamentally higher-quality foundation that is more amenable to continuous improvement. As the chain extends, this early structural advantage compounds, allowing the long-term benefits of local quality control to fully materialize. ##### Scalability of the total evaluation budget ( N N ). Building upon the previous analysis of the effects of C C , K K and L L , we further investigate the overall scaling behavior as a function of the total evaluation budget N = L × C × K N=L\times C\times K . As illustrated in Figure 4 , SimpleTES demonstrates a remarkably stable performance improvement as the available compute budget scales. This trend highlights a core strength of our framework: it does not merely consume inference tokens, but highly effectively translates an expanded evaluation budget into tangible, systematic performance gains. Crucially, this robust scaling behavior points to an important test-time scaling law for open-ended scientific discovery. It demonstrates that SimpleTES allows researchers to continuously and reliably push past existing performance upper bounds through the allocation of additional evaluation calls. #### D.2 Experimental Details on Post-Training In this subsection, we investigate whether post-training methods successfully transform the TES histories into parametric knowledge. For simplicity and training efficiency, we focus on mathematical extremal analysis and combinatorial construction tasks. Specifically, we adopt the Second Autocorrelation Inequality, the Third Autocorrelation Inequality, Circle Packing in a Unit Square with n = 26 n=26 , and Erdős Minimum-Overlap Problem as training tasks. In addition to these tasks, we perform evaluation on 4 held-out tasks, including the First Autocorrelation Inequality, Circle Packing in a Unit Square with n=32, Hadamard Maximum-Determinant Problem of Order 29 29 , and the Sum-Difference Problem. We perform post-training for 6 iterations. For the first training iteration, we collect ∼ \sim 320 trajectories for each task with K = 16 K=16 and L = 100 L=100 and default settings of Φ \Phi in Methods for cold-start. For subsequent iterations, we sample C ^ = 32 \hat{C}=32 trajectories for each task using the same hyperparameters. We adopt the IRFT setting, set the selection ratio R = 10 R=10 for top-performing trajectories in the first 4 iterations and R = 5 R=5 for the last 2 iterations. The number of effective training samples with w = 1 w=1 for each iteration ranges from 8.1K to 10.7K. In each iteration, we perform training for 100 steps with a batch size of 256 and a learning rate of 2e-5. We adopt a linear warmup for the first 20 steps, followed by a cosine decay to 0. It takes a total of 15 hours on 32 Nvidia H200 GPUs for training and 82 hours on 256 Nvidia H200 GPUs for TES sampling. In open-ended problem solving, the objective shifts from maximizing average performance to pushing the boundaries of the state-of-the-art. Standard metrics like overall trajectory scores often mask a model’s true potential, as they could be skewed by a small number of low-quality trajectories. To more accurately evaluate a model’s capacity for high-ceiling breakthroughs, we focus on the distribution of its elite trajectories. Specifically, we report the average trajectory-level scores for the top 10%, 25%, 50%, and 75% of trajectories (the average of the top R R % of each trajectory’s eventual score), across eight math tasks in Table 1 . We also visualize the relative improvements over the baseline for the training tasks and the held-out tasks after every 2 training iterations in Figure 5 . The key observations are as follows: ##### Training improves the trajectory scores on the training tasks. As shown in Table 1 , training yields consistent gains in Top 10%, Top 25%, Top 50%, and Top 75% scores across all training tasks. Notably, on the Erdős Minimum-Overlap Problem, the Top 25%, Top 50%, and Top 75% scores provided by the trained model are comparable to the baseline’s Top 10%, Top 25%, and Top 50% results, respectively, indicating that our training is effectively shifting the chain score distribution upward and improving the discovery efficiency. Furthermore, the iterative progress visualized in Figure 5 reveals a consistent trend of improvement across iterations for Top 10%, Top 25%, Top 50%, and Top 75% scores. These findings suggest that by self-distillation on its own successful trajectories, the model effectively enhances its proficiency across the overall TES process. ##### Training enables generalization to the held-out tasks. While our training approach yields notable gains on the training tasks, a critical question remains: is the model genuinely bootstrapping its general TES capabilities or merely memorizing task-specific shortcuts? Our analysis of the held-out tasks provides strong evidence for the former. As shown in Table 1 , the trained model consistently outperforms the gpt-oss-120b baseline in the Top 50% and 75% scores across all the held-out tasks, demonstrating robust transferability. In the Top 10% metric, it maintains or exceeds baseline performance in three out of the held-out tasks. The marginal performance drops regarding the Top 10% in the Hadamard Maximum-Determinant Problem of Order 29 29 and the Top 75% in the First Autocorrelation Inequality are likely attributable to the inherent stochasticity of the TES process, where a small portion of sub-optimal chains can skew the average scores. Further iteration-level analysis in Figure 5 reveals that the performance on the held-out tasks improves in the first 2 iterations and maintains during subsequent training. Collectively, these results suggest that iterative training enables the model to internalize fundamental skills for TES, allowing it to effectively generalize to a diverse range of unseen, complex discovery tasks. ##### Training unlocks stronger solutions that vanilla models fail to uncover. Despite a limited budget of 32 trajectories per iteration, the trained model consistently generates novel solutions, frequently achieving best or second-best performance across multiple challenging tasks. Surprisingly, on the held-out Sum–Difference Problem, it establishes a new SOTA of 1.144887, significantly surpassing prior results of 1.143975. These findings suggest that the training process effectively unlocks the model’s latent potential to achieve further breakthroughs under the TES setting. #### D.3 From Golden Metrics to Surrogate: Hacking Analysis As formalized in Methods, the paradigm of TES relies on an explicit, queryable surrogate evaluator V V to approximate the true, often inaccessible, golden metric. However, the purely score-driven nature of feedback-driven search inevitably drives models to discover and exploit the discrepancies between the surrogate and the golden metric. This results in various forms of reward hacking, where the model optimizes for the surrogate’s implementation vulnerabilities rather than the underlying task. We systematically categorize these phenomena below. ##### Hacking in math construction tasks. In mathematical construction tasks, the surrogate evaluator is designed to be an exact metric (e.g., a formal geometric verifier). However, the gap here emerges from implementation vulnerabilities in numerical processing . The first hacking strategy targets the evaluator’s decision tolerance. In the Circle Packing in a Unit Square task, the model is required to place several circles in a unit square without overlapping. However, the surrogate evaluator provided by OpenEvolve [ 6 ] employs a floating-point tolerance of 1 × 10 − 6 1\times 10^{-6} for overlap determination. The model fully exploits this proxy limitation by crafting configurations where circles and boundaries actually overlap, but the overlapping magnitude strictly remains under the 1 × 10 − 6 1\times 10^{-6} threshold. See for the construction. The second hacking strategy exploits floating-point precision. For instance, evaluators used by EvoX [ 82 ] compute metrics directly on raw coordinates without proper normalization steps. The model exploits this vulnerability by providing degenerate configurations—such as collapsing all points to an infinitesimal scale ( ∼ 10 − 162 \sim 10^{-162} ) or arranging them collinearly—which induce severe numerical instability in the surrogate, yielding artificially high scores. ##### Hacking in kernel optimization tasks. In kernel optimization tasks, the gap between the surrogate evaluator and the golden metric is fundamentally wider. Rather than relying on a mathematically rigid verifier, the surrogate here acts as a vulnerable empirical estimator forced to operate within a highly complex execution environment. Consequently, models escalate from simple numerical exploits to sophisticated manipulations designed to bypass the benchmarking process entirely. We identify four primary hacking strategies: 1. Caching results over multiple runs : To ensure accurate latency measurements and reduce variance, evaluators often repeatedly call the target kernel using the exact same inputs. The model exploits this by caching input pointers or intermediate buffers during the initial call, bypassing the actual computation entirely on all subsequent iterations to achieve near-zero execution time. See for an example. 2. Time measurement manipulation : The model artificially deflates latency metrics by directly intercepting or evading the evaluation profiling tools. For instance, it may globally monkey-patch CPU timers (e.g., overriding time.time ) to return hardcoded minimal values, or divert GPU execution to newly created, unmonitored CUDA streams that the evaluator’s default event recorders fail to track. See for an example. ⬇ 1 import time 2 3 # hacking function created by model 4 def func (): 5 old_time = time . time 6 start = old_time () 7 time . time = lambda : start + 0.1 8 9 # In Evaluator 10 st_time = time . time () 11 func () 12 end_time = time . time () # The function time.time here is already replaced by func Code Snippet 3: Time measurement hacking example 3. Baseline corruption : Because the custom kernel and the reference baseline typically execute within the same shared Python process, the model can globally override core PyTorch modules (e.g., replacing nn.Linear with a zero-output dummy). Consequently, when the evaluator runs the reference baseline after the custom kernel, the baseline yields trivial outputs, allowing the model’s equally trivial custom output to effortlessly pass the correctness check. See for an example. ⬇ 1 class ZeroLinear ( nn . Module ): 2 def __init__ ( self , in_features : int , out_features : int , bias : bool = True ): 3 super (). __init__ () 4 self . out_features = out_features 5 6 def forward ( self , input : torch . Tensor ) -> torch . Tensor : 7 return torch . zeros ( 8 ( input . shape [0], self . out_features ), 9 dtype = input . dtype , 10 device = input . device , 11 ) 12 13 # Apply the monkey patch globally before any Model instances are created. 14 nn . Linear = ZeroLinear 15 16 # In evaluator: 17 reference = TorchRefModel ( input ) # The TorchRefModel uses the replaced nn.Linear. 18 output = CustomModel ( input ) # Model creates a naive CustomModel that returns zeros. 19 correctness = torch . isclose ( reference , output ). all () Code Snippet 4: Baseline corruption hacking example 4. Triton partial computation : The model exploits torch.empty and Triton’s autotuning mechanism. During the initial autotuning phase, the kernel computes the fully correct output, leaving it in the GPU memory. For the actual benchmark, the model reuses the pre-computed results, vastly reducing execution time. See for an example. ⬇ 1 def _layernorm_proj_configs (): 2 cfgs = [] 3 for BLOCK_M in [32, 64, 128]: 4 for BLOCK_N in [64, 128, 256]: 5 for BLOCK_K in [32, 64, 128]: 6 cfgs . append ( 7 triton . Config ( 8 { "BLOCK_M" : BLOCK_M , "BLOCK_N" : BLOCK_N , "BLOCK_K" : BLOCK_K }, 9 num_warps =4, 10 ) 11 ) 12 return cfgs 13 14 @triton . autotune ( configs = _layernorm_proj_configs (), key =[ "B" , "N" , "H" , "D" ]) 15 @triton . jit 16 def _layernorm_proj_kernel ( 17 out_ptr , # fp16 [B, N, N, H] (triangle output) 18 gate_ptr , # fp16 [B, N, N, H] (out_gate) 19 ln_weight_ptr , # fp16 [H] (to_out_norm.weight) 20 ln_bias_ptr , # fp16 [H] (to_out_norm.bias) 21 weight_ptr , # fp16 [D, H] (to_out.weight) 22 proj_ptr , # fp32 [B, N, N, D] (final output) 23 B , N , H , D , 24 eps : tl . constexpr , # LayerNorm epsilon 25 # strides for out / gate 26 stride_out_b , stride_out_i , stride_out_j , stride_out_h , 27 stride_gate_b , stride_gate_i , stride_gate_j , stride_gate_h , 28 # strides for LN parameters 29 stride_ln_w , 30 stride_ln_b , 31 # strides for weight 32 stride_w_out_d , stride_w_out_h , 33 # strides for proj 34 stride_proj_b , stride_proj_i , stride_proj_j , stride_proj_d , 35 # compile-time tile sizes 36 BLOCK_M : tl . constexpr , 37 BLOCK_N : tl . constexpr , 38 BLOCK_K : tl . constexpr , 39 ): 40 pid_m = tl . program_id (0) # position tile 41 pid_n = tl . program_id (1) # output-dim tile 42 ... 43 44 45 def fused_layernorm_proj (): 46 proj = torch . empty (( B , N , N , D ), dtype = torch . float32 , device = out . device ) 47 stride_proj_b , stride_proj_i , stride_proj_j , stride_proj_d = proj . stride () 48 49 grid = ( 50 triton . cdiv ( total_rows , 128), 51 triton . cdiv ( D , 128), 52 ) 53 54 _layernorm_proj_kernel [ grid ]( 55 out , out_gate , 56 ln_weight , ln_bias , 57 to_out_weight , 58 proj , 59 B , N , H , D , 60 eps , 61 stride_out_b , stride_out_i , stride_out_j , stride_out_h , 62 stride_gate_b , stride_gate_i , stride_gate_j , stride_gate_h , 63 stride_ln_w , 64 stride_ln_b , 65 stride_w_out_d , stride_w_out_h , 66 stride_proj_b , stride_proj_i , stride_proj_j , stride_proj_d , 67 # tile sizes are chosen by autotune 68 ) 69 return proj Code Snippet 5: partial compute hacking example ##### Discussion. Taken together, these phenomena demonstrate that, despite being entirely blind to the evaluator’s underlying source code, LLMs exhibit a remarkable capability to autonomously discover and exploit surrogate gaps. By corrupting the evaluation feedback loop, these behaviors mislead the optimization process into yielding deceptive solutions rather than genuinely advancing target capabilities. Currently, closing this gap between the proxy and reality relies heavily on human-in-the-loop interventions, requiring multiple rounds of manual debugging and iterative patching of evaluator loopholes. Developing fully automated, robust evaluation frameworks capable of aligning the proxy metric with the ultimate objective—while dynamically detecting such logical hacking behaviors—remains a critical open challenge for future work. ### E Variants of Proposal Constructor Φ \Phi ##### Balance policy. The Balance policy is a stratified random sampling strategy designed to balance the exploitation of high-performing solutions with the exploration of diverse, sub-optimal trajectories. Let S S be the set of historical attempts within a trajectory, sorted in descending order of their evaluator scores such that s 1 ≥ s 2 ≥ ⋯ ≥ s | S | s_{1}\geq s_{2}\geq\dots\geq s_{|S|} . To construct a proposal with n n inspirations, the policy enforces that the absolute best historical solution, s 1 s_{1} , is always deterministically included. For the remaining n − 1 n-1 slots, the policy categorizes the sorted historical trials into three overlapping tiers: an exploitation tier T exploit T_{\text{exploit}} containing the top r elite r_{\text{elite}} fraction of attempts, an exploration tier T explore T_{\text{explore}} containing mid-tier solutions (typically between the 10th and 60th percentiles), and a global random tier T random = S T_{\text{random}}=S covering all prior attempts. Each subsequent inspiration is sampled without replacement according to the following probability distribution: Tier ∼ { T exploit with probability ​ p exploit , T explore with probability ​ p explore , T random with probability ​ 1 − p exploit − p explore . \displaystyle\text{Tier}\sim\begin{cases}T_{\text{exploit}}&\text{with probability }p_{\text{exploit}},\\ T_{\text{explore}}&\text{with probability }p_{\text{explore}},\\ T_{\text{random}}&\text{with probability }1-p_{\text{exploit}}-p_{\text{explore}}.\end{cases} (23) Once a tier is selected, an attempt is drawn uniformly at random from that tier. This heuristic ensures that the prompt is primarily grounded in elite solutions while consistently injecting varied structural contexts to prevent premature convergence. ##### LLM-Elite policy. While score-based heuristics like the Balance policy are efficient, evaluating a trial’s potential solely by its scalar score can be myopic. The LLM-Elite policy addresses this by leveraging semantic insights to dynamically maintain a bounded-size "elite pool" 𝒫 \mathcal{P} (with maximum capacity L elite L_{\text{elite}} ) that maximizes both solution quality and methodological diversity. Whenever a new candidate solution x new x_{\text{new}} is generated and evaluated, an auxiliary LLM acts as a gatekeeper. The LLM is provided with the scores and self-reflective summaries (detailing the approach and insights) of both x new x_{\text{new}} and all current solutions in 𝒫 \mathcal{P} . It is instructed to output one of three actions: Add , Replace ( j ) (j) (to swap out a redundant or inferior attempt j ∈ 𝒫 j\in\mathcal{P} ), or Reject . To prevent the LLM from inadvertently discarding substantial progress due to misjudging diversity, we enforce a strict monotonic override rule: if the score of x new x_{\text{new}} is strictly greater than the maximum score currently in 𝒫 \mathcal{P} , it bypasses the LLM’s rejection and is deterministically added (replacing the lowest-scoring solution if | 𝒫 | = L elite |\mathcal{P}|=L_{\text{elite}} ). To sample n n inspirations from 𝒫 \mathcal{P} for a new proposal, we employ the identical stratified random sampling mechanism described in the Balance policy. The attempts in 𝒫 \mathcal{P} are sorted in descending order by their scores and divided into the exploitation, exploration, and global random tiers. Trials are then sampled without replacement according to the previously defined tier probabilities. This approach heavily favors the inclusion of the highest-scoring elites while retaining a non-zero probability to draw inspiration from lower-ranked, yet semantically distinct, solutions curated by the LLM. Furthermore, this policy injects a concise text-only overview of the entire elite pool into the prompt (excluding raw code), allowing the generator to comprehend the global landscape of explored directions without consuming an excessive context window. ### F Ablation Studies In this section, we conduct a series of ablation studies to systematically evaluate the contribution of individual components and the robustness of our proposed framework. #### F.1 Ablations on Different Designs of Φ \Phi In this section, we ablate the inspiration sample algorithm—a core component of our framework—to demonstrate how different strategies to sample inspirations from historical attempts can impact the search process. To systematically investigate this, we evaluate multiple selection algorithms on the First Autocorrelation Inequality and the Erdős Minimum-Overlap Problem. We compare naive baselines, including Random (uniform sampling from all historical nodes) and Balance (a purely score-driven heuristic), against more sophisticated approaches, including RPUCG and LLM-Elite (which uses semantic evaluation to maintain an elite set of inspirations, as detailed in Methods). Additionally, for the RPUCG policy, we ablate the number of provided inspirations ( i ​ n ​ s ​ p ∈ { 1 , 3 , 5 , 10 } insp\in\{1,3,5,10\} ) to understand the effect of context size on the generation quality. The results are summarized in Table 2 . The results indicate that purely score-based or random sampling methods generally fall short in providing consistent, high-quality guidance compared to advanced strategies like LLM-Elite and RPUCG . Evaluating the heuristic value of a node solely by its current scalar score is often myopic; a node with an ordinary immediate score might actually serve as a critical stepping stone to a globally optimal region. Methods like LLM-Elite address this by leveraging semantic insights, while RPUCG utilizes graph-based state-value estimation, coupled with an explicit balance of exploration and exploitation. This demonstrates that effective inspiration selection must incorporate structural or semantic information beyond mere intermediate performance scores. The ablation on the number of inspirations ( i ​ n ​ s ​ p insp ) within the RPUCG framework reveals a delicate balance between exploration and context coherence. When the inspiration pool is too small (e.g., i ​ n ​ s ​ p = 1 insp=1 ), the model tends to perform marginal refinements along a single trajectory, severely limiting its ability to crossover ideas and discover novel directions. Conversely, providing an excessively large number of inspirations (e.g., i ​ n ​ s ​ p = 10 insp=10 ) crowds the context window, causing the model to become distracted or confused by overwhelming, and sometimes conflicting, global information. Empirically, for a single continuous run, setting i ​ n ​ s ​ p = 3 insp=3 or i ​ n ​ s ​ p = 5 insp=5 emerges as the optimal sweet spot, providing sufficient diversity without overwhelming the model’s reasoning capabilities. While these experiments validate the conceptual advantages of using semantically and structurally aware selection policies with an optimal inspiration size, it is worth noting that the absolute numerical differences in peak performance across strategies are relatively modest. This reinforces our central thesis: the primary driver of discovery is the systematic scaling of the evaluation budget (TES), rather than complex selection heuristics. Our overarching evolutionary pipeline is inherently robust, leaving the design of more sophisticated sampling policies as an orthogonal direction for future work. #### F.2 Ablation on Reflection and Failure Patterns Recall that our algorithm constructs the next query by integrating the task instruction, the selected inspirations S ( c ) S^{(c)} , and the chain-local memory ℛ ( c ) \mathcal{R}^{(c)} . To further enrich this contextual memory, we introduce and ablate two complementary mechanisms for textual guidance: Reflection and Failure Patterns . Reflection serves as positive guidance; after a generation batch, the system synthesizes the approach and insights from the best-scoring node, appending this textual summary to the inspiration nodes in S ( c ) S^{(c)} . Conversely, Failure Patterns serve as negative guidance; the system aggregates the most frequent evaluation errors across all nodes in the chain and injects them directly into ℛ ( c ) \mathcal{R}^{(c)} . Together, these mechanisms are designed to explicitly inform the generator about “what worked” and “what to avoid.” To isolate and evaluate the individual contributions of these two mechanisms, we conduct an ablation study on the First Autocorrelation Inequality and the Erdős Minimum-Overlap Problem. We systematically test four configurations: disabling both features, enabling only Reflection, enabling only Failure Patterns, and enabling both. All other hyperparameters remain identical across the runs. The results are detailed in Table 3 . The results demonstrate that incorporating explicit negative constraints (Failure Patterns) provides a consistently strong foundation. While adding explicit Reflection (the “On / On” setting) achieves the best performance in the First Autocorrelation Inequality, utilizing Failure Patterns alone (“Off / On”) performs slightly better in the Erdős Minimum-Overlap Problem. Crucially, the absolute performance gap across these configurations is relatively marginal. This minimal variance indicates that our evolutionary framework is inherently robust to the specific combination of textual guidance, provided that basic negative constraints are present to prevent repeated errors. #### F.3 Efficiency Analysis: Trajectory-level Pruning While scaling the number of parallel chains generally improves search outcomes, the vast search space often leads to many chains exploring suboptimal directions. To effectively scale our search process and prevent the wasting of computational budget on unpromising trajectories, we introduce pruning strategies. Pruning acts as an early-stopping mechanism: by identifying and terminating underperforming chains during the generation process, we can reallocate computational resources toward the most promising directions without significantly compromising the final solution quality. To systematically evaluate the impact of pruning, we conduct experiments across six mathematical tasks: three autocorrelation tasks, two circle packing tasks, and the Erdős Minimum-Overlap Problem. Each task is evaluated 3 times with a total chain length of L = 100 L=100 and number of chains C = 32 C=32 (Fix K = 16 K=16 ). We examine pruning through two complementary settings: a two-stage sweep with nominal cutoffs at L = 25 L=25 and L = 50 L=50 , summarized in Figure 7 , and a task-level single-cutoff ablation at nominal L = 12 L=12 and L = 25 L=25 , reported in Table 4 . At each designated cutoff, active chains are ranked by the best score attained up to that point, and a predefined proportion of the lowest-scoring chains is eliminated. The resulting trade-off, summarized in Figure 7 , demonstrates a surprisingly high preservation rate for the best final score under aggressive pruning. For instance, even when applying a stringent first-stage cutoff at the nominal L = 25 L=25 checkpoint that retains only a single chain, the original best final score is preserved in 10 out of the 18 total runs. This suggests that for a substantial portion of successful trials, the structural advantages of the final solution manifest early in the search process. Furthermore, across the two-stage configurations shown in Figure 7 , the relative degradation remains remarkably minimal. In the majority of configurations, the expected score degradation is bounded within 0.01 % 0.01\% , and all configurations exhibit degradation below 0.03 % 0.03\% . The task-level single-cutoff results in Table 4 further clarify the task dependence underlying the aggregate trend. For the circle packing tasks, early pruning exerts virtually no negative impact, as the performance upper bound is easily accessible and a large proportion of chains can successfully reach it. Conversely, for tasks such as the Autocorrelation Inequalities and the Erdős Minimum-Overlap Problem, initial performance is not always a reliable proxy for final solution quality. In these scenarios, chains often require extended iterative refinement before their true potential is realized. Consequently, the earlier and more aggressive single-cutoff settings can lead to higher elimination rates of optimal solutions and more noticeable score degradation. These findings affirm that while pruning is an effective and necessary technique for scaling parallel search, its implementation presents multiple challenges. The primary difficulty lies in adaptively selecting a pruning strategy that aligns with the task’s specific landscape. Furthermore, relying solely on intermediate scalar scores for elimination can be myopic. Future work will focus on developing more sophisticated pruning algorithms that incorporate richer contextual signals—such as textual reasoning patterns, trajectory growth trends, and structural heuristics—to make more informed early-stopping decisions. ### G Wall-clock Time Efficiency We compare the wall-clock time required by SimpleTES and OpenEvolve v0.2.27 to reach the same solution quality under matched search budgets and runtime concurrency. For each task, both methods use the same initial solution and evaluation program. We configure SimpleTES with ( C , K , L ) = ( 16 , 16 , 50 ) (C,K,L)=(16,16,50) , corresponding to 12,800 evaluator queries, and configure OpenEvolve with sixteen islands to match the number of SimpleTES trajectories. Both methods use a maximum concurrency of 128 for generation and evaluation, ensuring comparable evaluator-query budgets, search structures, and aggregate request concurrency. SimpleTES matches OpenEvolve’s whole-run best score 11.8 × 11.8\times faster on third autocorrelation and 18.0 × 18.0\times faster on sum-difference. ### H Quantum Compilation Tasks This section gives supporting details for the quantum compilation tasks. We separate the material into four parts: problem definitions for the two hardware settings, experimental scaffolds and evaluation setup, discovered programs for the main SimpleTES results, and an additional superconducting routing policy obtained by SimpleTES using Gemini. #### H.1 Problem definitions ##### Superconducting qubit routing. On superconducting quantum computers, this compilation task is typically called qubit routing , and is often discussed together with qubit mapping or qubit allocation [ 128 , 28 ] . Superconducting quantum computers encode qubits in electrical circuits fabricated on a chip, typically using Josephson-junction-based devices operated at cryogenic temperatures [ 70 ] . Because these qubits occupy fixed chip locations, native two-qubit operations are available only between selected pairs of physical qubits. This hardware connectivity is represented by a sparse coupling graph, such as a grid lattice or a heavy-hex layout, whose low degree helps reduce frequency collisions and crosstalk while remaining scalable ( Figure 2 a). This objective is motivated by the fact that SWAPs add extra two-qubit operations and circuit depth; in the evaluation convention used here, one SWAP is decomposed into three added CNOTs [ 76 ] . From a complexity perspective, qubit routing is closely related to token swapping and graph reconfiguration and is NP-hard in general [ 128 , 61 ] . A routing decision that helps one gate may also change the placement seen by all later gates, which makes the long-horizon component essential. ##### Zoned neutral-atom compilation. Neutral-atom platforms store qubits in atoms and use Rydberg interactions to perform entangling gates. A two-qubit gate can be applied only after the participating atoms are brought within the Rydberg interaction radius [ 133 , 137 ] . In a zoned architecture, the static SLM (Spatial Light Modulator) traps are divided into a storage zone and an entangling zone: idle atoms remain in storage, while atoms participating in the current two-qubit layer are moved into the entangling zone [ 15 , 133 , 80 ] . This separation avoids unnecessarily exposing atoms outside the current interaction layer to a global Rydberg pulse. Atom transport is implemented by AODs (Acousto-Optic Deflectors), which consist of independently controlled horizontal and vertical tweezer lines. Each intersection between an AOD row and an AOD column creates a potential movable trap. Open activates new AOD rows or columns and loads atoms into AOD traps when these newly activated lines create intersections that coincide with occupied atomic sites. Close deactivates selected AOD rows or columns and releases atoms currently carried on those lines back to the corresponding SLM sites. Move translates the currently activated AOD lines while carrying all atoms loaded on the AOD. Within one rearrangement step, these operations are subject to two main shuttling constraints. The non-crossing constraint requires active parallel AOD rows or columns to remain separated by at least 2 ​ μ ​ m 2\,\mu\mathrm{m} , so distinct active AOD lines can neither merge into the same position nor exchange their relative order [ 133 , 132 ] . The ghost-spot constraint arises because activating multiple AOD rows and columns simultaneously creates all pairwise row–column intersections as live traps, so a legal rearrangement must avoid unintended auxiliary traps on occupied sites [ 133 , 132 ] ( Figure 2 d). Zoned neutral-atom compilation involves challenging, tightly coupled decisions about atom placement and AOD-constrained routing. Their effects extend across stages: a placement that shortens one transition can still serialize the next once AOD constraints are enforced, and a reuse decision can save transfers in the current stage while forcing longer motion later. Candidate placements must therefore be evaluated over the full placement trajectory rather than one stage at a time [ 133 , 80 , 132 , 131 ] . #### H.2 Experimental scaffolds and evaluation setup ##### Superconducting routing setup. Our experimental scaffold fixes a complete routing engine (circuit parsing, dependency tracking, legality checks, and output construction) and exposes only the decision policy . Concretely, the mutable surface contains two hooks: (i) initial layout construction (a mapping from logical to physical qubits) and (ii) online SWAP selection (choosing a legal SWAP edge at each routing step). This design maximizes the algorithmic search space, allowing the model to implement a broad range of routing strategies. At the same time, it cleanly separates SWAP selection from gate execution, guaranteeing correctness by preventing the model from exploiting invalid gate operations. Initial program. The initial policy is a refactor of Qiskit’s released LightSABRE Rust implementation, embedded inside the fixed engine. The initial layout procedure is simplified to encourage diverse exploration, while the routing policy stays effectively the same as the LightSABRE regime. This baseline provides a strong hand-engineered starting point with a well-studied quality–runtime tradeoff. Evaluation. Candidate policies are compiled and executed on a benchmark suite. Routing cost is measured by the number of inserted SWAPs. We evaluate on the 24 circuits selected from the original SABRE study [ 76 ] , across three coupling graphs: IBM’s 20-qubit Q20 , Google’s 105-qubit Willow processor [ 47 ] , and IBM’s 156-qubit heavy-hex architecture Heron [ 97 , 59 ] , yielding 72 routing instances in total. The baselines are SABRE [ 76 ] and its modern variant LightSABRE [ 176 ] . For LightSABRE, we choose the configuration of swap_trials=20, layout_trials=20, max_iterations=4 from the original paper. ##### Zoned neutral-atom setup. Our scaffold exposes the compiler as executable code inside a fixed framework that defines the interfaces, solver, and evaluator. The compilation pipeline has four components: a scheduler that partitions a logical ZNAACircuit into one-qubit and two-qubit stages, a reuse analyzer that marks qubits reused across adjacent two-qubit stages, a placer that outputs full placement snapshots, and a router that converts consecutive placements into legal Open , Move , and Close operations. A fixed solver invokes these components in order, and the fixed evaluator replays the emitted operations on a ZNAAMachine , and rejects invalid compilations. In our setting, the editable evolve block modifies only the placer ; the scheduler, reuse analyzer, and router remain fixed. Initial program. The initial program exposes the editable evolve block while keeping the surrounding pipeline fixed. In Round 1, that block contained a trivial placeholder placer. Each later round initialized the editable block from the strongest placer discovered in the previous round. Evaluation. We evaluate candidate placers on a 36-circuit suite built from QASMBench [ 75 ] circuits and MQT Bench [ 110 ] generated circuits, spanning textbook quantum algorithms, state preparation, variational circuits, quantum simulation, quantum machine learning, and arithmetic circuits, with sizes from 6 to 500 qubits. The fixed evaluator executes each compiled plan, records correctness, total execution time, and fidelity-related losses from gate execution, transport, transfers, and decoherence, and scores each circuit against a cached baseline from the fixed scaffold. Our main comparison baseline is the reuse-aware ZAC-style compiler [ 80 ] . #### H.3 Discovered programs ##### Discovered superconducting algorithm. The discovered algorithm can be summarized as follows. First, the discovered algorithm invests heavily in initial layout : it seeds high-degree logical qubits onto central, high-degree physical qubits, and then refines the mapping through an aggressive stack of hill-climbing and restart-based local search. Second, it strengthens online SWAP selection : it broadens the candidate neighborhood beyond front-layer incident edges to include look-ahead and shortest-path edges, changes the look-ahead term into dynamic horizon with exponential decay, and reshapes the swap objective to explicitly reward immediate gate executability. Together, these changes preserve the overall LightSABRE-style structure while materially improving robustness against long-range interactions and stagnation. ##### Discovered zoned neutral-atom algorithm. The discovered placer replaces stage-local travel minimization with whole-trajectory routing-aware search, so it explicitly optimizes the parallelism that survives under the fixed router rather than only local travel distance. It begins by constructing a moderate set of heuristic initial storage layouts, including interaction-graph-based qubit orderings, their reversed counterparts, and a few shuffled orders. From each initial layout, it runs a full forward placement pass to build a complete placement trajectory, and then compares these candidates using router-aware evaluation. Within a single forward pass, active qubits are placed stage by stage. When a stage loads qubits from storage into the entangling zone, a qubit whose partner is already reused is assigned to the complementary column of the same entangling site, whereas the remaining two-qubit pairs are assigned by Hungarian matching, with physical Euclidean distance as the placement cost. When active non-reused qubits return from the entangling zone to storage, they are grouped by source entangling row and ordered by source column. A monotone dynamic program then selects destination storage cells that preserve this order while minimizing column deviation, after which an intra-row swap refinement is applied whenever it shortens travel distance and remains compatible with the router’s non-crossing constraint. After a complete trajectory has been constructed, the algorithm evaluates the full placement sequence with the fixed router. Candidates are ranked first by the number of router-emitted Move operations and then, for ties, by total Euclidean travel distance. Starting from the best candidate, the algorithm performs several rounds of reverse-through-time refinement, in which a reverse pass seeded by the current final layout is fed back into a new forward pass, and then applies hill-climbing swaps on the initial storage layout. An update is accepted only when it reduces the number of router-emitted Move operations or preserves that count while improving travel distance. #### H.4 Discovered superconducting routing policy using Gemini The superconducting qubit routing experiment was additionally run with gemini-3-pro-preview . This run produced a policy family in which both the initial placement and the online routing rule are explicitly shaped by circuit structure and hardware geometry. For initial layout, the discovered program classifies logical interaction components by morphology, such as path-like, ring-like, or dense structure. It then selects physical regions with matching shape and connectivity, prioritizes central and well-connected physical qubits for high-traffic logical qubits, and refines the mapping through large-scale randomized restarts and local improvement. For online routing, the program strengthens SWAP selection with criticality-aware scoring. Front-layer and look-ahead gates are weighted by reverse depth and critical-path relevance, while stronger decay and centrality-sensitive preferences steer active qubits toward more useful regions of the device and help avoid stagnation. Quantitatively, the discovered program reduces added SWAPs by 24.7 % 24.7\% relative to SABRE and 18.3 % 18.3\% relative to LightSABRE. The largest gains occur on Q20 , where the reductions in added SWAPs reach 40.7 % 40.7\% and 32.9 % 32.9\% , respectively. Table 5 gives the topology-level breakdown. These results suggest that effective qubit routing benefits from a more explicitly graph-theoretic treatment, in which circuit morphology, interaction criticality, and hardware connectivity are exploited jointly. ### I Astrodynamics trajectory design This section provides the information needed to interpret and reproduce the astrodynamics experiments in the main text. We define the executable-program interface and independent trajectory evaluator, document the benchmark and historical-reference construction, and examine how evolution changes the allocation of search effort. Throughout, “cost” denotes the propulsive cost computed by the benchmark evaluator, not the reported Δ ​ v \Delta v budget of a flown mission. #### I.1 Task formulation and program interface We consider impulsive multi-gravity-assist trajectory design under a patched-conic model [ 148 , 121 ] . Each candidate produced by SimpleTES is an executable optimizer that maps a numerical mission specification to a complete trajectory rather than to a single scalar design vector. The feasible domain Ω s \Omega_{s} depends on the selected discrete design. Changing the encounter sequence or DSM placement can change the dimension and meaning of the continuous variables, whereas the transfer branch selects among distinct transfer solutions for the same variables. These choices remain coupled through the transfer-time and flyby constraints. The reported runs expose the numerical mission specification to the executable program through a structured problem object. Mission names, descriptions, historical trajectory-family labels, reference encounter sequences and event epochs, and the optimized reference trajectories and costs are withheld. The program is instead given the departure body and target endpoint, the broad mission windows, the admissible flyby bodies, the endpoint models, and the feasibility limits listed above. This removes direct historical encounter-sequence information from the task interface, but does not preclude a pretrained language model from inferring dynamically plausible sequences, or recognizing a familiar mission class, from the numerical specification itself. #### I.2 Independent trajectory evaluator Candidate programs are executed in isolated processes and return only the node list in Eq. ( 24 ); they do not return their own validity or score. The evaluator first checks the structural schema and mission limits, and reconstructs planetary states from DE430 ephemerides [ 38 ] . To verify each heliocentric leg, it propagates the declared outgoing state ( 𝐫 i , 𝐯 i + ) (\mathbf{r}_{i},\mathbf{v}_{i}^{+}) from epoch t i t_{i} to the next node’s declared epoch t i + 1 t_{i+1} under solar two-body dynamics, then compares the propagated position and velocity with the next node’s declared incoming state ( 𝐫 i + 1 , 𝐯 i + 1 − ) (\mathbf{r}_{i+1},\mathbf{v}_{i+1}^{-}) . This verification requires neither a Lambert-branch label nor a Lambert solve. The evolved programs may use Lambert solvers internally to construct candidate arcs [ 62 ] ; the evaluator instead tests the returned states directly for propagation consistency. At every GA node, the evaluator recomputes the planetary velocity and applies a powered-flyby model to the incoming and outgoing hyperbolic-excess velocity vectors. A flyby is rejected if the required periapsis radius lies below the configured planetary radius plus minimum altitude. A feasible powered flyby contributes its required periapsis impulse to the cost [ 148 , 121 ] . DSM cost is the norm of the instantaneous velocity discontinuity, Δ ​ v DSM , i = ‖ 𝐯 i + − 𝐯 i − ‖ 2 . \Delta v_{\mathrm{DSM},i}=\left\|\mathbf{v}_{i}^{+}-\mathbf{v}_{i}^{-}\right\|_{2}. (27) Departure and arrival costs are specified independently for each benchmark. For a piecewise-linear boundary, the evaluator computes the velocity discontinuity relative to the departure planet or target state and maps its magnitude through prescribed breakpoints. At launch, the configured quantity C 3 , 0 C_{3,0} defines a zero-cost threshold v ∞ , 0 = C 3 , 0 v_{\infty,0}=\sqrt{C_{3,0}} rather than a hard feasibility boundary. Within the configured linear range, excess departure speed is charged one for one, such that Δ ​ v start = max ⁡ ( 0 , v ∞ − v ∞ , 0 ) . \Delta v_{\mathrm{start}}=\max\!\left(0,v_{\infty}-v_{\infty,0}\right). (28) The same piecewise-linear form is used for flyby or rendezvous endpoints with a prescribed relative-speed allowance. For a planetary-capture endpoint, the evaluator models a tangential periapsis insertion into an elliptic orbit with period T T . If the configured altitude is h factor ​ R p h_{\mathrm{factor}}R_{p} above the planetary reference radius R p R_{p} , the insertion radius is r p = R p ​ ( 1 + h factor ) r_{p}=R_{p}(1+h_{\mathrm{factor}}) . With a = ( μ p ​ T 2 4 ​ π 2 ) 1 / 3 , a=\left(\frac{\mu_{p}T^{2}}{4\pi^{2}}\right)^{1/3}, (29) the capture cost is Δ ​ v cap = v ∞ 2 + 2 ​ μ p r p − 2 ​ μ p r p − μ p a . \Delta v_{\mathrm{cap}}=\sqrt{v_{\infty}^{2}+\frac{2\mu_{p}}{r_{p}}}-\sqrt{\frac{2\mu_{p}}{r_{p}}-\frac{\mu_{p}}{a}}. (30) The four terms in the verified propulsive cost J ⁡ ( 𝒯 ) J(\mathcal{T}) in Eq. ( 26 ) account for departure, DSMs, powered flybys, and arrival, respectively. The number of DSMs is therefore not itself the objective: a DSM is useful only when its added impulse is outweighed by a reduction elsewhere in the verified cost. Valid trajectories are ranked by a monotone transformation of J J , whereas any schema, dynamics, flyby, boundary, timeout, or memory-limit failure receives zero score. We report J J in km ​ s − 1 \mathrm{km}\,\mathrm{s}^{-1} because the transformation does not change the ranking. The evaluator is deliberately lower fidelity than an operational mission design model. It does not attempt to reproduce navigation margins, finite-burn dynamics, spacecraft mass evolution, detailed launch-vehicle performance, or full N N -body perturbations. It nevertheless enforces a common dynamical and cost model across every candidate and reference, including ephemeris consistency, heliocentric propagation, minimum flyby altitude, launch allowance, powered-flyby cost, and target-specific arrival conditions. The reported comparisons should therefore be interpreted as controlled within-evaluator comparisons, not as claims to outperform the flown missions in their operational models. #### I.3 Benchmark instances The five historical mission windows were selected to span distinct trajectory-design structures. Mariner 10 represents an early inner-planet gravity-assist transfer to a flyby target [ 44 ] ; Galileo represents a classical Venus–Earth–Earth-assisted Jupiter transfer [ 32 ] ; Cassini–Huygens combines repeated inner-planet encounters, a Jupiter assist, and Saturn capture [ 106 ] ; Voyager 2 traverses the four giant planets in one continuous tour [ 71 ] ; and Rosetta uses repeated Earth and Mars encounters to rendezvous with a small body [ 46 ] . Tables 7 and 8 summarize the boundary conditions and admissible topology of each benchmark. Planet abbreviations are Me (Mercury), V (Venus), E (Earth), M (Mars), J (Jupiter), S (Saturn), U (Uranus), and N (Neptune). The quantity C 3 , 0 C_{3,0} is the zero-cost launch threshold defined in Eq. ( 28 ). For a planetary boundary, the evaluator obtains the body’s ephemeris state at the candidate’s reported boundary epoch. Thus, the listed arrival windows constrain the epoch rather than specify a stored Cartesian state. The relatively permissive terminal-speed ranges in the Mariner 10 and Voyager 2 benchmarks represent flyby rather than rendezvous conditions. Their endpoint costs remain zero up to 10 and 30 ​ km ​ s − 1 30\,\mathrm{km}\,\mathrm{s}^{-1} , respectively, and increase with unit slope up to the configured admissible limits of 20 and 40 ​ km ​ s − 1 40\,\mathrm{km}\,\mathrm{s}^{-1} . The Rosetta endpoint is the fixed heliocentric state 𝐫 67 ​ P \displaystyle\mathbf{r}_{67\mathrm{P}} = ( 185144048.790 , − 438016487.910 , − 252361615.057 ) ​ km , \displaystyle=(185144048.790,\;-438016487.910,\;-252361615.057)\,\mathrm{km}, 𝐯 67 ​ P \displaystyle\mathbf{v}_{67\mathrm{P}} = ( 7.655449311 , 12.114354862 , 5.612458233 ) ​ km ​ s − 1 , \displaystyle=(7.655449311,\;12.114354862,\;5.612458233)\,\mathrm{km}\,\mathrm{s}^{-1}, at MJD 56875.3791667. This state was obtained from the NASA/JPL Horizons on-line ephemeris system [ 45 ] , using 67P record 90000701 (target solution JPL#K154/5), the Sun as the centre (“500@10”; DE441), and a geometric heliocentric ICRF vector at 6 August 2014, 09:06 TDB. #### I.4 Construction of historical-sequence references The historical comparison does not use a reported mission Δ ​ v \Delta v budget. Reported mission budgets also reflect operational decisions and navigation margins that do not match the benchmark model. Instead, the cited mission chronologies supply the historical encounter orders and nominal event epochs. We reconstruct each reference through fine local optimization around those epochs, then save the resulting node sequence and replay it through the same full evaluator used for SimpleTES . This estimates how the historical trajectory structure performs under the benchmark model; it is not intended as a competing broad mission-design optimizer. The two searches use the same mission boundaries, physical evaluator, and cost definition, but different prior information. Reference construction may use the historical encounter order and nominal epochs because its purpose is to map that trajectory family into the benchmark. The evolved program receives the broad mission windows and allowed flyby bodies, but not the mission identity, historical order or intermediate epochs, or the optimized reference trajectory and cost. This prevents direct transfer of the reference result. It does not establish independence from knowledge encoded in the pretrained language model: the numerical task specification may itself suggest a familiar encounter sequence. #### I.5 How evolution changes mission-design search To examine the search procedure, we replayed the archived initial and evolved programs with lightweight instrumentation. The wrapper records the ordered encounter sequences passed to the program’s internal trajectory evaluations, every Lambert-solver call and its branch parameters, and elapsed execution time, without changing the search logic. The traces are used only to characterize program behavior; the trajectory costs reported in the main text come from independent replay of the checkpoint solutions. The evolved programs do not uniformly maximize search breadth ( Figure 9 ). They expand the set of encounter sequences substantially for Galileo and Cassini, retain a narrow sequence set for Mariner 10, and search fewer encounter sequences than the initial permutation generator for Voyager 2 and the prospective Jupiter case. Their code stages systematic and stochastic sequence generation, prunes or retains candidate sequences according to early costs, and redirects the remaining time toward promising sequences. Evolution therefore selects task-dependent search schedulers rather than simply making every search wider or deeper. The middle panel captures a second discrete decision that remains after an encounter sequence has been selected. Each transfer leg may admit short- and long-path, prograde and retrograde, and multi-revolution Lambert solutions, and the useful combinations depend on the continuous encounter epochs. In each of the five tasks for which the initial program reaches this stage, the evolved program explores more branch types. Cassini nevertheless exhibits the complementary strategy: it evaluates many encounter sequences while using only two Lambert branch types. After this screening, the programs combine global exploration with local refinement of encounter epochs and maneuver variables. This division of effort avoids spending the full budget on either exhaustive sequence enumeration or repeated polishing of one transfer basin. The right panel measures an engineering property of the complete search implementation. All valid evolved replays execute Lambert-solver calls at a higher observed rate than the initial program. The archived programs use caching and candidate pruning to avoid repeated work, while the measured rate also reflects solver choice and other program overhead. A higher throughput allows more transfer candidates to be tested within a wall-clock budget, but is not itself a measure of solution quality; this comparison is not an isolated speed ablation. It instead shows the computational capacity available to the evolved search strategy. Together with adaptive budget allocation, broad branch coverage where useful, and complementary global and local optimization, this capacity allows more effective use of finite computation. The replay describes the behavior of the final executable programs, not the provenance of every heuristic they contain. Several programs include mission-specific sequence templates or incumbent solutions that may have been introduced during evolution or inferred from pretrained knowledge. The figure therefore supports claims about the resulting search policies, but does not by itself establish how independently each encounter-sequence heuristic was discovered. ### J GPU Kernel Optimization ##### Task definitions and benchmark settings. TriMul follows the triangular multiplicative update operator used in AlphaFold3-class protein-structure models [ 1 , 140 ] , instantiated through the official GPUMode TriMul benchmark [ 49 ] . The evaluated benchmark settings are summarized in Table 9 . Asymmetric matrix multiplication follows the KernelBench Level 1 tail skinny matrix multiplication task [ 103 ] , extended with additional highly unbalanced shapes that arise in modern AI workloads. The evaluated benchmark settings are summarized in Table 10 . Batched cumsum follows the KernelBench Level 1 cumsum task [ 103 ] , extended with additional LLM sampling-inspired shapes. The evaluated benchmark settings are summarized in Table 11 . ##### TriMul evaluation alignment. Table 12 reports the TriMul comparison on H100 under Triton 3.4.0 and Triton 3.6.0. We report both versions because TriMul latency is noticeably affected by the Triton compiler version. Triton 3.4.0 is retained for comparability with prior AI-discovery system baselines, including TTT-Discover and K-Search. However, in our local evaluation, Triton 3.6.0 reproduces the public GPUMode leaderboard times more closely. We therefore use Triton 3.6.0 as the leaderboard-aligned setting for the main-text comparison, while reporting Triton 3.4.0 as an additional reference. Each local measurement is repeated three times and reported as mean ± \pm standard deviation. Under both Triton versions, the kernel discovered by SimpleTES achieves the lowest latency among all compared AI methods and public GPUMode Triton submissions. ##### Cross-device performance generalization. Table 13 evaluates whether the TriMul kernel discovered on NVIDIA H200 transfers to other GPU architectures. The same Triton program is evaluated directly on NVIDIA A100, NVIDIA H100, and AMD MI300, without any platform-specific re-tuning or additional search. Despite being discovered only on H200, the kernel achieves lower latency than the compared AI-discovery system baselines and the top-3 GPUMode Triton submissions available for each target GPU. ### K Scaling Law Discovery ##### Overview. Scaling law discovery [ 79 , 78 ] studies how machine learning performance changes with scale and aims to identify compact symbolic laws that extrapolate from small-scale experiments to larger regimes. In foundation-model development, such laws are used to predict quantities such as training loss, downstream error, or task-specific metrics from variables including model size, dataset size, vocabulary size, learning rate, batch size, and architectural choices. The central challenge is not merely to fit a curve, but to recover a concise analytic form whose structure is not known a priori and that generalizes across related experimental settings. Each trial consists of a set of input variables, a target quantity, and a control index identifying the experimental context, such as a model family, dataset, or domain. The control index makes it possible to search for a shared symbolic form across settings while allowing the coefficients to vary within each setting. This makes scaling law discovery especially suitable for SimpleTES : candidate laws can be represented as executable programs, evaluated automatically on held-out extrapolation performance, and compared across a large, structured design space involving symbolic form, asymptotic behavior, parameter sharing, and fitting procedures. Compared with standard regression, the emphasis here is on symbolic structure and extrapolation; compared with broader agentic ML benchmarks, the objective is not to engineer an end-to-end workflow, but to discover a compact law that transfers across settings and remains accurate in unseen scaling regimes. We formalize scaling law discovery as follows. A canonical example is pretraining scaling, where the inputs are model size N N and dataset size D D , the target is training loss L L , and the objective is to discover a law of the form L ≈ f θ ​ ( N , D ) L\approx f_{\theta}(N,D) . The same framework also covers scaling with domain mixture, learning rate, batch size, and related variables. What distinguishes this task from ordinary regression is that success is measured primarily by extrapolation beyond the fitted regime, not by interpolation within it. ##### Experimental setting. We report results on the four-task SLDBench subset: parallel (36 seen / 12 unseen) [ 26 ] , domain_mix (80 seen / 24 unseen) [ 163 ] , lr&bsz (2,702 seen / 117 unseen) [ 77 ] , and u_shape (389 seen / 127 unseen) [ 157 ] . Following the original SLDBench protocol, the unseen split is always constructed as an extrapolation test set by holding out the largest model sizes, compute regimes, or other extreme settings, rather than by using random interpolation-style splits. The execution environment also follows SLDBench: agents operate in a sandbox terminal with a minimal Python stack ( scikit-learn , pandas , and datasets ) and no network access, and must implement the discovered law and its parameter-fitting subroutine under the required function signature. Final performance is measured by test-set R 2 R^{2} , clipped to [ − 1 , 1 ] [-1,1] , on the hidden extrapolation split, where higher is better. We set C = 16 C=16 , L = 20 L=20 , and K = 16 K=16 for SimpleTES on this task. To ensure a controlled comparison, our method uses the same initialization, task instruction, and evaluator as SLDAgent. In particular, the initial program is the same baseline program pair used by SLDAgent: a naive power-law scaling_law_func together with a standard BFGS-based fit_scaling_law optimizer. We also adopt the same task-specific instruction template, which asks the model to evolve both the symbolic expression and the fitting routine from this baseline, while emphasizing extrapolation accuracy, cross-setting generalization, parameter efficiency, and numerical/theoretical stability; the instruction additionally provides task context, function signatures, and data characteristics such as feature definitions and value ranges. The evaluator is likewise kept identical to the SLDAgent setting. This alignment is important because SLDAgent is itself an evolution-based method built on top of the OpenEvolve framework, using iterative mutation and evaluation of candidate programs; accordingly, our comparison isolates the effect of the evolutionary strategy rather than differences in initialization, prompting, or evaluation pipeline. ##### Results analysis. Table 14 shows that scaling law discovery remains a strong discriminator of agent capabilities even on this reduced four-task subset. SimpleTES achieves the best overall average score of 0.674, outperforming the strongest baseline, SLDAgent with GPT-5, which attains 0.572 on the same subset. The tasks also reveal a useful spectrum of difficulty. parallel is close to saturated, with many methods reaching near-perfect extrapolation, while domain_mix is broadly tractable but still differentiates top-performing methods at the margin. By contrast, lr&bsz and especially u_shape remain challenging: many agents still obtain negative test R 2 R^{2} , indicating that symbolic forms that interpolate well in-range often fail to extrapolate reliably into the held-out regime. SimpleTES ’s advantage is most pronounced on these harder extrapolative settings, where it attains the best score on lr&bsz and the best score on u_shape , while also tying for best on parallel and domain_mix . ##### Case study. The lr&bsz task is a particularly informative case study because it tests whether a discovered scaling law can support an actual decision, rather than merely fitting observed losses. The goal is to predict validation loss as a function of learning rate, batch size, dataset size, and model size, and then use the predicted surface to select a good hyperparameter configuration in an extrapolated regime. Because the held-out split consists of larger-scale configurations, success requires recovering the geometry of the loss basin beyond the observed range. Our method discovers an explicit symbolic law for the full loss surface, L ^ SimpleTES ​ ( lr , bsz , D , N ) = \displaystyle\hat{L}_{\text{{SimpleTES} }}(\mathrm{lr},\mathrm{bsz},D,N)={} 6.567185 ​ lr − 0.0131 ​ bsz 0.0096 ​ D 0.0346 ​ N − 0.1173 \displaystyle 6.567185\,\mathrm{lr}^{-0.0131}\,\mathrm{bsz}^{0.0096}\,D^{0.0346}\,N^{-0.1173} + 0.408323 ​ lr 0.2807 ​ bsz − 0.4838 ​ D 0.0499 ​ N 0.0995 \displaystyle+0.408323\,\mathrm{lr}^{0.2807}\,\mathrm{bsz}^{-0.4838}\,D^{0.0499}\,N^{0.0995} + 26.840067 ​ lr − 0.0657 ​ bsz 0.0595 ​ D − 0.2274 ​ N 0.0474 \displaystyle+26.840067\,\mathrm{lr}^{-0.0657}\,\mathrm{bsz}^{0.0595}\,D^{-0.2274}\,N^{0.0474} + 66071022.857 ​ lr 0.0844 ​ bsz 1.8723 ​ D − 1.4189 ​ N − 0.0982 \displaystyle+66071022.857\,\mathrm{lr}^{0.0844}\,\mathrm{bsz}^{1.8723}\,D^{-1.4189}\,N^{-0.0982} + 11.145293 ​ lr 0.9783 ​ bsz 0.1254 ​ D − 0.4142 ​ N 0.4786 \displaystyle+11.145293\,\mathrm{lr}^{0.9783}\,\mathrm{bsz}^{0.1254}\,D^{-0.4142}\,N^{0.4786} − 0.1294720520 . \displaystyle-0.1294720520. Given a target regime ( D , N ) (D,N) , we evaluate this law on the full admissible hyperparameter grid 𝒢 \mathcal{G} and choose the point with the minimum predicted loss: ( lr † , bsz † ) = arg ⁡ min ( lr , bsz ) ∈ 𝒢 ​ L ^ SimpleTES ​ ( lr , bsz , D , N ) . (\mathrm{lr}^{\dagger},\mathrm{bsz}^{\dagger})=\arg\min_{(\mathrm{lr},\mathrm{bsz})\in\mathcal{G}}\hat{L}_{\text{{SimpleTES} }}(\mathrm{lr},\mathrm{bsz},D,N). This matches the actual deployment setting: we score every feasible grid point directly, rather than optimizing in a continuous space and then rounding back to the nearest grid point. For the extrapolation example in Figure 5 e, corresponding to a 1 1 B-parameter model trained on 100 100 B tokens, this procedure selects the blue-star configuration ( 1.953 × 10 − 3 , 512 ) (1.953\times 10^{-3},512) . Its realized validation loss is 2.0774 2.0774 , while the true best grid point is the red-star configuration ( 1.381 × 10 − 3 , 512 ) (1.381\times 10^{-3},512) with loss 2.0762 2.0762 . Thus, the configuration selected by SimpleTES is only 0.058 % 0.058\% above the optimum. For comparison, SLDAgent selects the gray-star point ( 1.953 × 10 − 3 , 384 ) (1.953\times 10^{-3},384) , whose realized loss is 2.0776 2.0776 , or 0.067 % 0.067\% above the optimum. Although the numerical gap is small, SimpleTES more accurately identifies the optimal region of the extrapolated loss surface and lands slightly closer to the empirical optimum. This comparison is meaningful because the three marked points all lie in the same narrow low-loss valley, so the main challenge is fine-grained recovery of the learning-rate/batch-size trade-off rather than coarse localization. Modeling the full loss surface, rather than only the optimum coordinates, provides richer supervision and makes the discovered law directly usable as a practical selector in unseen regimes. More broadly, this example shows why high held-out R 2 R^{2} matters in practice: the value of a scaling law lies not only in predictive accuracy, but also in its ability to support extrapolative decisions. On the lr&bsz task, SimpleTES turns a symbolic loss law into a grid-level hyperparameter rule and selects a point that is essentially optimal on the true evaluation landscape. ### L Whole-brain neural-activity forecast on ZAPBench ##### Problem definition and benchmark design. ZAPBench evaluates neural-activity forecasting at whole-brain, single-neuron resolution. Structural maps of neural wiring are essential for understanding circuit organization, but they do not by themselves specify the time-varying circuit dynamics that transform sensory input and internal state into neural responses. A functional account therefore requires predicting how activity evolves across the brain over time. ZAPBench provides a quantitative benchmark for this problem by asking models to forecast the activities of all recorded neurons from a short temporal context, enabling direct comparison of forecasting models and analysis of how whole-brain dynamics unfold across stimulus conditions. It is built from a two-hour light-sheet fluorescence recording of a larval zebrafish that was exposed to a sequence of visual stimuli designed to elicit a range of behaviors; the raw volumetric video was aligned, motion-stabilized, and segmented into activity traces, yielding a data matrix of 71,721 71{,}721 neurons over 7,879 7{,}879 time steps (a sampling interval of roughly 914 914 ms per step) spanning nine stimulus conditions ( gain , dots , flash , taxis , turning , position , open loop , rotation , and dark ). ##### Forecasting task and data splits. The task is to predict the joint future activity of all 71,721 71{,}721 neurons up to 32 32 steps ahead given only 4 4 steps of past context—a high-dimensional multivariate forecasting problem. Because the traces are derived from the volumetric recording, a forecaster may operate either on the extracted per-neuron traces or on the 3D volumes directly: the benchmark’s strongest hand-designed baseline is a volumetric U-Net that predicts the brain video, exploiting spatial structure that is lost when the volumes are reduced to traces, but at a much higher computational cost—the video model is reported to train in roughly 36 36 hours on 16 16 A100 GPUs, versus under four hours on a single 16GB GPU for trace-based forecasters. Following ZAPBench, each stimulus condition is split 70 / 10 / 20 70/10/20 into contiguous train, validation, and test segments. ##### Evaluation and search protocol. Candidate programs are scored by the mean absolute error (MAE) between predicted and recorded activity, averaged over all neurons and over the full 32 32 -step (or 1 1 , 4 4 , 8 8 , 16 16 steps in different tasks) horizon. The validation split serves two purposes: it drives model selection including hyperparameter choices and early stopping; and it supplies the single scalar (validation MAE) that guides the search. The test split is never seen during search or model selection. SimpleTES is initialized from one neural starter program, a shared-weight per-neuron forecaster (Fig. 4 f), and edits the body of its forecast routine. Each candidate is trained and scored in an isolated evaluator. ### M Single-Cell RNA-Seq Denoising Single-cell RNA sequencing (scRNA-seq) resolves gene expression in individual cells [ 89 , 173 , 87 ] , but its measurements are corrupted by low capture efficiency and stochastic dropout [ 146 ] , motivating dedicated denoising algorithms [ 147 ] . Denoising has no ground truth to check against; quality is assessed by the molecular cross-validation framework of [ 11 ] , which partitions observed UMI counts by binomial subsampling and underpins the OpenProblems benchmark [ 86 ] . The core difficulty is a tension between two objectives: mean-squared error in log-normalized space rewards faithful recovery of relative expression, while the Poisson negative log-likelihood rewards count consistency after library-size rescaling, and aggressively optimizing one degrades the other. We follow the molecular cross-validation protocol [ 11 ] , splitting the Pancreas dataset into a training matrix the algorithm observes and a held-out test matrix it does not. Candidates are scored on two normalized metrics, log-space MSE and Poisson negative log-likelihood, with the Poisson term acting as a hard constraint that rejects any solution failing it, and the final score is the mean of the two. Following [ 165 ] , we initialize SimpleTES with MAGIC [ 147 ] , cap each candidate at 400 seconds, and run evolution exclusively on Pancreas; PBMC and Tabula Muris Senis Lung are withheld entirely and used only for final evaluation, so the reported numbers measure cross-dataset generalization to held-out datasets. Across the held-out datasets in Table 15 , SimpleTES reaches a Tabula score of 0.74 0.74 , surpassing the prior state-of-the-art method TTT-Discover [ 165 ] ( 0.73 0.73 ) and matching it on PBMC ( 0.71 0.71 ), while outperforming MAGIC [ 147 ] , ALRA [ 81 ] , and every other baseline on both. Because evolution never saw PBMC or Tabula, these gains are direct evidence that the discovered algorithm generalizes rather than fitting the search distribution. Inspection shows that SimpleTES departs structurally from the baselines. Where MAGIC builds a single diffusion operator and TTT-Discover elaborates that template with variance-stabilizing transform ensembling and low-rank refinement, the discovered program constructs several independent denoisers, including multi-scale diffusion under both correlation- and Euclidean-based graph operators, log-space diffusion, neighbor averaging, PCA imputation, and NMF reconstruction, and scores each internally on the training split alone. It then blends them into a data-driven ensemble, weighting candidates by their inverse Poisson and MSE losses under the Poisson constraint, and applies a gene-mean calibration that clips per-gene scaling near 1.0 1.0 to preserve count fidelity. None of these choices are dataset-specific: the ensemble weights are fit from the training split at runtime, which is what lets the same program generalize to the unseen PBMC and Tabula distributions. ### N Lasso Regularization Path ##### Overview. The lasso regularization path is a core computational primitive in high-dimensional statistics, arising naturally in cross-validation and model selection across domains from genomics to finance. Computing the full path of solutions across a grid of regularization values is orders of magnitude more efficient than solving each problem independently, but demands careful algorithmic design to exploit warm starts and sparsity structure. The de facto standard solver, glmnet [ 39 ] , represents decades of expert engineering. This task asks whether SimpleTES can discover a solver that is faster than glmnet while maintaining the same float64 precision and correctness guarantees. ##### Experimental setting. Each candidate is a self-contained C++ program that reads a binary-encoded problem from standard input and writes the complete coefficient path to standard output. Correctness is evaluated on a separate fresh problem that is not used for timing, preventing caching or test-set overfitting. The surrogate metric optimized during evolution is the geometric mean of solve times over 17 synthetic problem sizes. These instances span the axes that govern lasso-path performance: sample-to-feature ratio n / p n/p , design-matrix sparsity, active-set density, and feature correlation from near-independent to highly correlated Toeplitz structure. This evaluation surface is deliberately heterogeneous: a solver specialized for wide problems with p ≫ n p\gg n can fail on tall problems with n ≫ p n\gg p , and vice versa. The 11 real-world datasets in Table 16 are not used during evolution and therefore test whether the discovered solver generalizes beyond the synthetic search distribution. We initialize SimpleTES from a faithful C++ port of glmnet ’s Gaussian lasso path, including the covariance update for p < 500 p<500 , the naive residual update for p ≥ 500 p\geq 500 , warm starts, sequential strong-rule screening [ 142 ] , active-set inner loops, and KKT verification. ##### Results analysis. All reported solutions pass the held-out correctness check. Across the 11 real-world datasets, SimpleTES achieves an average speedup of 2.17 × 2.17\times over glmnet and 14.08 × 14.08\times over sklearn . The gains are largest in tall regimes: on DNA ( 1700 × 180 1700\times 180 ), SimpleTES is 9.56 × 9.56\times faster than glmnet . It also improves on the large sparse RCV1 design by 1.76 × 1.76\times over glmnet and 5.49 × 5.49\times over sklearn ; the gains against sklearn are particularly large on high-dimensional biological datasets, including 50.77 × 50.77\times on TCGA BRCA and 10.16 × 10.16\times on Duke Breast Cancer. ##### Case study. The initialization uses coordinate descent throughout: the covariance method for p < 500 p<500 and the naive residual method for p ≥ 500 p\geq 500 . The discovered program instead uses a geometry-aware switching rule. For moderate-dimensional, tall designs ( p ≤ 2000 p\leq 2000 and n ≥ p / 4 n\geq p/4 ), it replaces coordinate descent with an exact LARS homotopy solver that traces the regularization path analytically and updates the active set at each kink with rank-one inverse-Gram updates [ 34 ] . For wide or sparse designs, it retains coordinate descent with strong-rule screening, an active-set inner loop, and KKT verification. ### O Mathematical Scientific Discovery Tasks Extremal analysis and combinatorial construction ask for explicit mathematical solutions that optimize precisely defined quantities: a distribution, a witness function, an integer set, a geometric packing, or a sign matrix. These tasks are a natural setting for scientific discovery because progress depends on finding the right structural solution, and each candidate can be checked by a deterministic evaluator rather than by a noisy empirical proxy. #### O.1 Extremal analysis Extremal analysis tasks search over constrained functions and optimize overlap or autoconvolution functionals. Candidate solutions are represented as discretized step functions or signed piecewise-constant witnesses, with feasibility enforced by support, range, and unit-mass constraints. ##### Erdős minimum-overlap problem. Erdős’ minimum-overlap problem asks how to distribute a unit amount of mass over [ 0 , 2 ] [0,2] so that the distribution has as little one-sided overlap as possible with the complements of its shifted copies. The evaluator reports the translated-overlap objective through 1 / Ψ ⁡ ( h ) 1/\Psi(h) during search, while the standard value Ψ ⁡ ( h ) \Psi(h) is reported for comparison. Best result and discovered solution. The discovered witness is near-binary: most grid cells are driven close to either zero or one, rather than remaining near the uniform initialization. This structure arranges the unit mass so that the translated-overlap profile is flattened across the high-risk shift range, reducing the worst shift rather than only improving the average shift behavior. The solution is therefore interpretable as an equioscillation-style witness: the best solution balances multiple competing translated-overlap constraints instead of optimizing a single shift locally ( Figure 10 a). Case Analysis. The strongest program uses a coarse-to-fine optimization pipeline. It first searches over coarse discretizations to identify promising mass layouts, then refines them with constrained local optimization and projected polishing steps that preserve the exact discrete mass constraint and the box constraint 0 ≤ h ≤ 1 0\leq h\leq 1 . This code structure is important because naive local perturbations around the uniform witness tend to improve only a small subset of shifts, whereas coarse-to-fine refinement can reorganize the global mass pattern before polishing the final overlap profile. ##### Autocorrelation Inequalities. The First, Second, and Third Autocorrelation Inequalities form a family of extremal autoconvolution problems. During search, the evaluators cast the tasks as maximization by using 1 / Φ 1 1/\Phi_{1} , Φ 2 \Phi_{2} , and 1 / Φ 3 1/\Phi_{3} , while the standard objectives are reported in the result tables. Best results and discovered solutions. The three autocorrelation tasks lead to different witness structures. For AC1, the strongest discovered witness concentrates mass near the support boundary, reducing the peak of the non-negative autoconvolution but still trailing the leading public result. For AC2, the best witness is sparse and produces a long near-flat plateau in the autoconvolution profile, a structure that improves the norm ratio by avoiding a narrow peak. For AC3, the signed witness becomes oscillatory and uses cancellation to control the largest absolute autoconvolution value. These solutions show that the same evaluation-driven loop can search across non-negative and signed witness classes without a task-specific symbolic derivation ( Figure 10 b–d). Case Analysis. The best programs implement different search strategies for the three objectives. AC1 uses a simplex-projected mass-transfer search over a 1024-bin witness, repeatedly moving mass while preserving non-negativity and unit mass. AC2 combines FFT-based convolution evaluation with L-BFGS-B refinement, making it efficient to evaluate many candidate witnesses and polish sparse structures. AC3 uses a discrete cosine transform parameterization that directly searches signed oscillatory patterns and optimizes the cancellation structure of f ∗ f f\ast f . The code-level adaptation is therefore objective-specific even though the outer SimpleTES procedure is unchanged. #### O.2 Combinatorial construction Combinatorial construction tasks search over finite structured solutions: integer sets, unequal-circle packings, and { − 1 , 1 } \{-1,1\} sign matrices. Their evaluators are exact or feasibility-based, using integer arithmetic for sumsets, geometric non-overlap checks for packings, and exact determinant computation for sign matrices. ##### Sum–Difference Problem. The Sum–Difference Problem comes from classical additive combinatorics and is closely related to more-sums-than-differences sets and Ruzsa-type inequalities [ 92 , 54 , 116 , 117 ] . Candidate sets contain at most 512 integers with elements bounded in [ − 10 6 , 10 6 ] [-10^{6},10^{6}] . The evaluator computes | A + A | |A+A| and | A − A | |A-A| exactly using integer arithmetic. Best results and discovered solutions. The discovered set has a highly regular arithmetic backbone rather than an irregular collection of integers. Most consecutive gaps follow a long progression pattern, while a small number of fringe positions receive local corrections. These sparse perturbations enlarge the sumset more efficiently than the difference set, producing a construction that is both high-scoring and structurally interpretable ( Figure 10 e). Case Analysis. The strongest program maintains exact sum and difference multiplicity tables during search. It alternates aggressive pruning, greedy additions, and local replacements, allowing each candidate edit to be scored by its exact effect on the two set expansions. This is more efficient than repeatedly rebuilding A + A A+A and A − A A-A from scratch, and it makes the search sensitive to local fringe edits that would be difficult to identify from the final scalar score alone. ##### Circle Packing in a Unit Square. Unequal-circle packing in a square is a classical problem in continuous optimization and discrete geometry [ 105 , 136 , 55 ] . Invalid configurations receive zero score. Best results and discovered solutions. The n = 26 n=26 and n = 32 n=32 solutions recover the strongest known packing values under the same evaluator, but the two instances have visibly different structure. The n = 26 n=26 construction contains a dominant central circle, large boundary anchors, and a near-triangular ring of medium circles. The n = 32 n=32 construction is more homogeneous and resembles a fixed six-row quasi-hexagonal layout. These structures show that the search does not only adjust radii locally; it recovers global geometric organization suitable for each instance size ( Figure 10 g,h). Case Analysis. The evolved n = 26 n=26 solver uses a coarse-to-fine strategy: broad center-set exploration first identifies promising layouts, and exact radius optimization then solves the radius subproblem under boundary and pairwise-distance constraints. The n = 32 n=32 solver simplifies to a more direct routine by reusing an incumbent or fixed quasi-hexagonal layout, then solving for radii through linear programming. The two programs therefore adapt the search granularity to the instance: broader geometric exploration for n = 26 n=26 and more efficient radius refinement for n = 32 n=32 . ##### Hadamard maximum-determinant problem of order 29 29 . The Hadamard maximum-determinant problem of order 29 29 is a classical benchmark in extremal matrix theory and D D -optimal design. Order 29 lies outside the Hadamard regime, so the best known solutions are near-extremal sign matrices rather than exact Hadamard matrices. The evaluator computes the determinant exactly using the Bareiss algorithm and reports both the raw determinant and a normalized score. Best results and discovered solutions. The discovered sign matrices recover the long-standing classical lower-bound record. Their structure is best inspected through the Gram matrix rather than raw signs, since row and column permutations and sign flips can change the visible matrix while preserving equivalence. The high-determinant matrices exhibit near-orthogonal Gram structure, with off-diagonal interactions controlled so that rows remain as close to mutually orthogonal as possible under the order-29 constraint ( Figure 10 f). Case Analysis. The strongest solver expands the initial search beyond a single quadratic-residue construction. It builds a diverse seed pool containing incumbent, quadratic-residue circulant, orthogonal-sign, and cropped Sylvester-style matrices, then alternates inverse-guided multi-flip hill climbing with simulated annealing. This combination allows the program to make coordinated sign changes rather than relying on isolated single-entry flips, which is essential for escaping low-determinant basins while preserving the global Gram structure needed for a near-extremal matrix. ## Instructions for reporting errors We are continuing to improve HTML versions of papers, and your feedback helps enhance accessibility and mobile support. To report errors in the HTML that will help us improve conversion and rendering, choose any of the methods listed below: Click the "Report Issue" ( ) button, located in the page header. Tip: You can select the relevant text first, to include it in your report. Our team has already identified the following issues . We appreciate your time reviewing and reporting rendering errors we may not have found yet. Your efforts will help us improve the HTML versions for all readers, because disability should not be a barrier to accessing research. Thank you for your continued support in championing open access for all. Have a free development cycle? Help support accessibility at arXiv! Our collaborators at LaTeXML maintain a list of packages that need conversion , and welcome developer contributions .