© The Authors 2024

1 Introduction

Cosmic rays (CRs) are energetic, ionised particles that are ubiquitously found in the interstellar medium (ISM). In the densest gas phases, they play a key role in triggering the chemistry and regulating the thermodynamical balance. In molecular gas, CRs collide with H₂ molecules and usually end up ionising them via the following reaction: $$\text{CR}+{\text{H}}_2\to {\text{H}}_2^{+}+{\text{e}}^{-}+\text{CR}.$$

Afterwards, the highly reactive ${\text{H}}_2^{+}$ immediately reacts with another H₂ molecule, producing the trihydrogen cation ${\text{H}}_3^{+}$. In doing so, CRs set the ionisation state of dense matter, where the photodissociating UV field is completely attenuated. This, in turn, strongly affects the dynamics of the dense gas. The electron fraction x(e⁻) determines the degree of coupling between the magnetic fields and the matter and it is therefore linked to the time scale for ambipolar diffusion. This corresponds to the drift between the neutral and the ionised flows, which is one of the proposed mechanisms behind the dissipation of the magnetic flux and enabling gravitational collapse (Mouschovias & Spitzer 1976).

With respect to chemical evolution in the ISM, ${\text{H}}_3^{+}$ is a pivotal species also since it drives the rich ion chemistry that is ultimately responsible for such phenomena as the production of CO and gas-phase water in cold regions. Furthermore, by reacting with deuterated hydrogen HD, ${\text{H}}_3^{+}$ is converted into H₂D⁺, which is the starting point of the deuteration process in the gas phase (see for instance Ceccarelli et al. 2014, and references therein).

In summary, CRs are crucial to the physics and chemistry of star-forming regions. The key question is how to measure their effects, particularly their ionisation rate per hydrogen molecule (ζ₂). The question of how to infer ζ₂ from typical observables has been the focus of several works dating back to the 1970s (see for instance Guelin et al. 1977; Wootten et al. 1979). In the diffuse medium (A_V ≲ 1 mag), ${\text{H}}_3^{+}$ can be directly observed in absorption towards bright infrared sources and its relatively simple kinetics can be solved to infer ζ₂. This is the approach followed for instance by McCall et al. (2003); Indriolo et al. (2007); Indriolo & McCall (2012), with typical results of the order of ζ₂ ≈ 10⁻¹⁶ s⁻¹.

At higher densities, the absorption lines of ${\text{H}}_3^{+}$ can no longer be detected. To address this situation, Guelin et al. (1977, 1982) proposed an analytic expression containing ζ₂, based on the kinetics of DCO⁺ and HCO⁺. This work was later expanded by Caselli et al. (1998, hereafter CWT98), who developed a system of two equations to first infer the ionisation fraction x(e⁻) and then ζ₂. Those authors ultimately used a more comprehensive chemical code to infer these parameters in a sample of dense cloud cores (see their Sect. 5.1 and Table 7), resulting in ζ₂ values on average lower than with the analytic expression. However, the analytic equations are still used nowadays (for instance Cabedo et al. 2023), due also to the fact that they depend on common observable quantities, such as¹ R_H = X(HCO⁺)/X(CO) and the deuteration level of HCO⁺, R_D = X(DCO⁺)/X(HCO⁺). It is important to note, however, that this approach was based on several assumptions, such as the abundance of HD. Furthermore, it was derived in specific conditions (for instance at a temperature of T_K = 10 K), which prevent its generalisation.

More recently, Bovino et al. (2020, hereafter BFL20) suggested a new analytical approach to determine ζ₂ in cold gas, based on the detection of H₂D⁺ (and, in particular, of its ortho spin state, oH₂D⁺, which can be observed from the ground). This method is, in turn, based on the formulation of Oka (2019) and the underlying idea is to infer the abundance of ${\text{H}}_3^{+}$ from that of its deuterated forms, together with the deuterium fraction measured from HCO⁺ isotopologues. It relies on fewer assumptions than CWT98, but it requires oH₂D⁺ data, which are observa-tionally more expensive in terms of observing time. Sabatini et al. (2020) applied this implementation to a sample of highmass star-forming clumps, using observations made with the Atacama Pathfinder EXperiment 12-m submillimeter telescope (APEX; Gusten et al. 2006) and the Institut de Radioastronomie Millimétrique (IRAM) 30m single-dish telescope. Their analysis yielded 7 × 10⁻¹⁸ < ζ₂/s⁻¹ < 6 × 10⁻¹⁷, in line with theoretical predictions (Padovani et al. 2018). Sabatini et al. (2023) used the BFL20 equation to estimate ζ₂ maps at high resolution in two high-mass clumps, using data from the Atacama Large Millimeter and sub-millimeter Array (ALMA), finding a remarkable agreement with the most recent predictions of cosmic-ray propagation and attenuation (Padovani et al. 2022).

Another possible approach is based on having an underlying chemical model to interpret the observational results, possibly using multiple tracers. This can be done either by comparing the abundances and column densities obtained from observations to the results of chemical codes (as done, for instance, by Caselli et al. 1998; Ceccarelli et al. 2014; Fontani et al. 2017; Favre et al. 2018) or employing radiative transfer analysis. In the latter case, used for instance in Redaelli et al. (2021b), one compares the synthetic spectra of several transitions based on the molecular abundances from the chemical modelling to the observed lines. This approach is more sophisticated, and can potentially constrain ζ₂ more accurately, but it also depends on the employed model and the detail of the physics included to study the given source. However, its applicability is limited by the considerable combined runtimes of a chemical code and radiative transfer.

In cases where evaluating ζ₂ is required in large surveys, an analytical expression presents a clear advantage in terms of applicability and use of computational resources. The equation from BFL20 was tested on several setups of physicochemical three-dimensional (3D) simulations varying, for instance, the collapse speed or the initial H₂ ortho-to-para ratio (OPR). However, a validation based on simulated observations was never performed. It is crucial to test whether (and how) the results are affected by observational effects depending, for instance, on the selected transitions or the line opacities. The same is true for the CWT98, which moreover was never tested, either on a simulation set or on synthetic observations.

In this paper, we aim to: (i) investigate the effects of telescope response and of several possible observational biases and (ii) compare two widely used equations to infer ζ₂ from observable tracers and look for the potential applicability limits of the formulas. For these reasons, we have selected two simulation setups of a dense core, where the ζ₂ is known by construction, and it varies in the typical range predicted by models at high densities (≈10⁻¹⁸−10⁻¹⁶ s⁻¹). These simulations are used as input for a radiative transfer code, producing synthetic spectra of the chemical tracers involved. These are then post-processed to simulate either single-dish or interferometric-like observations. From this point, the analysis follows the standard procedure to apply the analytical expressions: the datacubes are converted into column densities and abundances maps, on which the formulas can be evaluated. The resulting ζ₂ maps are compared to the actual input values of the ionisation rate.

The paper is organised as follows. Section 2 introduces the analytical methods put to the test. Section 3 illustrates the simulation setup, radiative transfer analysis, and the post-processing procedure to infer column density maps. The ζ₂ maps are computed and described in Sect. 4. First, we focus on the BFL20 formulation, performing several runs, starting from a reference one (where ζ₂ = 2.5 × 10⁻¹⁷ s⁻¹ and the evolutionary time is t = 100 kyr), and then exploring distinct post-processing types, radiative transfer methods, and a set of rotational lines at different evolutionary stages (Sect. 4.1). In Sect. 4.2, we test the CWT98 analytical approach as an alternative approach. Section 4.3 presents the comparison of the two methods on literature values of real observed objects. The results are discussed in Sect. 5. A final summary with concluding remarks is presented in Sect. 6.

2 Retriving ζ₂ using analytical expressions

The first analytical expression we aim to test is the one proposed by BFL20. In Appendix A, we present the detailed derivation of the final equation based on observable quantities (namely molecular column densities N_col), as: $${\zeta }_2=k_{\text{CO}}^{{\text{oH}}_3^{+}}\frac{1}{3}\times X(\text{CO})\times \frac{N_{\text{col}}\left({\text{oH}}_2{\text{D}}^{+}\right)}{R_{\text{D}}}\frac{1}{L},$$(1)

where the deuterium fraction of HCO⁺ is $$R_{\text{D}}=\frac{N_{\text{col}}\left({\text{DCO}}^{+}\right)}{68\times N_{\text{col}}\left({\text{H}}^{13}{\text{CO}}^{+}\right)},$$

and the CO abundance (with respect to H₂) is $$X(\text{CO})=\frac{557\times N_{\text{col}}\left({\text{C}}^{18}\text{O}\right)}{N\left({\text{H}}_2\right)}.$$

In Eq. (1), $k_{\text{CO}}^{{\text{oH}}_3^{+}}$ is the destruction rate of ${\text{oH}}_3^{+}$ by CO (see reaction rate k₇ in Appendix A), assumed here to be the main destruction path for ${\text{H}}_3^{+}$, which holds when the CO depletion factor is up to ≈100 assuming x(e⁻) ≈ 10⁻⁸ (see also the discussion in Appendix A for more details). Then, L is the path length over which the column densities are estimated. Note that for tracers that are typically optically thick (CO and HCO⁺ ), we use the corresponding optically thin isotopologues, assuming standard isotopic ratios for the local ISM: ¹⁸O/¹⁶O = 557 and ¹³C/¹²C = 68 (Wilson 1999). This choice is consistent with that made later in the radiative transfer analysis (see Sect. 3.3). We stress that Eq. (1) can be used only until when H₂D+ is the dominant deuterated species of ${\text{H}}_3^{+}$, namely, at the early stages of star-forming regions.

We aim to test the performance of the approach proposed by CWT98 as well, exploring its applicability and comparing it with the BFL20 approach. In these regards, we first used their equations (Eqs. (3) and (4) as they are described in Sect. 4.2). This is motivated by the fact that some papers have been using them in that exact formulation (although CWT98 already employed chemical modelling to interpret the observational results). The numerical constants in those equations, however, were derived at a temperature of 10 K, and without including ortho- and para-states of the involved species. To update the work of CWT98, in Appendix B, we followed the same approach but including the spin-state separation of H₂, ${\text{H}}_3^{+}$, and H₂D⁺. We also updated the reaction rates involved with the most recent values, and we keep their temperature dependency.

3 Simulations and post-processing

In this work, we performed a post-processing of 3D simulations of prestellar cores to produce the observables needed to test the aforementioned analytical methods. In the following subsections, we describe the set of simulations used in our test, the details of the radiative transfer and the type of post-processing performed afterwards. Finally, we discuss how we recover the molecular column densities and the gas total column density from the synthetic observations.

3.1 The simulation setup

We used a set of three-dimensional magneto-hydrodynamical (3D MHD) simulations of prestellar cores, with constant and variable cosmic-ray ionisation rates, obtained with the code GIZMO (Hopkins 2015). The simulations evolved an isothermal, turbulent, and magnetised Bonnor-Ebert sphere of 20 M_⊙ with a radius of 0.17 pc, resembling a collapsing prestellar core. After 100 kyr of evolution, we obtain a low-mass object, with a total mass of ~3 M_⊙ in the central 10 000 AU (corresponding to an average density of 〈ra(H₂)〉 ~ 10⁵ cm⁻³). The gas and dust temperatures are T_k = T_dust = 15 K. The size of the total simulation box is 0.6 pc. For the purposes of this work, we focussed on the central 0.3 pc containing the Bonnor-Ebert sphere.

The simulations include a state-of-the-art deuterated and spin-state chemical network, advanced in time alongside hydrodynamics. We refer to Bovino et al. (2019) for the complete description of the physical and chemical initial conditions. We highlight that the chemical code includes molecular depletion and thermal and cosmic-ray-induced desorption, but no further surface chemistry is considered. This assumption is justified since, at the low temperatures considered, the thermal desorption of any species is negligible, and the cosmic-ray-induced desorption timescale is longer than the collapse one (cf. Bovino et al. 2019). Hence, surface chemistry would have negligible effects on the final abundances of gas phase species. The initial H₂ ortho-to-para ratio is OPR = 0.1, consistent with the values obtained by large-scale simulations (Lupi et al. 2021).

