© The Authors 2025

1 Introduction

Dark clouds are commonly characterised by extreme conditions, such as temperatures below 20 K and high column densities that span 10²³−10²⁵ cm⁻² (Bergin & Tafalla 2007). The dense material is sufficient to completely shield the cloud from incoming photo-radiation, which leaves only cosmic rays (CRs) as a potential source of ionisation in the inner regions of the clouds. Cosmic rays are non-thermal charged particles whose origin depends on their energies, which extend from local sources such as stars (10⁹−10¹⁰ eV), galactic sources such as supernova remnants (10¹⁰−10¹⁵ eV), and extra-galactic sources such as active galactic nuclei (E > 10¹⁵ eV; see Padovani et al. 2020, as a review).

Within dark clouds, hydrogen is mainly in its molecular form, H₂, and it represents the major constituent, followed by CO and then other species. Because H₂ dominates the other species, CRs are expected to have a stronger impact on it and first form ${\mathrm{H}}_2^{+}$ and subsequently ${\mathrm{H}}_3^{+}$. After ${\mathrm{H}}_3^{+}$ is formed, the path to molecular formation is paved. It leads to the production of other molecules such as the OH⁺ and HCO⁺ ions (Dalgarno 2006; Indriolo & McCall 2013), the molecular form of atomic species such as O₂ and N₂ (Indriolo & McCall 2013), complex neutrals, H₂O and NH₃ (Gaches et al. 2022), and deuterated species such as H₂D⁺ (Caselli et al. 2019).

The latter is formed through the reaction of HD with ${\mathrm{H}}_3^{+}$ (e.g. Dalgarno & Lepp 1984). When it is formed, the path to deuteration begins. Because H₂D⁺ donates its deuterium atom to abundant neutral species such as CO, DCO⁺ is produced. At the same time, CO forms HCO⁺ via interaction with ${\mathrm{H}}_3^{+}$, $$\mathrm{C}\mathrm{O}+{\text{o-H}}_2{\mathrm{D}}^{+}⟶{\mathrm{D}\mathrm{C}\mathrm{O}}^{+}+{\mathrm{H}}_2,$$(1) $$\mathrm{C}\mathrm{O}+{\mathrm{H}}_3^{+}⟶{\mathrm{H}\mathrm{C}\mathrm{O}}^{+}+{\mathrm{H}}_2.$$(2)

The formation of N₂H⁺ is similar to that of HCO⁺, $${\mathrm{N}}_2+{\mathrm{H}}_3^{+}⟶{\mathrm{N}}_2{\mathrm{H}}^{+}+{\mathrm{H}}_2.$$(3)

Ionisation occurs through collisions between CRs and neutral molecules (H₂, N₂, and CO). This drives fast ion-neutral reactions in the environment. The number of ionisations of H₂ molecules (or H atoms) per unit time, referred to as the CR ionisation rate, therefore plays a fundamental role in determining the chemical state of star-forming regions. Throughout this paper, we indicate the ionisation rate for molecular hydrogen as ζ₂.

Several studies have attempted to determine the CR ionisation rate from observational tracers in different environments so far, including diffuse medium (Shaw et al. 2008; Indriolo & McCall 2012), high-mass star-forming regions (Sabatini et al. 2020, 2023), circumstellar discs (Ceccarelli et al. 2014), and low-mass dense cores (Caselli et al. 1998; Redaelli et al. 2021; Bialy et al. 2022). To do this, they relied on relations between the species abundances and ζ₂ (Caselli et al. 1998; Indriolo et al. 2007, 2012; Indriolo & McCall 2012; Bovino et al. 2020; Socci et al. 2024), or they explored the chemical evolution with low-dimensional (0D or 1D) models (Shaw et al. 2008; Shingledecker et al. 2016; Neufeld & Wolfire 2017; Redaelli et al. 2021). From a theoretical point of view, it is highly computationally expensive to propagate CRs in simulations alongside hydrodynamics because the particle velocities are high (close to the speed of light). This reduces the simulation time-step, the extended energy spectrum, and several uncertainties in the CR propagation process itself. Some studies have attempted to self-consistently couple CRs with magnetohydrodynamics only recently (MHD; Chan et al. 2019; Winner et al. 2019; Bustard & Zweibel 2021; Ogrodnik et al. 2021; Thomas et al. 2021), but with several limitations, and never alongside a proper nonequilibrium chemical evolution.

We explore the effect of a consistent CR propagation on the chemistry of prestellar cores simulated with the magnetohydrodynamics code GIZMO (Hopkins 2015; Hopkins & Raives 2016), which was previously reported by Bovino et al. (2019). The paper is organised as follows: In Section 2 we describe the numerical framework, that is, the CR propagator and the coupling of the CR ionisation rate to the chemical post-processing. In Section 3 we discuss our main results, and in Section 4 we draw our conclusions.

2 Numerical framework

Our numerical framework is composed of two parts, a CR propagation code that we describe in Section 2.1, and a chemical post-processing tool that we present in Section 2.2. The tool uses the output of the CR propagation code as input.

2.1 CR propagation scheme

The propagation scheme we considered aimed to mimic a continuous CR flux in which CRs are modelled as tracer particles (TPs) that propagate throughout the system (e.g. a molecular cloud). Because CRs are charged particles, they propagate along paths defined by the magnetic field lines. As CRs propagate through the cloud, they lose energy because of collisions with molecular gas. This results in an effective attenuation of the flux deeper in the cloud at increasingly higher column densities (Padovani & Gaches 2024).

The structure of the propagation scheme is as follows:

