{"id":"ad8b183a-c8be-4009-ad40-497da7393380","arxiv_id":"2412.09924","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":6,"one_line_summary":"In semi-global linear theory, thermal stratification bifurcates VSI surface modes into two branches, boosts their growth rates, and raises the radial-to-vertical kinetic energy ratio relative to isothermal disks.","lead":"This paper computes how vertical temperature stratification changes the vertical shear instability in planet-forming disks. Thermal stratification strengthens the instability's growth and shifts its energy toward radial motions, which changes how VSI turbulence is expected to appear in simulations and observations.","discovery_kind":"new_application","skeptic_critique":{"model":"deepseek-v4-flash","headline":"No-flow boundary at ζ=±5 sits exactly on the second-branch surface modes; without domain-size sensitivity tests, the claimed mode bifurcation may be a numerical artifact.","rationale":"The reader's weakest assumption correctly identifies the untested no-flow boundary at ζ = ±5 as the most load-bearing concern. The second-branch surface modes, which are a central new result, are localized at the boundary, so the boundary condition could be creating or significantly modifying them. The authors' convergence check addresses grid resolution, not domain-size independence, so the concern is not answered by the manuscript. The efficient-cooling/isothermal idealization is also important but is explicitly acknowledged and partially justified via the critical cooling-time argument, making it less decisive for the claimed bifurcation. I agree with the reader's conditional verdict: the paper is transparent and plausible, but the boundary sensitivity test is needed before the bifurcation claim can be fully accepted. Since my concern matches the reader's weakest assumption and the recommended conditional status remains appropriate, no verdict change is required.","tokens_in":14703,"tokens_out":7599,"duration_ms":82852,"concrete_test":"Recompute the eigenvalues and eigenfunctions of Eq. (29) for the n = 2 and n = 3 disks with the vertical domain extended to ζ = ±7.5, ±10, and ±15, using the same background profiles (extending the hydrostatic density, temperature, and shear beyond Z = 50 au) and the same no-flow boundary condition. Track the second-branch surface modes that peak at |ζ| ~ 5 in Fig. 6b,f. If their growth rates change by more than about 10% or if they disappear entirely as the domain grows, the bifurcation claim is not robust. As a complementary check, replace the no-flow condition with a sponge layer or an open boundary at ζ = ±5 and compare the mode structure.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central novel claim is that in thermally stratified disks the VSI surface modes bifurcate into two branches, with the second branch localized at |ζ| ~ 5 (Figs. 5 and 6) and attributed to strong shear near the disk surfaces. However, the eigenvalue problem is solved on the finite domain |ζ| ≤ 5 with the no-flow boundary condition ∂Π/∂ζ = 0 imposed at ζ = ±5. The eigenfunctions of the second-branch surface modes peak right at this boundary (Fig. 6b,f), so the boundary condition is not a distant, passive truncation but an active reflecting wall that can support or modify these modes. The authors verify spectral convergence by comparing N = 400 and N = 600 Chebyshev grids, but this only checks resolution for a fixed domain; it does not test sensitivity to the domain size. If the second-branch modes are boundary-dominated artifacts, their existence and growth rates would change when the boundaries are moved outward, undermining the bifurcation claim. A related internal inconsistency is that Eq. (22) defines ζ ≡ z/h0, which would place ζ = ±5 at Z ≈ 5 au (for R0 = 100 au, H0 = 10 au), whereas the text repeatedly states the boundaries are at Z = 50 au; clarifying this scaling is necessary to interpret the domain. The efficient-cooling/isothermal assumption is a recognized modeling limitation and is flagged by the authors, so it is less decisive than the untested boundary placement.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"This paper analyzes the linear stability of the vertical shear instability (VSI) in protoplanetary disks with vertical thermal stratification. The authors derive a semi-global eigenvalue equation (Eq. 29) for axisymmetric perturbations, adopting an isothermal equation of state and a background equilibrium with a vertically varying temperature profile. They solve the eigenvalue problem with Chebyshev collocation for isothermal (n=1) and stratified (n=2,3) disk models. They recover previously known surface and body modes in the isothermal limit, and for stratified disks they report a bifurcation of the surface modes into two branches, along with enhanced growth rates and an increased radial-to-vertical kinetic energy ratio. The paper concludes that VSI simulations will initially excite large-wavenumber surface modes, followed by small-wavenumber body modes, with the transition occurring earlier in more stratified disks.","tokens_in":15001,"tokens_out":14478,"duration_ms":136952,"significance":"If the findings are correct, the paper extends VSI linear theory to a more realistic disk temperature structure and makes concrete predictions about mode selection and kinetic energy partition that can be tested in upcoming simulations. The derivation of Eq. (29) is clear and reproduces the isothermal results of Nelson et al. (2013) and Barker & Latter (2015). The use of a public spectral code with a stated convergence check is a strength. However, the central novelty—the surface-mode bifurcation—depends on the placement of the computational boundary and on a vertical-coordinate scaling that is currently inconsistent between Eq. (22) and the text. These issues must be resolved before the main conclusions can be accepted.","major_comments":[{"comment":"The vertical coordinate scaling is inconsistent with the physical disk height. With H0 = 10 au and h0 = 0.1, the definition ζ ≡ z/h0, where z = Z/H0, gives ζ = Z/(h0H0) = Z/(1 au). The computational domain |ζ| ≤ 5 therefore corresponds to |Z| ≤ 5 au, not to |Z| = 50 au as stated repeatedly (e.g., 'strong shear at the disk surfaces (Z = 50 au)' in Section 5). This discrepancy changes the interpretation of the surface modes: the second-branch modes at |ζ| ≈ 5 would lie at Z ≈ 5 au, inside the stratified region and far from the temperature-transition height Zq = 3H = 30 au. The authors must either correct the scaling definition or the mapping to physical units and re-examine whether the reported mode structure and growth rates remain the same.","section":"4.1, Eq. (22); Section 5"},{"comment":"The existence and growth rates of the second-branch surface modes have not been shown to be independent of the numerical boundary. These modes are localized at |ζ| ≈ 5, which is exactly the position of the imposed no-flow boundary condition ∂Π/∂ζ = 0. The convergence check with N = 400 and 600 grid points only establishes resolution convergence for a fixed domain; it does not test the sensitivity to the domain size. If the boundary is moved outward (e.g., |ζ| ≤ 10 or 20), the second-branch modes could disappear or change growth rate substantially if they are supported by the reflecting wall rather than by the disk physics. The authors should repeat the eigenvalue calculations for larger domains and show that the bifurcation of surface modes, including the growth rates and eigenfunctions of the second branch, converges as the boundary recedes. Without this test, the central claim of a surface-mode bifurcation remains unverified.","section":"4.2.2, Figures 5 and 6"},{"comment":"The isothermal equation of state is justified by the critical cooling time criterion of Eq. (5), but the paper does not evaluate whether the actual cooling time in the modeled disks satisfies tcool ≲ tcrit. The discussion in Section 5 notes that the isothermal assumption 'may not be valid at large radii in PPDs,' which is a relevant caveat, but it leaves open the possibility that the reported growth-rate enhancements for n = 2 and 3 are overestimated if tcool exceeds tcrit in the regions where the new surface modes reside. The authors should provide a quantitative estimate of tcool for their disk parameters (or at least for a representative dust model) and state explicitly that the results apply only where efficient cooling holds.","section":"2, Eq. (5); Section 5"}],"minor_comments":[{"comment":"The title contains a typo: 'DICUSSION' should be 'DISCUSSION'.","section":"Section 5 title"},{"comment":"The statement 'surface modes disappear when the radial wavenumber k surpasses the threshold kcrit ∼ π/2' is internally inconsistent with the example k = π/5, which is smaller than π/2. The intended meaning is presumably that surface modes require k > kcrit; the wording should be corrected.","section":"4.2.1"},{"comment":"Equation (2) appears to have a rendering issue in the integral: 'Z Z 0' should be a definite integral from 0 to Z.","section":"2, Eq. (2)"},{"comment":"The sentence 'But, how does the ratio of the radial to vertical kinetic energy vary with the degree of thermal stratification.' is a fragment and should be completed or merged with the following sentence.","section":"4.2.3"},{"comment":"The critical wavenumbers kcrit for n = 2 and 3 are reported as π and 2π/3, but the method by which these values were determined is not described; please clarify the criterion used to define the disappearance of surface modes.","section":"4.2.2"},{"comment":"The x-axis tick labels in Figure 6 (e.g., '4 2 0 2 4') are ambiguous because the minus signs are not visible in the manuscript; please ensure the axis is clearly labeled with values such as -4, -2, 0, 2, 4.","section":"Figure 6"}],"recommendation":"major_revision","confidential_remarks":"The paper addresses a relevant problem and the linear framework is sound, but the scaling inconsistency and the untested boundary placement are serious enough that I could not recommend acceptance without revision. The authors should be asked to clarify the coordinate mapping and to add a domain-size convergence test. If the second-branch modes vanish or change when the boundary is moved, the paper's main claim would need to be substantially weakened. The subsequent simulation paper will likely provide a more decisive test of the predictions."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Main take: this is a competent, transparent semi-global VSI calculation for thermally stratified disks, with a clean derivation that reduces to earlier isothermal results. The new physics claim, the two-branch surface-mode structure, sits right on the numerical boundary and is not yet secured.\n\nWhat's new: it's the first semi-global VSI analysis with an irradiation-like T(Z) profile. Eq. (29) is the natural generalization of Nelson et al. (2013) and Barker & Latter (2015), and they verify convergence at N=400 and 600. The growth-rate enhancement with stratification and the RKE trends are plausible and testable, and no fitted parameter feeds the prediction. That's real credit.\n\nSoft spots: the no-flow boundary at ζ=±5 is not a distant truncation for the second branch—those eigenfunctions peak at |ζ|~5 (Fig. 6b,f). Spectral convergence at fixed domain does not test whether the mode survives a larger domain. Without that test, the bifurcation claim could be a boundary artifact. That's the load-bearing weakness. There's also a scaling inconsistency: Eq. (22) defines ζ≡z/h0, which for their parameters puts ζ=5 at Z=5 au, while the text repeatedly says Z=50 au. That needs fixing to know which modes are actually \"surface\" modes. And the kcrit sentence in §4.2.1 (\"disappear when k surpasses... illustrated for k=π/5\") is backwards as written.\n\nThe isothermal EOS is a real limitation, but they flag it and give a critical-cooling-time argument, so it's not a hidden flaw. The first branch at |ζ|~1.6 and the body modes are well away from the boundary and more robust.\n\nWho it's for: VSI and protoplanetary disk turbulence people, especially those designing simulations. It deserves a serious referee, but the referee should demand boundary-sensitivity tests and the scaling fix. I'd send it to review; I wouldn't cite the bifurcation claim until the boundary question is settled.","headline":"A clean linear VSI analysis for stratified disks whose central new claim—the surface-mode bifurcation—is undermined by untested boundary placement.","tokens_in":15576,"tokens_out":4936,"would_cite":false,"duration_ms":49765,"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":"Thermally stratified protoplanetary disks grow the vertical shear instability faster and with more radial kinetic energy.","keywords":["vertical shear instability","protoplanetary disks","thermal stratification","linear stability analysis","surface modes","body modes","vertical shear","disk turbulence"],"falsifier":"Recompute the eigenvalues of Eq. (29) with the vertical boundaries moved from $\\zeta = \\pm 5$ to $\\zeta = \\pm 7$ or replaced by radiative boundary conditions; if the second surface-mode branch's growth rate changes substantially or the branch disappears, the two-branch result is an artifact of the no-flow walls.","tokens_in":14474,"feed_emoji":"🪐","tokens_out":11164,"duration_ms":97763,"temperature":0.7,"pith_summary":"The paper tries to establish that the vertical shear instability (VSI), a purely hydrodynamic mechanism that can drive turbulence in protoplanetary disks, is substantially stronger in disks whose surfaces are hotter than their midplanes, as real irradiated disks are. Using a semi-global linear analysis (local in radius, global in height) of three disk models with atmosphere-to-midplane temperature ratios 1, 2, and 3, it finds that thermal stratification splits the previously known surface modes into two branches, raises their growth rates by up to a factor of about 2, and increases the share of kinetic energy in the radial direction. These changes matter because VSI turbulence regulates dust motions and planet formation, and because they give concrete predictions for what hydrodynamic simulations and molecular-line observations should see.","feed_headline":"Hotter disk surfaces speed up vertical shear instability","feed_subtitle":"Stratified disks grow the instability up to ~2x faster and with more radial motion, reshaping how turbulence stirs dust.","key_machinery":"The load-bearing object is a second-order linear ODE for the vertical structure of the density perturbation, Eq. (29): $$\\frac{\\$partial^{2}$ \\hat{\\Pi}}{\\partial \\$zeta^{2}$} + \\left(\\frac{\\partial \\bar{\\Pi}}{\\partial \\zeta} + \\frac{\\partial \\ln f}{\\partial \\zeta} + 2ik\\bar{q}\\right)\\frac{\\partial \\hat{\\Pi}}{\\partial \\zeta} - \\$sigma^{2}$ $k^{2}$ \\hat{\\Pi} = 0,$$ with no-flow boundary conditions at $\\zeta = \\pm 5$. Here $\\zeta$ is the scaled vertical coordinate, $f(\\zeta) = T(\\zeta)/T_{\\rm mid}$ is the vertical temperature profile, $\\bar{q}$ is the vertical shear parameter, $k$ is the scaled radial wavenumber, and $\\sigma$ is the complex eigenfrequency whose real part is the growth rate. This equation turns the VSI into a vertical eigenvalue problem: the eigenfunctions classify modes by where their vertical structure is concentrated, and the non-monotonic $\\bar{q}(\\zeta)$ produced by stratification is what splits the surface modes into two branches.","core_discovery":"The central claim is that thermal stratification changes both the spectrum and the character of VSI eigenmodes. In an isothermal disk the unstable modes separate into surface modes, confined to regions of strongest vertical shear, and body modes that extend across the disk. When the disk is vertically stratified, the shear profile $q = -R \\partial \\ln \\Omega / \\partial Z$ develops a local maximum away from the surfaces, and the surface modes bifurcate into two branches: one localized near that interior shear peak (around $Z = 16$ au for the $n = 2$ and 3 models at $R_0 = 100$ au) with the larger growth rate, and one near the disk surfaces with growth nearly unchanged from the isothermal case. Stratification also increases body-mode growth at small radial wavenumbers and raises the ratio of radial to vertical kinetic energy; for $k = \\pi/5$ the $n = 3$ disk's fundamental body mode has about three times the energy ratio of the $n = 1$ disk.","pith_inferences":["If the vertical boundaries were moved from $\\zeta = \\pm 5$ outward (or replaced by open conditions), the second surface-mode branch, localized at $|\\zeta| \\approx 5$, might shift or vanish, which would mean the two-branch bifurcation is partly a box effect.","Solving the same eigenvalue problem with finite cooling time and the Brunt–Väisälä frequency from Eq. (6) would test how much of the stratification-driven growth survives where the disk cools too slowly for the isothermal assumption.","The preference for small-$k_R$ body modes in stratified disks implies that numerical experiments need wide radial domains; periodic radial boxes could suppress exactly the modes stratification favors.","A larger radial-to-vertical energy ratio should make VSI turbulence in flared irradiated disks more anisotropic, so high-resolution molecular-line observations may be able to distinguish stratified from isothermal disk models without resolving individual modes."],"forward_implications":["In a stratified disk, the fastest-growing disturbances will be surface modes with large radial wavenumber $k_R$; body modes with small $k_R$ take over later, and the turnover happens earlier when the atmosphere is hotter.","The most unstable surface mode at $k = 5\\pi$ grows about 1.4 times ($n = 2$) and 2.1 times ($n = 3$) faster than its isothermal counterpart, so turbulence should develop quicker in more stratified disks.","Thermal stratification raises the radial-to-vertical kinetic-energy ratio $R_{\\rm KE}$; for the $k = \\pi/5$ fundamental mode the $n = 3$ disk reaches roughly three times the $n = 1$ value, which alters the expected velocity signature of the turbulence.","Applying the local theory to the density-weighted mean shear predicts growth rates, radial wavelengths, and turbulent amplitudes in the $n = 3$ disk about 5–15 times those in the $n = 1$ disk, making stratification a first-order control on VSI outcomes."],"supporting_citations":[{"why":"Supplies the semi-global reduced model and the original surface/body mode classification that this paper extends to stratified disks.","marker":"Nelson et al. (2013)"},{"why":"Provides the polytropic semi-global analysis and the analytic body-mode frequency formula used for comparison with the eigenfrequencies.","marker":"Barker & Latter (2015)"},{"why":"Gives the local dispersion relation and the maximum growth rate scaling $\\sim \\Omega q$ that motivate the mode expectations.","marker":"Latter & Papaloizou (2018)"},{"why":"Derives the critical cooling-time criterion $t_{\\rm cool} \\lesssim t_{\\rm crit}$ used to justify the isothermal assumption and bound its validity.","marker":"Lin & Youdin (2015)"},{"why":"Establishes the local VSI growth-rate scaling proportional to the disk aspect ratio and the basic instability criteria.","marker":"Urpin (2003)"},{"why":"Provides the midplane temperature model used to construct the stratified disk equilibria.","marker":"Barraza-Alfaro et al. (2021)"},{"why":"Supplies the pseudo-spectral eigensolver used to turn Eq. (29) into a matrix eigenvalue problem.","marker":"Burns et al. (2020)"}],"fun_headline_variants":["Stratified disks grow shear instability up to 2x faster","Thermal layering speeds vertical shear instability","Stratified protoplanetary disks boost VSI turbulence","Surface modes drive faster instability in stratified disks","Disk stratification accelerates VSI and radial motion"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The calculation assumes the gas cools infinitely fast (isothermal) and puts rigid walls at $\\zeta = \\pm 5$; the second surface-mode branch sits right at those walls, so its existence and growth rate may be set by the artificial boundary rather than by disk physics.","fun_headline_variants_meta":{"raw":{"variants":["Stratified disks grow shear instability up to 2x faster","Thermal layering speeds vertical shear instability","Stratified protoplanetary disks boost VSI turbulence","Surface modes drive faster instability in stratified disks","Disk stratification accelerates VSI and radial motion"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000287,"raw_usage":{"total_tokens":1686,"prompt_tokens":948,"completion_tokens":738,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":564,"completion_tokens_details":{"reasoning_tokens":664}},"tokens_in":564,"tokens_out":738,"duration_ms":7381,"temperature":1.0,"reasoning_tokens":664,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-11T16:33:52.544670+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Recompute the eigenvalues of Eq. (29) with the vertical boundaries moved from $\\zeta = \\pm 5$ to $\\zeta = \\pm 7$ or replaced by radiative boundary conditions; if the second surface-mode branch's growth rate changes substantially or the branch disappears, the two-branch result is an artifact of the no-flow walls.","supporting_citations":[],"review_version":1}