To investigate different ionisation states, we first used a simulation with constant ${\zeta }_2^{\text{s}}=2.5\times {10}^{-17}{\text{s}}^{-1}$, a value typically assumed for dense regions (this corresponds to the M1 case of Bovino et al. 2019). Throughout this paper, we use the superscript ‘s’ to denote the input ${\zeta }_2^{\text{s}}$ of the simulations, to avoid confusion with the ζ₂ retrieved from the synthetic observations. We performed an additional test on a simulation performed with constant ${\zeta }_2^{\text{s}}=2.5\times {10}^{-18}{\text{s}}^{-1}$, to expand the range of input ζ₂ explored and to be consistent with the tests performed by Bovino et al. (2020). We also simulated a variable ${\zeta }_2^{\text{s}}$ model, by postprocessing the same simulation by employing the framework developed by Ferrada-Chamorro et al. (2021) where ζ₂ is varied according to the density-dependent function reported in Ivlev et al. (2019). The details of the latter simulations will be presented in a forthcoming paper (Gaete-Espinoza et al., in prep.). Thus, this set of simulations covers the typical ζ₂ values predicted by the most recent models of CR propagation at high densities (for N(H₂) ≳ 10²² cm⁻² ζ₂ ≲ a few × 10⁻¹⁶ s⁻¹; cf. Padovani et al. 2022).

3.2 Description of tested runs

Since our goal is to compare the results of analytic expressions to infer ζ₂ from observables, we tested a variety of combinations of types of radiative transfer (assuming or deviating from the local-thermodynamic-equilibrium), of rotational transitions (ground state lines, or higher J transitions), and of post-processing (simulating single-dish or interferometric observations). Starting from a reference run, we modified one parameter at a time, to evaluate its effects. We performed a total of eight tests, reported in Table 1 and described in the following text. The details of how the radiative transfer is performed are discussed in Sect. 3.3, whilst Sect. 3.4 describes how the telescope response is simulated.

Runs 1 and 2 (the latter is considered the reference one throughout the paper) used the simulation with constant ${\zeta }_2^{\text{s}}=2.5\times {10}^{-17}{\text{s}}^{-1}$ at an evolutionary time of 50 and 100 kyr, respectively. The radiative transfer was performed in local-thermodynamic-equilibrium approximation (LTE), focussing on the molecular lines in the 215–370 GHz range. The postprocessing was single-dish-like, with a final beam size of 27″. Run 3 tested the low-${\zeta }_2^{\text{s}}$ case, where the only difference with respect to the reference run is the value ${\zeta }_2^{\text{s}}=2.5\times {10}^{-18}{\text{s}}^{-1}$. Run 4 explored the interferometric-like post-processing, (see Sect. 3.4.2). In runs 5 and 6, we modified the kind of radiative transfer, using a non-LTE approach. In particular, in run 5 we simulated again the high-J transition of C¹⁸O, DCO⁺, and H¹³CO⁺. In run 6, instead, we aimed to explore the effect of targeting the (1-0) transitions of C¹⁸O, H¹³CO⁺, and DCO⁺. This is beneficial to the observational studies that trace the 3 mm lines of these molecules, as in the pioneering work of Caselli et al. (1998). We also adopted a large beam size of ≈70″, to simulate poorly resolved observations, where the beam area is comparable to the source size. The final two runs adopted the LTE analysis and single-dish-like post-processing, performed on the simulation with variable ${\zeta }_2^{\text{s}}$ at 50 kyr (run 7) and 100 kyr (run 8), respectively.

3.3 Radiative transfer of MHD simulations

The total gas column density distribution and the molecular column densities are involved in the equations to infer ζ₂ (see Sect. 2). The radiative transfer of the dust and the molecular lines was performed with the POLARIS code² (Reissl et al. 2016; Brauer et al. 2017). Concerning the radiative transfer of the dust, we set the gas mean molecular weight to µ = 2.4 (Kauffmann et al. 2008) and the gas-to-dust mass ratio to 100 (Hildebrand 1983). We simulated the dust emission at a wavelength of 870 µm (corresponding to 345 GHz). From an observational perspective, this was the wavelength of the LABOCA instrument mounted on APEX, which was used to perform the all-sky survey ATLASGAL (Schuller et al. 2009). It is also close to the James Clerk Maxwell Telescope (JCMT) SCUBA II longer wavelength (850 µm). Finally, it represents the standard frequency for ALMA continuum observations in Band 7. We stress, however, how the choice of wavelength does not impact the results. For the dust model, we assumed pure silicate grains, with opacities taken from Laor & Draine (1993)³. The grain size distribution is a standard MRN (Mathis et al. 1977), with a power-law index of −3.5 between 5 nm and 0.25 µm. The grain density is 3.5 g cm⁻³, which is consistent with the value assumed in the simulations. These parameters are likely different from the real dust populations within dense cores, where for instance a mixture of carbonaceous and silicate grains is expected. However, our goal is not to reproduce the exact properties of a real dust population, but to be consistent in the various steps of the analysis, from the simulation to the radiative transfer.

For the molecular tracers, POLARIS needs as input the spec-troscopic description of the simulated transitions, which are summarised in Table 2. In the case of oH₂D⁺, the only line accessible from the ground is the (1_1,0 − 1_1,1) one at 372 GHz (Caselli et al. 2003). This can be observed, for instance, by ALMA in Band 7 and by APEX using, for instance, SEPIA345 or LAsMA. Concerning the other tracers (C¹⁸O, DCO⁺, and H¹³CO⁺), their transitions at 215–260 GHz are commonly observed. However, several studies have focussed on their lower-J transitions at 3 mm, hence, we also tested these lines in run 6. All molecular transitions were simulated over a total velocity range of 7.5 km s⁻¹ and a velocity resolution of 0.1 kms⁻¹. We assumed that the local standard-of-rest velocity of the source is 0 km ⁻¹.

POLARIS is able to perform different approaches of radiative transfer, including LTE and large velocity gradient (LVG). We performed six runs assuming LTE conditions, which allowed us to focus initially on the effect of the radiative transfer itself and of the response of the telescope on the inferred ζ₂ values. However, lines with high critical densities (n_c ≈ 10⁵−10⁶ cm⁻³), such as the high-J transitions of H¹³CO⁺ and DCO⁺, and oH₂D⁺, are likely to be sub-thermally excited. Two runs (nos. 5 and 6) hence explored a more realistic case, using the LVG approach.

POLARIS accepts a variety of grid types as input, in particular Voronoi grids. The simulations we consider were obtained with GIZMO, which samples the fluid using a set of discrete tracers representing a sort of cells with smoothed boundaries. In this respect, converting this volume discretisation to a Voronoi tessellation is the most natural and consistent choice, despite some differences existing between the two volume partition schemes⁴. We, hence, prepared the outputs of the simulations in the form of a Voronoi grid. In order to properly treat boundary cells, we added at the edges of the region of interest a set of virtual particles placed according to a cubic regular grid. We placed virtual particles at 1.5 times the simulated region size to avoid artefacts and we passed them to the SCIPY package (Virtanen et al. 2020) to build the Voronoi cells. We then cut the grid to match the original region, and we passed the information associated with every cell to POLARIS, including the IDs of the cell and its neighbours. In particular, the gas density, gas and dust temperatures, velocity field (3D), and the molecular mass fraction for each species (f_mol = ρ_mol/ρ_tot) were the input of the radiative transfer. We note that the chemical code does not treat oxygen or carbon fraction-ation. The mass fractions of C¹⁸O and H¹³CO⁺ were thereby derived from the main ones of the isotopologues (using the same standard isotopic ratios assumed in Sect. 2).

In all radiative transfer analyses, the grid size of the output maps or cubes produced by POLARIS was set to 256 pix × 256 pix. We aimed to produce synthetic observations both for a single dish-like and for an interferometer-like case, with the distance of the source set to 170 pc and 2kpc, respectively⁵. The final pixel and field-of-views (FoV) are 1″.4 and 6′ × 6′ (single-dish), and 0″.12 and 30″ × 30″ (interferometer-like).

3.4 Post-processing of the POLARIS output

The output of POLARIS consists of bi-dimensional maps (in Jy pix⁻¹), one for each wavelength for the continuum emission (a single one at 870 µm in our case) or one for each velocity channel set for the molecular lines. In the latter case, the first stage is to build the position-position-velocity datacube concatenating all the velocity channels. We now describe the approaches used to simulate a single-dish-like or interferometric response.

3.4.1 Single-dish analysis

In this case, we convolved the continuum maps and the molecular line datacubes to a specific beam size. For all the tests performed with the higher J transitions, we chose a beam size of 27″. This corresponds approximately to the APEX beam size at the lowest frequency in the sample, DCO⁺ (3−2) at 216 GHz. In the case of run 6, where we simulated the lower J = (1−0) lines (see Sect. 4.1.4 for more details), we selected a beam size of 70″. It is aimed at determining the effects of poorly resolved observations.

We introduced pixel by pixel in the data cubes and in the continuum fluxes some artificial Gaussian noise with zero mean and rms standard deviation. For the continuum maps, we used rms = 15mJybeam⁻¹ and rms = 100mJybeam⁻¹ for the cases at 27″ and 70″ of resolution, respectively. Concerning the line datacubes, we injected a noise with rms = 100 mK. For oH₂D⁺ in run 1, and all lines in run 5, this sensitivity is insufficient for significant detections (see Sect. 4.1.1 and 4.1.4 for more details). In these runs, we reduced the noise level to rms = 50 mK. These values are consistent with the typical rms of observational campaigns with APEX (cf. Sabatini et al. 2020). Run 3, performed with the lowest ${\zeta }_2^{\text{s}}$ value, present faint lines, and the noise level is reduced to rms = 1 mK (cf. Sect. 4.1.2).

3.4.2 Interferometer-like analysis

To simulate interferometer-like observations, we focussed on ALMA, which can cover the frequencies of the transitions analysed here, except for DCO⁺ (1−0). We hence used the output of POLARIS as input for the task SIMOBSERVE of CASA (version 6.4.3). Due to current limitations of SIMOBSERVE, it is not possible to add total power at the desired sensitivity. We hence simulated only the 12m and 7m-array observations. We chose the Cycle 8 configuration sets. We selected the integration times using the corresponding observing Tool (OT), setting the requested angular resolution to 1″ and the desired noise level to 100 mK and to 5 µJy for the continuum simulations. Table 3 summarises the integration times used in each run of SIMOBSERVE. Concerning the noise, we let the task construct the atmospheric model (using the option THERMAL_NOISE = TSYS-ATM).

The task SIMOBSERVE was called separately to simulate the 12m and the 7 m array observations. We then combined the output visibilities (using CONCAT), making sure that the relative weights are correct⁶. After that, we imaged the concatenated visibilities using the TCLEAN task. We used the MULTISCALE deconvolver (scales: [0, 5, 15]× pixel size), which is an appropriate choice in case of extended emission, such as in the simulated data. We selected a BRIGGS weighting, with ROBUST = 0.5. The noise threshold was set to 2σ. The final datacubes (or 2D images for continuum observations) have a FoV of 30″ × 30″, sampled with 250 pix × 250 pix.

The lack of total power observations leads to flux losses, due to the filtering of the large-scale emission which is particularly important for extended sources such as the core we simulate. Focussing on run 4, we estimate that between ≈20 and 65% of the flux in a 15″ area around the core is recovered, depending on the tracer. This is in line with simulations regarding the filter-out effect. For instance, Plunkett et al. (2023) found that up to 90% of the original flux can be lost in extended sources when single-dish data are not available.

3.5 Column density computation

The different approaches for computing the cosmic-ray ionisation rate depend on the column densities of the involved species.

To estimate them, we used the approach of Mangum & Shirley (2015): $$N_{\text{col }}=\frac{8\pi Q_{\text{rot}}\left(T_{\text{ex}}\right)v^3}{c^3{A}_{\text{ul}}{ɡ}_{\text{u}}}\times \frac{\exp \left(E_{\text{u}}/k_{\text{B}}{T}_{\text{ex}}\right)}{\exp \left(hv/k_{\text{B}}{T}_{\text{ex}}\right)-1}\times {{\displaystyle \int }}^{\text{}}{\tau }_v\text{d}V,$$(2)

where h, k_B, and c are the Planck constant, the Boltzmann constant, and the speed of light in vacuum; E_u is the upper-level energy, ɡ_u the upper-level multiplicity, v the line frequency, A_ul the Einstein coefficient for spontaneous emission, and Q_rot(T_ex) the partition function at the excitation temperature T_ex. We list the values used for the spectroscopic constant in Table 2. They are taken from the CDMS catalog⁷. The partition functions are from Giannetti et al. (2019) for oH₂D⁺, the CDMS catalogue for C¹⁸O and H¹³CO⁺, and Redaelli et al. (2019) for DCO⁺. To compute the partition function at the requested temperature, we interpolated the values linearly, when necessary. ∫ τ_vdV is the integral along the velocity axis of the optical depth τ_v computed channel-by-channel using (cf. Caselli et al. 2002) $${\tau }_v=-\ln \left[1-\frac{T_{\text{MB}}}{J_v\left(T_{\text{ex}}\right)-J_v\left(T_{\text{bg}}\right)}\right],$$(3)