• (i)

  TP generation: A specific number of TPs are created (determined by 2^p, with p set before the run) at a distance R from the cloud, and they are distributed randomly over a spherical surface. The initial direction of the TPs along the magnetic field lines is chosen to be towards the centre of the cloud, and all are assumed to start from a uniform isotropic distribution of pitch angles (α₀) over a hemisphere, where the pitch angle distribution ranges from 0 to π/2, sampled using N_α bins. For each initial value of the pitch angle, TPs are assigned an initial column density $$N_{\mathrm{e}\mathrm{f}\mathrm{f}}({\mathbf{x}}_0,{\alpha }_0)=\frac{2\times {10}^{21}{\mathrm{c}\mathrm{m}}^{-2}}{\cos {\alpha }_0},$$(4)

  where N_eff stands for the effective column density that is traversed by the TP, 2 × 10²¹ cm⁻² is the column density value at which hydrogen is mainly in its molecular form (Snow & McCall 2006), and x₀ is the initial position of the TP.

• (ii)

  Gas particle interpolation: In order to compute the column density that is traversed by CRs and then the CR ionisation rate, local gas properties at the position of the TPs are needed, in particular, the H₂ number density n(H₂) and the average magnetic field B, the latter setting the propagation direction. Because all these physical quantities vary throughout the cloud, we estimated them by means of kernel weighting, analogous to what is done for hydrodynamics (see Hopkins 2015). In particular, at the location x_k of the k-th TP, we determined the kernel size h_k that encompasses an effective number of gas-neighbours N_ngb = 32, defined by $$N_{\mathrm{n}\mathrm{g}\mathrm{b}}=\frac{4\pi }{3}{h}_k^3\sum _{j}W({\mathbf{x}}_j-{\mathbf{x}}_k,h_k),$$(5)

  where the sum is over all neighbours within h_k, x_j is the location of the jth neighbour, and the kernel function W(x_j − x_k, h_k) is defined by the standard cubic spline $$\overline{W}(q,h_k)=\frac{8}{\pi h_k^3}\{\begin{array}{ll}1-6q^2+6q^3,& 0\le q\le \frac{1}{2},\\ 2(1-q{)}^3,& \frac{1}{2}<q\le 1,\\ 0,& q>1,\end{array}$$(6)

  with q = |x_j − x_k|/h_k the normalised distance.

  The neighbour search is repeated at each step of the integration for all TPs and thus represents the highest computational cost of the scheme because Eq. (5) must be solved implicitly. To limit the computational cost and rapidly converge to the desired solution, we performed this operation using a Newton-Raphson iteration¹ (Springel 2005; Hopkins 2015). After the kernel size was determined, a last loop on the neighbours was performed to compute the smoothed magnetic field and the gas density, which were then used to update the column density.

• (iii)

  TPs propagation: The propagation of TPs is computed by means of a fourth-order Runge-Kutta integrator, where the step size for the kth TP is determined as Δℓ_k = 0.1Δx_k, with $\mathrm{\Delta }{x}_k=(4\pi /3{)}^{1/3}{N}_{\text{ngb\hspace{0.17em}}}^{-1/3}{h}_k$ the grid-equivalent inter-particle spacing, and motion occurs along the direction of the smoothed magnetic field. We note that the interpolated quantities were recomputed at every sub-step of the propagation by repeating (ii).

• (iv)

  Column density and CR ionisation rate update: Using the new location of the TP x′ and the new interpolated quantities n(x′) and B(x′), we updated the column density that is traversed by CRs as $$N_{\mathrm{e}\mathrm{f}\mathrm{f}}({\mathbf{x}}^{\mathrm{\prime }},{\alpha }_0)=N_{\mathrm{e}\mathrm{f}\mathrm{f}}(\mathbf{x},{\alpha }_{\mathbf{0}})+\frac{n\left({\mathbf{x}}^{\mathrm{\prime }}\right)\mathrm{\Delta }{\ell }_k}{\sqrt{1-B_k\left({\mathbf{x}}^{\mathrm{\prime }}\right){\sin }^2\left({\alpha }_0\right)/B_{k,0}}},$$(7)

  and the corresponding CR ionisation rate was re-estimated by integrating with the trapezoidal rule over all the pitch angles as $${\zeta }_2\left({\mathbf{x}}^{\mathrm{\prime }}\right)=\delta \left({\mathbf{x}}^{\mathrm{\prime }}\right){\int }_0^{\pi /2}f\left[N_{\mathrm{e}\mathrm{f}\mathrm{f}}({\mathbf{x}}^{\mathrm{\prime }},\alpha )\right]\sin \alpha \mathrm{d}\alpha ,$$(8)

  Fig. 1 Sketch of the CR propagation scheme based on the structure given in Section 2.1. The left panel describes the TP propagation through the domain. TPs (grey dots) are created from a spherical surface of radius R (dashed black circle) with radial inward trajectories (step i). TPs are then deactivated when the column density that is traversed exceeds 1029 cm−2 or when they leave the outer spherical surface at a distance Rout (solid black circle), as shown by the crossed grey dots. The right panel instead highlights the propagation of a single TP particle. The trajectory is shown by the curved solid blue arrow, and the magnetic field lines (B) are reported as dash-dotted grey lines. The pitch angle (α) between the TP velocity and the local magnetic field is shown as the dashed red line. The displacement of the TP between the previous (x) and the current (x′) location is defined by the step size Δℓ (step iii). To compute the effective column density of the TP, the average local density (n(x′)) and local magnetic field (B(x′)) were estimated from gas particles (pink dots) within the kernel size h (thin yellow arrows; step ii). Then, its value was used to update the CR ionisation rate (ζ2; step iv), which was finally deposited at the location of the original gas particle (thick yellow arrows; step v).

  where $\delta \left({\mathbf{x}}^{\mathrm{\prime }}\right)=\frac{B\left(x^{\mathrm{\prime }}\right)}{B_0}$ accounts for the increase in ζ₂ due to the growth in the magnetic flux, with B(x′) the magnetic field at the current position and B₀ the field at the TP injection, sin α accounts for the spherical distribution of the pitch angles (see Eq. (6) in Padovani et al. 2013) at the TP injection, and the f function converts N_eff into the corresponding ζ₂ according to the relations presented in Fig. F.1 of Padovani et al. (2018), where the three curves represent different levels of ζ₂ and account for the ionisation of H₂ by protons and primary and secondary electrons. The exact conversion defined by f is arbitrary and was chosen at the beginning of the calculation. We considered model ℋ (the average CR ionisation rate on diffuse media²) and model ℒ in particular, which comes from the data of the two Voyager spacecraft (see Padovani et al. 2022, for more details).

  In our set-up, the argument of the square root in Eq. (7) could become negative or equal to zero. When this occurs, the gyrocenter of the particle is mirrored backwards (e.g., Silsbee et al. 2018). It is therefore possible to define B_crit = B₀/ sin²(α₀) as the upper limit of the non-mirroring particles. The propagation scheme includes this effect by assuming that the reflected particles above B_crit are neglected in the evaluation of Eq. (8).

• (v)

  CR ionisation rate deposition: At the end of each iteration of the propagation algorithm, the ζ₂ value for each kth TP is deposited on the jth gas particles that contributed to its calculation, again weighted by the smoothing kernel. Because the same gas particles can contribute to the calculation of ζ₂ at different iterations, the actual value of ζ₂ for each gas particle was computed only at the end of the entire procedure as $${\zeta }_2^j=\frac{\sum _{k=1}^{N_p}W({\mathbf{x}}_j-{\mathbf{x}}_k,h_k){\zeta }_2\left({\mathbf{x}}_k\right)}{\sum _{k=1}^{N_p}W({\mathbf{x}}_j-{\mathbf{x}}_k,h_k)}.$$(9)

Steps (ii) to (iv) were repeated over all active TPs, where a TP was considered active until one of the following two conditions was fulfilled: a) the CR flux becomes negligible, which in our scheme corresponds to reaching a traversed column density N_eff > 10²⁹ cm⁻², or b) TPs leave the spherical region around the cloud, defined by an outer boundary spherical surface, which in our case was located at R_out = 1.02 max_j{r_j}, where r_j is the distance of the jth gas particle from the centre of the cloud. The calculation ended when all the TPs became inactive, namely, they completed their propagation throughout the cloud.

