REVIEW 3 major objections 5 minor 51 references
Impact of cosmic-ray propagation on the chemistry and ionisation fraction of dark clouds
T0 review · 3 major / 5 minor · reviewed 2026-08-06 · deepseek-v4-flash
Pith's one-line read Cosmic-ray propagation cuts core ionisation by an order of magnitude
desk verdict Useful, clearly described CR propagator for GIZMO; transport model is unvalidated against full CR solutions and the 'guarantee' in the abstract overstates, but the paper deserves serious review. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The engine is a tracer-particle cosmic-ray propagator. A population of tracer particles is released on a spherical surface around the cloud and stepped with a fourth-order Runge-Kutta integrator; at every step the local gas density and magnetic field are read from the same cubic-spline smoothing kernel used in the magnetohydrodynamics simulation, the magnetic field sets the propagation direction, and the traversed effective column density $N_{\mathrm{eff}}$ is accumulated with a mirroring factor $\sqrt{1 - B(x')\sin^2(\alpha_0)/B_0}$ so that particles whose pitch angle would mirror at $B_{\mathrm{crit}}$ are discarded. The ionisation rate is then obtained from the $f(N_{\mathrm{eff}})$ conversion curves of two interstellar cosmic-ray models and deposited back onto gas particles with kernel weights. The chemistry is evolved in post-processing with a 134-species, 4616-reaction network, and the observable side is carried by the H$_3^+$ column-density proxy $N(\mathrm{H}_3^+)_{\mathrm{proxy}} = N(\mathrm{o\text{-}H}_2\mathrm{D}^+)/(3 R_D)$ with $R_D = N(\mathrm{DCO}^+)/N(\mathrm{HCO}^+)$, leading to the electron-column-density formula $N(e^-)_{\mathrm{proxy}} = N(\mathrm{DCO}^+) + N(\mathrm{HCO}^+) + N(\mathrm{N}_2\mathrm{D}^+) + N(\mathrm{N}_2\mathrm{H}^+) + N(\mathrm{H}_3^+)_{\mathrm{proxy}}$.
What would settle it
A numerical experiment could settle the central claim: run the same magnetised prestellar core with a full cosmic-ray transport solution that includes pitch-angle scattering and energy losses, and compare the resulting ionisation-rate maps with the tracer-particle maps; a systematically different attenuation depth would bias every fitted relation. An observational check is to map the DCO$^+$/HCO$^+$ ratio across a high-mass clump at high resolution and compare the radial gradient with the gradient predicted from the propagated $\zeta_2$ fields; a flat ratio where the model predicts a steep drop would falsify the attenuation picture.
Extended reading notes
Core claim
The central claim is that cosmic-ray attenuation is a first-order chemical effect in dark clouds: when the cosmic-ray ionisation rate $\zeta_2$ is propagated self-consistently along magnetic field lines, it drops by more than an order of magnitude towards the core and differs by about an order of magnitude between the two adopted interstellar cosmic-ray spectra, models H and L. As a result, column densities of the deuterated ions DCO$^+$ and N$_2$D$^+$ deviate substantially from what a constant rate $\zeta_2 = 2.5 \times 10^{-17}\,\mathrm{s}^{-1}$ would give, with non-linear, time-dependent behaviour tied to deuterium fractionation chemistry. The paper further claims that $X(e^-)$ correlates linearly with $\zeta_2$, and that the electron column density can be recovered from the sum of the main observable positive ions plus an H$_3^+$ proxy, with an overall 65\% underestimate (a factor-of-two lower limit). Applying this formula to five high-mass clumps yields $X(e^-)$ between $1.8\times 10^{-8}$ and $2.7\times 10^{-8}$, in line with previous estimates for prestellar regions.
Load-bearing premise
Everything rests on the assumption that cosmic rays can be represented by tracer particles that follow the smoothed magnetic field, with no scattering, diffusion, or reacceleration, and that the conversion of traversed column density into ionisation rate is given by the chosen $f(N_{\mathrm{eff}})$ curves, a choice the paper itself calls arbitrary.
Editorial extensions
If this is right
- If the central claim is correct, cosmic-ray ionisation rates inside dense cores cannot be treated as constant: attenuation lowers $\zeta_2$ by over an order of magnitude toward the core, affecting chemical ages and deuteration histories.
- The analytical proxy formula gives a factor-of-two lower limit on $X(e^-)$ in high-mass clumps, placing observed values near $10^{-8}$ and providing a direct bridge from line observations to ionisation fraction.
- The power-law fits of $\zeta_2$ against $n(\mathrm{H}_2)$, with scatter of about $0.2\text{--}0.25$ dex, let three-dimensional simulations omit explicit cosmic-ray propagation and still assign a reasonable local ionisation rate from local density.
- Deuterated species are the most sensitive diagnostics of cosmic-ray propagation and of the assumed cosmic-ray spectrum, with DCO$^+$ and N$_2$D$^+$ column densities separating models H and L where non-deuterated ions barely do.
Reading between the lines
- If the proxy systematically underestimates electron abundance by about 65\%, observational estimates of $X(e^-)$ in dense star-forming clumps that rely on similar ion sums may be biased low, and correcting for missing ions or charged grains could raise them by up to a factor of two (our inference, not the paper's claim).
- The fitted $\zeta_2$--$n(\mathrm{H}_2)$ slopes are calibrated on one collapsing magnetised core and two cosmic-ray spectra; whether they hold in filaments, high-mass clumps, or different magnetic geometries is untested, and the scatter likely grows with line-of-sight averaging (our inference).
- Because the scheme omits cosmic-ray scattering, which generally lengthens particle paths, adding scattering would probably steepen the attenuation and make the chemistry even more sensitive to magnetic topology; this is a testable prediction for future transport-code comparisons (our inference).
- The same tracer-particle machinery could in principle accept time-dependent cosmic-ray spectra without changing the chemical post-processing, since the ionisation-rate maps are passed to the chemistry as inputs (our inference).
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper presents a numerical framework for computing the cosmic-ray ionisation rate (ζ2) in magnetohydrodynamic simulations of prestellar cores. Cosmic rays are treated as tracer particles that propagate along smoothed magnetic field lines, accumulate an effective column density (Eq. 7), and are converted to ζ2 via the f(Neff) curves of Padovani et al. (2018) for two CR spectra, models H and L (Eq. 8). The resulting ζ2 maps are used in post-processing with the KROME chemical network to study the abundances of HCO+, N2H+, DCO+, N2D+, and o-H2D+, comparing against a constant-ζ2 reference case. The authors derive empirical power-law fits between ζ2 and n(H2) (Eqs. 10–11) and between ζ2 and X(e−) (Eqs. 12–13), and propose an analytical formula (Eq. 16) to estimate the electron column density from observable ionic tracers, which they apply to five high-mass clumps, obtaining X(e−) ≈ 1–3 × 10−8.
Significance. If the underlying CR transport approximation is faithful, the framework offers a computationally inexpensive way to include spatially varying CR ionisation in astrochemical simulations, and the analytical proxy could be useful for interpreting observations. The paper is carefully structured, the convergence tests in Appendices A and B are systematic and show that the fiducial tracer-particle number is adequate, and the internal comparison between models H, L, and C is clear and internally consistent. The explicit functional forms in Eqs. (10)–(13) and (16) are falsifiable predictions that can be tested against future simulations and observations. However, the central result rests on a simplified propagation scheme that is not benchmarked against a full CR transport calculation, and the observational application inherits a known systematic bias that is not propagated into the quoted values.
major comments (3)
- [§2.1, Eqs. (7)–(8)] The CR propagation scheme is not validated against any reference solution for full CR transport. The model follows tracer particles along smoothed magnetic field lines with no scattering, pitch-angle diffusion, or reacceleration, and particles with B > B_crit are discarded from the integral in Eq. (8) rather than reflected. The authors themselves state in §2.1 (step iv) that 'the exact conversion defined by f is arbitrary.' Because the resulting ζ2 maps drive all subsequent results (Figs. 2, 3, and 5, and Eqs. 10–13), the paper's central claim that CR propagation significantly affects the chemistry is not yet demonstrated to be robust. A concrete test would be to compare the ζ2 distribution from the tracer-particle scheme with a Fokker–Planck or Monte Carlo CR transport calculation for the same MHD snapshots, at least in a simplified spherical or 1D geometry; without such a benchmark, the systematic bias from discarding mirrored particles remains unquantified.
- [§3.3, Eq. (16) and Fig. 7] The electron proxy in Eq. (16) is validated on the same simulations that are used to assess it, and the H3+ proxy in Eq. (14) is imported from a same-group calibration (Bovino et al. 2020). This in-sample validation does not establish the proxy's accuracy on independent data. More importantly, the text in §3.3 states that the proxy 'tends to underestimate the true electron abundance by 65% overall,' yet the abstract and §3.4 report X(e−) values from observations without applying or propagating this correction. The reader cannot tell whether the quoted observational values are intended as lower limits (as stated in §3.3) or as best estimates (as implied by the phrase 'reliable estimates' in the abstract). The paper should either apply the correction consistently or clearly report the factor-of-two lower-limit nature in the abstract and in the observational application.
- [§3.2, Eqs. (10)–(13)] The fitting relations between ζ2 and n(H2) and between ζ2 and X(e−) are derived from a single simulated prestellar core (one MHD realization, one set of initial conditions) at three epochs. The scatter quoted (σH = 0.25, σL = 0.22 dex) measures scatter within that single simulation, not variation across cloud masses, magnetic field strengths, or evolutionary states. The conclusions recommend these fits for use in 3D simulations without CR propagation, which is a strong claim based on one realization. The paper should either test the fits on additional cores with different properties or substantially qualify the generality of Eqs. (10)–(13).
minor comments (5)
- [Abstract] The phrase 'guarantees a reliable estimate' overstates the certainty given the acknowledged arbitrary choice of f in §2.1 and the lack of external validation; consider softening to 'provides a computationally efficient estimate.'
- [§2.1, Eq. (4)] The initial effective column density Neff(x0, α0) = 2×10^21 / cos(α0) cm^−2 is introduced without derivation or justification; please explain why this value and the 1/cos(α0) dependence are appropriate for the injection surface.
- [§3.1] The sentence 'In our run, this is the effect of the pitch angles that were not considered in Eq. (8), and it is related to the mirroring effect' is confusing; clarify whether the lowest ζ2 values arise from mirroring (discarding particles) or from the pitch-angle integration in Eq. (8).
- [Fig. 7 caption] 'lineal correlation' should be 'linear correlation.'
- [§3.4] The reported uncertainty of ~30% on X(e−) appears to propagate only the line-intensity uncertainties; please state explicitly whether the 65% systematic underestimate of the proxy is included in this uncertainty or treated separately.
Circularity Check
The CR-propagation/chemistry comparison is self-contained; the only partial circularity is the in-sample, same-group calibration of the electron-abundance proxy (Eqs. 14 and 16), which supports the observational X(e-) claim.
-
self citation load bearing
[Section 3.3, Eqs. (14) and (16); applied in Section 3.4]
"Considering that H3+ is not observable in the dense and cold regions of molecular clouds, we replaced it with the proxy provided by Bovino et al. (2020), N(H3+)proxy = 1/3 N(o-H2D+)/RD ... and we tested the reliability of the following formula that can also be applied to observational data: N(e-)proxy = N(DCO+) + N(HCO+) + N(N2D+) + N(N2H+) + N(H3+)proxy. (16) ... The analytical formula reported in Eq. (16) was validated by employing a simulation of a typical prestellar core..."
The 1/3 coefficient in Eq. (14) is imported from a prior article whose authors overlap with the present paper (Bovino et al. 2020; co-authors Bovino and Lupi here). The validation of Eq. (16) is carried out on the same MHD prestellar-core simulation and the same chemical network (Bovino et al. 2019; Ferrada-Chamorro et al. 2021) used for the rest of this study. Therefore the agreement shown in Fig. 7 and the claim in Section 3.4 that the formula 'provides reliable estimates' of X(e-) are in-sample confirmations of a same-group calibration rather than an independent prediction. The observational numbers themselves are external, so the circularity is partial; the paper's main comparison of CR propagation models H, L, and C does not depend on this proxy.
full rationale
The core derivation chain — tracer-particle CR transport (Eqs. 4–9), the ζ2 maps, and the H/L/C chemical comparison — is self-contained: ζ2 is computed from an explicitly adopted f(Neff) relation and the simulation's own kernel-weighted fields, and the deuterium-chemistry result is a model-to-model difference, not a fitted tautology. Equations (10)–(13) are presented transparently as fits to the simulated distributions, so their use as fitting formulas is not circular. The only load-bearing same-group element is the H3+ proxy used to build the electron-abundance estimator (Eqs. 14 and 16); because that proxy is imported from a citation with overlapping authorship and is validated on the same simulation family that produced it, the 'reliable estimate' wording for X(e-) in the abstract is partly in-sample. This does not invalidate the central CR-propagation claim, but it lowers the independent-evidence value of the observational X(e-) application.
Assumptions & free parameters
free parameters (7)
- initial effective column density =
2e21 cm^-2 / cos(alpha0)
- power-law slope and intercept, zeta2 vs n(H2), model H =
-0.22 and -15.06
- power-law slope and intercept, zeta2 vs n(H2), model L =
-0.12 and -16.45
- power-law slope and intercept, zeta2 vs X(e-), model H =
0.41 and -13.06
- power-law slope and intercept, zeta2 vs X(e-), model L =
0.32 and -14.53
- H3+ proxy coefficient in Eq. (14) =
1/3
- electron proxy systematic bias =
proxy is 65% lower than true (about 35% of true)
assumptions (5)
- domain assumption The M1 chemical network (134 species, 4616 reactions) from Bovino et al. (2019) accurately describes prestellar core chemistry.
- domain assumption The f(Neff) conversion curves from Padovani et al. (2018) models H and L are valid for the simulated cloud.
- ad hoc to paper CRs propagate along magnetic field lines without scattering or diffusion and are only attenuated by column density; reflected particles above B_crit are discarded.
- domain assumption Electron abundance equals the sum of the included positive ion abundances (charge neutrality, no charged grains).
- domain assumption The Bonnor-Ebert core simulation of Bovino et al. (2019) is a typical prestellar core.
Cite this review
Pith. "Pith review of Impact of cosmic-ray propagation on the chemistry and ionisation fraction of dark clouds." pith.science (2026). https://pith.science/paper/VGNVIRMX
@misc{pith2026250703832,
author = {Pith},
title = {Pith review of: Impact of cosmic-ray propagation on the chemistry and ionisation fraction of dark clouds},
year = {2026},
howpublished = {\url{https://pith.science/paper/VGNVIRMX}},
note = {Machine review of arXiv:2507.03832}
}
abstract
A proper modelling of the cosmic-ray ionisation rate within gas clouds is crucial to describe their chemical evolution accurately. However, this modelling is computationally demanding because it requires the propagation of cosmic rays throughout the cloud over time. We present a more efficient approach that simultaneously guarantees a reliable estimate of the cosmic-ray impact on the chemistry of prestellar cores. We introduce a numerical framework that mimics the cosmic-ray propagation within gas clouds and applies it to magnetohydrodynamic simulations performed with the code GIZMO. It simulates the cosmic-ray attenuation by computing the effective column density of H$_2$ that is traversed, which is estimated using the same kernel weighting approach as employed in the simulation. The obtained cosmic-ray ionisation rate is then used in post-processing to study the chemical evolution of the clouds. We found that cosmic-ray propagation affects deuterated and non-deuterated species significantly and that it depends on the assumed cosmic-ray spectrum. We explored correlations between the electron abundance, the cosmic-ray ionisation rate, and the abundance of the most relevant ions (HCO$^+$, N$_2$H$^+$, DCO$^+$, N$_2$D$^+$, and o-H$_2$D$^+$), with the purpose of finding simple expressions that link them. We provide an analytical formula to estimate the ionisation fraction, X(e$^-$), from observable tracers and applied it to existing observations of high-mass clumps. We obtained values of about 10$^{-8}$, which is in line with previous works and with expectations for dense clouds. We also provide a linear fit to calculate the cosmic-ray ionisation rate from the local H$_2$ density, which is to be employed in three-dimensional simulations that do not include cosmic-ray propagation.
Figures
Figures from the paper (3 more)
Reference graph
Works this paper leans on
-
[1]
Bergin, E. A., Plume, R., Williams, J. P., & Myers, P. C. 1999, ApJ, 512, 724
work page 1999
-
[2]
Bergin, E. A. & Tafalla, M. 2007, ARA&A, 45, 339
2007
- [3]
-
[4]
Bovino, S., Ferrada-Chamorro, S., Lupi, A., et al. 2019, ApJ, 887, 224
work page 2019
-
[5]
Bovino, S., Ferrada-Chamorro, S., Lupi, A., Schleicher, D. R. G., & Caselli, P. 2020, MNRAS, 495, L7
2020
-
[6]
& Zweibel, E
Bustard, C. & Zweibel, E. G. 2021, ApJ, 913, 106
2021
-
[7]
Caselli, P., Sipilä, O., & Harju, J. 2019, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 377, 20180401
work page 2019
-
[8]
M., Terzieva, R., & Herbst, E
Caselli, P., Walmsley, C. M., Terzieva, R., & Herbst, E. 1998, ApJ, 499, 234
1998
Show all 51 references
-
[9]
2014, ApJ, 790, L1
Ceccarelli, C., Dominik, C., López-Sepulcre, A., et al. 2014, ApJ, 790, L1
2014
-
[10]
K., Kereš, D., Hopkins, P
Chan, T. K., Kereš, D., Hopkins, P. F., et al. 2019, MNRAS, 488, 3716
2019
-
[11]
2006, Proceedings of the National Academy of Science, 103, 12269
Dalgarno, A. 2006, Proceedings of the National Academy of Science, 103, 12269
2006
-
[12]
& Lepp, S
Dalgarno, A. & Lepp, S. 1984, ApJ, 287, L47
1984
-
[13]
2021, MNRAS, 505, 3442
Ferrada-Chamorro, S., Lupi, A., & Bovino, S. 2021, MNRAS, 505, 3442
2021
-
[14]
Gaches, B. A. L., Bisbas, T. G., & Bialy, S. 2022, A&A, 658, A151
2022
-
[15]
2017, A&A, 603, A33
Giannetti, A., Leurini, S., Wyrowski, F., et al. 2017, A&A, 603, A33
2017
-
[16]
2014, A&A, 570, A65
Giannetti, A., Wyrowski, F., Brand, J., et al. 2014, A&A, 570, A65
2014
-
[17]
Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, MNRAS, 439, 2386 Güsten, R., Nyman, L. Å., Schilke, P., et al. 2006, A&A, 454, L13
2014
-
[18]
Hopkins, P. F. 2015, MNRAS, 450, 53
2015
-
[19]
Hopkins, P. F. & Raives, M. J. 2016, MNRAS, 455, 51
2016
-
[20]
R., Oka, T., & McCall, B
Indriolo, N., Geballe, T. R., Oka, T., & McCall, B. J. 2007, ApJ, 671, 1736
2007
-
[21]
& McCall, B
Indriolo, N. & McCall, B. J. 2012, ApJ, 745, 91
2012
-
[22]
& McCall, B
Indriolo, N. & McCall, B. J. 2013, Chemical Society Reviews, 42, 7763
2013
-
[23]
A., Gerin, M., et al
Indriolo, N., Neufeld, D. A., Gerin, M., et al. 2012, ApJ, 758, 83
2012
-
[24]
2022, ApJ, 939, 102
Li, S., Sanhueza, P., Lu, X., et al. 2022, ApJ, 939, 102
2022
-
[25]
& Bergin, E
Maret, S. & Bergin, E. A. 2007, ApJ, 664, 956
2007
-
[26]
2011, A&A, 526, A47
Maret, S., Hily-Blant, P., Pety, J., Bardeau, S., & Reynier, E. 2011, A&A, 526, A47
2011
-
[27]
2020, A&A, 634, A115
Miettinen, O. 2020, A&A, 634, A115
2020
-
[28]
Neufeld, D. A. & Wolfire, M. G. 2017, ApJ, 845, 163
2017
-
[29]
V ., Silsbee, K., et al
Obolentseva, M., Ivlev, A. V ., Silsbee, K., et al. 2024, ApJ, 973, 142
2024
-
[30]
A., Hanasz, M., & Wólta´nski, D
Ogrodnik, M. A., Hanasz, M., & Wólta´nski, D. 2021, ApJS, 253, 18
2021
-
[31]
2022, A&A, 658, A189
Padovani, M., Bialy, S., Galli, D., et al. 2022, A&A, 658, A189
2022
-
[32]
& Gaches, B
Padovani, M. & Gaches, B. 2024, in Astrochemical Modeling: Practical Aspects of Microphysics in Numerical Simulations. Edited by Stefano Bovino and Tommaso Grassi (Elsevier), 189–231
2024
-
[33]
2013, A&A, 560, A114
Padovani, M., Hennebelle, P., & Galli, D. 2013, A&A, 560, A114
2013
-
[34]
V ., Galli, D., & Caselli, P
Padovani, M., Ivlev, A. V ., Galli, D., & Caselli, P. 2018, A&A, 614, A111
2018
-
[35]
V ., Galli, D., et al
Padovani, M., Ivlev, A. V ., Galli, D., et al. 2020, Space Sci. Rev., 216, 29
2020
-
[36]
Patil, A., Huard, D., & Fonnesbeck, C. J. 2010, Journal of Statistical Software, 35, 1
2010
-
[37]
2024, A&A, 685, A67
Redaelli, E., Bovino, S., Lupi, A., et al. 2024, A&A, 685, A67
2024
-
[38]
2021, A&A, 656, A109
Redaelli, E., Sipilä, O., Padovani, M., et al. 2021, A&A, 656, A109
2021
-
[39]
2020, A&A, 644, A34
Sabatini, G., Bovino, S., Giannetti, A., et al. 2020, A&A, 644, A34
2020
-
[40]
2023, ApJ, 947, L18
Sabatini, G., Bovino, S., & Redaelli, E. 2023, ApJ, 947, L18
2023
-
[41]
2024, A&A, 692, A265
Sabatini, G., Bovino, S., Redaelli, E., et al. 2024, A&A, 692, A265
2024
-
[42]
M., Contreras, Y ., et al
Schuller, F., Menten, K. M., Contreras, Y ., et al. 2009, A&A, 504, 415
2009
-
[43]
J., Srianand, R., et al
Shaw, G., Ferland, G. J., Srianand, R., et al. 2008, ApJ, 675, 405
2008
-
[44]
N., Bergner, J
Shingledecker, C. N., Bergner, J. B., Le Gal, R., et al. 2016, ApJ, 830, 151
2016
-
[45]
V ., Padovani, M., & Caselli, P
Silsbee, K., Ivlev, A. V ., Padovani, M., & Caselli, P. 2018, ApJ, 863, 188
2018
-
[46]
Snow, T. P. & McCall, B. J. 2006, ARA&A, 44, 367
2006
-
[47]
2024, A&A, 687, A70
Socci, A., Sabatini, G., Padovani, M., Bovino, S., & Hacar, A. 2024, A&A, 687, A70
2024
-
[48]
2005, MNRAS, 364, 1105
Springel, V . 2005, MNRAS, 364, 1105
2005
-
[49]
2021, MNRAS, 503, 2242
Thomas, T., Pfrommer, C., & Pakmor, R. 2021, MNRAS, 503, 2242
2021
-
[50]
S., Wells, M
Urquhart, J. S., Wells, M. R. A., Pillai, T., et al. 2022, MNRAS, 510, 3389
2022
-
[51]
2019, MNRAS, 488, 2235 Article number, page 8 of 10 G
Winner, G., Pfrommer, C., Girichidis, P., & Pakmor, R. 2019, MNRAS, 488, 2235 Article number, page 8 of 10 G. Latrille et al.: Impact of cosmic-ray propagation on the chemistry and ionisation fraction of dark clouds Appendix A: Convergence tests Because of our choice of kernel...
2021
Reviewed August 6, 2026 · model on record in the stance chip above.
Discussion (0). Sign in to comment.