where T_MB is the line main beam temperature, J_v(T) is the equivalent Rayleigh-Jeans temperature at the line frequency⁸, and T_bg = 2.73 K is the background temperature. The obtained optical depth profiles are integrated along the velocity range [−1 : 1] km s⁻¹, which is large enough to include the line emission in all the synthetic cubes analysed for this work.

The excitation temperature for all the transitions is T_ex = T_k = 15 K when assuming LTE conditions. For the two non-LTE cases, we selected the excitation temperatures based on the critical densities of the analysed transitions and on available literature data⁹. The C¹⁸O first two rotational lines have relatively low critical densities (n_c ≈ 10³ cm⁻³), and it is hence reasonable to assume that they are thermalised by collisions with H₂, leading to T_ex = T_k = 15 K. The critical density of oH₂D⁺ is higher (n_c ≈ 10⁵ cm⁻³, Hugo et al. 2009), and the line is likely sub-thermally excited, leading to T_ex < T_k. We adopt T_ex = 10 K, which is frequently employed in the literature. For instance, Caselli et al. (2008) computed T_ex = 7–13 K in the envelope of protostellar cores that have gas temperatures of 10–15 K, close to that of our simulations; Friesen et al. (2014) adopted T_ex = 12 K; Redaelli etal. (2021a, 2022) used T_ex = 10 K. DCO⁺ and H¹³CO⁺ are isotopologues with similar critical densities, and it is reasonable to assume that the same rotational transitions share similar excitation temperatures. However, literature information about these are scarce. Using a full non-LTE modelling of the DCO⁺ lines in the well-known core L1544, Redaelli et al. (2019) found T_ex(1−0) = 5.7 K and T_ex(3−2) = 7.8 K. The core L1544 is, however, colder than our simulated cores. We, hence, used the online tool RADEX¹⁰ to confirm these values. Using n = 10⁵cm⁻³, T_k = 15 K, and N_col = 10¹²cm ⁻², we derived T_ex ≈ 8–12 K for the lower-J transitions and ≈5 K for the higher-J ones. We therefore set T_ex(1−0) = 10 K and T_ex(3−2) = 5.5 K for both isotopologues. In Appendix C we show that a 20% variation of these values does not affect our conclusions. The last column of Table 2 summarises the excitation temperature values used in the LVG analysis.

To estimate the uncertainties on the derived column density values (rmsN), we applied standard error propagation on Eq. (2), assuming that the frequency channels are independent and using the small-error approximation. Then, the uncertainty propagation leads to $$\begin{array}{c}{\text{rms}}_{\text{N}}=\frac{8\pi Q_{\text{rot}}\left(T_{\text{ex}}\right)v^3}{c^3{A}_{\text{ul}}{ɡ}_{\text{u}}}\times \frac{\exp \left(E_{\text{u}}/k_{\text{B}}{T}_{\text{ex}}\right)}{\exp \left(hv/k_{\text{B}}{T}_{\text{ex}}\right)-1}\times rms\\ \times \text{Δ}{V}_{\text{ch}}\sqrt{{\displaystyle \underset{k=c{h}_i}{\overset{c{h}_f}{{{\displaystyle \sum }}^{\text{}}}}}{\left(\frac{1}{J_v\left(T_{\text{ex}}\right)-J_v\left(T_{\text{bg}}\right)-T_{\text{MB}}^k}\right)}^2},\end{array}$$(4)

where ch_i and ch_f are channels corresponding to the velocity interval over which Eq. (3) is computed, $T_{\text{MB}}^k$ is the intensity of the k-th channel, and ΔV_ch is the channel width (in km s⁻¹).

To estimate the abundances, we derived the total gas column density map from the continuum map as $$N\left({\text{H}}_2\right)=f\times \frac{S_{\text{pix}}{D}^2}{B_v\left(T_{\text{dust}}\right){\kappa }_v{\mu }_{{\text{H}}_2}{m}_{\text{H}}{\text{Ω}}_{\text{pix}}},$$(5)

where f = 100 is the gas-to-dust mass ratio (Hildebrand 1983), D is the source distance, B_v(T_dust) is the Planck function at the dust temperature T_dust = 15 K, ${\mu }_{{\text{H}}_2}=2.8$ is the mean molecular weight per hydrogen molecule, m_H is the mass of the hydrogen atom, S_pix is the flux (in units of Jypix⁻¹), Ω_pix is the pixel size (in physical units), and κ_V is the dust opacity. For the latter, we used the output of POLARIS, which tabulates the opacities at the simulated wavelength: κ_{345 GHz} = 0.388 cm² g⁻¹. We estimated the uncertainties on the total column densities using Eq. (5), with the flux noise level of the continuum map as S_pix.

4 Resulting ζ₂ maps

We go on to apply the analytical expressions described in Sect. 2 to infer the ζ₂ maps. Uncertainties on derived quantities were computed pixel-per-pixel assuming standard error propagation calculated from the uncertainties on the column densities that are considered independent. We neglected, for instance, any source of uncertainty from the reaction rates. With the computed errors, we masked pixels where the signal-to-noise ratio is S/N ≤ 3.

4.1 The Bovino et al. (2020) method

To test the analytical method of Bovino et al. (2020, BFL20), we computed ζ₂ in the eight runs described in Table 1, varying the ζ₂ model (constant or variable), the type of radiative transfer (LTE or LVG), the post-processing method (single-dish or ALMA-like), and the frequency of the molecular lines. All the simulations are isothermal at 15 K, and, therefore, the rate coefficient in Eq. (1) is $k_{\text{CO}}^{{\text{oH}}_3^{+}}=2.3\times {10}^{-9}$ cm³ s⁻¹.

In the following subsections, we discuss in detail the results of each run, showing the obtained ζ₂ maps. In order to compare these values with the actual ones, we make use of the ζ₂ distributions in Fig. 1. Its panels show the histograms of the ionisation rate values extracted in the densest region of the core (namely where the H₂ column density is higher than 50% of its peak value), where the signal-to-noise ratio is higher. The medians of the distributions (vertical dashed lines) are directly compared with the actual value of 2.5 × 10⁻¹⁷ s⁻¹ in the runs with constant ${\zeta }_2^{\text{S}}$. In runs with variable ${\zeta }_2^{\text{S}}$, we compare pixel-by-pixel the ratio between actual and retrieved values, as discussed in Sect. 5.

Fig. 1 Histogram of the distribution of ζ2 (in units of 10−17 s−1) of the runs as labelled at the top of each panel. We represent the core region where N(H2) > 0.5 × max(N(H2)) (namely the lowest contour in the panels of Figs. 2, 3, and 4). The median value, together with the median uncertainty, is reported in the top-right corner of each panel and shown as a vertical dashed line with a grey-shaded area. Cf. Table 1.

4.1.1 Runs 1–2: constant ${\zeta }_2^{\text{S}}$ and single-dish like analysis

The first runs we tested have constant ${\zeta }_2^{\text{S}}=2.5\times {10}^{-17}$ s^-1, coupled with LTE radiative transfer of the high-J transitions and single-dish-like analysis. We explored two evolutionary stages, 50 kyr (run 1) and 100 kyr (run 2, the reference run). The resulting ζ₂ maps are shown in Fig. 2, top row. We employed Eq. (1) with L = 0.3 pc, which represents the path length (along the line of sight) over which the column densities are integrated. In our case, this corresponds to the size of the simulated box (0.3 pc × 0.3 pc × 0.3 pc). A detailed discussion on how to choose L and the derived uncertainty is presented in Sect. 4.3 and 5.

The obtained ζ₂ values span the range ≈(2–8) × 10⁻¹⁷ s⁻¹, with medians in the denser parts of the core of ζ₂ = (3.4 ± 0.6) × 10⁻¹⁷ s⁻¹ (50kyr) and ζ₂ = (4.0 ± 0.4) × 10⁻¹⁷ s⁻¹ (100 kyr), as shown in Fig. 1. These values should be compared with the actual ${\zeta }_2^{\text{s}}=2.5\times {10}^{-17}$ s⁻¹. We conclude that, in these runs, the BFL20 reproduces the ${\zeta }_2^{\text{S}}$ within a factor 1.5–1.6 on average.

Concerning the morphology of the retrieved ζ₂ maps, the one at 50 kyr shows a smaller scatter around the median value than the run at 100 kyr (see Fig. 1), mainly because we can infer ζ₂ only for positions where N(H₂) > 4.8 × 10²² cm⁻². This is because at this early stage, the oH₂D+ abundance is at most X(oH₂D⁺) = 7 × 10⁻¹¹, producing a line peak intensity of 0.5 K¹¹. For comparison, at 100kyr, the oH₂D⁺ abundance reaches 4 × 10⁻¹⁰, and the transition is as bright as 3 K. This limits the area where ζ₂ is computed with S/N > 3 in run 1.

The histogram from run 2 spans a larger range of values than run 1 and presents a tail at higher values, because the retrieved ζ₂ map presents an increase with increasing distance from the core centre, especially in the northern and western directions (cf. top-right panel of Fig. 2). A further enhancement up to ζ₂ ~ 9 × 10⁻¹⁷ s⁻¹ is visible in the south-eastern part of the source (note that this does not affect the histogram, which focuses on the high H₂ column density region to improve the S/N). In Sect. 5, we discuss more in detail the implication of the spatial trends seen in the ζ₂ maps.

4.1.2 Run 3: low ${\zeta }_2^{\text{S}}$ and single-dish like analysis

To further expand the range of input ${\zeta }_2^{\text{S}}$ values explored and to be consistent with the tests performed by Bovino et al. (2020), we analysed an additional simulation where the ζ₂ is kept constant on the value 2.5 × 10⁻¹⁸ s⁻¹ (low ${\zeta }_2^{\text{S}}$ case). We considered the evolutionary time of 100 kyr. The radiative transfer was performed as described in Sect. 3.3, adopting LTE conditions and focusing on the high-J transitions. The post-processing was single-dishlike, with a convolution beam size of 27″. The setup, hence, is identical to the reference run no. 2 except for the input ${\zeta }_2^{\text{S}}$ value and the injected noise level. With this ${\zeta }_2^{\text{S}}$ value, the deuteration process is slow, and the abundances of deuterated species (H₂D⁺, DCO⁺) at 100 kyr are orders of magnitudes lower than in the tests with ${\zeta }_2^{\text{S}}={10}^{-17}-{10}^{-16}{\text{s}}^{-1}$. We reduced the simulated noise level to rms = 1 mK in the post-processing, to compute the column density of all species significantly.

Using the BFL20 method, we obtained the map shown in the bottom-left panel of Fig. 2. The corresponding histogram of the distribution of values is shown in Fig. 1. The computed values are in the range (2–8) × 10⁻¹⁸ s⁻¹, and the median value of (4.8 ± 0.8) × 10⁻¹⁸ s⁻¹ is less than a factor of two higher than the input ${\zeta }_2^{\text{S}}$.

Fig. 2 Resulting ζ2 maps obtained with Eq. (1) in runs 1 to 4. The run ID is reported in the top-left corner, and the key parameters are included at the top of each panel. Run 2 is taken as a reference throughout the rest of this work. The white contours show the 50, 70, and 90% of the N(H2) peaks, which are: 6.67 × 1022 cm−2 (core in single-dish like analysis at 100 kyr); 6.86 × 1022 cm−2 (core in single-dish-like analysis at 50 kyr); 2.51 × 1022 cm−2 (core in ALMAlike analysis at 100 kyr). The beam size and scalebar are shown in the bottom corners of each panel. Note that we show a zoom-in of the central 0.15 pc.

4.1.3 Run 4: ALMA-like analysis

Here, we focus on the case with constant ${\zeta }_2^{\text{S}}=2.5\times {10}^{-17}{\text{s}}^{-1}$ and an ALMA-like post-processing as described in Sect. 3.4.2. It is important to discuss the chosen value of L. In the single-dish-like analysis, the integrated intensity (or optical depth) is computed along the whole simulated line of sight, namely along the whole length of the simulation box (0.3 pc). In the ALMA-like analysis, on the other hand, this is not the case. Once the SIMOBSERVE task is run, the interferometer acts as a low spatial-frequency filter, and the emission over scales larger than the so-called maximum recoverable scale (θ_mrs) is filtered out. We hence select the θ_mrs = 15″ that the ALMA OT predicts for the chosen antenna configuration in the oH₂D⁺ (1_1,0 – 1_1,1) setup. This corresponds to L = 0.15 pc.

The resulting ζ₂ map is shown in the bottom-right panel of Fig. 2. The histogram of run 4 (top-right panel of Fig. 1) is asymmetric, with a global maximum at low ζ₂ values (1 × 10⁻¹⁷ s⁻¹), and a tail up to 7 × 10⁻¹⁷ s⁻¹. This is due to the spatial gradient seen in the bottom-right panel of Fig. 2, where ζ₂ increases as N(H₂) decreases. In the central part of the core, the actual ${\zeta }_2^{\text{S}}$ is well recovered, as confirmed by comparing the median value 〈ζ₂〉 = (2.4 ± 0.4) × 10⁻¹⁷ s⁻¹ with ${\zeta }_2^{\text{S}}=2.5\times {10}^{-17}{\text{s}}^{-1}$.