The final output of our propagation algorithm, a sketch of which is shown in Fig. 1, consisted of a 3D distribution of ζ₂ generated from the magnetised cloud at a specific time during its evolution. By taking different time steps from the simulation, it is then possible to obtain the evolution of the ζ₂ distribution over time as the gas evolves and to use it as an input to the chemical post-processing, as explained in the next section.

2.2 Chemical post-processing

The output of the CR propagator can be used to chemically evolve the cloud with time-dependent chemistry over the timescale of the simulation. The calculations were performed via the post-processing tool developed by Ferrada-Chamorro et al. (2021), which was coupled to the thermo-chemistry library KROME (Grassi et al. 2014). The chemical network we employed includes 134 species and 4616 reactions (for more details, see the ‘M1’ case and Section 2 of Bovino et al. 2019) and was evolved throughout the entire history of the cloud. We obtained the ζ₂ distribution at each timestep, considering models ℋ, ℒ, and an additional reference case (ℒ), where the CR ionisation rate was kept constant to a typical value of ζ₂ = 2.5 × 10⁻¹⁷ s⁻¹ (hereafter ζ_c). We used ζ_c to compare the impact of a consistent propagation of CRs on several species (HCO⁺, N₂H⁺, DCO⁺, N₂D⁺, and o-H₂D⁺) that are commonly used as tracers for ζ₂. We first determined the value of ζ₂ every three snapshots of the simulation, corresponding to a time interval of 1.5 kyr. This was done to decrease the computational costs. Following Ferrada-Chamorro et al. (2021), we only considered 10% of all gas particles and interpolated the chemical abundances for the remaining ones. The subset particles were randomly selected according to the procedure and the convergence tests reported by Ferrada-Chamorro et al. (2021), which guarantee a proper ratio of the core and background particles.

Fig. 2 Time evolution of the density-weighted ζ2 for models ℋ (top panels) and ℒ (bottom panels) in a sphere of 0.5 pc radius, inside which CRs are propagated. The green arrows correspond to the x-y projection of the magnetic field lines. The two different models yield a difference of one order of magnitude within the high-density region.

3 Results

3.1 Distribution of the CR ionisation rate

Before investigating the chemical evolution in the core, we analysed the distribution of ζ₂ resulting from our propagator. The code was applied to each snapshot of the simulation within a spherical region of radius 0.5 pc centred on the highest-density region of the cloud.

Figure 2 shows density-weighted maps of ζ₂ for ℋ (top panels) and ℒ (bottom panels) at three different times: 30.7, 95.0, and 153.3 kyr. The x–y projection of the magnetic field lines is overplotted as green arrows. ζ₂ differs by approximately one order of magnitude between the two models and spans a range of 0.6 dex. The interval below 10^−16.2 s⁻¹ (ℋ) and 10^−17.1 s⁻¹ (ℒ) covers the high-density region, as expected by considering the CR attenuation, and the highest CR ionisation rate is found in the background region where CRs are almost unattenuated. Strong twisting and bending of the magnetic field lines is also observed, especially in the highest-density regions, where the line density (hence the magnetic field intensity) also increases. This complex pattern results in a significant increase in the effective column density that is traversed by CRs, and hence, in a stronger attenuation.

