{"id":"8574e77e-da73-453c-9f0f-0dbe6f1a356d","arxiv_id":"2506.15230","paper_version":2,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":6,"one_line_summary":"Single-star MESA models with a mass-loss-driven down-stream phase imply R136a2 and R136a3 had initial masses of at least 500 solar masses, above the still redward-evolving R136a1.","lead":"The paper fits new rotating single-star models to the three most massive known stars in the R136 cluster and argues that the two less luminous ones were actually born heavier, above 500 solar masses. This matters because it would push the empirical upper mass limit for stars upward and change predictions for supernovae and black holes.","discovery_kind":"new_application","skeptic_critique":{"model":"deepseek-v4-flash","headline":"MCMC for R136a2/a3 excludes the observed mass-loss rates that the down-stream models overpredict by ~10x, so the headline initial masses rest on an incomplete confrontation.","rationale":"The reader's weakest-assumption flag identified the mass-loss switch and the exclusion of two constraints; my concern sharpens this into a specific, testable inconsistency: the models that fit L, Teff, and Y for R136a2/a3 overpredict their observed mass-loss rates by about an order of magnitude, and the MCMC excludes exactly that observable (Appendix C.2). This is the most load-bearing issue because the down-stream interpretation is not an independent inference from helium abundance alone; it is produced by the same high WR-type mass-loss rates that are in tension with the observed winds. If the observed Mdot is included, the posterior may no longer select the down-stream branch. The paper is otherwise careful: the grid is large, the MESA models and Zenodo reproduction package are provided, the MCMC methodology is transparent, and the authors openly report the flat chi-square above 500 Msun and the mass-loss discrepancy. These strengths support a conditional acceptance, not rejection. My proposed test directly settles whether the central claim survives a complete confrontation with available data. I therefore keep the reader's CONDITIONAL verdict unchanged, with the condition made more specific: the headline initial masses for R136a2/a3 should not be treated as settled until the Mdot constraint is included or the WR-rate calibration is revised and the MCMC rerun.","tokens_in":42084,"tokens_out":5027,"duration_ms":57415,"concrete_test":"Re-run the MCMC for R136a2/a3 with log Mdot from Brands et al. (2022) included as a fourth constraint (Gaussian likelihood with log Mdot = -4.48 ± 0.12 and -4.64 ± 0.08, respectively), alongside log L, log Teff, and log Ysurf, using the same interpolated grids. If the posterior moves off the down-stream branch or the initial-mass distribution shifts below ~500 Msun, the central claim fails under a full confrontation. As a complementary check, repeat the analysis on a grid built with the Nugis & Lamers (2002) WR rates or with the Sander & Vink (2020) rates scaled down to match the observed Mdot, and determine whether any down-stream solutions with Mini ≥ 500 Msun survive.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central inference that R136a2 and R136a3 are in a mass-loss-driven 'down-stream' phase from ≳500 Msun progenitors rests on models whose predicted current mass-loss rates are roughly an order of magnitude higher than the spectroscopically measured rates for these same stars, and the MCMC likelihood deliberately omits log Mdot (and surface N) because including them worsens the fit (Appendix C.2, Section 4.6.1). This is not a generic uncertainty: the down-stream branch is defined by the switch to Sander & Vink (2020) WR-type rates at X<0.4 (Section 2.6.2), and the same switch that drives the rapid mass drop and mass turnover (Figure 3) also makes the models overpredict the observed mass-loss rates of R136a2/a3 by roughly an order of magnitude. If the observed rates are taken at face value, these stars cannot presently be in the high-Mdot WR phase assumed by the down-stream models, and the helium enrichment must be explained by another channel (e.g., rotational mixing or binarity), removing the need for Mini≳500 Msun. The paper explicitly acknowledges this tension but defers it to future calibration of the WR rates and switch condition. Furthermore, the factor-of-two and Nugis & Lamers (2002) wind tests are not propagated through the MCMC, so the sensitivity of the headline initial masses to the WR-rate normalization is unquantified. Because the headline masses for a2/a3 are derived from a likelihood that excludes a directly relevant, contradictory observable, the 'rigorous confrontation' is incomplete as it stands.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"This paper presents rotating single-star MESA evolutionary models for very massive stars (100–800 Msun) at LMC metallicity, adopting an empirical wind-momentum–luminosity relation calibrated from R136 stars for the optically-thin regime and switching to Sander & Vink (2020) Wolf-Rayet-type rates when the surface hydrogen mass fraction drops below 0.4. The authors build interpolated grids and use an MCMC analysis to fit the luminosity, effective temperature, and surface helium abundance of the three WNh stars R136a1/a2/a3, inferring an initial mass of 346±41 Msun for R136a1 and ≳500 Msun for both R136a2 and R136a3, which they interpret as being in a 'down-stream' phase of rapid mass loss and decreasing luminosity after the mass-turnover.","tokens_in":42417,"tokens_out":6226,"duration_ms":58641,"significance":"If the central inference holds, the paper would establish that R136a2 and R136a3, currently less luminous than R136a1, were initially more massive and are now in a mass-loss-dominated evolutionary stage, with implications for the stellar upper-mass limit, the interpretation of WNh stars, and the formation of massive black holes and supernovae. The paper has notable strengths: the full model grid, inlists, and analysis scripts are openly available on Zenodo; the parameter space in mass, rotation, angular-momentum transport, and MLT++ is explored systematically; the mass-turnover phenomenon is clearly identified; and the sensitivity to alternative wind prescriptions is tested for at least one branch. The work also produces falsifiable predictions, such as Type Ib/c supernova progenitors, no pair-instability supernovae, no GRB progenitors, and a measurable luminosity decline that can be checked with future monitoring.","major_comments":[{"comment":"The MCMC likelihood is constructed from log L, log Teff, and log Ysurf only; log Mdot and surface nitrogen are deliberately excluded because including them worsens the fit (Appendix C.2). As acknowledged in §4.6.1, the posterior mass-loss rates for R136a2 and R136a3 are about an order of magnitude higher than the values spectroscopically determined by Brands et al. (2022). This is not a peripheral tension: the down-stream phase into which the paper places R136a2/a3 is defined by the switch to Sander & Vink (2020) WR-type rates at X<0.4 (§2.6.2), and the same switch drives both the rapid mass drop (Figure 3) and the factor-of-ten overprediction. If the measured mass-loss rates are correct, these stars cannot presently be in the high-Mdot WR phase, and the surface helium enrichment must be produced by another mechanism, removing the need for Mini>500 Msun. I request that the authors either include log Mdot in the likelihood with realistic systematic uncertainties, or provide a dedicated sensitivity test showing how the inferred initial masses shift when the WR rates are reduced by the factor needed to match the observed rates.","section":"§4.6.1 and Appendix C.2"},{"comment":"The fiducial optically-thin mass-loss rates are the empirical WLR calibration of Brands et al. (2022), and both the fitted observables (L, Teff, Ysurf) and the mass-loss rates used for comparison in Table 2 are also from Brands et al. (2022). This means the models are calibrated to the same system they are then used to fit, reducing the independence of the evolutionary inference. In particular, the WLR parameters x and log D0 were derived including R136a1-a3, so the early main-sequence mass loss is effectively tuned to reproduce these stars. The paper should quantify how the inferred initial masses and ages change when an alternative WLR (e.g., from a sample excluding the WNh stars, or from the theoretical Vink et al. 2001 rates) is adopted.","section":"§2.6.1 and §4.1"},{"comment":"The sensitivity tests with halved mass-loss rates and with Nugis & Lamers (2002) WR rates are run for a single branch (TS, 300 km/s, MLT++) and are not propagated through the MCMC. This matters because the inferred initial masses for R136a2/a3 are the central quantitative claim, and the WR-rate normalization is estimated to be uncertain by an order of magnitude (§2.6.2). Section 3.2 shows that halving the mass-loss rates roughly doubles the maximum pre-supernova mass (from 24 to 67 Msun), so the ≳500 Msun lower limits are unlikely to be robust to the WR-rate uncertainty. I ask that the reduced chi^2 projection onto initial mass (Appendix C.5) be recomputed for the alternative wind branches, or that an explicit argument be given why the MCMC grids need not include them.","section":"§2.6.4 and §3.2"}],"minor_comments":[{"comment":"Radcliff136 should be written as 'Radcliffe 136'.","section":"§1"},{"comment":"Taylor-Spruit should be 'Tayler-Spruit' to match the notation used elsewhere in the paper.","section":"Figure 4 caption"},{"comment":"mesaSoftware Development Kit should be 'MESA Software Development Kit'.","section":"§2.9"},{"comment":"The statement that the mass-loss and surface nitrogen constraints 'cancel out' is qualitative; please show the individual chi^2 distributions or a corner plot with these quantities included so the reader can assess the trade-off.","section":"Appendix C.2"}],"recommendation":"major_revision","confidential_remarks":"The paper is honest about its limitations, but the gap between the headline claim (Mini≳500 Msun) and the excluded mass-loss constraints is substantial. I recommend major revision: the authors should be encouraged to include the mass-loss rate in the likelihood or to provide a quantitative sensitivity analysis, and the circularity associated with the WLR calibration deserves a prominent caveat. Otherwise the manuscript is well-structured, the models are described in detail, and the data availability is exemplary for the field."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The paper deserves a serious referee, but the headline inference is more fragile than the abstract lets on. The genuinely new step is using surface helium to break the mass-turnover degeneracy and propose that R136a2 and R136a3, currently fainter than R136a1, were initially more massive and are now in the mass-loss-driven down-stream phase. That is a real contribution, even though the down-stream trajectory itself was already in prior work.\n\nWhat the paper does well: it ships a large grid (608 MESA models), a dense interpolation, and a Zenodo reproduction package. The MCMC is documented and the authors are unusually candid about its limitations—they state that the χ² is flat above 500 M_sun for the two stars, and they explicitly report the order-of-magnitude uncertainty in the Sander & Vink WR rates. The factor-of-two and Nugis & Lamers tests are useful sensitivity checks, and the discussion of AM transport effects on GRB and PISN progenitors is a solid byproduct.\n\nThe soft spots are in proportion. The stress-test concern is not a strawman: Appendix C.2 says the likelihood excludes log Ṁ and surface N because including them worsens the fit, and Section 4.6.1 notes the down-stream models predict mass-loss rates for R136a2/a3 about ten times higher than the spectroscopically derived values. That discrepancy matters because the down-stream branch is defined by the switch to WR-type rates at X<0.4. If the observed rates are taken at face value, these stars are not in the high-Ṁ phase, and the helium enrichment would need rotational mixing or binarity, removing the need for ≳500 M_sun progenitors. The authors acknowledge this tension but defer it. The circularity concern—WLR calibration from the same stars—is real but minor-to-moderate; the factor-of-two tests help, though they are not propagated through the MCMC.\n\nBottom line: the paper is honest, reproducible, and important enough to referee, but the central claim is a plausible lower limit, not a determination. A referee should ask the authors to propagate the WR-rate normalization and switch condition through the MCMC and to report the fit with and without the mass-loss and nitrogen constraints. I would send it to review, and I would bring it to a reading group for the discussion.","headline":"A careful, reproducible grid of VMS models whose headline ≳500 M_sun initial masses for R136a2/a3 depend on a mass-loss switch that disagrees with those stars' observed wind rates by an order of magnitude.","tokens_in":43025,"tokens_out":3539,"would_cite":true,"duration_ms":34438,"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":"Single-star evolution can explain the three most massive stars known, with R136a2 and R136a3 born heavier than record-holder R136a1.","keywords":["very massive stars","stellar evolution","mass loss","Wolf-Rayet stars","R136","30 Doradus","helium enrichment","supernova progenitors"],"falsifier":"Show that R136a2 or R136a3 has a surface helium mass fraction below about 0.50, or measure its luminosity holding steady over a decade rather than dropping by roughly 1 dex per Myr, and the down-stream, initially-heavier-than-R136a1 interpretation is ruled out.","tokens_in":41794,"feed_emoji":"⭐","tokens_out":5571,"duration_ms":55118,"temperature":0.7,"pith_summary":"The paper argues that single-star evolution, with no binary merger, can explain the three most massive stars known: R136a1, R136a2, and R136a3 in the LMC's R136 cluster. Its key move is a mass-loss switch: once a very massive star's surface hydrogen drops below 40 percent, Wolf-Rayet-type winds strip it so fast that its luminosity and mass plummet, so the initially heaviest star can end up lighter than a star born lighter. This resolves the helium conundrum -- the currently brightest star is the least helium-rich -- by making R136a2 and R136a3 the initially most massive, now in the 'down-stream' phase. Fitting luminosity, effective temperature, and surface helium with a Markov-Chain Monte Carlo analysis over interpolated grids yields an initial mass around $346 \\pm 41\\,M_\\odot$ for R136a1 and $\\gtrsim 500\\,M_\\odot$ for the other two. If right, the record-breaking stars need not be exotic merger products, but their winds reshape their fate, precluding pair-instability supernovae and leading to Type Ib/c if they explode.","feed_headline":"R136a2 and R136a3 likely started above 500 solar masses","feed_subtitle":"Stellar-wind stripping can make the heaviest newborns fade fastest, solving a helium puzzle among the most massive stars known.","key_machinery":"The load-bearing mechanism is a two-regime mass-loss prescription. While surface hydrogen $X>0.4$, winds follow the empirically calibrated wind-momentum luminosity relation; below $X=0.4$, the models switch to Wolf-Rayet-type rates, which are so strong that the star's mass and luminosity plunge on a short timescale. This switch creates a 'down-stream' trajectory and the mass turnover that lets initial masses be reconstructed only with added constraints like surface helium. The comparison machinery is an interpolated grid of 20 million points with a Markov-Chain Monte Carlo fit to three observables: luminosity, effective temperature, and surface helium mass fraction.","core_discovery":"The central discovery is the mass turnover: at LMC metallicity, an initially 800 $M_\\odot$ model loses mass so rapidly once Wolf-Rayet-type winds begin that it becomes less massive than an initially 100 $M_\\odot$ model within the main sequence. Applied to R136, the turnover means the current luminosity ranking is not the birth ranking. The authors propose R136a2 and R136a3 are past the turnover, having shed hundreds of solar masses and revealing near-core helium, while R136a1 is younger and still evolving redward. From simultaneous fits to luminosity, effective temperature, and surface helium, they estimate $M_\\mathrm{ini}=346\\pm41\\,M_\\odot$ for R136a1 and $\\gtrsim500\\,M_\\odot$ for R136a2 and R136a3, with ages 1.02, 1.45, and 1.48 Myr. The fiducial grid predicts no pair-instability and no gamma-ray-burst progenitors, and any supernova would be Type Ib/c.","pith_inferences":["If R136a2 and R136a3 really started above 500 $M_\\odot$, their formation in the local universe strains standard massive-star formation theories; testing that would require studying the molecular-cloud conditions that could concentrate hundreds of solar masses in one core.","The paper's own factor-of-two wind test shows the outcome is brittle: reducing all mass-loss rates by half roughly doubles pre-supernova masses and can restore pulsational pair-instability, so a future empirical calibration of Wolf-Rayet wind rates at low metallicity would be a sharper test than any new spectrum of R136 itself.","A cleaner observational discriminator is the helium abundance of R136a1, the most disputed measurement; if a re-analysis pushed it above about 0.6, R136a1 too could be in the down-stream phase, which would change the inferred age and initial mass."],"forward_implications":["At LMC metallicity, all stars born above about 250 $M_\\odot$ should become helium-rich main-sequence stars, regardless of initial spin or internal coupling scheme.","The mass turnover means present-day masses in young clusters cannot be read back as initial masses; by the time the first supernova explodes, no model star exceeds roughly 130 $M_\\odot$ in the fiducial grid.","R136a2 and R136a3 should currently be losing luminosity at about 1 dex per Myr while staying hot, a change of roughly 25 ppm per decade in visual magnitude.","The fiducial models predict no pair-instability supernovae and no long gamma-ray-burst progenitors at LMC metallicity; the most massive pre-supernova stars have carbon-oxygen cores below the pair-instability limit.","If these stars explode, the supernovae should be Type Ib/c, and the stars should not form exceptionally massive black holes."],"supporting_citations":[{"why":"Supplies the spectroscopic observables (luminosity, effective temperature, surface helium) and the wind-momentum luminosity calibration used for the optically-thin wind regime.","marker":"Brands et al. (2022)"},{"why":"Provides the Wolf-Rayet-type mass-loss rates that trigger the rapid mass and luminosity drop in the down-stream phase.","marker":"Sander & Vink (2020)"},{"why":"Established the WNh classification and current masses for the R136 stars, defining the objects the models must explain.","marker":"Crowther et al. (2010)"},{"why":"Provides independent measurements of luminosity, effective temperature, surface helium, and mass-loss rates that serve as alternative constraints.","marker":"Bestenlehner et al. (2020)"},{"why":"Earlier VMS evolution models at LMC metallicity that this work extends and compares against, especially on rotational and mass-loss treatment.","marker":"Köhler et al. (2015)"},{"why":"Identified the down-stream or vertical evolution trajectory in non-rotating models, which this paper places in a rotating, calibrated context.","marker":"Higgins et al. (2022)"},{"why":"Scrutinized the optically-thin to optically-thick wind transition, framing the switch condition used in the models.","marker":"Sabhahit et al. (2022)"}],"fun_headline_variants":["R136a1 wasn't heaviest at birth: two others topped 500","Mass turnover: R136's heaviest were not born heaviest","Helium hints: R136a2, a3 began above 500 solar masses","Stellar winds reverse R136's mass order"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The inference that R136a2 and R136a3 began above 500 solar masses rests on switching to Wolf-Rayet mass-loss rates once surface hydrogen drops below 40 percent, rates that are order-of-magnitude uncertain and calibrated for hydrogen-free stars; if that switch happens later or weaker, the mass turnover and inferred birth masses change.","fun_headline_variants_meta":{"raw":{"variants":["R136a1 wasn't heaviest at birth: two others topped 500","Mass turnover: R136's heaviest were not born heaviest","Helium hints: R136a2, a3 began above 500 solar masses","Stellar winds reverse R136's mass order"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000594,"raw_usage":{"total_tokens":2904,"prompt_tokens":1186,"completion_tokens":1718,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":802,"completion_tokens_details":{"reasoning_tokens":1640}},"tokens_in":802,"tokens_out":1718,"duration_ms":14293,"temperature":1.0,"reasoning_tokens":1640,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-15T19:39:46.783010+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Show that R136a2 or R136a3 has a surface helium mass fraction below about 0.50, or measure its luminosity holding steady over a decade rather than dropping by roughly 1 dex per Myr, and the down-stream, initially-heavier-than-R136a1 interpretation is ruled out.","supporting_citations":[],"review_version":1}