4.1.4 Runs 5–6: LVG radiative transfer and low-J transition

Runs 5 and 6 explored the effects of the type of radiative transfer performed and of the rotational levels of the lines used in the analysis. Both runs used the LVG option in POLARIS. Run 5 employed the DCO⁺, H¹³CO⁺, and DCO⁺ transitions at frequencies ≈215–260 GHz. The line intensities are generally lower than in the corresponding LTE calculation. The change is the smallest for the C¹⁸O (2−1) line (≈15%), which is expected as this transition is thermalised. On the contrary, the DCO⁺ and H¹³CO⁺ (3−2) fluxes are reduced by a factor of up to 3 and 10, respectively. In fact, due to their high critical densities, these transitions are subthermally excited. This required reducing the simulated noise in this run to 50 mK, to improve the final S/N. Table 2 reports the excitation temperatures used to compute the column densities, following what is stated in Sect. 3.5.

The resulting ζ₂ map is shown in the central panel of Fig. 3, and it spans values in the range (0.5–2.5) × 10⁻¹⁷ s⁻¹. The median in the high-density part of the core (see Fig. 1) is 〈ζ₂〉 = (1.1 ± 0.2) × 10⁻¹⁷ s⁻¹, hence a factor 2.3 smaller than the actual ${\zeta }_2^{\text{S}}$ value of the simulations. The ζ₂ distribution is asymmetric, as a consequence of the increasing trend of ζ₂ with decreasing N(H₂) that has also been noted in the previous runs at 100 kyr.

In run 6, we explored the effect of targeting the low-J transitions of C¹⁸O, H¹³CO⁺, and DCO⁺, as these lines are often targeted by spectroscopic surveys at 3 mm. We adopt a large beam size of ≈70″, to simulate unresolved observations. The resulting ζ₂ map, shown in Fig. 3, presents the flattest distributions of values, with the smallest scatter around the median (see Fig. 1). This happens because the area where we recover the ζ₂ map is comparable to the beam size, and hence any spatial trend is smoothed out by the poor resolution. The resulting median 〈ζ₂〉 = (3.3 ± 0.7) × 10⁻¹⁷ s⁻¹ agree with the actual ${\zeta }_2^{\text{S}}=2.5\times {10}^{-17}{\text{s}}^{-1}$ within a factor of 1.3.

Fig. 3 Summary of runs performed in LVG approximation, compared with the reference case (run 2, left panel). The central panel refers to the analysis performed on the high-J transitions for DCO+, H13CO+, and C18O, whilst the low-J ones are used in the map in the right panel. The white contours are the same as in Fig. 2. Note that the colour scale is the same across all the panels. The beam size and scalebar are shown in the bottom left and right corners of each panel.

Fig. 4 Maps of the density-averaged ${\zeta }_2^{\text{S}}$ in the simulation with variable cosmic rays (top row) at 50 kyr (left) and 100 kyr (right). Bottom row is the same as Fig. 2, but for runs 7 and 8. Note: the colour scale is the same within rows, but it changes from the top panels to the bottom ones.

4.1.5 Runs 7–8: variable ${\zeta }_2^{\text{S}}$

Runs 7 and 8 employed a variable ${\zeta }_2^{\text{S}}$, as described in Sect. 3.1. In this case, it is not straightforward to compare the resulting ζ₂ maps with the simulation value, which is a three-dimensional, spatially-dependent quantity. For this comparison, we computed the line-of-sight density-averaged ${\zeta }_2^{\text{S}}$ at timesteps 50 and 100 kyr. The results are shown in the top row of Fig. 4. The maps show that the cosmic-ray ionisation rate decreases from ≈2 × 10⁻¹⁶ s⁻¹ at low densities down to ≈1.5 × 10⁻¹⁶ s⁻¹ in the core’s centre, with a variation of 25%. Moreover, the average value in these simulations is almost one order of magnitude larger than in those with constant ζ₂, offering us the chance to test a high ζ₂ case.

The maps of ζ₂ computed using Eq. (1) are shown in Fig. 4 (bottom panels). In general, they tend to underestimate the simulation values, and the disagreement is larger at the earlier timestep. Overall, our results underestimate the actual ζ₂ of a factor ≈3. Concerning the decrease of ζ₂ with increasing total column density, we note that at 50 kyr the extension of the retrieved ζ₂ map corresponds to only a few beams, and no clear spatial trend is seen. In the later timestep, a positive gradient is visible from the core’s centre to the westernmost outskirts, but no symmetric trend is visible in the other directions. A localised enhancement (ζ₂ = (7–9) × 10⁻¹⁷ s⁻¹) is seen in the south-east corner of the source, with no clear counterpart in the actual ${\zeta }_2^{\text{S}}$ map, where instead a local decrease is visible in this area. We conclude that the resulting ζ₂ map does not reproduce the morphology of the actual one, as we further discuss in Sect. 5.2, but provides an accurate average estimate of ${\zeta }_2^{\text{S}}$. In addition, we note that the ${\zeta }_2^{\text{S}}$ gradient in the original simulations is smaller than the intrinsic error of the analytical formula.

Fig. 5 Maps of the key quantities employed by CWT98: DCO+ deuteration level (top-left panel), CO depletion factor (top-right), electronic fraction (bottom-left), and ζ2 (bottom-right). Note that the colour scales for the last two quantities are in units of 10−6 and 10−13 s−1, respectively. These maps assume the molecular column density and total gas column density derived in run 2. The mean uncertainties are 2.6 × 10−4 (RD), 0.3 (fD), 3 × 10−6 $\left(x_{{\text{e}}^{-}}\right)$, and 1.5 × 10−13 s−1 (ζ2). The beam size and scalebar are shown in the bottom-left panel. The contours show N(H2), as in Fig. 2.

4.2 Using the analytical method from Caselli et al. (1998)

We now focus on the analytical method proposed by Caselli et al. (1998, CWT98 hereon) to test the limitations already discussed in Caselli et al. (2002), providing robust evidence via an accurate methodology. The CWT98 approach has the advantage of depending on commonly observed tracers. We note that the equations depend on the H₂ volume density, which we estimate as n(H₂) = N(H₂)/L, where L is set on the same value used for the BFL20 method for each run.

4.2.1 The reference run 2

For the sake of observational applicability, we tested the behaviour of the original equations (Eqs. (3) and (4) of CWT98). Here, we present the results for the reference case (run 2). In Fig. 5, we show the relevant required quantities, in particular, the deuteration level of HCO⁺ (top left panel), and the CO depletion factor $$f_{\text{D}}=\frac{X^{\text{st}}(\text{CO})}{X(\text{CO})},$$

where X^st(CO) = 1.2 × 10⁻⁴ is the CO standard abundance. The HCO⁺ deuteration level peaks at R_D = (6.9 ± 0.3) × 10⁻³ towards the core’s centre, where the CO depletion reaches f_D = (15.8 ± 0.6). Hence, the deuteration level is smaller than the values spanned by the cores of CWT98, but it fulfils the requirement R_D < 0.023f_D under which the equations can be applied. We further discuss this in Sect. 5. The derived values for the electronic fraction, shown in the bottom-left panel, are in the range (1 – 10) × 10⁻⁶. These are overestimated by more than two orders of magnitude compared to the original simulations, where x(e⁻) ≈ a few × 10⁻⁸. This error propagates to the resulting ζ₂ map (bottom-right panel of Fig. 5). The equation overestimates the original ${\zeta }_2^{\text{S}}$ value by more than three orders of magnitudes, especially at the core’s outskirts.

The original method made several simplifications and assumptions, such as, for example, the reaction rates at constant temperature (10 K), and the lack of ortho- and para-state separation. Furthermore, several reaction rates have been updated since then. We have, therefore, derived the equations again using the same theoretical approach of CWT98, but with the formalism of BFL20, to show that the large overestimates produced by the method are not ascribed to these parameters but rather to the approximations made to obtain the formula. The derivation is illustrated in Appendix B. We then computed the electronic fraction and the cosmic-ray ionisation rate using the updated set of Eq. (B.6). We set the HD abundance to X(HD) = 1.5 × 10⁻⁵ (Kong et al. 2015), and the para-fraction of ${\text{H}}_3^{+}$ to f_para = 0.7. The latter is consistent with estimates of this parameter in diffuse clouds (see for instance, Crabtree & McCall 2012, and references therein). These values have also been verified in the simulation, and they agree within less than a factor of two (see also Lupi et al. 2021).

The resulting maps are shown in Fig. 6. Towards the core’s centre, the x(e⁻) values are 15–20% lower than those derived with the original equations, but still strongly overestimated. As a consequence, in this area, the new estimates for ζ₂ are a factor of ≈2 lower than those from the original derivation, but we still find ζ₂ ≈ (1.5 ± 0.2) × 10⁻¹³s⁻¹, namely more than three orders of magnitude higher than the actual value.

4.2.2 CWT98 results in all tested runs

We now describe the behaviour of the CWT98 method applied to all the remaining runs. We adopt the new formulation of the method, described in Appendix B. The histograms of the resulting maps, focussing on the central part of the core, are presented in Fig. 7. The most notable feature is that in all tests the retrieved values overestimate the actual ones. The median values in the runs performed on the simulation using ${\zeta }_2^{\text{s}}=\cos t=2.5\times {10}^{-17}{\text{s}}^{-1}$ range from (4.8 ± 1.4) × 10^–15 s⁻¹ (run 5) to (9.7 ± 1.5) × 10⁻¹³ s⁻¹ (run 1), namely an overestimation of two to four orders of magnitude. In the low ζ₂ case (run 3), we obtain 〈ζ₂〉 = (1.41 ± 0.07) × 10⁻¹¹ s⁻¹ (overestimated by almost seven orders of magnitude). The two tests performed on simulations with ${\zeta }_2^{\text{s}}=\mathrm{var}\approx {10}^{-16}{\text{s}}^{-1}$ result in (1.68 ± 0.11) × 10⁻¹³ s⁻¹ (run 7) and (3.1 ± 0.2) × 10⁻¹⁴ s⁻¹ (run 8), again over-estimating the actual ζ₂ by two to three orders of magnitude.

Fig. 6 Electronic fraction (top panel) and ζ2 map (bottom panel) from the reference run 2, using the updated formulation of CWT98. The contours show H2 column density as in Fig. 2. The mean uncertainties are 3 × 10–6 for x(e–) and 0.8 × 10–13 s–1 for ζ2. The beam size and scalebar are shown in the bottom corners of the lower panel.

4.3 Comparison of the methods on real observations

We now compare the two analytical methods on real observations of prestellar cores found in the literature. We found three sources for which all the needed observables are available: L1544, L183, and TMC-1C. The literature values of the required quantities are listed in Table 4. For these cores, the H₂ volume density is well characterised, and we used this quantity directly in the CWT98 Eq. (B.6). For the BFL20 method, a discussion on the parameter L is required. The physical meaning of L is the length of the path on the line-of-sight along which the column densities are computed; in other words, L is the depth of the emitting source. This, in the simulations, is known by construction. In our setup, we cut a (0.3 pc)³ subcube in the initial larger simulation box (0.6 pc) which is initialised with molecular gas; thus, it is full of CO, for instance. Emitting gas, therefore, is found along the whole box, and the choice of L to be equal to the size of the subcube is well justified. If we were to cut a smaller subcube, L should be adjusted accordingly, since the depth of the emitting gas would also be reduced (see also Appendix D for further details). This does not apply to real observations. Real cores are finite, and the length of the emitting gas is limited along the length of the line of sight. Thus, L has to be computed as the source size along the line of sight, which however is unknown. For isolated prestellar cores, we propose to compute L by considering the 20% isocontour of the N(H₂) peak. Using this prescription, the sizes of the three analysed cores are 0.20–0.35 pc. The final ζ₂ values are also reported in Table 4.

The ζ₂ values computed with BFL20 are in the range (0.71.7) × 10⁻¹⁷s⁻¹, whilst with CWT98 we obtain (0.2–2.0) × 10⁻¹⁴ s⁻¹. Note that there are many uncertainties in the analysis. For instance, the literature values have been computed with data from different telescopes (hence at different resolutions). Furthermore, we assume T = 10 K, whilst some of the sources might be colder (cf. Caselli et al. 2008). However, the uncertainties on these quantities cannot explain a difference of three orders of magnitude between the two methods. These examples, hence, confirm that CWT98 tends to produce overestimated results compared to BFL20.