The time evolution of the ζ₂ distribution for the two models in a 0.2 × 0.2 × 0.2 pc³ box centred at the peak of the high-density region is shown in Fig. 3. Both distributions show that regardless of the initial ζ₂ computed at injection, the attenuation reduces the initial ionisation rate by more than one order of magnitude.

In the intervals 10^−16.4 to 10^−15.8 s⁻¹ (ℋ) and 10^−17.3 to 10^−16.7 s⁻¹ (ℒ), the ζ₂ distribution does not vary over time, but for the low- and high-end tails. In particular, the low-ζ₂ tail, which is associated with the innermost region of the cloud (namely only about 1% of the total mass) simply fluctuates over time without any clear evolutionary trend. When the lowest ionisation rates per distribution, ~10^−17.1 (ℋ) and ~10⁻¹⁸ s⁻¹ (ℒ), are compared with Fig. F.1 of Padovani et al. (2018), these ionisation rate values correspond to column densities of ~10^25–26 cm⁻². 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 and not to a high column density accumulated by TPs.

Fig. 3 Time evolution of the distribution of ζ2 for the gas particles in the simulation for models ℋ (top) and ℒ (bottom). The distributions are determined from a box of 0.2 × 0.2 × 0.2 pc3, which only includes particles from the simulated cloud, i.e., with identical mass. The rightmost edge of the distribution shrinks over time as the gas particles within the box become denser. The low-end tail of the distribution expands towards lower ionisation rates, which is consistent with the evolution of the gas density.

3.2 ζ₂ dependence on local properties

Because the self-consistent propagation of CRs in numerical simulations is expensive, effective models that relate the gas density to ζ₂ are often employed. With our framework, we have the unique ability to determine more accurate semi-analytic prescriptions and how these prescriptions might depend on the CR model that is considered. The attenuation of ζ₂ is expected to significantly depend on the gas density, under the assumption that high volume densities correspond higher column densities (Sabatini et al. 2023; Socci et al. 2024). In Fig. 4 we explore this dependence by showing the correlation between the CR ionisation rate and the number density of H₂, where the viridis histogram corresponds to the ℋ model and the magma histogram to the ℒ model. The trend is clearly negative in both models, which we fit with power laws. The power law is shown as solid and dashed red lines and given by the following expressions: $${\log }_{10}{\zeta }_2^{\mathcal{H}}=-0.22{\log }_{10}n({\mathrm{H}}_2)-15.06,$$(10) $${\log }_{10}{\zeta }_2^{\mathcal{L}}=-0.12{\log }_{10}n({\mathrm{H}}_2)-16.45,$$(11)

for models ℋ and ℒ respectively. Due to the 3D structure of the system, the different lines of sight result in a slightly different CR ionisation rate, with a scatter of σ^ℋ = 0.25 and σ^ℒ = 0.22 for the two models (shown as shaded grey areas), calculated as the root mean square (RMS) between the actual ζ₂ and our best fits.

Fig. 4 Correlation between ζ2 and n(H2) at 153.3 kyr. The viridis and magma color maps represent the ℋ and ℒ models, respectively. The solid and dashed red lines correspond to our best-fit relations, and the shaded grey areas show the intrinsic scatter discussed in Section 3.2.

3.3 Ion chemistry and physical correlations

We present in Fig. 5 the radial profiles obtained from column density maps of the ions HCO⁺, N₂H⁺, DCO⁺, N₂D⁺, and o-H₂D⁺, together with ζ₂ for the inner 0.1 pc of the simulated prestellar core. The profiles display the ℋ and ℒ models and ℒ at time 153.3 kyr. The column densities were calculated by integrating the gas cloud along a line of sight that was aligned with the z-axis over a depth of 0.2 pc. In order to exclude background gas outside the initial Bonnor-Ebert sphere representing our core, we removed gas particles with number densities below log[n/(cm⁻³)] < −3.78, with n calculated as ρ(2.3m_p)⁻¹ with m_p the proton mass. The radial profiles were then computed by averaging the integrated quantities in circular annuli. The radial profiles (panel f) clearly show the effect of the CR attenuation in the central part (highest density). This brings the ionisation rate from model ℋ close to ζ_c. As a consequence, for deuterated species in particular (panels b, c, and e), the column densities of the analysed tracers approach each other. The ℒ model shows more marked differences, in particular, for deuterated species such as DCO⁺ and N₂D⁺. This is due to the efficiency of the formation channels. In particular, after checking the reaction fluxes in the different cases, we found that reactions that are directly dependent on ζ₂ and on (multiple) deuterated species (e.g. H₂D⁺, D₂H⁺, and ${\mathrm{D}}_2^{+}$) become much more relevant for high ζ₂ values by boosting N₂D⁺ and even more DCO⁺, as a consequence of the CO + N₂D⁺ → N₂ + DCO⁺ reaction. Overall, the effect of the different CR ionisation rate spectra is highly non-linear because of the attenuation effects and the different chemical evolution of these species over time.

