{"id":"3028bbdf-a616-4207-acb2-ea507f4bf1bc","arxiv_id":"2505.24841","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":2,"one_line_summary":"A Python rewrite of STARBURST99 with new low-metallicity and very massive star models predicts a 0.3 dex boost in HI ionising flux when the upper mass limit is raised from 120 to 300 solar masses.","lead":"pySTARBURST99 is a new Python version of the widely used STARBURST99 code for modeling the light of star-forming galaxies. It adds new stellar models for low-metallicity and very massive stars, and predicts that very massive stars increase hydrogen-ionizing light by 0.3 dex in the first 2 million years.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The VMS ionising-flux increase hinges on unsampled FASTWIND spectra in the 300–500 Msol regime and may be artificially amplified by extrapolation; a direct check of input coverage is needed.","rationale":"The reader's weakest assumption—missing Eddington-enhanced mass loss in VMS evolution—is a genuine, author-acknowledged physical limitation (Sect. 2.1). It could shift VMS temperatures and lifetimes, and therefore the 0.3 dex result. After rereading the manuscript, I find the spectral-coverage concern to be more immediately load-bearing and more tractable to test: the integrated SED is built by assigning each isochrone point the nearest (or interpolated) FASTWIND model, and the paper never demonstrates that the extended grid contains the VMS region where ionising photons are produced. The 33% grid extension is asserted but never quantified in T_eff–log g coordinates. If the VMS tail is extrapolated, the headline ionising-flux increase is not a prediction of the models but a numerical artifact. The paper's own verification of pySTARBURST99 against FORTRAN (Sect. 3.1) shows up to 300% local flux differences at identical output times, so the claimed 'essentially no difference' is somewhat overstated—but this is not the central weak point. The strongest checks would be to (1) verify grid coverage of the VMS tracks, (2) re-run the 1–2 Myr HI ionising flux after computing spectra at any missing VMS grid points, and (3) optionally recompute with a VMS mass-loss prescription that includes Eddington-enhanced rates to bracket the physical uncertainty. None of these undermine the code-release contribution, but they should be resolved before the 0.3 dex result is used as a headline prediction.","tokens_in":32931,"tokens_out":1758,"duration_ms":15696,"concrete_test":"Inspect the delivered FASTWIND grid files (or the GitHub repository) and map the 180, 250, 300, and 500 Msol GENEC tracks onto the FASTWIND grid in the (T_eff, log g, luminosity) plane. If any VMS evolutionary point lies outside the convex hull of the grid, recompute the SEDs by adding one or more FASTWIND models at those coordinates. If the resulting 0.3 dex HI ionising flux increase at 1–2 Myr changes by more than ~0.05 dex, the headline result is not robust.","verdict_should_be":"CONDITIONAL","load_bearing_attack":"The central claim that extending the upper IMF mass limit from 120 to 300 Msol raises HI ionising flux by ~0.3 dex in the first 2 Myr depends on the FASTWIND grid actually covering the VMS parameter-space region where those stars contribute ionising photons. The paper says the grid was extended by 33% 'to account for evolutionary predictions of VMS' (Sect. 2.2), but it never states the maximum effective temperature, luminosity, or log g of the VMS grid points, nor how the spectrum for a 300 Msol star is obtained when the nearest grid model lies outside the star's parameters. The VMS models (Martinet et al. 2023) reach 300 Msol at Z=0.006 and Z=0.0 (Table 1), and the 300 Msol FASTWIND models are critical for the 0.3 dex result. If those spectra are extrapolated from lower-mass grid points rather than computed at the VMS parameters, the ionising flux boost could be an artifact. Additionally, the Z0 VMS track at 300 Msol reaches very high temperatures; Fig. 2 claims coverage 'at all metallicities,' but the Z0 grid is stated to use Z=1e-6 for numerical reasons (Sect. 2.2), so the Z0 VMS spectra are not truly metal-free. This does not invalidate the porting work, but it makes the headline numerical result less secure than presented.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper presents pySTARBURST99, a Python port of the STARBURST99 population synthesis code, and combines it with new GENEC evolutionary tracks (rotating and non-rotating, metallicities from Z=0.02 down to Z=0.0, and initial masses up to 300-500 Msol) and a new grid of FASTWIND synthetic spectra. The authors verify pySTARBURST99 against the FORTRAN version for the older GENEC/WMBASIC inputs, then use the new inputs to predict SEDs, HI/HeI/HeII ionising fluxes, bolometric luminosities, wind powers, H-alpha equivalent widths, and UV beta-slopes. The headline result is an increase in HI ionising flux of about 0.3 dex in the first 2 Myr when the upper IMF mass limit is raised from 120 to 300 Msol. The code and model grids are publicly available, with tabulated predictions in the appendix.","tokens_in":33167,"tokens_out":6177,"duration_ms":63859,"significance":"If the VMS-related predictions are robust, this paper provides a valuable community resource: a modern, open Python implementation of a widely used code, newly consistent low-metallicity evolution+atmosphere grids, and a clear set of falsifiable predictions for extreme star-forming populations. The explicit cross-checks against STARBURST99 in Section 3.1 and Appendix C, the tabulated output values, and the public code release are concrete strengths. The main significance risk is that the headline 0.3 dex ionising-flux boost depends on VMS evolutionary tracks that lack Eddington-enhanced mass loss and on a FASTWIND grid whose VMS coverage is not quantitatively documented; both issues are fixable with additional analysis and should be addressed before the predictions are used for quantitative inference.","major_comments":[{"comment":"The paper states in §2.1 that the VMS models of Martinet et al. (2023) 'do not contain a general increase in mass-loss rate for VMS on the main sequence', despite the physical expectation of Eddington-enhanced mass loss. The headline result of a ~0.3 dex increase in HI ionising flux when the upper mass limit is raised from 120 to 300 Msol (Abstract; §3.2; Table 2) depends on the temperatures, luminosities, and lifetimes of 180-300 Msol stars. If VMS mass loss is underestimated, these stars would be cooler or have shorter lifetimes, and the predicted ionising flux boost would change. Please quantify the sensitivity, for example by recomputing the isochrones with an enhanced mass-loss prescription or by applying a bracketing multiplicative factor to the VMS mass-loss rates and rerunning the synthesis. Without this, the 0.3 dex result rests on a known missing physical ingredient.","section":"§2.1, §3.2, Table 2"},{"comment":"The extension of the FASTWIND grid for VMS is described only as a '33% increase in the size of the model grid' and is illustrated in Fig. 2. The paper does not provide the maximum effective temperature, luminosity, or surface gravity of the added grid points, nor does it state whether the 180-300 Msol GENEC tracks lie inside the grid or require extrapolation. Since the M300 columns of Table 2 and Fig. 7 rely on these spectra, the 0.3 dex HI ionising flux increase could be an artifact of extrapolation rather than a physical prediction. Please add a table of the grid parameter ranges (Teff, log g, mass-loss rate) for each metallicity, and show that the VMS tracks are covered by the grid or explain how interpolation/extrapolation is performed.","section":"§2.2, Fig. 2"},{"comment":"The verification against STARBURST99 shows flux differences up to 300% at wavelengths <1000 A when both codes are evaluated at 1.01 Myr, and the agreement is recovered only by comparing the pySTARBURST99 output at 1.04 Myr. The text attributes this to the precision of the time increment, but the native time-step and its effect on the Appendix C comparisons are not quantified. Please state the time resolution of the isochrone outputs and confirm that the same 0.03 Myr offset (or similar) accounts for the residuals in all the Appendix C comparisons; otherwise the claim that pySTARBURST99 'faithfully reproduces' STARBURST99 is not fully established.","section":"§3.1, Appendix C"}],"minor_comments":[{"comment":"The caption reads 'WMBASICspectra at Z=0.2', which appears to be a typo for Z=0.02.","section":"Fig. 23"},{"comment":"The sentence 'Available GENECevolutionary models for initial masses from from 1 to 120M⊙' contains a duplicated 'from'.","section":"§2.1"},{"comment":"Please justify the choice of Z=10^-5 for the zero-metallicity terminal wind speed and state how the wind-power predictions depend on this assumed value; currently the Z0 wind-power entries in Table 7 are mostly empty, so the reader cannot assess the impact.","section":"§3.3, Table 7"},{"comment":"The footnote that Z0 FASTWIND models are computed with Z=10^-6 should be referenced explicitly when discussing the Z0 columns of Table 2, since the abstract and §3.2 describe these as zero-metallicity predictions; the paper's argument that Q(H) is insensitive to input metallicity for fixed stellar parameters mitigates this, but the statement should be made at the point of the Z0 predictions.","section":"§2.2, §3.2"}],"recommendation":"major_revision","confidential_remarks":"The manuscript is within the journal's scope and the code release is a clear community service. The main risk is that the headline VMS ionising-flux result is presented before the underlying model limitations are quantified; the revision should add grid coverage tables and mass-loss sensitivity tests. I would not reject on this basis, as the issues appear fixable with additional analysis within the manuscript's scope."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Bottom line: this is a solid, useful code update rather than a breakthrough science result. The Python port and the new FASTWIND grid matched to the GENEC metallicities are real work the community will want to use, and the paper is refreshingly honest about what it does not include: Eddington-enhanced mass loss for VMS, X-rays and clumping in the atmospheres, tailored mass-loss rates, and observational validation. The headline 0.3 dex HI ionising flux boost from extending the IMF to 300 Msun is a new quantitative prediction, but it inherits the uncertainties of the VMS tracks and the FASTWIND coverage, and both need scrutiny.\n\nWhat it does well: the verification against the FORTRAN STARBURST99 is comprehensive (solar, LMC, SMC, IZw18, Z0, with and without rotation, in the appendices) and mostly convincing. The 300% short-wavelength flux differences at 1.01 Myr look alarming, but the paper explains them as time-step offset; at 1.04 Myr the agreement is within 5%. Still, 'essentially no difference' overstates what Figure 3 shows, and a referee should ask for residuals at matched times. The sample output tables are a nice practical touch.\n\nThe soft spots, in descending order of severity:\n\n1. The VMS tracks (Martinet et al. 2023) deliberately omit a general Eddington-enhanced mass-loss rate on the main sequence, stated in Sect. 2.1. If real VMS lose mass faster, their Teff, luminosity, and lifetimes change, which directly affects the 0.3 dex HI boost. This is a genuine limitation, not a hidden one, but it means the headline number is a model-dependent prediction.\n\n2. The stress-test worry about the FASTWIND grid is real and unaddressed. The paper states a 33% grid extension and points to Figure 2, but never publishes the grid boundaries (max Teff, log g, luminosity) or states whether the 300-500 Msun spectra are computed at the VMS parameters or extrapolated from lower-mass grid points. The Z0 spectra use Z=1e-6 for numerical reasons, which is fair, but the claim of coverage at all metallicities needs the actual grid table. Without that, the VMS ionising-flux boost cannot be independently checked.\n\n3. No comparison to observations. The authors explicitly defer this to a future paper. That is acceptable for a code paper, but it means the new predictions are unvalidated.\n\n4. The code links are not direct URLs and there are no version identifiers, so reproducing the exact results from the manuscript alone is not possible. Minor but fixable.\n\nOn citations: the heavy GENEC co-authorship is not a problem; the predictions are outputs of published tracks and new atmosphere grids, not fitted to the target results.\n\nOverall: the porting is careful, the model choices are defensible, and the limitations are stated. This deserves peer review and, after the grid metadata and versioning are fixed, publication in a specialist journal. I would not cite the 0.3 dex number as established, but I would cite the code and the new model grid.","headline":"A useful, honest update to a standard population synthesis code, with a headline VMS ionising-flux boost that is plausible but rests on model physics and grid coverage the paper does not fully document.","tokens_in":33825,"tokens_out":4384,"would_cite":true,"duration_ms":38629,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"The updated population-synthesis code pySTARBURST99 reproduces STARBURST99 and predicts that allowing stars up to 300 solar masses raises early hydrogen-ionising flux by 0.3 dex (about a factor of two).","keywords":["population synthesis","starburst galaxies","very massive stars","ionising flux","stellar evolution","model atmospheres","UV spectral slope","Python code"],"falsifier":"Re-run the same synthesis with a very-massive-star grid that includes enhanced main-sequence mass loss, or measure the ionising flux of a young star-forming region whose upper mass limit is independently known; if the predicted difference between 120 and 300 solar-mass upper limits vanishes, or observations rule out the boost, the headline claim is wrong.","tokens_in":32673,"feed_emoji":"🌟","tokens_out":13956,"duration_ms":139994,"temperature":0.7,"pith_summary":"pySTARBURST99 claims to be a faithful, more capable successor to STARBURST99: it reproduces the old code's spectral outputs when given the same stellar inputs, and it extends the models to lower metallicities, rotating stars, and very massive stars (initial masses above roughly 100 solar masses) up to $300$–$500\\,M_\\odot$. The paper combines new GENEC evolutionary tracks with a new grid of FASTWIND model-atmosphere spectra, then uses the updated code to predict the properties of young starburst galaxies: ionising fluxes, spectral energy distributions (SEDs), bolometric luminosity, wind power, hydrogen-line equivalent widths, and UV $\\beta$ slopes. Its headline result is that raising the upper mass limit from $120$ to $300\\,M_\\odot$ raises the hydrogen-ionising flux by 0.3 dex during the first 2 Myr, an effect that matters for interpreting the earliest, most luminous phases of star-forming galaxies and for reionisation-era predictions. The paper also finds that metallicity barely changes early ionising flux but strongly boosts late-time flux at low metallicity, and that rotation keeps ionising fluxes higher for longer.","feed_headline":"Very massive stars double early ionising flux","feed_subtitle":"pySTARBURST99 shows that allowing 300-solar-mass stars nearly doubles the hydrogen-ionising light of young starbursts.","key_machinery":"The engine is isochrone synthesis: GENEC evolutionary tracks are interpolated track-to-track to arbitrary mass resolution and turned into isochrones, and each stellar position is assigned a spectrum from a new FASTWIND model-atmosphere grid built to cover the extended parameter space, including very massive stars. This grid replaces the WMBASIC low-resolution SED library and is matched in metallicity to the evolutionary tracks ($Z=0.0$, $0.0004$, $0.002$, $0.006$, $0.014$, $0.02$). The Python port uses SciPy and NumPy interpolation so that runtimes stay comparable to the FORTRAN version; the new very-massive-star tracks extend to $300\\,M_\\odot$ at low metallicity and $500\\,M_\\odot$ at solar metallicity, and their early hot phase is what produces the 0.3 dex ionising flux boost.","core_discovery":"The central claim is that a modernised population synthesis code can both reproduce a well-tested legacy tool and extend it into new physical territory. Concretely, pySTARBURST99 produces SEDs that agree with STARBURST99 to within a few percent once small time-step and interpolation differences are accounted for, giving the authors confidence to adopt new GENEC tracks (including rotation and stars up to $300$–$500\\,M_\\odot$) and a new FASTWIND spectral library. Using these inputs, the code predicts that extending the initial-mass upper limit from $120$ to $300\\,M_\\odot$ increases the H I ionising flux by 0.3 dex in the first 2 Myr, that this boost disappears once the very massive stars die out within about 3 Myr, that metallicity has little effect on early H I ionising flux (0.015 dex across $Z=0.02$ to $0.0$) but lower metallicity raises later H I flux by about 1 dex, and that rotating models maintain higher ionising fluxes after 2 Myr. Similar behaviour holds for He I and He II ionising fluxes, bolometric luminosity, and wind momentum, while the H-$\\alpha$ equivalent width and UV $\\beta$ slope show more complex dependence on very massive stars.","pith_inferences":["Beyond the paper's stated results, the missing general increase in main-sequence mass loss for very massive stars is the main open lever: if such mass loss is strong, the temperatures and lifetimes of $180$–$300\\,M_\\odot$ stars change, and the 0.3 dex early ionising flux boost could move substantially; a direct test would be to rerun the isochrones with an enhanced mass-loss prescription.","The paper matched FASTWIND metallicities to the evolutionary tracks, so users comparing older WMBASIC-based models (computed at $Z=0.02$ for solar) with the new ones should attribute part of any flux difference to the change in atmosphere metallicity rather than to stellar evolution alone.","Because binary interactions are excluded, late-time ionising flux predictions are likely lower bounds for real young populations; the paper itself notes that binary and stripped-star channels can raise ionising flux at ages beyond 10 Myr by about an order of magnitude, so combining pySTARBURST99 with binary population synthesis is a natural next test.","For users, a practical consequence is that the current Python release covers low-resolution SEDs and derived quantities but not high-resolution UV line profiles, so line-profile work should wait for the planned high-resolution spectral library."],"forward_implications":["Raising the upper mass limit from $120$ to $300\\,M_\\odot$ raises H I ionising flux by 0.3 dex in the first 2 Myr, after which the flux returns to ordinary levels once very massive stars disappear within about 3 Myr.","Metallicity has almost no effect on H I ionising flux before 2 Myr (0.015 dex from $Z=0.02$ to $0.0$), but after about 3 Myr lower metallicity raises the flux by roughly 1 dex, with zero-metallicity populations highest.","Rotating populations keep H I, He I, and He II ionising fluxes, bolometric luminosity, and wind momentum higher for longer than non-rotating populations, roughly from 3 to 10 Myr.","Including very massive stars boosts early wind momentum by about 0.43 to 0.47 dex, while overall wind momentum decreases toward low metallicity.","The H-alpha equivalent width rises with very massive stars at first but dips between about 1.6 and 2.2 Myr because those stars cool sharply, then recovers when Wolf-Rayet stars appear."],"supporting_citations":[{"why":"Provides the baseline STARBURST99 implementation of GENEC tracks that pySTARBURST99 must reproduce and extend.","marker":"Leitherer et al. 2014"},{"why":"Supplies the solar and SMC GENEC evolutionary tracks used for standard populations.","marker":"Ekström et al. 2012"},{"why":"Supplies the SMC-metallicity GENEC tracks and the mass-loss metallicity scaling used at $Z=0.002$.","marker":"Georgy et al. 2013"},{"why":"Supplies the zero-metallicity GENEC tracks that produce the extreme late-time ionising flux.","marker":"Murphy et al. 2021"},{"why":"Supplies the very massive star tracks up to $300$–$500\\,M_\\odot$ that drive the 0.3 dex ionising flux boost.","marker":"Martinet et al. 2023"},{"why":"Introduces the FASTWIND atmosphere code used to generate the new synthetic spectral grid.","marker":"Santolaya-Rey et al. 1997"},{"why":"Provides the systematic WMBASIC-versus-FASTWIND comparison that supports replacing the old spectral library.","marker":"Puls et al. 2005"},{"why":"Defines the WMBASIC OB spectral library that the new FASTWIND grid replaces.","marker":"Leitherer et al. 2010"},{"why":"Provides the main-sequence mass-loss prescription used in the evolutionary tracks and wind-power predictions.","marker":"Vink et al. 2001"}],"fun_headline_variants":["pySTARBURST99: very massive stars double early ionising flux","New pySTARBURST99 predicts double ionising flux from very massive stars","Python port of STARBURST99 nearly doubles early ionising flux","300-solar-mass stars double ionising flux in first 2 Myr","pySTARBURST99 extends to 500 Msol stars and doubles ionising flux"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that the evolutionary tracks for stars above 120 solar masses describe those stars correctly, even though the models do not include a general increase in mass loss for such stars while they steadily burn hydrogen in their cores; if these stars shed mass much faster than assumed, their temperatures, lifetimes, and luminosities would change, and the paper's headline 0.3 dex boost in ionising flux would change with them.","fun_headline_variants_meta":{"raw":{"variants":["pySTARBURST99: very massive stars double early ionising flux","New pySTARBURST99 predicts double ionising flux from very massive stars","Python port of STARBURST99 nearly doubles early ionising flux","300-solar-mass stars double ionising flux in first 2 Myr","pySTARBURST99 extends to 500 Msol stars and doubles ionising flux"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001667,"raw_usage":{"total_tokens":6721,"prompt_tokens":1161,"completion_tokens":5560,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":777,"completion_tokens_details":{"reasoning_tokens":5453}},"tokens_in":777,"tokens_out":5560,"duration_ms":41396,"temperature":1.0,"reasoning_tokens":5453,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-07T12:12:37.102257+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Re-run the same synthesis with a very-massive-star grid that includes enhanced main-sequence mass loss, or measure the ionising flux of a young star-forming region whose upper mass limit is independently known; if the predicted difference between 120 and 300 solar-mass upper limits vanishes, or observations rule out the boost, the headline claim is wrong.","supporting_citations":[],"review_version":1}