The actual ζ₂ value in these objects is not known, however Redaelli et al. (2021b) derived 3.0 × 10⁻¹⁷ s⁻¹ in L1544, and Fuente et al. (2019) found ζ₂ = (5−18) × 10⁻¹⁷ s⁻¹ in the translucent cloud associated with TMC1. Both papers used extensive modelling of the sources, coupling chemical models with radiative transfer analysis on a large set of molecular tracers. Pagani et al. (2009) explored the chemistry and structure of L183, assuming ionisation rates in the range (0.1−10) × 10⁻¹⁶ s⁻¹ and, even though a definite best-value for this parameter is not given, the authors used ζ₂ = 2 × 10⁻¹⁷ s⁻¹ in their most detailed modelling. Furthermore, the most recent models of CR propagation in the dense gas predict ζ₂ values of at most a few 10⁻¹⁶s⁻¹, unless a local source of CR re-acceleration is present (cf. Padovani et al. 2016, 2018, 2022). It is safe to assume, in summary, that values as high as ζ₂ = 10⁻¹⁴ s⁻¹ are excluded for these quiescent and dense cores.

5 Discussion

5.1 Results and limitations of the methods

In Fig. 8, we summarise the results obtained with the BFL20 method in the eight runs of Sect. 4.1. In the case of constant ${\zeta }_2^{\text{s}}$, we show the average ratio between the derived values and the reference ${\zeta }_2^s=2.5\times {10}^{-17}{\text{s}}^{-1}$ (runs 1, 2, and 4 to 6) or ${\zeta }_2^s=2.5\times {10}^{-18}{\text{s}}^{-1}$ (run 3). For runs 7 and 8, where the input ${\zeta }_2^s$ is variable, we show the median of the ratio between the derived values and the input values (see top panels of Fig. 4). Error bars are computed as three times the median uncertainties over the pixels considered to evaluate the median. As for the histograms in Fig. 1, we consider only the densest part of the core.

Figure 8 shows that the retrieved values are within a factor of 1.5-3 from the actual ones. The offset is not constant, nor systematic. In cases with low and constant ${\zeta }_2^s$, we see that Eq. (1) overestimates the input value by a factor of up to 1.5 (except run 5). On the other hand, for the two runs with variable ${\zeta }_2^s$ (runs 7 and 8), the resulting ζ₂ maps tend to underestimate the actual values by a factor of ≈3. Overall, the BFL20 formula represents a robust and reliable method to estimate the order of magnitude of ζ₂ in dense regions. As for similar analytical methods, even if the BFL20 method shows to be accurate within a factor of 2–3, several aspects should be taken into account when it is applied to actual observations. The first one is that the analytic expression depends on the column density of four molecular tracers. If any of these is affected by a systematic error, this will propagate to the ζ₂ estimation. Column densities strongly depend on the chosen excitation temperature values. By performing an LTE analysis, we have initially avoided this problem, fixing the T_ex for all the molecular tracers to the constant gas temperature. In the LVG runs, we selected T_ex looking for literature references and checking the selected values with non-LTE tools (such as RADEX). Indirectly, this work hence provides good estimates of the T_ex of several commonly observed transitions, in the considered density (n ≈ 10⁵ cm⁻³) and gas temperature (T_K = 15 K) regimes. In reality, the problem of choosing the correct T_ex has no straightforward solution, especially in the case of subthermally-excited lines. We strongly suggest (when possible) using multiple lines of the same tracer, which allows us to constrain their excitation conditions better.

When using optically-thin isotopologues to infer the total column densities of molecular tracers, particular caution has to be paid to the assumed isotopic ratios, since fractionation processes can lead to significant variation from the elemental iso-topic ratios (see, for instance, the discussion of Colzi et al. 2020 on ¹²C/¹³C). In this work, we avoided this problem by assuming consistent isotopic ratios throughout the postprocessing analysis (cf. Sects. 3.1 and 3.3).

Another crucial parameter in Eq. (1) is the scale length L employed to obtain the column densities, as discussed in Sect. 4. In our analysis, L is known by construction from the size of the simulation box, since this is a subcube cut out of a larger simulation initialised with molecular gas. We justify this choice further in Appendix D. For observed cores, their size (and, in particular, their depths) are in general unknown. We suggest as a prescription to use the contour where N(H₂) is higher than 20% of its peak value to estimate L. In the three real objects we analysed, this approach leads to ζ₂ values in agreement with estimations of this quantity derived with completely different methods. Note that applying this rule to the reference run 2 results in L = 0.15 pc, hence overestimating ζ₂ by a factor ≈2, still within the uncertainties of the method. This prescription should not be used on interferometric data with a maximum recoverable scale smaller than the actual source size, as this leads to emission filtering. In these cases, we speculate that θ_mrs is a proper scale to estimate L. This uncertainty on L also affects the approach of CWT98, where however this is mitigated if the gas volume density is known from other observables.

We highlight how the core fraction where we are able to infer ζ₂ is different in each case (see Figs. 2, 3, and 4). This is driven, in particular, by the column density of oH₂D+^. At earlier times, the abundance of this species is lower, hence the extension of the core where it is detected is smaller. This highlights the main limitation of this method, which relies on the detection of oH₂D⁺(1_1.0 – 1_1.1) (see also Sect. 5.3 for details).

On the other hand, the analytic method of CWT98 overestimates ζ₂ by several orders of magnitudes in the cases explored in this work, where it should not be employed (see also the discussion in Sect. 3 of Caselli et al. 2002). This is due to the overestimation of the electron fraction caused by the neglected kinetics of H₂D⁺ employed to derive Eq. (B.2), in particular, the reactions producing the doubly- and triply-deuterated forms of H⁺, and other destruction channels involving neutrals. Concerning the importance of D₂H⁺ and ${\text{D}}_3^{+}$, Caselli et al. (2008) developed analytical equations where the deuteration level of HCO⁺ is expressed in terms of all the deuterated forms of H⁺.

We note that when ζ₂ is derived in the reference case with the analytic expression of CWT98 (10⁻¹⁴−10⁻¹³ s^-1), its value is one to two orders of magnitude higher than in the original paper (where ζ₂ = 10⁻¹⁶−10⁻¹⁴ s⁻¹ was found when assuming ƒ_D = 5). This is because the simulated core presents relatively low levels of deuteration and a high depletion factor. In Appendix B, we show how the CWT98 formula leads to increasingly high ζ₂ values when the deuteration is low, and the CO depletion is high (cf. Fig. B.2). In the reference case (run 2), the core centre presents R_D ≈ 0.4–0.7% and ƒ_D ≈ 10–15. These are not the typical values observed in CWT98, where most sources present a deuteration fraction of a few per cent (see also Butner et al. 1995; Williams et al. 1998). Note that, given the lack of information available back then on the CO depletion, ƒ_D > 5 was not included in their analysis (catastrophic CO freeze-out was first measured 1 yr later, Caselli et al. 1999). Furthermore, the assumed temperature T = 15 K is higher than typical gas temperatures observed in low-mass prestellar cores (Bergin et al. 2006; Crapsi et al. 2007). However, in run 8, with variable ζ₂ at 100 kyr, the core’s centre presents R_D ≈ 2% and ƒ_D ≈ 3–4, closer to the properties of the objects analysed by CWT98. The retrieved ζ₂ is (3.1 ± 0.2) × 10⁻¹⁴ s⁻¹ (see Fig. 7), overestimating the actual ζ₂ by two orders of magnitudes, confirming the limitation of this approach.

The relatively low deuteration level of the reference case is due to the assumed initial conditions, in particular the initial H₂ OPR. By using constant-density one-zone models, we foundthat when OPR = 10⁻³ (as reported in the dense and evolved prestellar cores, for instance Kong et al. 2015), the HCO⁺ deuteration level increases by about one order of magnitude, and the derived ζ₂ values decrease to ζ₂ = 10⁻¹⁶−10⁻¹⁵ s⁻¹, still a factor of 10-100 more than the actual value. These tests show that the CWT98 analytic method has a marked dependence on the initial OPR, conversely to BFL20, which is relatively unaffected by this parameter, as already discussed in Bovino et al. (2020). The scope of this work is to retrieve ζ2 under the typical observational conditions while reproducing the exact physical details of a specific observed object is beyond our aims. The latter was done, for example, by Bovino et al. (2021), where the simulation was designed to reproduce six observed cores in Ophiucus.

Fig. 7 Same as Fig. 1, but the ionisation rate is estimated using the new formulation of the CWT98 method presented in Appendix B, summarised in Eq. (B.6). Note that the ζ2 values are normalised to 10−13 s−1 in all panels. The median (± median uncertainty) is reported in the top-right corner and shown with the vertical dashed line and shaded area in each panel. The reference run is labelled with "Ref" in the bottom-right corner.

Fig. 8 Ratio between ζ2 retrieved in each run (x-axis labels), and the actual value ${\zeta }_2^s$. The ratios are median computed in each map, in the same regions used to produce the histograms of Fig. 1. The shaded areas show a variation of a factor of 2 (darker shade) or 3 (lighter shade). The dashed line represents where the ratio is equal to 1. The error bars are 3σ median uncertainties.

5.2 Morphology of the resulting ζ₂ maps

The obtained ζ₂ maps allow us to comment also on the morphology of the retrieved ionisation rate. By looking at Figs. 2 and 3, it is clear that for most of the runs where ζ2 is computed in an extended part of the source (namely when N(H₂) ≲50% of its peak), ζ2 presents a positive gradient with increasing distance from the core’s centre. This is also seen in the histograms in Fig. 1, where asymmetric tails towards the high ζ2 values are seen in runs 2, 4, and 5. Since in these runs, the actual ${\zeta }_2^s$ is constant, these spatial trends are not real. These considerations holds also for CWT98 (see Figs. 5 and 6). The single-dish test performed with the larger beam size (run 6) presents the flattest distribution and the smallest scatter around the median because the core is not spatially resolved. On the contrary, in the tests with variable ${\zeta }_2^s$, the retrieved maps do not show the gradient with the increasing density present in the simulation (see Sect. 4.1.5 for more details). We conclude that apparent spatial trends should not be trusted, and averages across the densest regions of the source should be instead considered to express the resulting ζ2.

5.3 Observational feasibility

This work aims to provide observers with reliable methods to estimate ζ2 in real sources. It is hence important to assess the observability of the proposed tracers. The continuum observations are not particularly challenging and do not represent the limiting aspect of this endeavour. Conversely, the feasibility of the molecular line observations requires a more detailed discussion.

We first focus on single-dish facilities. As reported in Table 2, the H¹³CO⁺ (3–2), DCO⁺ (3–2), ad C¹⁸O (2–1) lines can be observed by APEX, namely with the nFLASH23O instrument. Using its observing time calculator¹², and standard input values, we compute that on-source times of 1.0–1.4 h are sufficient to reach rms = 100 mK in a 3′ × 3′ FoV and a 0.1 km s^–1 channel. The corresponding (1-0) transitions are covered, for instance, by the EMIR receiver mounted on the IRAM 30m telescope. The time estimator¹³ predicts that, in winter, 1 h of on-source time is enough to reach a sensitivity of 100 mK, with a 3′ × 3′ FoV and using a 0.1 km s^–1 resolution.

The oH₂D⁺(1₁₀ – 1_1,1)transition is the most challenging in terms of sensitivity, due to its high frequency when compared to the other transitions analysed here. The same requirements in terms of spectral resolution, sky coverage, and sensitivity made for the other lines would lead to 25 h of on-source time using the multi-beam LAsMA receiver mounted on APEX. However, by relaxing the requirements to a FoV of 2′ × 2′ (still able to cover the portion of the core where N(H₂) ≳ 2.5 × 10²² cm^–2) and downgrading the resolution to 0.15 km s⁻¹, the on-source time reduces to 7.5 h, which is manageable in a small project. The downgraded spectral resolution does not impact the computation of the column density, since the lines are still resolved by at least 3 channels. The sensitivity required by the low ζ₂ case in run 3(1 mK), on the other hand, is currently well beyond the capabilities of existing single-dish facilities, even in the case of single-point observations at the centre of the source.

The observing times required with ALMA are listed in Table 3. Considering that several lines can be observed simultaneously, these observations appear feasible. We conclude that in most cases the observations required to compute ζ2 are feasable, as proved by the increasing number of campaigns recently reported (cf. Giannetti et al. 2019; Sabatini et al. 2020; Redaelli et al. 2021a, 2022).

6 Summary and conclusions

In this work, we tested two analytical methods to retrieve the cosmic-ray ionization rate ζ₂ in dense gas. We used synthetic molecular and continuum data, produced via radiative transfer analysis on a set of 3D simulations that include the chemistry of the involved molecular tracers. This allows us to evaluate with accuracy the loss of information (and then the accuracy of the method) when simulating realistic observations.