We show in Fig. 6 the correlation between the abundance of electrons X(e⁻)and ζ₂ over time for the two models with variable CR ionisation rates. The abundance was calculated from column densities obtained by integrating the gas cloud properties along the z-axis of the line of sight, as X(mol) = N(mol)/N(H₂). The results for the ℋ model are shown using the outermost color map, and those for model ℒ are reported using the innermost color map. Overall, the two quantities show a clear linear correlation that remains stable over time. The distributions were fit with a power law of the form $${\log }_{10}{\zeta }_2^{\mathcal{H}}=0.41{\log }_{10}X({\mathrm{e}}^{-})-13.06,$$(12) $${\log }_{10}{\zeta }_2^{\mathcal{L}}=0.32{\log }_{10}X({\mathrm{e}}^{-})-14.53.$$(13)

At first approximation, a linear correlation is also expected for the electron abundance (or column density) and the sum of the main positive ions (DCO⁺, HCO⁺, N₂H⁺, N₂D⁺, and ${\mathrm{H}}_3^{+}$). Considering that ${\mathrm{H}}_3^{+}$ 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{\left({\mathrm{H}}_3^{+}\right)}_{\text{proxy\hspace{0.17em}}}=\frac{1}{3}\frac{N\left({\text{o-H}}_2{\mathrm{D}}^{+}\right)}{R_{\mathrm{D}}},$$(14)

where $$R_{\mathrm{D}}=\frac{N({\mathrm{D}\mathrm{C}\mathrm{O}}^{+})}{N({\mathrm{H}\mathrm{C}\mathrm{O}}^{+})},$$(15)

and we tested the reliability of the following formula that can also be applied to observational data: $$\begin{array}{rl}N({\mathrm{e}}^{-}{)}_{\text{proxy\hspace{0.17em}}}=& N({\mathrm{D}\mathrm{C}\mathrm{O}}^{+})+N({\mathrm{H}\mathrm{C}\mathrm{O}}^{+})+N({\mathrm{N}}_2{\mathrm{D}}^{+})\\ & +N({\mathrm{N}}_2{\mathrm{H}}^{+})+N({\mathrm{H}}_3^{+}{)}_{\text{proxy}}.\end{array}$$(16)

We show in Fig. 7 the distributions of ${\mathrm{H}}_3^{+}$ and e⁻ against their proxies for models ℋ and ℒ averaged over time. The upper panel shows a very tight linear correlation between $N({\mathrm{H}}_3^{+})$ and $N{\left({\mathrm{H}}_3^{+}\right)}_{\text{proxy}}$ up to ~10¹³−10¹⁴ cm⁻² for both models. A deviation then appears, while the relation still remains linear, but becomes steeper. This region corresponds to the central 10 000 au, and the deviation is likely due to the neglected isotopologues in the estimation of ${\mathrm{H}}_3^{+}$ (see Bovino et al. 2019, for a discussion). A similar behavior is also observed for the electron column density (bottom panel of Fig. 7).

Our proxy tends to underestimate the true electron abundance by 65% overall. This was determined as the RMS error of the ratio of the distribution of N(e⁻)_proxy and the actual N(e⁻) because only a limited number of ions was included in the formula and the approximations made. The analytical formula reported in Eq. (16) was validated by employing a simulation of a typical prestellar core, specifically, on sub-parsec scales, but it also holds on larger scales and has a wider and more general applicability within the approximation we introduced so far. First, in the sum we neglected the contribution of ions that might become important (e.g. H₃O⁺ among others) and the charged grains that might be relevant in star-forming regions, which particularly affects the total negative charge. In addition, when we probed ${\mathrm{H}}_3^{+}$ through its first deuterated product (H₂D⁺), we neglected other isotopologues that become important when the core central densities increase (see also Bovino et al. 2020; Redaelli et al. 2024). Finally, we worked with averaged quantities (column densities) to estimate local physical variables. For these reasons, the estimates of the ionisation fraction we provide represent a robust (within a factor of two) lower limit of the real value.

Fig. 5 Radial profiles of the column densities of HCO+(a), DCO+ (b), o-H2D+(c), N2H+(d), and N2D+(e) ions and the ζ2 distribution (f) at 153.3 kyr. The ζ2 radial profile is calculated in circular annuli from the average along the line of sight (as shown in Fig. 2). The CR ionisation models ℋ, ℒ, and ℒ are represented by solid, dashed, and dotted black lines, respectively.

Fig. 6 Correlation between X(e−)and ζ2 obtained from the 2D maps (line-of-sight averages) for models ℋ (outermost color scale) and ℒ (innermost color scale), combining different time snapshots (t = 30.7, 95.0, and 153.3 kyr). The solid and dashed lines correspond to power-law fits. The cross located at log10(ζc) = −16.6 represents the median of the electron abundance (X(e−)) for the ℒ model, and the error bars are the 1σ of the data.

Fig. 7 Pixel-by-pixel 2D histogram distributions for the ℋ (viridis) and ℒ (magma) model by combining different times (t = 30.7, 95.0, and 153.3 kyr). The upper panel correlates the ${\mathrm{H}}_3^{+}$ column density obtained from KROME and ${\mathrm{H}}_3^{+}$ calculated as Eq. (14). The bottom panel shows the correlation between the e− column density and its proxy calculated by Eq. (16). The solid black line represents a 1:1 straight line and serves as an indicator for the lineal correlation formed by the 2D histogram distributions. The grey shadow indicates a region spanning one order of magnitude (−0.5 and +0.5) around the 1:1 black line.

3.4 Ionisation fraction in observed clumps

