{"id":"630a5b87-b183-49f9-a46e-083793d8fa4a","arxiv_id":"2505.04681","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":3,"one_line_summary":"Symbolic regression on FIRE-2 simulations yields analytic star formation rate equations involving gas surface density, stellar surface density, and gas velocity dispersion that outperform classical Kennicutt-Schmidt laws on test data.","lead":"This paper applies symbolic regression to the FIRE-2 galaxy simulations to discover compact equations that predict star formation rates from local gas and stellar properties. The resulting equations beat standard analytic star formation laws on held-out test data and may provide new subgrid recipes for galaxy simulations.","discovery_kind":"new_application","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Test-set outperformance may be inflated by random pixel-level splitting: training and test pixels share galaxies/snapshots with strong spatial correlations, so the reported R^2 gains over analytic models do not yet demonstrate generalization; leave-one-galaxy-out cross-validation is needed.","rationale":"The most load-bearing part of the paper is the test-set evaluation: all conclusions about outperformance and physical convergence rest on R^2 and sigma_plane from a 20% pixel-level test split. Because the split is random within the same seven galaxies and their snapshots, the test pixels are not independent of the training pixels; the interstellar medium is spatially correlated, and separate snapshots of the same galaxy are dynamically related. Symbolic regression is a flexible procedure that can fit autocorrelated structure and then appear accurate on a correlated test set. The model-selection step, which selects equations with training loss below the analytic baseline, also means the reported test R^2 does not validate the chosen equation form against independent data. Thus the headline claim that PySR equations outperform all analytic models on FIRE-2 is not yet established for unseen galaxies or independent ISM patches. The 100 Myr convergence to log Sigma_SFR = log Sigma_gas + (log sigma_gas,z + log Sigma_star)^0.86 - 6.1 could be a real equilibrium relation, but it could also be a selection artifact on one correlated sample; holding out whole galaxies is the decisive check. I do not see this as a fatal flaw: the pipeline is clearly described, the synthetic test in Appendix D shows the machinery can recover a known relation and ignore a noise feature, and the authors are appropriately cautious about stochasticity at 10 Myr. Still, the test-set independence issue is the single most load-bearing weakness. The reader's weakest assumption, about FIRE-2 fidelity and expressibility of the local variables, is related but not identical; the decisive issue is not only what physics is in the simulations but whether the reported test metric measures generalization at all. My recommended action is therefore to keep the conditional acceptance, with the added requirement that galaxy-level cross-validation be reported.","tokens_in":31670,"tokens_out":8497,"duration_ms":89790,"concrete_test":"Leave-one-galaxy-out test: train the full pipeline (XGBoost+SHAP feature selection, PySR search, and selection rule) on pixels from six FIRE-2 galaxies, then evaluate the selected equations on all pixels of the held-out galaxy; repeat for each of the seven galaxies. Calibrate the analytic baselines on the training fold only, not the test fold. Record test R^2, sigma_plane, and the form of the best equation for both 10 Myr and 100 Myr targets. If the mean held-out R^2 for the 100 Myr equation is close to the random-split value (about 0.55) and remains above the best analytic baseline (about 0.42), the outperformance survives; if it drops toward or below baseline, or if the equation form changes strongly with the held-out galaxy, the central claim of a robust, universal top-down relation is not supported.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim rests on the test-set comparison in Sections 3.1-3.2: PySR equations are said to outperform six analytic/XGBoost baselines and to converge to a single 100 Myr scaling relation. But the train/test split in Section 2.1.1 is made pixel-by-pixel over the seven FIRE-2 galaxies and their snapshots, not by galaxy or snapshot. Pixels in the test set share the same galaxy, snapshot, and local ISM environment with training pixels; star-forming regions are strongly spatially autocorrelated on scales well above the 750 pc pixel size. The fitted equations, and especially the model-selection step that chooses equations with training loss below the analytic baseline, can therefore exploit information from the very structures that also appear in the test set. The reported R^2 gains are small for 10 Myr (0.413 vs 0.402 over SAF20) and about 0.13 for 100 Myr (0.555 vs 0.420 over FG13), but without cluster-level out-of-sample evaluation there is no way to know whether these gains reflect physical generalization or in-sample spatial correlation. Separately, the analytic baselines have free parameters fit to the test set (Section 3.1), which complicates the comparison even though it does not favor the SR models. The convergence of the 100 Myr equations to a single form may be a real equilibrium relation, but it could also be a selection artifact on one correlated sample; a held-out-galaxy test is the decisive check.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"This manuscript applies the symbolic regression package PySR to 750 pc pixel maps from seven FIRE-2 Latte galaxies, aiming to discover closed-form expressions for the star formation rate surface density averaged over 10 Myr and 100 Myr. The input features are reduced from 76 candidates to eight via XGBoost and SHAP analysis supplemented by physical arguments; the outputs are evaluated on a held-out 20% test set and compared with XGBoost, the Kennicutt-Schmidt K98 relation, three bottom-up models (KDM12, SFK15, SAF20), and one top-down model (FG13). The paper reports that the found equations outperform all analytic models on the test set, that the 100 Myr equations converge to a scaling relation of the form log Sigma_SFR ~ log Sigma_gas + (log sigma_gas,z + log Sigma_*)^0.86 - 6.1, and interprets this as support for a top-down, feedback-regulated origin of the Kennicutt-Schmidt relation. An appendix presents a synthetic recovery test based on FG13 as validation of the pipeline.","tokens_in":31975,"tokens_out":5886,"duration_ms":58527,"significance":"If the central claims hold, this is a useful contribution: it provides interpretable, data-driven parameterizations of star formation from a modern high-resolution simulation suite, with potential applications as subgrid recipes in large-volume simulations and semi-analytic models. The pipeline itself is transferable to other datasets. The manuscript has real strengths: the use of a held-out test set, a synthetic recovery test in Appendix D demonstrating that the method can rediscover a known relation while ignoring a pure-noise variable, an explicit complexity/loss trade-off analysis, and a candid discussion of caveats in Section 4.1. However, the significance of the headline claims depends on whether the reported test-set gains reflect genuine generalization; the current validation strategy does not yet establish that.","major_comments":[{"comment":"The central claim that the PySR equations \"statistically outperform\" the analytic models rests on a random pixel-level 80/20 split, but pixels drawn from the same galaxy and snapshot are strongly spatially correlated at 750 pc. The test set is therefore not independent of the training set, and the modest margins (10 Myr: 0.413 vs 0.402 for SAF20; 100 Myr: 0.555 vs 0.420 for FG13) could be inflated by this leakage. The synthetic FG13 recovery test in Appendix D is encouraging but does not address this issue, since its synthetic data inherit the same pixel-level structure. I request a leave-one-galaxy-out or leave-one-snapshot-out cross-validation, with per-fold R^2, sigma_plane, and found equations reported; without it, the generalization claim in the abstract and Section 5 is not established.","section":"Section 2.1.1 and Table 3"},{"comment":"The model-selection step is not fully out-of-sample. In Section 3.2 the baseline loss is computed on the test set, while the selection criterion applied to PySR candidates is the final training loss; this leaks test information into the choice of the top four equations. In addition, Section 3.1 states that the analytic baselines have their normalizations or efficiencies fitted to the test set. Fitting the baselines to the test set is conservative in that it favors the analytic models, but using the test set to set the selection threshold weakens the claim that the reported test-set R^2 values are a clean out-of-sample comparison. Please use a separate validation set, or an internal cross-validation loop, for both the baseline loss and the selection of equations, and reserve the test set for the final comparison only.","section":"Sections 3.1 and 3.2"},{"comment":"No uncertainties are reported for R^2, sigma_plane, or the fitted exponents and constants. At the 10 Myr timescale the difference between the best PySR equation and SAF20 is 0.413 vs 0.402, which is too small to interpret without error bars. Please add resampling-based uncertainties (for example, bootstrap over pixels or over snapshots) or repeated PySR runs with different random seeds, and report the scatter of the coefficients across those runs. This is needed to support the quantitative claim of statistical outperformance over the analytic models.","section":"Tables 3 and 4"}],"minor_comments":[{"comment":"The equations use \"log\" without specifying the base, while expressions such as e^{log f_gas} imply a natural logarithm and the text states training was performed in log10 space; please define the convention and verify that all equations are internally consistent.","section":"Table 4 and Figures 6-9"},{"comment":"Equation 5 uses alpha_i but the surrounding text only lists quantile values 0.1, 0.25, 0.5, 0.75, and 0.9; please state explicitly that alpha_i takes these values and that N_alpha = 5.","section":"Equation 5"},{"comment":"The caption refers to \"Equation 2.1.2\" for the quantile loss, which should be Equation 5.","section":"Figure 4 caption"},{"comment":"The pixel selection criteria are described only qualitatively as \"both recent star formation and gas,\" with a Sigma_gas cut shown in Figure 3; please give the exact thresholds and the final pixel counts after selection.","section":"Section 2.2 and Figure 3"},{"comment":"There are typographical errors, including \"for that for that\" in Section 2.2 and \"casual\" for \"causal\" in Section 4.1.2; a careful proofread is needed.","section":"Throughout"}],"recommendation":"major_revision","confidential_remarks":"The manuscript presents a genuinely useful pipeline and an honest discussion of limitations, and the synthetic recovery test is a strong piece of validation. The main obstacle is the validation strategy: the pixel-level split and test-informed model selection undermine the headline generalization claim. I believe this is fixable with galaxy/snapshot-level cross-validation and a clean validation/test separation, which is why I recommend major revision rather than rejection. I would also encourage the authors to add uncertainties to the R^2 and to the fitted coefficients, since the 10 Myr comparison is close."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The headline: this is a legitimate first attempt at symbolic regression for star formation parameterization, and the 100 Myr result is a tidy scaling relation that looks physically sensible. The paper is not a flashy claim; it's an honest, careful application of a known ML tool to FIRE-2, with a synthetic recovery test that actually gives the method some credibility. But the central quantitative claim — that the discovered equations beat analytic models on a held-out test set — is more optimistic than it reads, because the train/test split is pixel-level rather than by galaxy or snapshot. Strong spatial autocorrelation means the test set shares structures with training data, so the R² gains don't yet demonstrate generalization to new galaxies.\n\nWhat's genuinely new: PySR has been applied to cosmology and black hole scaling before, but not to spatially resolved star formation in galaxy simulations. The custom loss combining MSE and quantile terms is a reasonable choice, and the model-selection criterion — require training loss below the best analytic baseline — is clever and physically motivated. The synthetic FG13 test (Appendix D) is a real strength: it shows the pipeline can recover a known relation and reject a noise variable. The 100 Myr convergence across complexities is striking and consistent with top-down/feedback-regulated ideas. The authors also do themselves credit by listing caveats: poor stochasticity at 10 Myr, coarse 750 pc resolution, and the minor contribution of the quantile loss.\n\nWhere I'd push back: the pixel-level split is the load-bearing weakness. Leave-one-galaxy-out (or at least snapshot-based splitting) would be decisive. The R² values are modest (0.41 and 0.55) with no uncertainties, so the outperformance over SAF20 and FG13 is small relative to likely noise. They also exclude non-star-forming pixels and apply a Σgas cut, which is fine for a subgrid model but not for a universal law. The analytic baselines are fitted to the test set, which actually favors them, so the fact that SR still wins is some evidence, but it doesn't fix the lack of error bars.\n\nWho's it for: anyone thinking about subgrid SF prescriptions in simulations or SAMs, and anyone using symbolic regression in astrophysics. It deserves a serious peer review, not a desk reject, but I'd send it back for a cluster-level validation before acceptance. If that holds, this becomes a citable recipe; if not, the paper still works as a method demonstration.","headline":"A well-constructed symbolic regression study that recovers a top-down-like SF law at 100 Myr, but the test-set validation is more optimistic than it looks because the split is pixel-level.","tokens_in":32519,"tokens_out":3474,"would_cite":true,"duration_ms":35856,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"Symbolic regression on kiloparsec galaxy maps finds star formation equations that beat all tested analytic models, with the 100 Myr fits converging to a single scaling relation.","keywords":["star formation","symbolic regression","Kennicutt-Schmidt relation","FIRE-2 simulations","interstellar medium","galaxy evolution","scaling relations","machine learning"],"falsifier":"Take the converged 100 Myr equation and apply it to 750 pc pixels of observed galaxies, or of an independent simulation suite with a different feedback prescription, and compare predicted versus measured $\\Sigma_{\\rm SFR}$; if the relation fails on galaxies with different star formation histories while still fitting FIRE-2, the convergence is a property of the training data rather than a universal star formation law.","tokens_in":31477,"feed_emoji":"🌌","tokens_out":12693,"duration_ms":105979,"temperature":0.7,"pith_summary":"This paper sets out to replace the hand-assumed star formation recipes used in galaxy simulations with equations that are discovered from data. It applies symbolic regression to 750 pc pixel maps of seven simulated Milky Way-mass galaxies and searches for the formula that best predicts the star formation rate surface density, $\\Sigma_{\\rm SFR}$, from eight local variables. The discovered equations beat every analytic star formation model the authors compare against, and at the 100 Myr averaging timescale the top equations all converge to nearly the same scaling relation in $\\Sigma_{\\rm gas}$, $\\Sigma_*$, and $\\sigma_{\\rm gas,z}$. The authors read this convergence as evidence that the Kennicutt–Schmidt relation is governed top-down, by feedback-regulated disk equilibrium, and that replacing empirical subgrid recipes with such an equation could improve large-volume galaxy simulations and semi-analytic models.","feed_headline":"Symbolic regression beats analytic star formation models","feed_subtitle":"The 100 Myr fits collapse into one scaling relation that links gas, stars, and turbulence.","key_machinery":"The central machinery is symbolic regression, a genetic algorithm that searches over equations by representing them as expression trees and evolving them through mutation and crossover, returning a Pareto front that trades equation complexity against fit. The paper runs the search in logarithmic space, using a loss that combines mean squared error with a quantile loss meant to preserve scatter, and selects only equations whose training loss is lower than the best analytic star formation model. The load-bearing output is the converged 100 Myr scaling relation, whose persistent structure — the combination of $\\log \\Sigma_{\\rm gas}$ with $(\\log \\sigma_{\\rm gas,z} + \\log \\Sigma_*)^{0.86}$ — is the object that carries the paper's physical interpretation.","core_discovery":"The central claim is that a closed-form, data-driven description of star formation exists at kiloparsec scales and that the search finds it. On the 100 Myr averaging timescale, every selected equation converges to $$\\log \\Sigma_{\\rm SFR} \\approx \\log \\Sigma_{\\rm gas} + (\\log \\sigma_{\\rm gas,z} + \\log \\Sigma_*)^{0.86} - 6.1,$$ with only small variations in the exponents and constant across the selected equations. On the 10 Myr timescale, the equations are less stable and capture less variance, which the paper attributes to stochastic feedback that cannot be represented as a deterministic function of instantaneous local variables. Because $\\Sigma_*$ and $\\sigma_{\\rm gas,z}$ appear alongside $\\Sigma_{\\rm gas}$ in the converged relation, the paper interprets the result as favoring 'top-down' feedback-regulated models over 'bottom-up' local-cloud models of the Kennicutt–Schmidt relation, and notes that the same variables appear in the dynamical-equilibrium pressure used by those theories.","pith_inferences":["A natural test the paper leaves open is whether the same converged 100 Myr equation reappears when the pipeline is trained on an independent hydrodynamical simulation with different feedback physics; if it does, the relation is likely a property of real star-forming disks rather than of FIRE-2's subgrid model.","The authors note the quantile-loss term contributed only 10 to 15 percent of the total loss; a distribution-matching loss such as KL divergence is an obvious extension that could fix the observed underestimation of scatter in the feature-space planes.","The paper's own resolution caveat points to a sharp experiment: rerun on 200 pc or finer maps; if the equation's form shifts toward the bottom-up variables (free-fall time, Mach number, virial parameter), the claimed top-down dominance is resolution-dependent.","Because the relation includes the stellar surface density as a proxy for disk potential, it may also connect to observed dynamical-equilibrium pressure correlations in nearby galaxy surveys; testing on those maps would tell whether the equation extrapolates beyond simulated galaxies."],"forward_implications":["The 100 Myr equation can be used directly as a subgrid star formation recipe in large-volume cosmological simulations, replacing the empirical Kennicutt–Schmidt power law, and on FIRE-2-like test pixels it should outperform all analytic recipes examined here.","If the convergence is physical rather than a fitting artifact, the resolved Kennicutt–Schmidt scatter seen in real galaxies is partly determined by $\\Sigma_*$ and $\\sigma_{\\rm gas,z}$, so observations that ignore those variables will miss real information.","Simulations that do not resolve the vertical structure of the disk cannot supply $\\sigma_{\\rm gas,z}$ and $\\Sigma_*$ reliably, so the discovered equation imposes a resolution requirement on the data it is applied to.","The instability of the 10 Myr equations implies that a deterministic, instantaneous star formation law is not the right target for stochastic short-timescale star formation; averaging over roughly a cloud lifetime is needed before the relation stabilizes."],"supporting_citations":[{"why":"Supplies the symbolic regression package whose genetic algorithm searches for candidate equations and whose score ranks them by the loss-complexity trade-off.","marker":"M. Cranmer 2023"},{"why":"Supplies the XGBoost regressor used to rank the 76 candidate variables down to 8 and to set the $R^2$ upper bound for equation-based models.","marker":"T. Chen & C. Guestrin 2016"},{"why":"Describes the FIRE-2 star formation and feedback model that produced the simulations used as training data.","marker":"P. F. Hopkins et al. 2018"},{"why":"Provides the public data release of the FIRE-2 Latte galaxies from which the pixel maps are drawn.","marker":"A. Wetzel et al. 2023"},{"why":"Defines the 750 pc face-on mapping and the pixel-level extraction of gas, stellar, and velocity-dispersion variables.","marker":"M. E. Orr et al. 2020"},{"why":"Justifies excluding redshift as a feature by finding no significant redshift evolution in the resolved Kennicutt–Schmidt relation.","marker":"M. E. Orr et al. 2018"},{"why":"Provides the top-down feedback-regulated model used as the 100 Myr baseline and as the synthetic dataset validating equation recovery.","marker":"C.-A. Faucher-Giguère et al. 2013"},{"why":"Provides the bottom-up multi-freefall model used as the 10 Myr baseline that selected equations must beat.","marker":"D. M. Salim et al. 2020"},{"why":"Defines the empirical Kennicutt–Schmidt relation whose plane anchors the comparison and whose normalization is the K98 baseline.","marker":"R. C. Kennicutt 1998"}],"fun_headline_variants":["Symbolic regression uncovers star formation scaling law","AI finds simple formula for star formation rate","Star formation law emerges from machine learning","New star formation equation derived from simulations","Star formation scaling relation discovered by AI"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The entire derivation assumes the FIRE-2 simulations are a faithful, unbiased picture of real star formation and that, at 750 parsec resolution, the star formation rate in a pixel is fully determined by the eight local quantities measured in that same pixel; if star formation history, sub-cloud physics, or the pixels excluded for not having recent star formation carry information, every fitted equation is biased.","fun_headline_variants_meta":{"raw":{"variants":["Symbolic regression uncovers star formation scaling law","AI finds simple formula for star formation rate","Star formation law emerges from machine learning","New star formation equation derived from simulations","Star formation scaling relation discovered by AI"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.00143,"raw_usage":{"total_tokens":5825,"prompt_tokens":1061,"completion_tokens":4764,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":677,"completion_tokens_details":{"reasoning_tokens":4700}},"tokens_in":677,"tokens_out":4764,"duration_ms":29771,"temperature":1.0,"reasoning_tokens":4700,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-15T23:23:05.928569+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Take the converged 100 Myr equation and apply it to 750 pc pixels of observed galaxies, or of an independent simulation suite with a different feedback prescription, and compare predicted versus measured $\\Sigma_{\\rm SFR}$; if the relation fails on galaxies with different star formation histories while still fitting FIRE-2, the convergence is a property of the training data rather than a universal star formation law.","supporting_citations":[],"review_version":1}