The method of Caselli et al. (1998) has several limitations by construction, such as R_D < 0.029 × ƒ_D (to avoid x(e) < 0 in the new formulation derived in Appendix B). Furthermore, this analytical approach strongly depends on the H₂ initial OPR. This limits its applicability, especially when the OPR is reset to much higher values than those in cold cores by conditions such as shocks or protostellar outflows. In our reference case, this method overestimates ${\zeta }_2^s$ by up to four orders of magnitude. In particular, in tests where the deuteration level is a few %, hence similar to what is observed in several prestellar cores, the Caselli et al. (1998) method overestimates by two orders of magnitude the actual ${\zeta }_2^s$.

On the contrary, the method of Bovino et al. (2020) is generally accurate (within a factor of 2-3) in retrieving the actual ${\zeta }_2^s$. Its main limitation is linked to the level of total deuteration, since at late evolutionary stages or at very high densities (n ≳ 10⁷ cm^–3) H₂D⁺ is converted into doubly and triply deuterated forms, and it is not a reliable tracer of the total ${\text{H}}_3^{+}$ abundance anymore. This leads to underestimating the actual ${\zeta }_2^s$, as already pointed out in the original paper (Bovino et al. 2020).

As a direct example of the application of the two formulae on observational data, we explored three well-known prestellar objects, with recent literature data on the quantities involved in the calculation. We showed that the ζ₂ values obtained with BFL20 are in overall good agreement with estimations of the same quantities obtained with non-analytical methods. The results of CWT98 are two to three orders of magnitude higher, as seen also in the tests on the simulations. We highlight, however, that to establish the methodology proposed by Bovino et al. (2020), a statistical sample of observed cores and a proper comparison with theoretical models of CR propagation are needed.

To conclude, we have discussed the feasibility of the observations necessary to use two commonly employed analytical methods to retrieve ζ₂. Despite the observational challenges, they are accessible with currently available radio facilities. When the physical structure of a source is well known, coupling a chemical code with radiative transfer using multiple tracers could be employed to infer the cosmic-ray ionisation rate, even if its results might be affected by the parameters’ degeneracy. For all the other sources (when this approach is not a viable option), the method of Bovino et al. (2020) is a model-independent and reliable analytical method to investigate ζ₂ in dense regions.

Acknowledgements

The authors acknowledge the referee’s comments that led to the manuscript’s improvement. E.R. and P.C. acknowledge the support of the Max Planck Society. S.B. is financially supported by ANID Fondecyt Regular (project #1220033), and the ANID BASAL projects ACE210002 and FB210003. A.L. acknowledges funding from MIUR under the grant PRIN 2017-MB8AEZ. G.S. acknowledges the projects PRIN-MUR 2020 MUR BEYOND-2p ("Astro-chemistry beyond the second period elements", Prot. 2020AFB3FX) and INAF-Minigrant 2023 TRIESTE ("TRacing the chemIcal hEritage of our originS: from proTostars to planEts"). The authors acknowledge the Kultrun Astronomy Hybrid Cluster for providing HPC resources that have contributed to the research results reported in this paper.

Appendix A Derivation of Bovino et al. (2020) (BFL20)

Here, we follow the derivation of Eq. (1). The main reactions in our framework, considering the different isomers and isotopologues (but ${\text{D}}_3^{+}$) for the formation of HCO⁺ and DCO⁺ are: $$\begin{array}{c}{\text{oH}}_2{\text{D}}^{+}+\text{CO}\stackrel{k_1}{\to }{\text{DCO}}^{+}+{\text{oH}}_2,\\ {\text{pH}}_2{\text{D}}^{+}+\text{CO}\stackrel{k_2}{\to }{\text{DCO}}^{+}+{\text{pH}}_2,\\ {\text{pD}}_2{\text{H}}^{+}+\text{CO}\stackrel{k_3}{\to }{\text{DCO}}^{+}+\text{HD},\\ {\text{oD}}_2{\text{H}}^{+}+\text{CO}\stackrel{k_4}{\to }{\text{DCO}}^{+}+\text{HD},\\ {\text{pH}}_3^{+}+\text{CO}\stackrel{k_5}{\to }{\text{HCO}}^{+}+{\text{pH}}_2,\\ {\text{pH}}_3^{+}+\text{CO}\stackrel{k_6}{\to }{\text{HCO}}^{+}+{\text{oH}}_2,\\ \begin{array}{l}{\text{oH}}_3^{+}+\text{CO}\stackrel{k_7}{\to }{\text{HCO}}^{+}+{\text{oH}}_2,\\ \begin{array}{r}\hfill {\text{pD}}_2{\text{H}}^{+}+\text{CO}\stackrel{k_8}{\to }{\text{HCO}}^{+}+{\text{pD}}_2,\\ \hfill {\text{pH}}_3^{+}+\text{CO}\stackrel{k_6}{\to }{\text{HCO}}^{+}+{\text{oH}}_2,\\ \hfill {\text{oH}}_3^{+}+\text{CO}\stackrel{k_7}{\to }{\text{HCO}}^{+}+{\text{oH}}_2,\\ \hfill {\text{pD}}_2{\text{H}}^{+}+\text{CO}\stackrel{k_8}{\to }{\text{HCO}}^{+}+{\text{pD}}_2,\\ \hfill {\text{oD}}_2{\text{H}}^{+}+\text{CO}\stackrel{k_9}{\to }{\text{HCO}}^{+}+{\text{oD}}_2,\\ \hfill {\text{pH}}_2{\text{D}}^{+}+\text{CO}\stackrel{k_{10}}{\to }{\text{HCO}}^{+}+\text{HD},\\ \hfill {\text{oH}}_2{\text{D}}^{+}+\text{CO}\stackrel{k_{11}}{\to }{\text{HCO}}^{+}+\text{HD}.\end{array}\end{array}\end{array}$$

For the destruction, we consider only dissociative recombinations, as per: $$\begin{array}{cc}& {\text{DCO}}^{+}+{\text{e}}^{-}\stackrel{{\beta }_1}{\to }\text{D}+\text{CO},\\ & {\text{HCO}}^{+}+{\text{e}}^{-}\stackrel{{\beta }_2}{\to }\text{H}+\text{CO}.\end{array}$$

The kinetic equations are expressed as: $$\frac{dn\left({\text{DCO}}^{+}\right)}{dt}=k_{1}n(\text{CO})n\left({\text{oH}}_2{\text{D}}^{+}\right)+k_{2}n(\text{CO})n\left({\text{pH}}_2{\text{D}}^{+}\right)+k_{3}n(\text{CO})n\left({\text{pD}}_2{\text{H}}^{+}\right)+k_{4}n(\text{CO})n\left({\text{oD}}_2{\text{H}}^{+}\right)-{\beta }_{1}n\left({\text{DCO}}^{+}\right)n\left({\text{e}}^{-}\right),$$

and $$\begin{array}{cc}\frac{dn\left({\text{HCO}}^{+}\right)}{dt}=& \left(k_5+k_6\right)n(\text{CO})n\left({\text{pH}}_3^{+}\right)+k_{7}n(\text{CO})n\left({\text{oH}}_3^{+}\right)+k_{8}n(\text{CO})n\left({\text{pD}}_2{\text{H}}^{+}\right)+k_{9}n(\text{CO})n\left({\text{oD}}_2{\text{H}}^{+}\right)+\\ & +k_{10}n(\text{CO})n\left({\text{pH}}_2{\text{D}}^{+}\right)+k_{11}n(\text{CO})n\left({\text{oH}}_2{\text{D}}^{+}\right)-{\beta }_{2}n\left({\text{HCO}}^{+}\right)n\left({\text{e}}^{-}\right).\end{array}$$

Assuming steady-state¹⁴ and taking the ratio between the two equations we obtain $$R_D=\frac{n\left({\text{DCO}}^{+}\right)}{n\left({\text{HCO}}^{+}\right)}=\left(\frac{{\beta }_2}{{\beta }_1}\right)\times \frac{k_{1}n(\text{CO})n\left({\text{oH}}_2{\text{D}}^{+}\right)+k_{2}n(\text{CO})n\left({\text{pH}}_2{\text{D}}^{+}\right)+k_{3}n(\text{CO})n\left({\text{pD}}_2{\text{H}}^{+}\right)+k_{4}n(\text{CO})n\left({\text{oD}}_2{\text{H}}^{+}\right)}{\left(k_5+k_6\right)n(\text{CO})n\left({\text{pH}}_3^{+}\right)+k_{7}n(\text{CO})n\left({\text{oH}}_3^{+}\right)+k_{8}n(\text{CO})n\left({\text{pD}}_2{\text{H}}^{+}\right)+k_{9}n(\text{CO})n\left({\text{oD}}_2{\text{H}}^{+}\right)+k_{10}n(\text{CO})n\left({\text{pH}}_2{\text{D}}^{+}\right)+k_{11}n(\text{CO})n\left({\text{oH}}_2{\text{D}}^{+}\right)}.$$

Using the following relations among the reaction rates $$\begin{array}{l}\begin{array}{cc}& k_1=k_2,\\ & k_3=k_4,\\ & k_5+k_6=k_7,\\ & k_8=k_9,\\ & k_{10}=k_{11},\end{array}\\ {\beta }_2={\beta }_1.\end{array}$$

The final equation is expressed as (noting that n(CO) and n(e⁻) cancel each other out): $$R_D=\frac{n\left({\text{DCO}}^{+}\right)}{n\left({\text{HCO}}^{+}\right)}=\frac{k_1\left[n\left({\text{oH}}_2{\text{D}}^{+}\right)+n\left({\text{pH}}_2{\text{D}}^{+}\right)\right]+k_3\left[n\left({\text{pD}}_2{\text{H}}^{+}\right)+n\left({\text{oD}}_2{\text{H}}^{+}\right)\right]}{k_7\left[n\left({\text{oH}}_3^{+}\right)+n\left({\text{pH}}_3^{+}\right)\right]+k_8\left[n\left({\text{pD}}_2{\text{H}}^{+}\right)+n\left({\text{oD}}_2{\text{H}}^{+}\right)\right]+k_{10}\left[n\left({\text{pH}}_2{\text{D}}^{+}\right)+n\left({\text{oH}}_2{\text{D}}^{+}\right)\right]}.$$(A.1)

In order to simplify Eq. (A.1), we can exploit further relations between the reaction rates, mainly linked to their branching ratios: k₁/k₇ = 1/3, k₁/k_l0 = 1/2, k₈/k₃ = 1/2, and k₃/k₇ = 2/3. Moreover, we neglect the correction for para and ortho species that cannot be observed, the contribution from the doubly deuterated isotopologue, and the formation channel of HCO⁺ via H₂D⁺, arriving at¹⁵ $$n\left({\text{H}}_3^{+}\right)=\frac{1}{3}\frac{n\left({\text{oH}}_2{\text{D}}^{+}\right)}{R_D}.$$(A.2)

By inserting Eq. (A.2) in (see for instance Oka 2019, and also the derivation in Appendix B) as follows: $${\zeta }_2=\frac{k_{7}n\left({\text{H}}_3^{+}\right)n(\text{CO})}{n\left({\text{H}}_2\right)}\text{, }$$(A.3)

we obtain the final formula for the cosmic-ray ionisation rate $${\zeta }_2=\frac{1}{3}\frac{k_{7}n\left({\text{oH}}_2{\text{D}}^{+}\right)}{R_D}\frac{n(\text{CO})}{n\left({\text{H}}_2\right)}.$$(A.4)

This equation is valid as long as oH₂D⁺ is the dominant deuterated form of ${\text{H}}_3^{+}$ and when the reaction with CO is more important than dissociative recombination in the destruction of H⁺₃. The first limitation implies that when deuteration levels become higher and oH₂D⁺ is converted into its doubly and triply deuterated isotopologues, Eq. (A.4) cannot be used anymore. Concerning the destruction pathways, we can investigate at which CO abundance (as a function of the electronic fraction) its reaction with H₃⁺ dominates over the dissociative recombination (see the right-hand side of Eq. (B.5)). We find that this holds for f_D < 9 × l0⁻⁷/x(e⁻) (assuming f_para = 0.7, see also Appendix B). For x(e) = 10⁻⁸, close to the values found in the reference run, the reaction with CO is dominant if f_D ≲ 90, which is verified in our simulations. However, at very high densities, when f_D > 100, this assumption might not hold anymore.

By introducing average quantities integrated over the path L along the line of sight, Eq. (A.4) finally becomes $${\zeta }_2=k_7\frac{1}{3}\times X(\text{CO})\times \frac{N_{\text{col}}\left({\text{oH}}_2{\text{D}}^{+}\right)}{R_{\text{D}}}\frac{1}{L}.$$(A.5)

We note that $k_7=k_{\text{CO}}^{{\text{oH}}_3^{+}}$ is given in the main text.

Appendix B Derivation of Caselli et al. (1998) (CWT98)

We now illustrate the derivation of the equations used in Caselli et al. (1998), which in turn are based on previous works (for instance Guelin et al. 1977, 1982). In particular, we aim to follow the same approach as those papers, including this time the ortho and para separation for all the involved species.

The first part of the equations is the same as illustrated in Appendix A, and involves balancing the destruction and formation pathways of HCO⁺ and DCO⁺ , arriving at (see Eq. (A.1)) $$R_D=\frac{1}{3}\frac{n\left({\text{H}}_2{\text{D}}^{+}\right)}{n\left({\text{H}}_3^{+}\right)},$$(B.1)

where we neglect the reactions involving doubly and triply deuterated ${\text{H}}_3^{+}$ and reactions 10 and 11). In this case, however, we look for a way to express the ratio $n\left({\text{H}}_2{\text{D}}^{+}\right)/n\left({\text{H}}_3^{+}\right)$ as a function of the densities of CO, HD, and of the electron fraction. To this aim, we have to write the reactions involved in the formation and destruction of H₂D⁺ (in its ortho and para states). For the formation pathways, we have: $$\begin{array}{cc}& {\text{pH}}_3^{+}+\text{HD}\stackrel{k_{12}}{\to }{\text{pH}}_2{\text{D}}^{+}+{\text{pH}}_2,\\ & {\text{pH}}_3^{+}+\text{HD}\stackrel{k_{13}}{\to }{\text{pH}}_2{\text{D}}^{+}+{\text{oH}}_2,\\ & {\text{pH}}_3^{+}+\text{HD}\stackrel{k_{14}}{\to }{\text{oH}}_2{\text{D}}^{+}+{\text{pH}}_2,\\ & \begin{array}{l}{\text{pH}}_3^{+}+\text{HD}\stackrel{k_{15}}{\to }{\text{oH}}_2{\text{D}}^{+}+{\text{oH}}_2,\\ \begin{array}{cc}& {\text{oH}}_3^{+}+\text{HD}\stackrel{k_{16}}{\to }{\text{pH}}_2{\text{D}}^{+}+{\text{oH}}_2,\\ & {\text{oH}}_3^{+}+\text{HD}\stackrel{k_{17}}{\to }{\text{oH}}_2{\text{D}}^{+}+{\text{oH}}_2,\\ & {\text{oH}}_3^{+}+\text{HD}\stackrel{k_{18}}{\to }{\text{oH}}_2{\text{D}}^{+}+{\text{pH}}_2.\end{array}\end{array}\end{array}$$