To show an application of the analytical formula presented in Eq. (16), we employed recent available data for five high-mass clumps at different evolutionary stages: namely G13.18+0.06, G14.11-0.57, G14.63-0.58, G15.72-0.59, and G19.88-0.54. These sources were selected from the APEX³ Telescope Large Area Survey of the Galaxy (ATLASGAL; Schuller et al. 2009), where they have reliable estimates of kinematic distances, masses, luminosities, and dust temperatures (e.g. Urquhart et al. 2022), employed in this derivation.

We used the data presented by Sabatini et al. (2020, 2024) for all tracers in Eq. (16). Sabatini et al. (2020) observed o-H₂D⁺ (1₁₀−1₁₁), N₂H⁺ (4–3), DCO⁺ (1–0), and H¹³CO⁺ (1–0) with the APEX and IRAM-30m telescopes. Additional N₂D⁺(3–2) data were obtained with the nFLASH receiver mounted at APEX through two projects (project-IDs: C-0105.F-9715C-2020 and C-0107.F-9711-2021; PIs: S. Bovino and G. Sabatini, respectively). We refer to Sabatini et al. (2020, 2024) for a full description of the observations and the data analysis.

Considering the recent updates of dust temperatures for ATLASGAL clumps (Urquhart et al. 2022), and following Sabatini et al. (2020), we derived the column densities using MCWeeds (Giannetti et al. 2017), which is an interface between Weeds (Maret et al. 2011) and the Bayesian statistical models of PyMC (Patil et al. 2010), and we assumed T_ex = T_dust equal to 20.3 K (G13.18+0.06), 15.8 K (G14.11-0.57), 19.1 K (G14.63-0.58), 12.1 K (G15.72-0.59), and 20.3 K (G19.88-0.54). The final column densities are in the range ~(0.8–3.2) × 10¹² cm⁻² for o-H₂D⁺, ~(3.7–5.9) × 10¹¹ cm⁻² for N₂D⁺, ~(0.2–3.0) × 10¹³ cm⁻² for N₂H⁺, ~(0.6–4.0) × 10¹³ cm⁻² for DCO⁺, and ~(1.6–6.3) × 10¹³ cm⁻² for H¹³CO⁺. This matches the typical abundances relative to H₂ in high-mass star-forming clumps (e.g. Miettinen 2020; Sabatini et al. 2020; Li et al. 2022). The average error on the computed column densities is ≲20%. We calculated N(HCO⁺) from N(H¹³CO⁺) taking into account a ¹²C/¹³C ratio that changes as a function of the galactocentric distance (D_GC) of each source; that is, [¹²C]/¹³C] = ${6.1}_{-1.8}^{+1.1}$ D_GC + ${14.3}_{-7.2}^{+7.7}$ (Giannetti et al. 2014). When applying Eq. (16) to this range of column densities, we obtained N(e⁻)in the range ~(1.1–3.9) ± 0.5 × 10¹⁵ cm⁻², corresponding to X(e⁻)from 1.8 × 10⁻⁸ to 2.7 × 10⁻⁸. We note that a typical uncertainty of ~30% is associated with X(e⁻), which results from propagating the uncertainty on the integrated line intensity.

Regardless of the scale difference between observations (~1 pc) and simulation (~0.1 pc), HCO⁺ is the dominating ion in both cases, together with ${\mathrm{H}}_3^{+}$ (and its isotopologues) that is necessary to recover the electron abundance in the central region. This confirms the reliability of the method even when it is applied on different scales. The simulations and observations we analysed clearly represent specific physical cases. We expect the formula to work also in situations in which other ions dominate the overall sum, in particular, cases in which small scales are properly resolved by high-resolution observations (e.g. using the Atacama Large Millimeter Array; ALMA) that might unveil emission and a contribution from more compact tracers.

4 Conclusions

We have introduced a CR propagation code for accurately estimating the CR ionisation rate in MHD simulations. We coupled our code to the outputs of MHD simulations performed with GIZMO by implementing an interpolation scheme of gas properties consistent with the smoothing kernel procedure in the simulation. The propagator maps CRs as a flux of TPs for which the effective column density that is traversed following the magnetic field lines is determined by means of a Runge-Kutta integrator. The CR ionisation rate at the location of each TP was then estimated according to the ℋ and ℒ models by Padovani et al. (2018) and was deposited onto the gas tracers in the original simulation in a self-consistent way. The results of the propagator were then used to investigate the evolution of several ions whose abundances are directly affected by CRs.

We tested the convergence of the code and found the optimal number of TPs (see Appendix A) that balance accuracy and performance (10 240). We therefore demonstrated that this number is adequate enough to generate a reliable characterisation of the CR ionisation rate of the object because the differences from the most expensive run (32 times more TPs) are minimal and mostly affect the background region.

After the optimal number of TPs was found, we analysed the time evolution of the chemistry in a collapsing magnetised prestellar core from Bovino et al. (2019) by employing the KROME framework in post-processing as described by Ferrada-Chamorro et al. (2021). We focussed on three different CR models, two variable cases (ℋ and ℒ in Padovani et al. 2018) evolved by means of the propagation scheme, and a reference constant case with ζ₂ = 2.5 × 10⁻¹⁷ s⁻¹, with the aim to assess the CR propagation impact on the chemistry.

We found that a proper CR propagation scheme can produce significant differences (with a non-linear behaviour) for key chemical tracers, in particular, on deuterated species such as DCO⁺ and N₂D⁺. We also analysed the correlations between the electron abundance X(e⁻)and ζ₂ and between the electron column density and the column density of the most relevant positive ions. We found that clear correlations exist in both cases, and the two CR ionisation models ℋ and ℒ only yield a different normalisation. In particular, we provided an analytical function to estimate the electron abundance in star-forming regions with a standard deviation of 0.6, and we showed that when applied to observational data, it provides reliable estimates. For early stages high-mass clumps, we reported values around 10⁻⁸ that agree with previous estimates of prestellar regions (e.g. Caselli et al. 1998; Bergin et al. 1999; Maret & Bergin 2007). We analysed the correlation between ζ₂ and n(H₂) and determined the best-fit expressions, with an uncertainty of ~25% that could be employed as an approximate estimator for ζ₂ in numerical simulations without a self-consistent CR propagation scheme.

