{"id":"d9dc19d1-ec61-4ac2-9606-fe2e0876cdc3","arxiv_id":"2505.02600","paper_version":3,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":6,"one_line_summary":"Stationary envelope models with nuclear burning are validated against MESA and used to map the heat flow between the burning envelope and the crust of accreting neutron stars.","lead":"This paper builds accretion-rate-dependent boundary conditions that let long-term neutron star cooling simulations include heat from thermonuclear burning in the envelope. The resulting maps show when and how much of that heat flows into the crust for persistent and transient accretors.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The MESA validation is run in non-relativistic mode, but the NSCool runs use relativistic boundary tables; the unquantified GR shift in the Tb-Lb relation could change Qb by hundreds of keV.","rationale":"The paper is a serious and transparent contribution: it provides the MESA inlists, an online Lb table, and uses public codes, and it openly lists several limitations. The reader's conditional verdict is reasonable. I focused on a concern that is closely related to but more specific than the reader's convergence worry. The MESA comparison in Section 2 is deliberately non-relativistic, while the boundary conditions used in Sections 3-5 are relativistic. Since Appendix A shows that the omitted GR factors enter the envelope equations at order e^Lambda, and since the central quantitative conclusion depends on the precise location of the Lb=0 crossing, this mismatch deserves a direct check. If the check shows a small shift, the central claim stands; if it shows a large shift, the Qb maps and the 'few hundred keV' bound would need revision. This does not change my agreement with the conditional verdict, because the paper already invites further checks and the concern is testable with its own tools. The reader's weaker-convergence issue is also real: Figure 3's caption states that T at 1e7 g/cm3 is still increasing, so the burst-averaged validation is not fully converged. I did not elevate that to the main attack because the paper's own Figure 4 shows agreement near Lb=0, the physically important case, and because the GR-consistency gap is more sharply defined and more directly connected to the boundary conditions actually used in the NSCool runs.","tokens_in":31081,"tokens_out":8131,"duration_ms":106462,"concrete_test":"Using the published stationary envelope code, generate Tb-Lb tables for Mdot = 1e-10, 1e-9, and 1e-8 Msun/yr with Tb in the Lb<0 region, once including GR and once with the same non-relativistic equations used in the MESA comparison. Quantify Delta Lb/Lb and the shift of the Lb=0 crossing temperature. If the crossing shifts by more than about 1% in Tb, or if Delta Qb exceeds about 100 keV, recompute the NSCool frames with fast neutrino cooling and Qimp=0 using the non-relativistic table; if no new model reaches Qb <= -1 MeV/baryon, the shallow-heating conclusion survives. The online table of Lb values makes this check straightforward.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Section 2 validates the stationary envelope models against MESA only after 'we adapted our stationary envelope code to be non-relativistic as well,' because MESA's GR correction factors are incomplete (Appendix A). The boundary conditions actually implemented in NSCool (Section 3) 'include general relativistic effects that were artificially excluded' from that comparison. Nowhere is the relativistic versus non-relativistic Tb-Lb relation compared. For a 1.4 Msun, 11.56 km star, the metric factors are e^Phi ~ 0.76 and e^Lambda ~ 1.3, and Appendix A shows the missing factors enter the temperature-gradient and luminosity equations at order e^Lambda (Eqs. A4 and A9). The sign and magnitude of Qb are residual differences between envelope nuclear luminosity and crust cooling; even a 20-30% shift in Lb(Tb) can move the Lb=0 crossing temperature and change Qb by hundreds of keV. Thus the central bound 'at most a few hundreds keVs' is not directly covered by the MESA comparison. The authors do flag the first-hours transient limitation and the Qimp fudge factor, but not this validation gap.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper develops a practical method for including thermonuclear burning in the envelope of accreting neutron stars in long-term thermal evolution simulations. The authors compute stationary accreting envelope models with a 380-isotope nuclear network, convert them into mass-accretion-rate-dependent T_b–L_b boundary conditions at ρ_b = 10^7 g cm⁻³, and validate the stationary approximation against time-dependent MESA envelope models, with the caveat that the MESA runs are not fully relaxed. These boundary conditions are then implemented in NSCool and used to explore persistent and transient accretion over a large parameter space spanning neutrino luminosity, accretion rate, shallow heating strength and depth, and crust impurity parameter. The central result is that heat can flow from the envelope into the crust (Q_b < 0), but the leakage reaches at most a few hundred keV per accreted baryon in the most extreme cases, so it cannot by itself explain the MeV-scale shallow heating inferred from observations. The paper also provides useful by-products, including tables of boundary conditions and a discussion of the implications for X-ray burst quenching.","tokens_in":31245,"tokens_out":5286,"duration_ms":69415,"significance":"If the central quantitative bound survives scrutiny, this is a valuable negative result: it separates envelope thermonuclear burning from the class of mechanisms that could produce the observationally required shallow heating, and it provides a reusable tool for future transient-source modeling. The paper's strengths are its transparency, the large parameter scan (over 7,000 models), the independent MESA comparison, the public availability of the NSCool framework and MESA inlists, and the fact that Q_b is an emergent output rather than a fitted quantity. Several limitations are explicitly disclosed by the authors, including incomplete convergence of the MESA runs and the inability of stationary envelope models to describe the first hours of an outburst. The main concern is a verification gap: the validation was performed in a non-relativistic code setup while the production boundary conditions include general-relativistic effects.","major_comments":[{"comment":"The validation of the stationary envelope models is carried out in non-relativistic mode, while the boundary conditions actually used in NSCool include general-relativistic effects. Section 2 states that the stationary code was adapted to be non-relativistic to match MESA, and Section 3 states that the production envelope models include GR effects 'artificially excluded' from the comparison. Appendix A shows that MESA's c_grav implementation misses factors of e^Λ in the temperature-gradient and luminosity equations (Eqs. A4 and A9). The manuscript, however, never compares the relativistic and non-relativistic T_b–L_b relations. Because Q_b is a residual difference between envelope nuclear luminosity and crust thermal flux, a shift in L_b(T_b) of the order implied by e^Λ ≈ 1.3 in a 1.4 M_⊙, 11.56 km model can move the L_b = 0 crossing and change Q_b by hundreds of keV. This gap directly affects the central claim in Section 6 that negative Q_b reaches 'at most a few hundreds keVs'. I request a quantitative comparison of the GR and non-GR T_b–L_b curves (or an equivalent sensitivity analysis) and an assessment of its effect on Figures 6 and 8.","section":"Sections 2, 3, and Appendix A"},{"comment":"The MESA validation is not a true steady-state test. The text states that 'reaching the stationary states would require much more computing time', and the caption to Figure 3 notes that T at ρ = 10^7 g cm⁻³ is still increasing with time. The comparison therefore uses burst averages that are still evolving toward the stationary curves. Figure 4, with the boundary placed at ρ_b = 10^7 g cm⁻³, is reassuring, particularly near L_b = 0, but it is not a quantitative convergence test. Since the stationarity assumption is load-bearing for the method, I ask the authors to estimate the remaining drift (for example, by extrapolating the trend in time or by quoting an implied uncertainty in the derived T_b–L_b relation and in Q_b) so that the possible bias in the boundary conditions can be judged.","section":"Section 2, Figures 3 and 4"}],"minor_comments":[{"comment":"The unit 'Mev' should be 'MeV' throughout the figure labels and captions (e.g., 'Q_b [Mev baryon−1]').","section":"Figures 6, 8, 10–14"},{"comment":"The caption says 'suppressed deep crustal heating according to the model of Haensel & Zdunik (2008)', but the text of Appendix B and the referenced model describe the suppressed heating case of Gusakov & Chugunov (2020); the caption should be corrected.","section":"Figure 10 caption"},{"comment":"The notation eT_8 for the redshifted core temperature is introduced without defining the tilde; use a single, consistently defined symbol such as ~T_8 throughout the text and equations.","section":"Equation (2) and Section 4"},{"comment":"The text states that for a cold crust the inward energy can be 'very significant, actually comparable to the amount of shallow heating', while Section 6 later explains that the self-consistent models reduce this to a few hundred keV. A forward reference or caveat at the first mention would prevent the reader from over-interpreting Figure 5.","section":"Section 3, right panel of Figure 5"},{"comment":"The reference 'Page, Garibay, Nava-Callejas, & Cavecchi 2025' has no journal or publication status; it should be updated or marked as in preparation, especially since it is cited for numerical values of neutrino luminosities.","section":"References"}],"recommendation":"major_revision","confidential_remarks":"The GR verification gap is the main reason for my recommendation. I see no sign of circularity or fabrication in the modeling chain; the issue is a missing numerical comparison between the physics used for validation and the physics used in production. If the authors can show that the relativistic and non-relativistic T_b–L_b relations differ by much less than the spread in Figures 6 and 8, the paper would be suitable for acceptance. I also encourage them to be explicit about the uncertainty from incomplete MESA relaxation, since the current text relies on the reader accepting that the trend in Figure 4 has converged sufficiently."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Dear colleague,\n\nHere's my take on Nava-Callejas, Page, and Cavecchi (arXiv:2505.02600). The paper should be engaged seriously: it gives the community an accretion-rate-dependent Tb-Lb boundary condition that includes thermonuclear burning, and it shows that in self-consistent long-term models the envelope can deliver at most a few hundred keV per baryon into the crust, too little to explain the MeV shallow heating inferred from cooling curves. The authors do not oversell; they explicitly note the earlier watershed results from Hanawa & Fujimoto, Miralda-Escude, Zdunik, and Dohi. What is new is the systematic construction of boundary tables from a 380-isotope network, the MESA burst-average validation, and the wide parameter survey with NSCool across accretion rates, neutrino luminosities, impurity, and shallow-heating choices. That is a useful package, and the public inlist on Zenodo makes it reproducible.\n\nThe soft spots are real but not fatal. The MESA runs do not reach steady state, and the authors admit this; the comparison at rho_b = 10^7 g/cm^3 shows convergence of burst averages to stationary curves, but the temperature is still rising. The more specific worry from the stress test is the GR versus non-relativistic mismatch: the MESA validation is done with the stationary code in non-relativistic mode because of MESA's incomplete GR corrections, while the actual NSCool boundary tables include GR. The authors never compare the two versions of the Tb-Lb relation. Since the metric factors are of order 25-30%, and Qb is a difference between envelope nuclear luminosity and crust cooling, it is plausible the Lb = 0 crossing shifts by tens of percent, which could move Qb by hundreds of keV. That does not undermine the headline claim that envelope leakage still falls short of MeV-level shallow heating, but it does mean the 'at most a few hundred keV' bound is not tightly calibrated. A direct relativistic versus non-relativistic plot of Tb-Lb for a couple of accretion rates would be the obvious addition.\n\nI would bring this to a reading group and would cite it if I worked on NS cooling. It deserves a serious referee; the referee should push for a quantification of the GR shift and a statement of the uncertainty on Qb, but the methodological core is sound.","headline":"A useful and reproducible envelope boundary-condition package that supports the conclusion that thermonuclear leakage can't explain shallow heating, though the GR validation gap leaves the quantitative bound a bit fuzzy.","tokens_in":31890,"tokens_out":2412,"would_cite":true,"duration_ms":29041,"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":"Thermonuclear burning in accreting neutron star envelopes can deposit at most a few hundred keV per baryon into the crust, far below the MeV-scale shallow heating that observations require.","keywords":["accreting neutron stars","thermonuclear burning","shallow heating","crust thermal evolution","X-ray bursts","neutrino cooling","envelope boundary conditions","quiescent cooling"],"falsifier":"Run a time-dependent envelope simulation at a low accretion rate for many burst cycles until the temperature at the boundary density stops rising, then compare the burst-averaged Q_b with the stationary prediction; if the averages fall systematically below the stationary curves, the few-hundred-keV ceiling is an artifact of incomplete convergence rather than a physical limit.","tokens_in":2028,"feed_emoji":"🔥","tokens_out":2530,"duration_ms":137382,"temperature":0.7,"pith_summary":"This paper asks whether the thermonuclear burning of accreted hydrogen and helium on a neutron star's surface can heat the star's interior, and how much. The authors build a set of stationary envelope models that include a full nuclear network and turn them into boundary conditions for long-term simulations of the whole star, validating the stationary approximation against time-dependent burst simulations. Evolving these models for a million years across a broad grid of accretion rates, crust conductivities, shallow-heating depths, and neutrino cooling strengths, they find that heat does sometimes flow from the burning envelope into the crust. In the most extreme cases that inward flow amounts to at most a few hundred keV per accreted baryon, far less than the one to three MeV per baryon of shallow heating that observations of cooling neutron stars typically demand. The paper concludes that thermonuclear leakage is a real but quantitatively limited heat source, active mainly at low accretion rates and in stars whose cores cool rapidly through neutrinos.","feed_headline":"Neutron star envelopes leak too little heat to explain shallow heating","feed_subtitle":"Self-consistent long-term models cap envelope-to-crust leakage at a few hundred keV per baryon.","key_machinery":"The load-bearing object is the stationary accreting envelope model used as a thermal-structural boundary condition at density rho_b = $10^{7}$ g $cm^{-3}$. Each model integrates the envelope inward from the photosphere with a 380-isotope nuclear network, so the resulting relation between boundary temperature T_b and boundary luminosity L_b automatically contains the heat produced by hydrogen, helium, and rp-process burning; the key feature is a temperature inversion in which the luminosity at the base, L_b = Mdot Q_b / m_u, can become negative, meaning thermonuclear energy leaks into the crust. The stationary maps are checked against burst averages of a time-dependent envelope code and then imposed as outer boundary conditions on a neutron star cooling code, letting the star relax to steady state under constant or periodically repeating accretion. The machinery's work is to translate a large parameter space, including accretion rate, outburst duty cycle, crust impurity parameter, shallow-heating strength and depth, and slow versus fast neutrino cooling, into one number: Q_b at the crust-envelope interface.","core_discovery":"The central claim is that a self-consistent treatment of accretion, envelope nuclear burning, and crust/core thermal evolution limits the envelope's contribution to interior heating. Using stationary accreting-envelope models at a boundary density of rho_b = $10^{7}$ g $cm^{-3}$ that include hydrogen, helium, and rp-process burning, the authors construct accretion-rate-dependent relations between the boundary temperature T_b and the boundary luminosity L_b, validated by comparing burst averages from a time-dependent envelope code. When these relations are imposed on a full neutron star cooling code run to steady state, values of Q_b (the heat per accreted baryon crossing the crust-envelope interface) that are negative, meaning heat flows from the envelope into the crust, reach only a few hundred keV per baryon in the most extreme cases. This is far below the 1 to 3 MeV per baryon shallow heating that observations require, so envelope thermonuclear burning cannot by itself explain shallow heating. The sign and magnitude of Q_b are controlled by the competition between core neutrino cooling, which cools the interior and favors inward flow, and accretion rate plus shallow heating, which warm the crust and favor outward flow.","pith_inferences":["Editorial inference: if the stationary-average picture survives true steady-state checks, the missing MeV-scale shallow heat must be sought in non-thermonuclear mechanisms such as interface shear, compositional settling, or deep-crust nuclear reactions beyond the standard deep crustal heating.","Editorial inference: applying the same boundary-condition machinery to helium-rich or metal-poor accreted compositions would shift the burning layers and could move the inward-flux ceiling either up or down.","Editorial inference: the agreement between burst averages and stationary curves suggests the time-averaged envelope structure is set mainly by the mean accretion rate, so stochastic burst-to-burst variability should not change Q_b by more than the few-hundred-keV ceiling; a dedicated long-duration burst simulation could test this.","Testable extension: the published boundary-condition tables can be plugged into independent thermal evolution codes to see whether the inward-flux regions shift when different equations of state or neutrino emission models are adopted."],"forward_implications":["Observed MeV-scale shallow heating in quiescent neutron stars must come from mechanisms other than envelope thermonuclear leakage, since the inward flux is capped near a few hundred keV per baryon.","Inward envelope-to-crust heat flow is mostly confined to low mass accretion rates, roughly below 3 x 10^-9 solar masses per year, and to cores with strong fast-neutrino cooling.","Transiently accreting sources, with colder cores from quiescent periods, are more likely to show heat flowing from the envelope into the crust than persistently accreting sources.","The inward heat flow is self-limiting: the layers just below the envelope heat up within about a day, quenching the temperature inversion and reducing the flux.","The computed maps of crust-to-envelope heat flow provide a diagnostic that can be inverted with observed cooling curves and burst-quenching accretion rates to constrain crust impurity and core neutrino luminosity."],"supporting_citations":[{"why":"Supplies the stationary envelope code and the 380-isotope nuclear network that generate the boundary conditions used throughout the paper.","marker":"Nava-Callejas et al. 2024a"},{"why":"Provides the time-dependent stellar evolution code whose burst simulations are used to validate the stationary envelope approximation.","marker":"Paxton et al. 2011"},{"why":"Provides the deep crustal heating profile, about 1.9 MeV per baryon, used as a heating source in the long-term cooling calculations.","marker":"Haensel & Zdunik 2008"},{"why":"Established the envelope boundary-condition method that this work generalizes to accreting, burning envelopes.","marker":"Gudmundsson et al. 1983"},{"why":"Defined the observational shallow-heating requirement that the paper's results are measured against.","marker":"Brown & Cumming 2009"},{"why":"First identified the temperature inversion with heat flowing from the burning envelope into the interior, the phenomenon modeled here self-consistently.","marker":"Hanawa & Fujimoto 1984"},{"why":"Supplies the normalization ranges for slow and fast neutrino luminosities scanned in the parameter study.","marker":"Ofengeim et al. 2017"},{"why":"Established the long-term steady-state thermal balance that justifies evolving the models for a million years.","marker":"Colpi et al. 2001"},{"why":"Provides the electron-ion and impurity scattering frequencies underlying the crust thermal conductivity treatment.","marker":"Yakovlev & Urpin 1980"}],"fun_headline_variants":["Envelope burning cannot explain shallow neutron star heating","Envelope-to-crust heat capped below shallow heating needs","Accreting envelope leaks too little heat to match observations","Thermonuclear envelope heating falls far short of observations","Neutron star envelope burning yields tiny crustal heat input"],"cache_read_input_tokens":33920,"weakest_assumption_plain":"The entire quantitative map rests on treating the envelope's time-averaged structure as a stationary sequence of models: the time-dependent burst simulations used for validation never reached a true steady state, and the temperature at the boundary density was still rising when those runs stopped.","fun_headline_variants_meta":{"raw":{"variants":["Envelope burning cannot explain shallow neutron star heating","Envelope-to-crust heat capped below shallow heating needs","Accreting envelope leaks too little heat to match observations","Thermonuclear envelope heating falls far short of observations","Neutron star envelope burning yields tiny crustal heat input"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000224,"raw_usage":{"total_tokens":1427,"prompt_tokens":880,"completion_tokens":547,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":496,"completion_tokens_details":{"reasoning_tokens":468}},"tokens_in":496,"tokens_out":547,"duration_ms":6111,"temperature":1.0,"reasoning_tokens":468,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-16T00:47:40.103585+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run a time-dependent envelope simulation at a low accretion rate for many burst cycles until the temperature at the boundary density stops rising, then compare the burst-averaged Q_b with the stationary prediction; if the averages fall systematically below the stationary curves, the few-hundred-keV ceiling is an artifact of incomplete convergence rather than a physical limit.","supporting_citations":[{"cited_title":"D., Fortin , M., Haensel , P., Yakovlev , D","cited_arxiv_id":null,"evidence_quote":"Supplies the normalization ranges for slow and fast neutrino luminosities scanned in the parameter study."},{"cited_title":"2001, , 548, L175, 10.1086/319107","cited_arxiv_id":null,"evidence_quote":"Established the long-term steady-state thermal balance that justifies evolving the models for a million years."}],"review_version":1}