The destruction pathways, instead, involve the reactions with CO (with rates k₁, k₂, k₁₀, and k₁₁, see above)¹⁶ and the following dissociative recombinations: $$\begin{array}{cc}& {\text{oH}}_2{\text{D}}^{+}+{\text{e}}^{-}\stackrel{{\beta }_3}{\to }\text{H}+\text{H}+\text{D},\\ & {\text{oH}}_2{\text{D}}^{+}+{\text{e}}^{-}\stackrel{{\beta }_4}{\to }\text{D}+{\text{oH}}_2,\\ & \begin{array}{l}{\text{oH}}_2{\text{D}}^{+}+{\text{e}}^{-}\stackrel{{\beta }_5}{\to }\text{H}+\text{H},\\ \begin{array}{cc}& {\text{pH}}_2{\text{D}}^{+}+{\text{e}}^{-}\stackrel{{\beta }_6}{\to }\text{H}+\text{H}+\text{D},\\ & {\text{pH}}_2{\text{D}}^{+}+{\text{e}}^{-}\stackrel{{\beta }_7}{\to }\text{D}+{\text{pH}}_2,\\ & {\text{pH}}_2{\text{D}}^{+}+{\text{e}}^{-}\stackrel{{\beta }_8}{\to }\text{H}+\text{HD}.\end{array}\end{array}\end{array}$$

We note that we have neglected all the reactions involving doubly and triply deuterated H⁺, so as to remain consistent with the simplification done to obtain Eq. (B.1). For several of the involved reaction rates, it is possible to show that $$\begin{array}{ll}\hfill & k_{12}+k_{13}+k_{14}+k_{15}\approx k_{16}+k_{17}+k_{18},\hfill \\ \hfill & {\beta }_3\approx {\beta }_6,\hfill \\ \hfill & {\beta }_4\approx {\beta }_7,\hfill \\ \hfill & {\beta }_5\approx {\beta }_8.\hfill \end{array}$$

These relations do not hold exactly, but we will show that in the temperature range here considered the agreement is reasonably good. The first relation is reported in the left panel of Fig. B.1. For temperatures ≲ 25 K, the discrepancy is lower than 10%, and in the range 10 − 15 K, the difference is 1 − 3%. We hence assume that equality holds. For the various rates of dissociative recombination (β₃ to β₈), the difference at 15K is ≈6%, but it quickly rises above 25% outside the range 10 – 20K. We hence suggest extreme caution in using these and the following relations outside this temperature range. However, without these assumptions, it is in practice impossible to properly re-derive and upgrade the formula proposed by CWT98.

We can now write the kinetic equations for the para and ortho species separately as: $$\begin{array}{ll}\frac{dn\left({\text{oH}}_2{\text{D}}^{+}\right)}{dt}\hfill & =k_{14}n\left({\text{pH}}_3^{+}\right)n(\text{HD})+k_{15}n\left({\text{pH}}_3^{+}\right)n(\text{HD})+k_{17}n\left({\text{oH}}_3^{+}\right)n(\text{HD})+k_{18}n\left({\text{oH}}_3^{+}\right)n(\text{HD})\hfill \\ \hfill & -{\text{k}}_{1}n\left({\text{oH}}_2{\text{D}}^{+}\right)n(\text{CO})-k_{11}n\left({\text{oH}}_2{\text{D}}^{+}\right)n(\text{CO})-{\beta }_{3}n\left({\text{oH}}_2{\text{D}}^{+}\right)n\left({\text{e}}^{-}\right)-{\beta }_{4}n\left({\text{oH}}_2{\text{D}}^{+}\right)n\left({\text{e}}^{-}\right)-{\beta }_{5}n\left({\text{oH}}_2{\text{D}}^{+}\right)n\left({\text{e}}^{-}\right),\hfill \end{array}$$

and $$\begin{array}{ll}\frac{dn\left({\text{pH}}_2{\text{D}}^{+}\right)}{dt}\hfill & =k_{12}n\left({\text{pH}}_3^{+}\right)n(\text{HD})+k_{13}n\left({\text{pH}}_3^{+}\right)n(\text{HD})+k_{16}n\left({\text{oH}}_3^{+}\right)n(\text{HD})-k_{2}n\left({\text{pH}}_2{\text{D}}^{+}\right)n(\text{CO})+\hfill \\ \hfill & -k_{10}n\left({\text{pH}}_2{\text{D}}^{+}\right)n(\text{CO})-{\beta }_{6}n\left({\text{pH}}_2{\text{D}}^{+}\right)n\left({\text{e}}^{-}\right)-{\beta }_{7}n\left({\text{pH}}_2{\text{D}}^{+}\right)n(\text{e})-{\beta }_{8}n\left({\text{pH}}_2{\text{D}}^{+}\right)n\left({\text{e}}^{-}\right).\hfill \end{array}$$

Assuming the steady state, we can re-write the two equations above as: $$\begin{array}{l}\hfill \left(k_{14}+k_{15}\right)n\left({\text{pH}}_3^{+}\right)n(\text{HD})+\left(k_{17}+k_{18}\right)n\left({\text{oH}}_3^{+}\right)n(\text{HD})=\left(k_1+k_{11}\right)n\left({\text{oH}}_2{\text{D}}^{+}\right)n(\text{CO})+\left({\beta }_3+{\beta }_4+{\beta }_5\right)n\left({\text{oH}}_2{\text{D}}^{+}\right)n\left({\text{e}}^{-}\right),\\ \hfill \left(k_{12}+k_{13}\right)n\left({\text{pH}}_3^{+}\right)n(\text{HD})+k_{16}n\left({\text{oH}}_3^{+}\right)n(\text{HD})=\left(k_2+k_{10}\right)n\left({\text{pH}}_2{\text{D}}^{+}\right)n(\text{CO})+\left({\beta }_6+{\beta }_7+{\beta }_8\right)n\left({\text{pH}}_2{\text{D}}^{+}\right)n\left({\text{e}}^{-}\right).\end{array}$$

By summing the two equations above and exploiting the relations between the reaction rates, we arrive at: $$\left(k_{12}+k_{13}+k_{14}+k_{15}\right)n\left({\text{H}}_3^{+}\right)n(\text{HD})=\left(k_1+k_{11}\right)n\left({\text{H}}_2{\text{D}}^{+}\right)n(\text{CO})+\left({\beta }_3+{\beta }_4+{\beta }_5\right)n\left({\text{H}}_2{\text{D}}^{+}\right)n\left({\text{e}}^{-}\right),$$

which allows us to rewrite Eq. (B.1) as: $$R_{\text{D}}=\frac{1}{3}\frac{n\left({\text{H}}_2{\text{D}}^{+}\right)}{n\left({\text{H}}_3^{+}\right)}=\frac{1}{3}\frac{\left(k_{12}+k_{13}+k_{14}+k_{15}\right)n(\text{HD})}{\left(k_1+k_{11}\right)n(\text{CO})+\left({\beta }_3+{\beta }_4+{\beta }_5\right)n\left({\text{e}}^{-}\right)}.$$(B.2)

Fig. B.1 Ratio between the sum of the rates of the reactions between ${\text{pH}}_3^{+}$ and HD and those involving ${\text{oH}}_3^{+}$ and HD, as a function of temperature. Right panel: the ratio between the total rate of dissociative recombination of ${\text{pH}}_3^{+}$ and that of ${\text{oH}}_3^{+}$, as a function of temperature. In both panels, the horizontal dashed lines correspond to unity.

The next step is to express the quantity R_H = n(HCO⁺)/n(CO) as a function of the cosmic ray ionisation rate and the electronic fraction. First, we solve the kinetic equation for HCO⁺ in steady-state, again neglecting all terms containing D₂H⁺ and ${\text{D}}_3^{+}$ (see above): $$R_{\text{H}}=\frac{n\left({\text{HCO}}^{+}\right)}{n(\text{CO})}=\frac{k_{7}n\left({\text{H}}_3^{+}\right)}{{\beta }_{2}n\left({\text{e}}^{-}\right)}.$$(B.3)

To find an expression for n(H⁺), we solve its kinetics. Its total formation rate is æ₂n(H₂). The destruction pathways instead are separated in the ortho and para species, and involve the reaction with CO (reactions 5 to 7) and the following dissociative recombinations: $$\begin{array}{cc}& {\text{oH}}_3^{+}+{\text{e}}^{-}\stackrel{{\beta }_9}{\to }\text{H}+\text{H}+\text{H},\\ & {\text{oH}}_3^{+}+{\text{e}}^{-}\stackrel{{\beta }_{10}}{\to }\text{H}+{\text{oH}}_2,\\ & \begin{array}{l}{\text{pH}}_3^{+}+{\text{e}}^{-}\stackrel{{\beta }_{11}}{\to }\text{H}+\text{H}+\text{H},\\ \begin{array}{cc}& {\text{pH}}_3^{+}+{\text{e}}^{-}\stackrel{{\beta }_{12}}{\to }\text{H}+{\text{pH}}_2,\\ & {\text{pH}}_3^{+}+{\text{e}}^{-}\stackrel{{\beta }_{13}}{\to }\text{H}+{\text{oH}}_2.\end{array}\end{array}\end{array}$$

The destruction rates of the two species are then: $$\begin{array}{l}{\frac{dn\left({\text{oH}}_3^{+}\right)}{dt}|}_{\text{destr }}=-\left[k_{7}n\left({\text{oH}}_3^{+}\right)n(\text{CO})+\left({\beta }_9+{\beta }_{10}\right)\left({\text{oH}}_3^{+}\right)n\left({\text{e}}^{-}\right)\right],\\ {\frac{dn\left({\text{pH}}_3^{+}\right)}{dt}|}_{\text{destr }}=-\left[\left(k_5+k_6\right)n\left({\text{pH}}_3^{+}\right)n(\text{CO})+\left({\beta }_{11}+{\beta }_{12}+{\beta }_{13}\right)\left({\text{pH}}_3^{+}\right)n\left({\text{e}}^{-}\right)\right].\end{array}$$