Finally, the increasing capabilities of modern astronomical facilities such as ALMA represent a viable method for testing our predictions. Additional high-sensitivity and high-resolution observations of the molecular tracers involved in Eq. (16) are needed to reveal possible fluctuations in electron abundance on core scales. In this context, intensive observational campaigns such as the recently accepted ALMA large program UNveiling the Initial Conditions of high-mass star-formation (UNIC)⁴, which targets ten massive clumps with potentially more than 130 embedded cores, offer a unique opportunity to extend the results reported in this work under a variety of physical and chemical conditions in the near future.

Acknowledgements

GL gratefully acknowledges the financial support of ANID-Subdirección de Capital Humano/Magíster Nacional/2024-22241873 and Millenium Nucleus NCN23_002 (TITANs). SB acknowledges BASAL Centro de Astrofisica y Tecnologias Afines (CATA), project number AFB-17002. GS acknowledges the PRIN-MUR 2020 BEYOND-2p (Astrochemistry beyond the Second period elements, Prot. 2020AFB3FX), the project ASI-Astrobiologia 2023 MIGLIORA (Modeling Chemical Complexity, F83C23000800005), the INAF Mini-grant 2023 TRIESTE (“TRacing the chemIcal hEritage of our originS: from proTostars to planEts”; PI: G. Sabatini), and the National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.1, Call for tender No. 104 published on 2.2.2022 by the Italian Ministry of University and Research (MUR), funded by the European Union – NextGenerationEU-Project Title 2022JC2Y93 Chemical Origins: linking the fossil composition of the Solar System with the chemistry of protoplanetary disks – CUP J53D23001600006 – Grant Assignment Decree No. 962 adopted on 30.06.2023 by the Italian Ministry of MUR.

Appendix A Convergence tests

Because of our choice of kernel size, the scheme does not guarantee that all gas particles in the original simulation are used for the calculation. When this happens, ζ₂ is set to zero, which would yield an inaccurate chemical evolution. In order to avoid this issue, there are two possible solutions. The computationally cheapest one, that however we have not considered in this work, is similar in spirit to the approach in Ferrada-Chamorro et al. (2021), and interpolates ζ₂ using the values of the neighbouring gas particles. The other possibility, that we have employed here, as it guarantees a much higher accuracy, requires a sufficiently large number of TPs, so that no regions of the cloud are left uncovered during the TP propagation.

In order to identify the minimum number of particles required to reach convergence, that is our ‘fiducial’ case, we explored different values of TPs in a suite of test runs, starting from 5120 TPs and gradually increasing their number by a factor of two, as reported in Table A.1.

In Fig. A.1, we show the distribution of ζ₂ per particle for the ℋ case (left-hand panels), and for the ℒ case (right-hand panel), comparing each run in Table A.1 with our fiducial case (X2). Most of the particles are located on the rightmost side of the distribution, with a long tail extending to low ζ₂ values. This reflects that gas particles with large ζ₂ are located near the edge of the cloud where the attenuation is not strong yet, before the steep decrease begins as we move towards the centre of the cloud, which corresponds to the leftmost tail where the attenuation is the strongest. In general, we find a reasonable agreement in the shape of the distribution among the runs, which is consistent with the fact that the distribution is mainly determined by the physical properties of the cloud, which does not vary. The only differences we find can thus be attributed to the number of TPs used, which affects the smoothing procedure and the integration step. On one hand, X1 exhibits an overestimation at the edges of the distribution, because of the poor spatial sampling that is more sensitive to outliers. On the other hand, the increase of TPs slowly carves the high-ζ₂ region, thanks to a better spatial sampling which removes the outliers, but does not significantly change the low-ζ₂ tail, which instead shows stochastic variations due to the small number of gas particles in the central region. This comparison implies that the X2 case represents the best compromise between accuracy and efficiency for the chemical analysis we are interested in. For sake of completeness, in Appendix B we analyse the relative error in the line-of-sight projection of ζ₂ between our fiducial run and the other five runs.

Fig. A.1 ζ2 distribution for ℋ (left) and ℒ (right) at 153.3 kyr in a 0.2 × 0.2 × 0.2 pc3 box. In each row, we compared one of the runs in Table A.1 (in cyan) with our fiducial case (X2, in pink), with the matching regions in purple. fpart represents the fraction of particles.

Appendix B Relative error of expensive runs

Fig. B.1 shows the relative difference between ζ₂ obtained in our fiducial case, X2, and the other runs employing different numbers of TPs, averaged along a line-of-sight aligned with the z-axis. The largest and smallest relative error among the pixels are reported in the bottom left corner of every panel, corresponding to X1, X4, X8, X16, and X32 from top to bottom, respectively. We clearly observe larger deviations for the X1 case, whereas the other cases oscillate around 10–20% almost independently of resolution.

To further check the robustness of our calculations, we determined the error within a spherical region of radius 0.065 pc (identified by the green dashed circle), comparable to the typical size of star-forming filaments. The relative error in this case drops by about a factor of 2, reaching at most 3–13% for both ionisation models. The use of 10240 particles (our base) seems therefore adequate, as higher resolution runs do not significantly improve the results, but at the same time require a much higher computational cost. Obviously, on the largest scales explored here, the errors are larger, but this is expected as we start to encompass more background gas as well, which has lower column densities, hence is more easily affected by variations in ζ₂.

