{"id":"43fb977e-a2a4-431b-8215-e4b1b9fa53d0","arxiv_id":"2608.13058","paper_version":1,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":5,"one_line_summary":"With a more accurate cooling model, protostellar discs fragment or form spirals under a different parameter range than earlier simulations suggested, including fragmentation in compact discs and stability up to 0.4 stellar masses.","lead":"This paper uses improved radiative cooling in simulations to map when young protostellar discs become gravitationally unstable and fragment. It finds that compact 50 au discs can fragment, that massive discs can stay stable, and that grand-design spirals are rare, which changes expectations for early planet formation.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Equation 8's self-gravitating scale-height correction has the wrong Q3D→∞ limit as printed, so the central cooling correction and the Fig. 4 parameter-space shift are not yet established.","rationale":"The reader identified the modified Lombardi cooling approximation as the weakest assumption, and I agree that the central claim depends on it. My stress-test sharpens this into a concrete, internally checkable problem: Eq. 8 as printed does not have the correct Q3D→∞ limit for a self-gravitating scale-height correction. This is not a disagreement with the broader consensus about GI; it is a potential internal inconsistency in the formula that sets the column density and hence the cooling rate. If Eq. 8 is a typographical error, the manuscript does not state the correct version, and the quantitative boundaries in Fig. 4 are not reproducible as written. If Eq. 8 is actually what was used, then the method contains a 25% scale-height bias in the opposite direction from the stated purpose, and the claimed shift in parameter space could be partly an artefact of that bias. Either way, a re-derivation of Eq. 8 and a targeted rerun of a boundary simulation would settle whether the concern lands. I keep the reader's CONDITIONAL verdict rather than moving to REJECT because the concern is specific and testable, and the paper's resolution and initial-condition checks provide some independent support for the qualitative picture. The condition should be: verify or correct Eq. 8 before the stability boundaries are relied upon.","tokens_in":22291,"tokens_out":10579,"duration_ms":111122,"concrete_test":"Independently re-derive H0/H* from vertical hydrostatic equilibrium, dP/dz = -ρ(Ω²z + 4πG∫ρ dz'), using Q3D as defined in Eq. 9, and evaluate the result at Q3D = 0.3, 1, 5, and ∞. If the derived formula differs from Eq. 8 (in particular, if it tends to 1 rather than √(π/2)), correct Eq. 8 in the code and rerun at least one boundary simulation from Fig. 4 whose classification is load-bearing, e.g. M* = 0.5 M⊙, R_out = 50 au, M_d/M* = 0.4. If the rerun changes from fragmentation to no fragmentation (or vice versa), the headline 50 au fragmentation claim and the corresponding stability boundaries are not robust to the scale-height correction.","verdict_should_be":"CONDITIONAL","load_bearing_attack":"The central claim—that the GI parameter space shifts when optical depths are estimated with the modified Lombardi method—rests directly on the self-gravitating scale-height correction in Eq. 8. As printed in §2.1, Eq. 8 gives H0/H* = √(π/2)/√(1 + 1/(Q3D√(π/2))). In the weak-self-gravity limit Q3D→∞, this tends to √(π/2) ≈ 1.25, not 1. The text introduces Eq. 8 by saying that in a self-gravitating disc H0 is reduced compared with the non-self-gravitating H* = c_s/Ω_K, and the standard vertical-hydrostatic result for the same Q3D definition (Eq. 9) is H/H* = (1 + 1/Q3D)^-1/2, which tends to 1. As written, Eq. 8 therefore contradicts its stated purpose and overestimates the scale height by 25% in the passive-disc limit. Because the modified Lombardi column density is obtained by averaging this H0 in inverse quadrature (Eq. 7), the pseudo-mean optical depth (Eq. 10), the radiative cooling/heating rates (Eqs. 11–12), and ultimately the fragmentation/stability classifications in Fig. 4 all inherit this factor. The 50 au fragmentation threshold near M_d/M_* ≈ 0.4 and the 100–200 au stability limits are set in exactly the optically thin outer regions where a 25% column-density error is most consequential. The method is not benchmarked against full radiative transfer here, and the prior validation is in a self-cited paper; but the more specific problem is that the formula as printed has the wrong asymptotic behaviour, so the quantitative boundaries in Fig. 4 are conditional on a correction that has not been shown to be correct.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"This paper revisits the parameter space for gravitational instability (GI) in young protostellar discs using SPH simulations with four approximate radiative cooling methods. The authors introduce and promote a 'modified Lombardi' cooling approximation that adjusts the local column density estimate for the reduced scale height of a self-gravitating disc. After comparing the four methods on selected cases, they run a grid of over 60 simulations with stellar irradiation for stellar masses 0.1-1.0 Msun, disc outer radii 50/100/200 au, and disc-to-star mass ratios up to ~1.4. The central claims are that compact 50 au discs can fragment for Md/M* above about 0.4, that 100-200 au discs can remain stable up to Md/M* of about 0.4 and sometimes above 1, and that large-scale grand-design spirals are uncommon, forming mainly around the lowest-mass stars or in the most compact discs. The paper also compares its outcomes with earlier Stamatellos-method results and with Haworth et al. (2020), and includes appendices testing initial-condition sensitivity and resolution.","tokens_in":22645,"tokens_out":7398,"duration_ms":81723,"significance":"If the quantitative parameter-space shifts are correct, the paper would meaningfully change expectations for where and when GI-driven fragmentation and spiral structure operate in embedded protostellar discs, with direct implications for planet formation pathways and for interpreting observations of young discs. The method comparison in Section 3.1 and the appendices on initial conditions and resolution are useful and provide checks that are not always present in parameter studies of this kind. The public availability of the modified code is a further strength. However, the central quantitative conclusions rest on the modified Lombardi scale-height correction, and one load-bearing formula in the manuscript has an incorrect asymptotic limit as printed. Because the column-density estimate feeds directly into the optical depth, cooling/heating rates, and the Fig. 4 stability boundaries, the claims in their present form are not fully established.","major_comments":[{"comment":"As printed, Eq. (8) gives H0/H* = sqrt(pi/2) / sqrt(1 + 1/(Q3D sqrt(pi/2))), which tends to sqrt(pi/2) ~ 1.25 as Q3D -> infinity. The text states that H0 is the scale height of a self-gravitating disc and H* = c_s/Omega_K is the non-self-gravitating scale height, so the physically required limit is H0/H* -> 1; the usual vertical-hydrostatic result for this Q3D definition is H/H* = (1 + 1/Q3D)^(-1/2). The printed expression therefore contradicts its stated purpose and introduces a 25% error in the scale height in exactly the weakly self-gravitating limit. Since the modified Lombardi column density is obtained by averaging this H0 in inverse quadrature (Eq. 7), the pseudo-mean optical depth (Eq. 10), the radiative cooling and heating rates (Eqs. 11-12), and ultimately the fragmentation/stability classifications in Fig. 4 all inherit this offset. Please correct Eq. (8), and state explicitly whether the simulations used the printed expression or the corrected one; if the code used the corrected expression, a text-only fix and a statement to that effect would resolve this point.","section":"Section 2.1, Eq. (8)"},{"comment":"The prefactor in Eq. (19) is inconsistent with the accompanying text. The text says the factor 1/sqrt(pi/2) is from H/H* = sqrt(pi/2), but Eq. (19) divides by sqrt(pi/2), which would produce H/H* = 1/sqrt(pi/2) rather than sqrt(pi/2). This affects the normalization of the initial sound-speed profile and therefore the initial Toomre Q of the discs. Please correct either the equation or the explanation, and confirm that the initial-condition robustness tests in Appendix A remain valid under the corrected normalization.","section":"Section 2.3, Eq. (19)"},{"comment":"The main quantitative outcomes, such as the Md/M* ~ 0.4 fragmentation threshold for 50 au discs and the stability limits for 100-200 au discs, are presented as sharp parameter-space boundaries in Fig. 4, but the classification into axisymmetric, faint spiral, extended spiral, and fragmentation is visual rather than defined by reproducible quantitative criteria. Given that Fig. 3 shows the disc evolution is highly sensitive to the thermal state and optical depth, the boundaries in Fig. 4 have an unquantified classification uncertainty. Please provide explicit criteria (for example, clump formation and survival thresholds, or spiral amplitude/contrast measures) and, if possible, an indication of how the quoted Md/M* thresholds shift when those criteria are varied by reasonable amounts.","section":"Section 3.2 and Fig. 4"},{"comment":"The modified Lombardi approximation is the basis for the paper's central conclusions, but it is not benchmarked against a full radiative-transfer calculation in this work; the validation is in the self-cited Young et al. (2024), and Section 4.4 itself lists geometries, shadowing, and dust-density variations for which the method is not suitable. A direct comparison of the modified Lombardi method against ray-tracing or Monte Carlo radiative transfer for at least one representative self-gravitating disc, even in the optically thin outer regions that set the fragmentation boundaries, would substantially increase confidence that the Fig. 4 parameter-space shift is physical rather than an artifact of the approximation.","section":"Section 2.1 and Section 4.4"}],"minor_comments":[{"comment":"The symbol u_i is described as internal energy in some equations and as internal energy density in Eq. (21); please make the notation consistent throughout.","section":"Section 2.1, Eq. (11) and Eq. (21)"},{"comment":"The black and red symbols in Fig. 4 are difficult to distinguish in the grayscale version of the manuscript; a legend with distinct marker sizes or shapes would help, since the comparison between the modified Lombardi and Stamatellos outcomes is central to the paper.","section":"Fig. 4"},{"comment":"The text says the simulations were evolved to at least 10 ORPs, but the Fig. 4 symbols do not indicate the simulation duration or the time at which the final classification was made; including this information would clarify cases where structure develops after several ORPs.","section":"Section 3.2, paragraph before Fig. 4"},{"comment":"The comparison with Haworth et al. (2020) is useful, but it is shown only for 50 au and 200 au discs; extending it to 100 au, where the new stability thresholds are a key result, would make the comparison more complete.","section":"Appendix C"},{"comment":"The statement that 'all codes that employ FLD' cannot model shadowing is correct, but the subsequent sentence about ray tracing should also mention that the present method's neglect of dust settling is shared with recent ray-tracing work; consider making this comparison explicit.","section":"Section 4.4"}],"recommendation":"major_revision","confidential_remarks":"The central issue is the incorrect asymptotic limit in Eq. (8). If this is simply a typo and the code implements the physically correct expression, the paper's quantitative conclusions may be unaffected, and the revision could be text-only plus a clarifying statement. If the code actually implements the printed expression, the claimed thresholds in Fig. 4 would need to be re-evaluated. I recommend asking the authors to state precisely which expression was implemented and to add a benchmark against full radiative transfer for at least one case. The self-citation to Young et al. (2024) is not by itself disqualifying, but an independent check would materially increase confidence in the method's accuracy."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"First, the one thing you should know: this is a systematic, carefully executed parameter study that plausibly revises the parameter space for gravitational instability in young protostellar discs. The catch is that the central cooling correction, Eq. 8, has the wrong asymptotic limit as printed, so the quantitative boundaries are not yet established without a fix.\n\nThe paper compares four radiative cooling approximations in SPH and then runs a grid of over 60 simulations covering stellar masses 0.1–1 M⊙, disc radii 50–200 au, and disc-to-star mass ratios up to 1.4. The reported differences from Haworth et al. (2020) and Cadman et al. (2020a) — 50 au discs can fragment for M_d/M_* ≳ 0.4, larger discs stay stable to ≳0.4 M_*, and grand-design spirals are rare — are physically plausible and potentially important for early planet formation. The paper also does useful methodological work: the comparison of the four cooling methods in Section 3.1 is clear, and the initial-condition and resolution tests (Appendices A and B) give real assurance.\n\nThe main issue is the scale-height correction at the heart of the modified Lombardi method. Equation 8 gives H0/H* = sqrt(π/2)/sqrt(1 + 1/(Q3D sqrt(π/2))). In the weak-self-gravity limit Q3D → ∞, that tends to sqrt(π/2) ≈ 1.25, not 1. The text says H0 should be reduced relative to H*, and the standard vertical-hydrostatic result, their Eq. 9, tends to 1. So either the equation is a typo, or the implementation overestimates the scale height by 25% in the outer, optically thin regions where the Fig. 4 boundaries are most sensitive. Since the column density, optical depth, and cooling rates all inherit that factor from Eq. 7, the fragmentation/stability thresholds are conditional on a correction.\n\nOther soft spots are minor by comparison. The method is not benchmarked against full radiative transfer in this paper; prior validation is in self-cited Young et al. (2024). The spiral classification is qualitative, the Fig. 4 boundaries have no error bars, and 'long-lived' means 20 outer rotation periods (a few thousand years). Data aren't released yet either.\n\nOverall, I think this is a serious paper with a likely-correct central shift, but it needs a major revision to fix Eq. 8 and ideally add at least one direct benchmark against full RT. I would send it to referees, and if the fix confirms the implementation, I'd be happy to cite it.","headline":"A serious, wide-ranging GI parameter study with a plausible shift in the fragmentation boundary, but the core cooling equation (Eq. 8) has a wrong asymptotic limit as printed, so the quantitative results need a fix before I'd trust them.","tokens_in":23224,"tokens_out":6871,"would_cite":true,"duration_ms":66446,"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":"With a more accurate radiative cooling treatment, compact protostellar discs can fragment at moderate mass ratios while extended discs stay stable at much higher masses than earlier simulations found.","keywords":["gravitational instability","protostellar discs","radiative cooling approximation","disc fragmentation","spiral structure","stellar irradiation","smoothed particle hydrodynamics","Toomre Q"],"falsifier":"Take a 50 au disc around a 0.5 $M_\\odot$ star at $M_d/M_* = 0.5$ and run it with full ray-tracing radiative transfer rather than the modified Lombardi cooling; if it stays axisymmetric while the modified-Lombardi version fragments (or vice versa), or if the mid-plane optical-depth profiles disagree, the claimed parameter-space shift does not survive.","tokens_in":22062,"feed_emoji":"🪐","tokens_out":5788,"duration_ms":53861,"temperature":0.7,"pith_summary":"The paper is trying to establish where gravitational instability actually operates in young protostellar discs once radiative cooling is treated more accurately. Earlier simulations used an optical-depth estimate that systematically overestimated the column density, making discs look cooler and more unstable than they are. With the corrected method, compact 50 au discs can fragment at moderate disc-to-star mass ratios, while large discs remain stable at higher masses than previously thought, and grand-design spiral arms become uncommon. This matters because it changes where and how directly collapsed planets can form, and what observers should expect to see in the youngest discs.","feed_headline":"A cooling-model fix moves disc fragmentation inward","feed_subtitle":"Compact 50-au discs can fragment at high disc-to-star mass; big discs stay stable longer.","key_machinery":"The central object is the modified Lombardi radiative cooling approximation: the pseudo-mean column density $\\bar{\\Sigma}_i$ above a gas parcel is estimated from the local pressure gradient (Lombardi's method) and then averaged in inverse quadrature with a scale height $H_0$ reduced by self-gravity via Eq. (8), $H_0/H_* = \\sqrt{\\pi/2}/\\sqrt{1+1/(Q_{3\\rm D}\\sqrt{\\pi/2})}$. This lowers the estimated optical depth and therefore raises the mid-plane temperature from stellar irradiation, making discs harder to destabilise than the older Stamatellos method. Coupled to flux-limited diffusion, this method sets cooling and heating rates that respond to local structure, and the paper uses it to redraw the fragmentation and spiral-formation boundaries.","core_discovery":"The paper argues that with the modified Lombardi radiative cooling approximation — an estimate of the column density above each gas parcel that accounts for the reduced scale height of a self-gravitating disc — the parameter space for gravitational instability changes materially. In simulations spanning 0.1 to 1 $M_\\odot$ host stars and discs of 50, 100, and 200 au, discs remain stable to disc-to-star mass ratios $M_d/M_* \\gtrsim 0.4$, and for low-mass stars extended discs can exceed $M_d/M_* = 1$ without fragmenting. The same models show fragmentation inside 50 au discs at $M_d/M_* \\gtrsim 0.4$, placing GI-driven planet formation at 20–40 au rather than only beyond roughly 70–100 au. Large-scale, grand-design spiral arms form only for $M_* \\lesssim 0.3\\,M_\\odot$ or in the most compact discs; long-lived spirals tend to be faint, flocculent, and hard to observe.","pith_inferences":["If the modified Lombardi estimate of optical depth is right, then constant-$\\beta_{\\rm cool}$ disc simulations — and the fragmentation thresholds derived from them — should be recalibrated, since the effective cooling slope in the outer disc is much steeper than a constant value.","The shift of fragmentation inward to 20–40 au predicts that young discs around low-luminosity stars are the best targets for detecting GI-driven signatures, because only there are extended spirals long-lived enough to observe.","A direct head-to-head against full ray-tracing radiative transfer for the same initial conditions, with and without dust settling, would test whether the remaining column-density approximation, not the underlying thermodynamics, is what moves the fragmentation boundary.","The stable massive discs the paper finds would be the ones to host grain-growth-assisted planet formation; a long self-regulated phase gives dust time to grow and drift into the spiral arms."],"forward_implications":["Fragmentation can occur in 50 au discs when $M_d/M_* \\gtrsim 0.4$, so direct-collapse planet formation is not confined to very extended discs; fragments appear at 20–30 au in the most compact cases.","Discs around young stars can remain gravitationally stable at $M_d/M_*$ at least 0.4, and above 1 for low-mass stars, retaining a large reservoir of solid material for planet formation.","Large-scale grand-design spirals are rare under this cooling treatment; the typical long-lived GI signature is a faint, flocculent, low-contrast spiral, which may be hard to detect.","GI can be active ($Q_{\\rm min} \\lesssim 1.5$) in discs whose spirals are so faint that no clear spiral is observable, which could explain the scarcity of observed GI spirals.","Simulations should be run for at least 10 outer rotation periods before declaring a disc stable, since some discs develop spirals or fragments only after several ORPs."],"supporting_citations":[{"why":"Introduces the pseudo-cloud optical-depth cooling approximation from which the column-density estimates derive; this is the earlier-method baseline.","marker":"Stamatellos et al. (2007)"},{"why":"Provides the pressure-gradient column-density estimate that the modified Lombardi method builds on.","marker":"Lombardi et al. (2015)"},{"why":"Introduces the modified Lombardi self-gravitating scale-height correction (Eq. 8) that the paper treats as its more accurate cooling method.","marker":"Young et al. (2024)"},{"why":"Adds flux-limited diffusion radiative transfer in the hybrid scheme used in every simulation.","marker":"Forgan et al. (2009)"},{"why":"Defines $\\beta_{\\rm cool}$ and the cooling-time threshold that links rapid cooling to fragmentation.","marker":"Gammie (2001)"},{"why":"Supplies the prior parameter-space results and MIST stellar luminosities against which the paper's outcomes are compared.","marker":"Haworth et al. (2020)"},{"why":"Another prior parameter study whose spiral and fragmentation boundaries are redrawn by the new cooling treatment.","marker":"Cadman et al. (2020a)"},{"why":"Demonstrates that constant $\\beta_{\\rm cool}$ is unrealistic and that earlier column-density estimates overestimate optical depth, motivating the method comparison.","marker":"Mercer et al. (2018)"},{"why":"Provides the phantom SPH code in which all simulations and the modified cooling implementation are run.","marker":"Price et al. (2018)"}],"fun_headline_variants":["Improved cooling model brings disc fragmentation inward","Cooling-model fix shifts disc fragmentation to 20–40 au","Fragmentation closer to star: GI planet formation at 20–40 au","Discs can be 40% stellar mass without fragmenting"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"All conclusions rest on the modified Lombardi approximation — especially the self-gravitating scale-height correction of Eq. (8) — being a faithful substitute for full radiative transfer in self-gravitating discs; the paper validates it against earlier approximations but not against full radiative transfer.","fun_headline_variants_meta":{"raw":{"variants":["Improved cooling model brings disc fragmentation inward","Cooling-model fix shifts disc fragmentation to 20–40 au","Fragmentation closer to star: GI planet formation at 20–40 au","Discs can be 40% stellar mass without fragmenting"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000743,"raw_usage":{"total_tokens":3368,"prompt_tokens":1054,"completion_tokens":2314,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":670,"completion_tokens_details":{"reasoning_tokens":2243}},"tokens_in":670,"tokens_out":2314,"duration_ms":17751,"temperature":1.0,"reasoning_tokens":2243,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-15T17:39:43.556709+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Take a 50 au disc around a 0.5 $M_\\odot$ star at $M_d/M_* = 0.5$ and run it with full ray-tracing radiative transfer rather than the modified Lombardi cooling; if it stays axisymmetric while the modified-Lombardi version fragments (or vice versa), or if the mid-plane optical-depth profiles disagree, the claimed parameter-space shift does not survive.","supporting_citations":[],"review_version":1}