From these equations, we can then compute the destruction rate for the total ${\text{H}}_3^{+}$ density, and set it equal to the total formation rate, obtaining: $${\zeta }_{2}n\left({\text{H}}_2\right)=k_{7}n\left({\text{H}}_3^{+}\right)n(\text{CO})+n\left({\text{e}}^{-}\right)\left[\left({\beta }_{11}+{\beta }_{12}+{\beta }_{13}\right)n\left({\text{pH}}_3^{+}\right)+\left({\beta }_9+{\beta }_{10}\right)\left({\text{oH}}_3^{+}\right)\right].$$(B.4)

To further simplify Eq. (B.4), we focus on the dissociative recombination rates of ortho- and para-${\text{H}}_3^{+}$. Their ratio is shown in the right panel of Fig. B.1, where one can see how at low temperatures (T ≲ 20 K), the reaction rates of ${\text{pH}}_3^{+}$ is about one order of magnitude higher than that of ${\text{oH}}_3^{+}$. We will hence neglect the second term, and introduce the para fraction $f_{\text{para }}={\text{pH}}_3^{+}/{\text{H}}_3^{+}$, to write: $${\zeta }_{2}n\left({\text{H}}_2\right)=k_{7}n\left({\text{H}}_3^{+}\right)n(\text{CO})+n\left({\text{e}}^{-}\right){\beta }_{{\text{pH}}_3^{+}}n\left({\text{H}}_3^{+}\right)f_{\text{para }},$$(B.5)

where ${\beta }_{{\text{pH}}_3^{+}}={\beta }_{11}+{\beta }_{12}+{\beta }_{13}$ Equation (B.5) can be solved for $n\left({\text{H}}_3^{+}\right)$, and then inserted in Eq. (B.3). The system of equations to infer the electron fraction and ζ₂ becomes: $$\begin{array}{cc}x\left({\text{e}}^{-}\right)& =\frac{1}{{\beta }_3+{\beta }_4+{\beta }_5}\left[\frac{1}{3}\frac{\left({\text{k}}_{12}+{\text{k}}_{13}+{\text{k}}_{14}+{\text{k}}_{15}\right)X(\text{HD})}{R_{\text{D}}}-\left({\text{k}}_1+{\text{k}}_{11}\right)X(\text{CO})\right],\hfill \\ \hfill {\zeta }_2& =\left[k_{7}X(\text{CO})+{\beta }_{{\text{pH}}_3^{+}}{f}_{\text{para}}x\left({\text{e}}^{-}\right)\right]\frac{{\beta }_2{R}_{\text{H}}x\left({\text{e}}^{-}\right)n\left({\text{H}}_2\right)}{k_7},\hfill \end{array}$$(B.6)

where now quantities are expressed in terms of abundances, rather than volume densities.

Equations (B.6) have a mathematical limitation, in that for certain combinations of R_D and X(CO) (or, equivalently, of f_D) the first one yields negative values for the electron fraction. By computing the reaction rates at 15 K, and assuming X(HD) = 1.5 × 10⁻⁵ and f_para = 0.7, we find that the threshold is R_D = 0.029 × f_D. Note a small variation to the original limitation of R_D = 0.023 × f_D. For deuteration levels higher than this limit, the equation cannot be applied. At 10 K, the new rates of Eqs. (B.6) 10 – 50% are lower than the original equations (3 and 4) of CWT98. As a result, the updated equations lead to electron fractions lower by 20% and ζ₂ values lower by 50% than the original derivation.

Figure B.2 shows the dependency of ζ₂ (normalised by the quantity R_H × n(H₂)) as a function of the deuterium fraction and CO depletion factor. One can see that, for decreasing deuteration levels, the quantity ζ₂/(R_H × n(H₂)) increases by several orders of magnitude. Since ζ₂ depends linearly on R_H and n(H₂), this translates into an equal increase also of this quantity. The plot also shows that for the deuteration values observed in dense prestellar cores (R_D = 0.01 – 0.1), there is a strong dependency on the depletion factor, up to f_D ≈ 10 – 20.

Fig. B.2 ζ2 values (divided by RH × n(H2)) obtained with Eq. (B.6) at 15 K, as a function of RD and fD. The bottom-right corner is missing because it violates the condition RD < 0.029 × fD. The plot shows that the quantity ζ2/(RH × n(H2)) increases by up to four orders of magnitude when the deuteration fraction decreases from 0.1 to 10−3 and that the CO depletion factor plays a significant role at RD typical for dense gas (0.01 – 0.1).

Appendix C T_ex values for DCO⁺ and H¹³CO⁺ in LVG analysis

In Sect. 3.5 we discussed the choice of T_ex values for each transition processed with the LVG radiative transfer. Whilst those for C¹⁸O and oH₂D⁺ are well documented in the literature, this is not the case for DCO⁺ and H¹³CO⁺ transitions. We have corroborated the values we chose by using the online tool RADEX. However, in this Appendix, we investigate how a variation of the T_ex values of ≈20% affect the inferred ζ₂ values. The results are shown in Fig. C.1, where, in the upper panel, we analyse run 5 that uses the DCO⁺ and H¹³CO⁺ (3-2) lines, and in the lower panel we present run 6, which instead focuses on the first rotational transitions.

In run 5, we explore a T_ex variation of 1K. The resulting ζ₂ values change by ≈25%. Note that the variation is stronger when the T_ex is lowered, due to the T_ex dependence of Eq. 2. For run 6, the values are changed by 2K, exploring the range 8 – 12 K, which leads to a smaller variation of the inferred ζ₂ maps (less than 10%). In all cases, the derived ζ₂ values do not change significantly compared to the uncertainties, and the median values are still less than a factor of three from the actual input ζ₂ of the simulation.

Appendix D The choice of L in single-dish-like runs

To apply the analytical equations in the case of single-dish-like post-processing of the simulations, we set L = 0.3 pc, corresponding to the size of the analysed cut-out from the original simulation box (0.6 × 0.6 × 0.6 pc³). To support this choice and to show that this is the relevant physical quantity to take into account, we perform a further additional test. We analyse the same simulation as in the reference run n. 2, but this time we halve the size of the cut-out box before the radiative transfer. Hence, L = 0.15 pc. The subsequent steps in terms of radiative transfer and post-processing are identical to those followed in run 2.

Fig. C.1 Histograms of the derived ζ2 values with the BFL20 method in run 5 (upper panel, high-J transitions) and run 6 (lower panel, low-J transitions). The different colours show different assumptions for the excitation temperature of DCO+ and H13CO+ (assumed to be equal). In particular, the reference values used in the main text are shown in red (hence these data are the same presented in Fig. 1), whilst the blue/green curves show a positive/negative variation of ≈20% of that value, respectively (labelled in the top-right corner of each panel). The median values and uncertainties are shown with the vertical dashed lines and shaded areas and with coloured text in each panel.

Our simulation consists of a box filled with molecular gas. As a consequence, most of the analysed quantities change when L = 0.15 pc. The oH₂D⁺ column density is the only one that is not affected by this change since this molecule is highly concentrated in the densest part of the core and has a high critical density. H¹³CO⁺ and C¹⁸O suffer the largest change, and their retrieved column densities decrease by a factor of up to 2 – 3. In the simulation, these two molecules are abundant everywhere, and hence cutting a smaller portion of the simulation box affects significantly their total density on the line of sight. The DCO⁺ column density decreases by 10 – 20%. The total gas density N(H₂) derived from the dust thermal emission decreases marginally (~ 5% or less). Because of these changes, R_D increases and X(CO) decreases, however not at the same rate. We compare the ζ₂ maps obtained in the two tests with distinct L computed with the BFL20 method in Fig. D.1 (left and middle panel). The distributions of values are shown as histograms in the right panel. The maps are similar both in morphology and in the range of values. The histograms confirm these conclusions. The distributions are comparable and the median values (shown with vertical, dashed lines) are consistent with each other: 〈ζ₂〉 = (4.0 ± 0.7) × 10⁻¹⁷ s⁻¹ (L = 0.3 pc) and 〈ζ₂〉 = (3.7 ± 0.8) × 10⁻¹⁷ s⁻¹ (L = 0.15 pc). If we used the old value L = 0.3 pc in the new test with a smaller box, the blue histogram would shift to the left (towards lower values) by a factor of two.

This test suggests that using the size of the cut-out box to estimate L in these runs is an appropriate choice. We stress again that this is a consequence of the particular simulation we are investigating, which represents a rather dense medium where most of the molecules of interest are abundant in the entire box. This is not the case for isolated cores such as the ones tested in Sect. 4.3 of the manuscript, for which the prescription based on the N(H₂) isocontour is appropriate.

Fig. D.1 Comparison of the resulting ζ2 maps (in units of 10−17 s−1) on the two runs with distinct box sizes L but otherwise identical, using the BFL20 method. Panel a) shows the run with L = 0.3 pc (namely the reference run 2), and panel b) shows the results with L = 0.15 pc. The colorbar is kept fixed to ease the comparison. The histogram distributions are compared in panel c), labelled in the top-right corner. The median values (median uncertainties) are shown with the vertical dashed lines (shaded areas).

References

All Tables

All Figures

Fig. 1 Histogram of the distribution of ζ2 (in units of 10−17 s−1) of the runs as labelled at the top of each panel. We represent the core region where N(H2) > 0.5 × max(N(H2)) (namely the lowest contour in the panels of Figs. 2, 3, and 4). The median value, together with the median uncertainty, is reported in the top-right corner of each panel and shown as a vertical dashed line with a grey-shaded area. Cf. Table 1.

In the text

Fig. 2 Resulting ζ2 maps obtained with Eq. (1) in runs 1 to 4. The run ID is reported in the top-left corner, and the key parameters are included at the top of each panel. Run 2 is taken as a reference throughout the rest of this work. The white contours show the 50, 70, and 90% of the N(H2) peaks, which are: 6.67 × 1022 cm−2 (core in single-dish like analysis at 100 kyr); 6.86 × 1022 cm−2 (core in single-dish-like analysis at 50 kyr); 2.51 × 1022 cm−2 (core in ALMAlike analysis at 100 kyr). The beam size and scalebar are shown in the bottom corners of each panel. Note that we show a zoom-in of the central 0.15 pc.

In the text

Fig. 3 Summary of runs performed in LVG approximation, compared with the reference case (run 2, left panel). The central panel refers to the analysis performed on the high-J transitions for DCO+, H13CO+, and C18O, whilst the low-J ones are used in the map in the right panel. The white contours are the same as in Fig. 2. Note that the colour scale is the same across all the panels. The beam size and scalebar are shown in the bottom left and right corners of each panel.

In the text

	Fig. 4 Maps of the density-averaged ${\zeta }_2^{\text{S}}$ in the simulation with variable cosmic rays (top row) at 50 kyr (left) and 100 kyr (right). Bottom row is the same as Fig. 2, but for runs 7 and 8. Note: the colour scale is the same within rows, but it changes from the top panels to the bottom ones.
In the text

Fig. 5 Maps of the key quantities employed by CWT98: DCO+ deuteration level (top-left panel), CO depletion factor (top-right), electronic fraction (bottom-left), and ζ2 (bottom-right). Note that the colour scales for the last two quantities are in units of 10−6 and 10−13 s−1, respectively. These maps assume the molecular column density and total gas column density derived in run 2. The mean uncertainties are 2.6 × 10−4 (RD), 0.3 (fD), 3 × 10−6 $\left(x_{{\text{e}}^{-}}\right)$, and 1.5 × 10−13 s−1 (ζ2). The beam size and scalebar are shown in the bottom-left panel. The contours show N(H2), as in Fig. 2.

In the text

	Fig. 6 Electronic fraction (top panel) and ζ2 map (bottom panel) from the reference run 2, using the updated formulation of CWT98. The contours show H2 column density as in Fig. 2. The mean uncertainties are 3 × 10–6 for x(e–) and 0.8 × 10–13 s–1 for ζ2. The beam size and scalebar are shown in the bottom corners of the lower panel.
In the text

Fig. 7 Same as Fig. 1, but the ionisation rate is estimated using the new formulation of the CWT98 method presented in Appendix B, summarised in Eq. (B.6). Note that the ζ2 values are normalised to 10−13 s−1 in all panels. The median (± median uncertainty) is reported in the top-right corner and shown with the vertical dashed line and shaded area in each panel. The reference run is labelled with "Ref" in the bottom-right corner.

In the text

	Fig. 8 Ratio between ζ2 retrieved in each run (x-axis labels), and the actual value ${\zeta }_2^s$. The ratios are median computed in each map, in the same regions used to produce the histograms of Fig. 1. The shaded areas show a variation of a factor of 2 (darker shade) or 3 (lighter shade). The dashed line represents where the ratio is equal to 1. The error bars are 3σ median uncertainties.
In the text