Fig. B.1 Relative error (ϵ) for ℋ (left) and ℒ (right) within a 0.5 × 0.5 × 0.2 pc3 at 153.3 kyr. We chose this box size to take into account the 0.5 pc of injection and a depth length similar to the filament scale (0.2 pc) to avoid the inclusion of background particles. The errors are calculated between X2 case and the rest. The panels are, from top to bottom: X1, X4, X8, X16 and X32. The bluer and redder indicate underestimation and overestimation of X2 relative to the others. The maximum and minimum errors for each comparison are written in the bottom-left corner. All green spheres have 0.065 pc radii.

References

All Tables

All Figures

Fig. 1 Sketch of the CR propagation scheme based on the structure given in Section 2.1. The left panel describes the TP propagation through the domain. TPs (grey dots) are created from a spherical surface of radius R (dashed black circle) with radial inward trajectories (step i). TPs are then deactivated when the column density that is traversed exceeds 1029 cm−2 or when they leave the outer spherical surface at a distance Rout (solid black circle), as shown by the crossed grey dots. The right panel instead highlights the propagation of a single TP particle. The trajectory is shown by the curved solid blue arrow, and the magnetic field lines (B) are reported as dash-dotted grey lines. The pitch angle (α) between the TP velocity and the local magnetic field is shown as the dashed red line. The displacement of the TP between the previous (x) and the current (x′) location is defined by the step size Δℓ (step iii). To compute the effective column density of the TP, the average local density (n(x′)) and local magnetic field (B(x′)) were estimated from gas particles (pink dots) within the kernel size h (thin yellow arrows; step ii). Then, its value was used to update the CR ionisation rate (ζ2; step iv), which was finally deposited at the location of the original gas particle (thick yellow arrows; step v).

In the text

	Fig. 2 Time evolution of the density-weighted ζ2 for models ℋ (top panels) and ℒ (bottom panels) in a sphere of 0.5 pc radius, inside which CRs are propagated. The green arrows correspond to the x-y projection of the magnetic field lines. The two different models yield a difference of one order of magnitude within the high-density region.
In the text

Fig. 3 Time evolution of the distribution of ζ2 for the gas particles in the simulation for models ℋ (top) and ℒ (bottom). The distributions are determined from a box of 0.2 × 0.2 × 0.2 pc3, which only includes particles from the simulated cloud, i.e., with identical mass. The rightmost edge of the distribution shrinks over time as the gas particles within the box become denser. The low-end tail of the distribution expands towards lower ionisation rates, which is consistent with the evolution of the gas density.

In the text

	Fig. 4 Correlation between ζ2 and n(H2) at 153.3 kyr. The viridis and magma color maps represent the ℋ and ℒ models, respectively. The solid and dashed red lines correspond to our best-fit relations, and the shaded grey areas show the intrinsic scatter discussed in Section 3.2.
In the text

	Fig. 5 Radial profiles of the column densities of HCO+(a), DCO+ (b), o-H2D+(c), N2H+(d), and N2D+(e) ions and the ζ2 distribution (f) at 153.3 kyr. The ζ2 radial profile is calculated in circular annuli from the average along the line of sight (as shown in Fig. 2). The CR ionisation models ℋ, ℒ, and ℒ are represented by solid, dashed, and dotted black lines, respectively.
In the text

Fig. 6 Correlation between X(e−)and ζ2 obtained from the 2D maps (line-of-sight averages) for models ℋ (outermost color scale) and ℒ (innermost color scale), combining different time snapshots (t = 30.7, 95.0, and 153.3 kyr). The solid and dashed lines correspond to power-law fits. The cross located at log10(ζc) = −16.6 represents the median of the electron abundance (X(e−)) for the ℒ model, and the error bars are the 1σ of the data.

In the text

Fig. 7 Pixel-by-pixel 2D histogram distributions for the ℋ (viridis) and ℒ (magma) model by combining different times (t = 30.7, 95.0, and 153.3 kyr). The upper panel correlates the ${\mathrm{H}}_3^{+}$ column density obtained from KROME and ${\mathrm{H}}_3^{+}$ calculated as Eq. (14). The bottom panel shows the correlation between the e− column density and its proxy calculated by Eq. (16). The solid black line represents a 1:1 straight line and serves as an indicator for the lineal correlation formed by the 2D histogram distributions. The grey shadow indicates a region spanning one order of magnitude (−0.5 and +0.5) around the 1:1 black line.

In the text

	Fig. A.1 ζ2 distribution for ℋ (left) and ℒ (right) at 153.3 kyr in a 0.2 × 0.2 × 0.2 pc3 box. In each row, we compared one of the runs in Table A.1 (in cyan) with our fiducial case (X2, in pink), with the matching regions in purple. fpart represents the fraction of particles.
In the text

Fig. B.1 Relative error (ϵ) for ℋ (left) and ℒ (right) within a 0.5 × 0.5 × 0.2 pc3 at 153.3 kyr. We chose this box size to take into account the 0.5 pc of injection and a depth length similar to the filament scale (0.2 pc) to avoid the inclusion of background particles. The errors are calculated between X2 case and the rest. The panels are, from top to bottom: X1, X4, X8, X16 and X32. The bluer and redder indicate underestimation and overestimation of X2 relative to the others. The maximum and minimum errors for each comparison are written in the bottom-left corner. All green spheres have 0.065 pc radii.

In the text