© E. Redaelli et al. 2021

1 Introduction

Cosmic rays (CRs) are ubiquitous in the interstellar medium (ISM), and play a leading role in determining its ionisation degree. In the denser regions of molecular clouds, where the column density is N ≳ 10²¹ cm⁻² (corresponding to visual extinction in magnitude A_V ≳ 1), ultraviolet (UV) photons of the interstellar radiation field are efficiently absorbed and cannot penetrate. In these physical conditions, CRs become the main ionising agent of the gas, with important consequences. First of all, they affect the chemistry by producing ions. In particular, the first reaction is the ionisation of hydrogen molecules: $${\text{CR+H}}_{\text{2}}\to {\text{H}}_{\text{2}}^{\text{+}}+{\text{e}}^{-}+\text{CR}.$$(1)

The ionised ${\text{H}}_2^{+}$ quickly reacts with another hydrogen molecule, producing the fundamental trihydrogen cation (${\text{H}}_3^{+}$). In turn, ${\text{H}}_3^{+}$ is the starting point for the chain of reactions between charged and neutral species, producing several other key molecules, such as HCO⁺, via a reaction with CO. Indirectly, CRs also regulate another important chemical process in the cold phases of the ISM: deuteration. Arguably the most important deuterated species in the gas phase is H₂D⁺, which is produced by isotopic exchange reactions between ${\text{H}}_3^{+}$ and HD (Dalgarno & Lepp 1984). Moreover, CRs produce He⁺, which can contribute to the destruction of CO, liberating C atoms in the gas phase which can then form carbon-chain molecules (see e.g. Ruffle et al. 1999), and to the formation of N⁺ from N₂, initiating the reaction chain leading to ammonia formation.

Furthermore, by determining the ionisation fraction of the ISM, CRs also affect its dynamical evolution. The ISM is threaded by magneticfields (see e.g. the results of Planck Collaboration Int. XXXV 2016 for the magnetic properties of molecular clouds), which exert their influence on charged particles directly. These transfer momentum to the neutral component through drag forces (collisions). The ionisation fraction therefore determines the degree of coupling between the gas and the magnetic field. In particular, the fraction of free electrons determines the timescale of ambipolar diffusion, that is to say the process of drift of neutral matter across the field lines which is the proposed mechanism that allows the gravitational collapse of magnetised cores (Mouschovias & Spitzer 1976).

Despite its importance, observational estimates of the CR ionisation rate are still scarce, uncertain, and scattered over a wide range of values. The most direct and reliable way to determine the ionisation rate is represented by direct observations of ${\text{H}}_3^{+}$. These, however, are possible only in diffuse gas via infrared spectroscopy. Lacking a permanent electric dipole, in fact, ${\text{H}}_3^{+}$ is not observable in emission in the interstellar gas at high densities. In dense clouds, van der Tak & van Dishoeck (2000) used ${\text{H}}_3^{+}$ absorption lines in combination with H¹³CO⁺ rotational lines towards background massive protostars, and they determined a CR ionisation rate of ≈ 3 × 10⁻¹⁷ s⁻¹. More recently, Indriolo & McCall (2012) – expanding the sample from Indriolo et al. (2007) – reported 21 detections towards diffuse clouds, determining CR ionisation rate per hydrogen molecule (ζ₂) in the range [1.7−10.6] × 10⁻¹⁶ s⁻¹. Neufeld & Wolfire (2017) derived a CR ionisation rate per hydrogen atom (ζ₁) of the order of a few 10⁻¹⁶ s⁻¹ in the Galactic disc, using ${\text{H}}_3^{+}$ and other ionised tracers. Their data support a steep dependency of ζ₂ on the gas column density (ζ₂ ∝ N⁻¹).

In regions so dense that no background source is visible at the near infrared wavelengths of the ${\text{H}}_3^{+}$ lines, no ${\text{H}}_3^{+}$ absorption can be observed. In these conditions, a commonly used approach is to compare the abundances of molecular tracers sensitive to the ionisation fraction (such as HCO⁺) with the results of chemical models (see e.g. Black et al. 1978). This method was later used for instance by Caselli et al. (1998) towards several pre-stellar and proto-stellar cores, using observations of DCO⁺, HCO⁺, and CO. The results span the range 10⁻¹⁸−10⁻¹⁶ s⁻¹. Bovino et al. (2020) proposed a new analytical method to assess ζ₂ using observations of CO, of the deuterium fraction (for instance from HCO⁺ and DCO⁺), and of ortho-H₂D⁺. The method was applied by Sabatini et al. (2020) to a sample of massive clumps, finding ζ₂ ≈ [0.7−6] × 10⁻¹⁷ s⁻¹.

From the theoretical point of view, one of the first studies on CR ionisation was performed by Spitzer & Tomasko (1968), who computed the CR ionisation rate of atomic hydrogen, yielding ζ₁ ≈ 6.8 × 10⁻¹⁸ s⁻¹. This can be transformed into ζ₂ using the relation ζ₂ = (2.3∕1.5) × ζ₁ (Glassgold & Langer 1974), obtaining the value ζ₂ ≈ 10⁻¹⁷ s⁻¹ which is still considered the standard value for molecular clouds. However, the assumption that ζ₂ is independenton the ISM density is incorrect, as also shown by the different observational results. Cosmic rays lose energy in their interactions with the ambient medium, which modifies their spectrum. These losses hence cause an attenuation of ζ₂, which decreases with gas column density (N). This processhas been studied in detail by Padovani et al. (2009), who computed the relation between the CR ionisation rate and the gas column density using: $${\zeta }_2(N)=4\pi [1+\Phi (N)]{\sum }_k{\displaystyle {\int }_{E_{\text{ion}}}^{\infty }{j}_k}(E,N){\sigma }_k^{\text{ion}}(E)\text{d}E,$$(2)

where k indicates the primary CR species (e.g. protons, electrons, positrons, and heavier nuclei), E_ion = 15.44 eV is the energy threshold for H₂ ionisation, ${\sigma }_k^{\text{ion}}(E)$ is the ionisation cross section of H₂ by the species k (Rudd et al. 1992; Kim et al. 2000), and j_k(E, N) is its local energy spectrum, which depends on N due to energy losses. Secondary CR electrons, produced due to the gas ionisation by the primary CR species, lead to additional H₂ ionisation. Their contribution is described by Φ(N), which is the ratio between the secondary to primary H₂ ionisation rates (see Ivlev et al. 2021). The local CR spectrum is commonly calculated using the continuous slowing-down approximation (see e.g. Padovani et al. 2009), which assumes that CRs propagate freely along the local magnetic field lines. As was shown by Silsbee et al. (2018), the results are generally applicable to arbitrary magnetic field configurations, including converging field lines: as long as the field strength increases monotonically along the field line, the effects of magnetic mirroring and focussing practically cancel each other. The column density N in this case is measured along the field lines.

Padovani et al. (2018) (hereafter PI18) further expanded the work of Padovani et al. (2009) considering the regime of higher column densities, and using the most recent data from the Voyager missions to constrain the CR interstellar flux (Cummings et al. 2016; Stone et al. 2019). Similarly to Ivlev et al. (2015), PI18 considered two models for the interstellar flux of protons: the first, called the ‘high’ model, is characterised by ζ₂ ≳ 10⁻¹⁶ s⁻¹ up to N ≈ a few × 10²³ cm⁻², and reproduces the high ζ₂ measured in the diffuse medium. The ‘low’ model, which instead implies ζ₂ ≈ a few 10⁻¹⁷ s⁻¹, is the one thatfits the Voyager 1 and 2 data. It is important to notice that the spectrum derived by the Voyager missions could represent a lower limit, due to (i) persisting contamination from CRs from inside the heliosphere (Scherer et al. 2008), and to (ii) the possible influence of local sources, which make the observed spectrum deviate from the average Galactic one (Phan et al. 2021).

In this work, we investigate the CR ionisation rate in the prototypical pre-stellar core L1544, combining CR propagation models fromPI18 with a state-of-the-art chemical model, and comparing the simulated spectra of several rotational lines of molecular ions with recent observations. L1544 is situated in the Taurus molecular cloud at a distance of ≈ 135 pc (Schlafly et al. 2014), and it shows bright continuum emission at millimetre and sub-mllimetre wavelengths (Ward-Thompson et al. 1999; Doty et al. 2004; Spezzano et al. 2016). Its physical structure has been modelled by means of a quasi-equilibrium Bonnor-Ebert sphere using N₂H⁺ and carbon monoxide spectra (e.g. Keto et al. 2015). This model predicts low temperatures (T ≈6 K) and high volume densities (n ≳ 10⁶ cm⁻³) within the central 2000 AU, in agreement with observations (Crapsi et al. 2007; Chacón-Tanarro et al. 2017; Caselli et al. 2019).

In a recent work (Redaelli et al. 2019, hearafter RB19), we combined this physical model with a gas-grain chemical model to investigate the deuteration maps of HCO⁺ and N₂H⁺ (diazenylium) across the source. We now perform a similar analysis, focussing this time on ζ₂. In Sect. 2 we present the chemical and physical models used, and the resulting abundance profiles are discussed in Sect. 3. The analysis of the radiative transfer results for the different cases are reported in Sect. 4. Section 5 summarises the main findings and conclusions of this work.

2 Model description and observations

2.1 Physical model

The physical model of L1544 consists of a one-dimensional, quasi-equilibrium Bonnor-Ebert sphere, developed by Keto et al. (2015) fitting C¹⁸O, H₂O, and N₂H⁺ lines. It comprises the gas temperature, dust temperature, volume density, and infall velocity profile, and it has been successfully used to model several molecules at the dust peak of L1544 (Bizzocchi et al. 2013; Redaelli et al. 2018, 2019). In the latest paper, we modelled the N₂H⁺ and HCO⁺ isotopologues using a modified infall velocity profile (scaled up by a factor 1.75), which was originally introduced by Bizzocchi et al. (2013) to reproduce the N₂H⁺ (1–0) transition with a constant abundance profile. We realised that this enhancement of the infall velocity leads to incorrect spectral profiles (e.g. strong double-peak profiles in N₂D⁺), which are inconsistent with the observations. For this work we therefore used the original velocity profile of Keto et al. (2015). The volume density profile and the corresponding integrated column density profile are shown in the left and central panels of Fig. 1.

Fig. 1 Left panel: H2 volume density profile (n) in thephysical model of L1544. Central panel: column density profile. Right panel: ζ2 profiles for the $\mathcal{F}$ (dashed bluecurve), $\mathcal{L}$ (solid red), and $\mathcal{H}$ (dashed-dotted green) models. The profiles for the models $\mathcal{L}$ and $\mathcal{H}$ PI18 have been derived integrating the volume density profile of the core model to obtain the column density profile, and then using the relation between ζ2 and N.

2.2 Chemical models

The physical model of L1544 has been coupled to the gas-grain chemical model presented in Sipilä et al. (2015a,b, 2019). The core model is divided into concentric shells, and the evolution of the chemistry is then solved in each shell separately. The combination of the result from all the shells at a common time step gives the radial profile of the molecular abundance. The radial grid consists of 35 points covering the whole core size, assumed to be 0.32 pc (see Fig. 1), with a resolution that ranges from 3.5 × 10⁻³ pc in the central parts to 2 × 10⁻² pc in the outskirts of the core model. The initial conditions are summarised in Table 1. The dust population consists of grains of 0.1 μm of radius and a grain material density of ρ_g = 2.5 g cm⁻³. The effect of the embedding cloud was simulated by assuming an external visual extinction of A_V = 2 mag¹, which attenuates the external UV field. We tested four different time steps: t = 10⁵, 5 × 10⁵, 10⁶, and 5 × 10⁶ yr.

The chemical code takes as input the CR ionisation-rate ζ₂, which can be either constant, or can have a radial dependence. For the purposes of this paper, we have tested three models: the first model, called the fiducial model ($\mathcal{F}$), has a constant value ζ₂ = 1.3 × 10⁻¹⁷ s⁻¹; the second one reproduces the high model ($\mathcal{H}$) from PI18, whilst the third one corresponds to the low model ($\mathcal{L}$) from PI18. The first model is similar to that used in RB19. In order to implement models $\mathcal{L}$ and $\mathcal{H}$, we have used the parametrisation presented in Appendix F of PI18 to compute the ζ₂ profile from the column density profile. The radial profiles of the three ζ₂ models are presented in the right panel of Fig. 1. In models $\mathcal{L}$ and $\mathcal{H}$, the ζ₂ profile is computed assuming a column density value corresponding to A_V = 2 at the edge of the core. The three models are well separated in the range of ζ₂ values that they span, and can be ordered by increasing ζ₂ as $\mathcal{F}$, $\mathcal{L}$, and then $\mathcal{H}$.

In our approach, we assume for simplicity that CRs propagate radially along the field lines in a spherically symmetric geometry, but we make no assumption on the field geometry (except the one that the field strength does not have multiple minima along the field lines, see Sect. 1). Accordingly, we compute the column density for CR attenuation by integrating radially the density from the cloud’s boundary to the point of interest. This is the minimum possible column density that characterises the CRs attenuation at any interior point of the cloud. In general, the attenuation of CRs propagating along magnetic field lines will be larger than this minimum value, by a factor of the order of a few times that depends on the actual field geometry (see Padovani et al. 2013 for an accurate calculation of the column density in a realistic cloud’s magnetic field and density profile). We note that at small scales (< 10⁻³ pc = 200 AU), where the field lines can become entangled and present complex morphology with local minima of the field strength (see again Padovani et al. 2013), the magnetic mirroring and focussing may no longer balance each other, and then calculating the local ionisation rate becomes a nontrivial problem (Silsbee et al. 2018). However, here we investigate processes occurring at substantially larger physical scales (>10⁻² pc = 2000 AU), and therefore we can safely employ the results of PI18.

We stress here that the models $\mathcal{F}$, $\mathcal{L}$, and $\mathcal{H}$ differ in the shape of theζ₂ profile, whilst the rest of the physical structure (in terms of velocity, temperature, and density) remains the same. This approach is not self-consistent, since the ionisation rate affects the temperature of the gas. A complete physical model should be computed self-consistently, which is however beyond the scope of the present paper. We do not expect this effect to be significant in the dense gas traced by the targeted molecular transitions, as the gas at these densities (n ≳ 10⁵ cm⁻³) is mainly cooled by dust. This cooling is efficient enough that we do not expect important changes (see also Ivlev et al. 2019). We however discuss this point in more detail in Appendix A. Our chemical network does not include the chemistry of oxygen isotopes. The abundance profiles for HC¹⁸O⁺, which is the observed isotopologue, are derived diving those of HCO⁺ by the isotopic ratio ¹⁶O∕¹⁸O = 557 (Wilson 1999).

Fig. 2 Top panels: abundance profiles for the four analysed species (from left to right: N2H+, N2D+, HC18O+, DCO+) at t = 106 yr. The colours refer to different models used for the CR ionisation rate: model $\mathcal{F}$ (blue), model $\mathcal{L}$ (red), and model $\mathcal{H}$ (green). Bottom panels: the green and red solid curves show the ratio of model $\mathcal{H}$ and model $\mathcal{L}$ to model $\mathcal{F}$, respectively. The dashed, horizontal lines represent the inverse ratios of the corrective factors fcorr used to reproduce the observed lines, as reported in Sect. 4.1 and Table 2. To give an example, the radiative transfer of DCO+ lines requires fcorr = 0.9 in model $\mathcal{L}$, and fcorr = 1.5 in model $\mathcal{F}$; therefore thered solid curve $\mathcal{L}$ / $\mathcal{F}$ in the rightmost panel has to be compared with the red dashed line at ${(0.9/1.5)}^{-1}=1.67$.

2.3 Observations

The abundance profiles derived from our chemical network can be used to obtain synthetic spectra to be compared with observations. For this, we have used the data presented in RB19. Specifically, we have used single-dish data obtained with the Institut de Radioastronomie Millimétrique (IRAM) 30 m telescope of: N₂H⁺ (1–0) and (3–2); N₂D⁺ (1–0), (2–1), and (3–2); DCO⁺ (1–0), (2–1), and (3–2); HC¹⁸O⁺ (1–0). The latest is preferred to the optically thick main isotopologue. All the spectra were obtained as single pointings taken at the millimetre dust peak of L1544 (RA(J2000) = 05^h04′17″.21, Dec(J2000) = 25°10′42″.8). We refer to RB19 for a complete description of the observational data. We remark however the unique features of this dataset, which make them particularly suitable for the purposes of our work. First of all, with the exception of the DCO⁺ (1–0) transition, all the data have been acquired with good spectral resolution (ΔV = 0.02−0.08 km s⁻¹), crucial to correctly recover the complex line profiles arising from a combination of molecular depletion, self-absorption, and kinematic effects. Furthermore, the different rotational transitions of a given species have distinct critical densities, and hence trace different depths into the core. For instance, at 10 K the critical densities of the N₂H⁺ (1–0) and (3–2) transitions are n_c = 6 × 10⁴ cm⁻³ and n_c = 1.4 × 10⁶ cm⁻³, respectively (Shirley 2015). This fact, combined with the improved angular resolution at higher frequencies, allows to test the radial profiles of these molecules reliably, even though the spatial resolution of the single-dish data is limited for the low-J transitions.

3 Results

Figure 2 shows the simulated abundance profiles for the different species of interest at t = 10⁶ yr. We show in Appendix B the plots for the remaining time steps. From the comparison of the models, we can derive some general conclusions. As the evolutionary time increases, N₂H⁺ isotopologues show stronger and stronger depletion towards the core’s centre. The abundance peak increases, but shifts outwards, whilst the abundance at smaller radii decreases. This effect is visible also for DCO⁺ and HC¹⁸O⁺, but it is less evident because these molecules are already significantly depleted early on (t ≲10⁵ yr) at radii smaller than0.04 pc.

The effect of varying ζ₂ is species-dependent. HCO⁺ (and as a consequence HC¹⁸O⁺) presents a clear, monotonic effect at every time step: the higher the value of ζ₂, the higher the overall abundance profile, due to the increase of the ionisation fraction. This trend is partially visible also for DCO⁺, but the differences in the abundance profiles of the various models are smaller. The behaviour of N₂H⁺ isotopologues is more complex, and time-dependent. In general, we note that the higher the overall ζ₂, the stronger is the depletion, and the earlier this phenomenon starts showing up. Model $\mathcal{H}$ predicts a strong drop in the abundance of N₂H⁺ and N₂D⁺ starting at t = 10⁶ yr. At t = 5 × 10⁶ yr, also model $\mathcal{L}$ (which has an intermediate ζ₂ value between models $\mathcal{F}$ and $\mathcal{H}$) predicts an enhanced decrease of the abundances. For N₂D⁺ in particular, the high model tends to always predict the lowest abundances with respect to the other two models.

4 Analysis and discussion

4.1 Non-LTE modelling at t = 10⁶ yr

In order to produce the synthetic spectra to be compared with the observations, we used the non-local-thermodynamic-equilibrium (non-LTE) radiative transfer code MOLLIE (Keto 1990; Keto et al. 2004). MOLLIE requires several inputs: the source’s physical structure; the molecular abundance profile, obtained from the chemical models, which can be modified by a corrective factor f_corr ²; the collisional coefficients of the molecular transitions (see RB19 for references); the turbulent line broadening that is summed in quadrature to the thermal one. For this we adopted σ_turb = 0.075 km s⁻¹, following RB19. The code, which is particularly suited to model crowded hyperfine structures prone to selective photon trapping, requires long computational times to converge (of the order of hours using 20 cores for N₂H⁺ and N₂D⁺), and a full parameter-space exploration is not feasible, in particular to determine f_corr. As in RB19, we therefore made several tests to determine by trial-and-error the f_corr value that matches the observed line intensities for each species.

Before performing the ray-tracing to compute the synthetic spectra, MOLLIE interpolates the 1D physical model and abundance profile to a 3D array (mantaining the spherical symmetry), which consists of three nested grids that become finer going towards the core centre. Going from the largest to the smallest grid, each of them covers half of the radial range with respect to the previous one, doubling the spatial resolution. The finest grid has a cell size of 1.7 × 10⁻³ pc, which corresponds to ≈2.6″ at the distance of L1544. We conclude that the models produced by MOLLIE have much higher resolution than the observations. In order to allow for a proper comparison with the observed data, we then convolved the results from MOLLIE to the corresponding beam size for each transition, and we extracted the spectra at the centre of the models.

The results of the chemical simulations tend to over- or under-estimate the observed fluxes, in particular when one tries to model several transitions of different species at the same time. We base the comparison among different models mainly on three parameters:

• 1.

  the corrective f_corr required to reproduce the observed fluxes; the closer these are to unity, the better the model performs;

• 2.

  the agreement with the line profiles (in term of asymmetries, double-peak features, self-absorption, and intensity ratios of the hyperfine components, if present);

• 3.

  the correct prediction of the intensity ratios of different rotational transitions for the same species. This is particularly important, because due to the different excitation conditions of the distinct rotational lines, their intensity ratios probe the level of depletion of the species at high densities.

In a first test, we focussed on the evolutionary time of t = 10⁶ yr, which in RB19 provided the best fit for all the molecules of interest. We hence ran MOLLIE with the abundance profiles at this time step, performing a few (2−5) tests to find the f_corr value that provides a good overall match to the observed fluxes and line shapes for all the observed transitions of each molecule, based in particularon the parameters listed above at points 2 and 3. In Appendix C we validate this approach to determine f_corr through the comparison with a more rigorous, quantitative approach, showing that the two methods produce consistent results within the uncertainties. The best-fit results at t = 10⁶ yr are shown in Figs. 3–6, and they are discussed in detail in the following.

N₂H⁺

None of the models at t = 10⁶ yr is able to reproduce the intensity ratio of the two observed rotational transitions for N₂H⁺. Models $\mathcal{F}$ and $\mathcal{L}$ tend to overestimate the flux of the (3–2) transition with respect to that of the (1–0), whilst the opposite is true for model $\mathcal{H}$. This last model, furthermore, requires an enhancement of 400% of the overall abundance profile in order to match the observed line intensities. The abundance obtained with model $\mathcal{F}$ needs to be increased by f_corr = 1.5, whilst model $\mathcal{L}$ predicts an abundance profile that agrees reasonably well with the observations, despite underestimating by 20% the intensity of the main component of the (1–0) line, and overestimating by 9% that of the (3–2) line.

N₂D⁺

All models tend to underestimate the fluxes of the three N₂D⁺ transitions, similarly to what was found in RB19. Model $\mathcal{H}$ is characterised by the strongest underestimation, and the N₂D⁺ abundance must be enhanced of a factor of 20 in order to match the observed (1–0) line. Still the higher J-transitions appear underestimated. Model $\mathcal{L}$ requires the smallest f_corr factor, and the synthetic line profiles are in good agreement with the observed ones, with the exception of the (1–0) transition, which is underestimated by ≈15%.

HC¹⁸O⁺

Also for this molecule, model $\mathcal{L}$ provides thebest agreement in terms of flux. None of the models however is able to reproduce the symmetric double-peakfeatures seen in the spectrum, which may depend on the kinematic structure of the source (Ferrer et al., in prep.).

DCO⁺

This molecule is the one that requires the smallest f_corr factors among the targeted species, and is the one for which the intensity ratios of the different J-transitions are better reproduced. Model $\mathcal{H}$ has to be reduced by a factor of ≈3 to reproduce the observed fluxes, whereas model $\mathcal{F}$ has to be enhanced by f_corr = 1.5. For model $\mathcal{L}$, a modificationof 10% is enough to obtain a good agreement with the observations.

Based on these results, we conclude first of all that the model $\mathcal{H}$, with the highest ζ₂ profile, is excluded by our data. In particular, it severely underestimates the abundance of N₂D⁺, and to a lesser extent also that of N₂H⁺, whilst it overestimates the abundances of HCO⁺ isotopologues. Models $\mathcal{F}$ and $\mathcal{L}$ are instead in better agreement with the observations The latter appears better than model $\mathcal{F}$, since it requires the smallest f_corr corrective factors for three out of the four targeted species. It is in fact remarkable that a single chemical model is able to reproduce simultaneously several transitions of N₂H⁺, HC¹⁸O⁺, and DCO⁺ (with a correction of only 10%), at the same time underestimating the fluxes of the N₂D⁺ lines by the smallest factor.

The bottom panels of Fig. 2 show the ratio between the models ($\mathcal{L}$ / $\mathcal{F}$ and $\mathcal{H}$ / $\mathcal{F}$), compared with the inverse ratio of the corresponding f_corr values used to obtain the best-fit solutions. Qualitatively, the range of radii where the model ratios are similar to the f_corr ratios gives us an idea of which parts of the source contribute the most to the observed transitions. For most of the combinations of models and species, it can be seen that this agreement is found for radii smaller than a few × 10⁻² pc, corresponding to volume densities higher than ≈10⁵ cm⁻³ (see left panel of Fig. 1). This is consistent with the high critical densities of the analysed transitions, which are in the range from a few 10⁴to ≈10⁶ cm⁻³.

Fig. 3 Spectra obtained with the radiative transfer code (model, in red) overlaid to the observations (black) for N2 H+ (1–0) and (3–2) lines. The columns refer to the three models: $\mathcal{F}$ in the left panels, $\mathcal{L}$ in the central panels, and $\mathcal{H}$ in the right ones. The corrective factors fcorr used to modify the whole abundance profile to match the observed fluxes is indicated above the top panel of each column. All models are at t = 106 yr.

4.2 Testing different evolutionary stages

The previous subsection has shown that model $\mathcal{L}$ provides the best fit to the observations at the evolutionary time of one million years. We can however test whether different times produce better agreements. To do so, we have run the non-LTE radiative transfer on model $\mathcal{L}$ at three further time steps: t = 10⁵ yr, t = 5 × 10⁵ yr, and t = 5 × 10⁶ yr. These abundance profiles are shown in Appendix B (together with the ones for models $\mathcal{F}$ and $\mathcal{H}$). As previously described, we tested a few f_corr values in order to match the observed fluxes.

There are two main ways in which the changes in abundance as a function of time affect the simulated line profiles. First of all, the chemical evolution will lead to an overall increase or decrease of the total molecular column density, which modifies the predicted fluxes. Secondly, as time goes on, depletion (due to freeze out onto the dust grains, either of the targeted molecule or of its precursors) becomes more significant. This changes the intensity ratios of the different transitions, since the higher-J transition (having higher n_c) will be less excited with respect to the (1–0) lines. This effect was already observed in RB19, and it stresses again how crucial it is to have access to distinct transitions with different critical density of the same molecular tracer to constrain the molecular abundance profile.

Concerning N₂H⁺, no time step significantly improves the agreement with the observations. In particular, the model at t = 10⁵ yr does not match the intensity ratios of the (3–2) to the (1–0) line and of the different hyperfine components in the (1–0) transition. The model at t = 5 × 10⁶ yr underestimates the observed peak intensity by more than 50% for both transitions. The models at t = 5 × 10⁵ yr and t = 10⁶ yr are similar in terms of the ratio of the two rotational lines intensities, with the model at t = 10⁶ yr providing a slightly better agreement.

The lines of N₂D⁺, on the other hand, are better reproduced when adopting an earlier evolutionary time (t =5 × 10⁵ yr), and f_corr = 3. In order to better understand the evolution of the N₂D⁺ profile, we show in Fig. 7 the abundance profiles at four time steps for model $\mathcal{L}$. As the chemistry evolves, the peak increases in value and it shifts outwards, whilst at the centre the abundance decreases. In Sect. 4.1we showed how at t = 10⁶ yr the line intensity ratios are well reproduced, but the abundance profile has to be increased by a factor f_corr = 4. The earliest and the latest time steps perform worse because they predict too little or too strong depletion, respectively. In the former case (t = 10⁵ yr) the fluxes of the (2–1) and (3–2) transitions are overestimated with respect to the (1–0) line, and the other way around for the latter case (t = 5 × 10⁶ yr). The intermediate time step t = 5 × 10⁵ yr presents instead the best combination of total molecular density and depletion level.

For DCO⁺ and HC¹⁸O⁺, the model at t = 10⁶ yr requires the smallest corrective factors to reproduce the observed lines, and it is hence considered the best-fit model. We highlight that using distinct time steps to reproduce different molecular tracers does not represent a contradiction, due to the fact that our physical model is static and the chemistry is evolved in post-processing. The physical age of the source is linked to the dynamical timescale, and not to the chemical one. Moreover, distinct molecules are sensitive to different physical conditions. Species such as N₂D⁺, which is formed only when CO is depleted from the gas phase and temperatures are low, trace mainly the very dense, inner parts of the source, where the dynamical timescale is shorter (since it depends inversely on square root of the volume density). This could explain why N₂D⁺ is better reproduced at an earlier chemical timescale with respect to the other species.

A better performance of the two remaining ζ₂ models can be excluded based on the following arguments. The main problems of model $\mathcal{H}$ are the severe underestimation of N₂D⁺ abundance (by more than one order of magnitude) and a large overestimation of HC¹⁸O⁺ abundance (by a factor of 4). At any time step, model $\mathcal{H}$ presents the highest abundance of HCO⁺ (and thus of HC¹⁸O⁺), which will always be overestimated. Concerning N₂D⁺, model $\mathcal{H}$ always predicts the lowest abundance, with the exception of the earliest evolutionary time (t =10⁵ yr, see Fig. B.1). This time step, however, presents almost no signs of molecular depletion, and we have shown how this kind of models is not able to reproduce the intensity ratios of the distinct J transitions.

The main problems of model $\mathcal{F}$ at t = 10⁶ yr are the underestimation of HC¹⁸O⁺ (by a factor of 2.6) and of N₂D⁺ (by a factor of 5). The former problem will not be solved at any time step, since the HCO⁺ abundance is always the lowest in the fiducial model with respect to the other models. Concerning N₂D⁺, we tested the same time-step that provides the best fit for the model $\mathcal{L}$: t = 5 × 10⁵ yr and f_corr = 3. We find that with t his abundance profile the (3–2) line is reproduced, but the (2–1) and the (1–0) transition are overestimated by ≈ 30%. We conclude that changing the evolutionary stage does not provide better agreement with the observations for both model $\mathcal{H}$ and $\mathcal{F}$, with respect to model $\mathcal{L}$.

Fig. 4 As Fig. 3, but for N2D+.

Fig. 5 As Fig. 3, but for HC18O+.

Fig. 6 As Fig. 3, but for DCO+.

Fig. 7 Abundance profiles for N2D+ in model $\mathcal{L}$ at four different time steps, from t = 105 yr to t = 5 × 106 yr. One can see that as the chemistry evolves, N2D+ is affected by stronger depletion, with a higher contrast between the central abundance and the peak value.

4.3 The electron fraction

Since CRs are the main ionising source in dense gas, they directly determine the abundance of free electrons, or electron fraction (x(e) = n(e)∕n). As explained in Sect. 1, this parameter regulates the degree of coupling between the gas and the magnetic fields, and it is linked to the timescale for ambipolar diffusion. The results of the chemical modelling include also this parameter, which is shown in Fig. 8 for the three models $\mathcal{F}$, $\mathcal{L}$, and $\mathcal{H}$. We focus on the model results at t = 10⁶ yr, but we highlight that, in the inner parts of the core (r < 0.05 pc), x(e) values change only by a few percent when varying the time step.

Figure 8 shows that the electron fraction increases with increasing ζ₂, as expected due to the nature of the ionisations. At the core centre, the value obtained in model $\mathcal{L}$ is x(e)= 10⁻⁹. This result is in agreement with the one from Caselli et al. (2002), who used a more simplified chemical model on a subset of the molecular transitions presented in this work to estimate this quantity.

Fig. 8 x(e) profiles in model $\mathcal{F}$ (dashed bluecurve), $\mathcal{L}$ (solid red), and $\mathcal{H}$ (dashed-dotted green). All the models are shown at t = 106 yr.

4.4 Simulating the cloud evolution

As a further test on the fiducial model with constant ζ₂ = 1.3 × 10⁻¹⁷ s⁻¹, we changed the initial conditions of the chemical model. We simulated first a cloud-like evolution for t_CE = 10⁶ yr. For this initial simulation, we used the following physical conditions: n = 10³ cm⁻³, T =15 K, and external extinction A_V = 2 mag. The final abundances of this model were then taken as initial conditions for the simulation of the core. The time step values t refer to the time elapsed after the end of the initial cloud simulation (i.e. we reset t = 0 at the time t_CE).

In Fig. 9 we show the final abundance profiles for the species of interest from the two models, with and without cloud evolution. N₂H⁺ and N₂D⁺ show the strongest variations. The model with cloud evolution presents stronger depletion towards the centre at any time step. For HCO⁺ isotopologues, the differences are smaller, especially for HCO⁺. In the case of DCO⁺, the molecular abundance is generally smaller in the model with cloud evolution, with the exception of the earliest time step (t =10⁵ yr).

As a result of these differences, the cloud-evolution model requires in general higher corrective factors f_corr to reproduce the observed fluxes. We report in Appendix D the complete set of comparisons between the synthetic and observed spectra at t = 10⁶ yr. We discuss here the case of N₂D⁺, which, as shown also in the previous subsections, is sensitive to the level of depletion due to the high critical densities of the (2–1) and (3–2) rotational lines. We report the results of MOLLIE for this model, multiplied by f_corr = 5, in Fig. 10, compared with the best-fit solution found for model $\mathcal{F}$. Whilst the first rotational transition is well reproduced both in its line profile and in the intensity ratios of the distinct hyperfine components, the (2–1) and (3–2) transitions are underestimated by 25 and 45% of the peak intensity, respectively. These results for N₂D⁺, together with the fact that the corrective factors f_corr are higher also for the other species (see Appendix D), suggest that the initial conditions obtained with the initial cloud evolution model are not consistent with our observations.

We point out that there are significant uncertainties involved in modelling an initial cloud-like evolution, associated for instance with choosing the appropriate physical conditions for the cloud stage or selecting the value of t_CE. Furthermore, neglecting the dynamical evolution from cloud to core is also a source of uncertainty.

4.5 Constant ζ₂ = 3 × 10⁻¹⁷ s⁻¹ model

The models so far explored are characterised by average ζ₂ values that are significantly different: ζ₂ = const = 1.3 × 10⁻¹⁷ s⁻¹ in model $\mathcal{F}$, ⟨ ζ₂ ⟩ = 3 × 10⁻¹⁷ s⁻¹ in model $\mathcal{L}$, and ⟨ ζ₂ ⟩ = 2 × 10⁻¹⁶ s⁻¹ in model $\mathcal{H}$. The last two values refer to gas column densities in the range N = [0.2−5] × 10²² cm⁻² (see Fig. 1), obtained by integrating the volume density profile from Keto et al. (2015). It is hence important to verify whether our observations are truly sensitive to the CR attenuation towards the core centre, or if they can only distinguish among different average values of ζ₂.

In order to test this, we ran another chemical model, identical to the fiducial one but adopting a constant CR ionisation rate equal to ζ₂ = 3 × 10⁻¹⁷ s⁻¹ (the average value of model $\mathcal{L}$). The resulting abundance profiles are shown in Fig. 11, in comparison to those of model $\mathcal{L}$. The two models produce very similar profiles, in particular at the two time steps that provide the best-fit (t =5 × 10⁵ yr and t = 10⁶ yr). The median relative differences at these time steps are 6−9% for N₂H⁺, 8− 10% for N₂D⁺, 16− 17% for HCO⁺, and 10−13% for DCO⁺. As a consequence, also the synthetic line profiles, shown in Appendix E, appear similar.

Given the uncertainties of this work (e.g. the assumed physical model, or the uncertainties relative to the chemical network, such as the rate coefficients), we cannot make stringent conclusions on these two models. In particular, we lack the spatial resolution to conclusively assess if the molecular line data support the model of CR attenuation in L1544. We highlight that, in fact, the angular resolutions of the low-J transitions are in the 27−36″ range, which correspond to 0.018−0.024 pc (3600−4900 AU) at the distance of L1544.

Neufeld & Wolfire (2017) used observations of several ions (and in particular of ${\text{H}}_3^{+}$) to derive the CR ionisation rate in a sample of diffuse clouds. Their data supported a stronger attenuation of the CRs with increasing density with respect to the models of PI18, with a power-law dependence following ζ₂ ∝ N^−0.92 (see their Fig. 6). The reader may note that in the range of column densities spanned by the physical model, the model $\mathcal{L}$ predicts a shallow dependency on the column density, approximately following ζ₂ ∝ N^−0.2. We tested the steeper relation derived from Neufeld & Wolfire (2017) in our chemical model³. The resulting abundances underestimate the observed transitions with respect to model $\mathcal{L}$. The reason for this behaviour is that the relation found by Neufeld & Wolfire (2017), when applied to the range of densities of L1544, implies that for radii smaller than 0.01 pc the CR ionisation rate is ζ₂ < 3 × 10⁻¹⁷ s⁻¹. As a consequence, in the central regions – where the targeted transitions are excited, see Sect. 4.1 – the lower CR ionisation rate leads in general to lower predicted abundances, in particular for HC¹⁸O⁺ and DCO⁺, for which the dependence on ζ₂ is monotonic, as discussed in Sect. 3.

Fig. 9 Top panels: abundance profiles obtained with the model with constant ζ2 = 1.3 × 10−17 s−1, with no initial cloud evolution (model $\mathcal{F}$, solid curves) and with a cloud evolution of tCE = 106 yr to set the initial conditions (dashed curves). The different colours refer to the four time steps investigated in this work (labelled in the first panel, in year). Bottom panels: ratios between the abundance profiles of the model with and without cloud evolution. The colour code is the same as for the upper panels.

Fig. 10 Best-fit model for N2D+ lines obtained with the fiducial model with initial cloud evolution and ζ2 = 1.3 × 10−17 s−1, at t = 106 yr (red curves),overlaid to the observations (black). The used fcorr value is indicated at the top. The modelled lines produced using model $\mathcal{F}$ with no cloud evolution at t = 106 yr (and fcorr = 5) are shown in green. It is clear that the model with cloud evolution underestimates the (2–1) and (3–2) lines.

5 Summary and conclusions

In this work we studied the effect of the CR ionisation rate on the results of gas-grain chemical simulations, focussing on four species: N₂H⁺, N₂D⁺, HC¹⁸O⁺, and DCO⁺. The predicted abundance profiles have been used to perform a full non-LTE radiative transfer analysis, and the synthetic spectra were compared with recent, high-sensitivity spectra of the molecules at the dust peak of the prototypical pre-stellar core L1544.

The first important conclusion is that the model with high ζ₂ (corresponding to the high model of PI18), characterised by an average value of ⟨ ζ₂ ⟩ = 2 × 10⁻¹⁶ s⁻¹, is excluded by our observations. In particular, this model strongly underestimates the N₂D⁺ abundance, while it overestimates by a factor of 3−4 the intensities of HC¹⁸O⁺ and DCO⁺ lines.

The model $\mathcal{L}$ from PI18 is the one that provides the overall best agreement with the observations. It is remarkable that, with the exception of N₂D⁺ – which is underestimated by all the models here considered – the other targeted molecules are reproduced within 10%. The fiducial model, with standard ζ₂ = 1.3 × 10⁻¹⁷ s⁻¹, performs worse than model $\mathcal{L}$, and requires corrections in the range f_corr = 1.5−5 at t = 10⁶ yr to reproduce the observed rotational lines. Changing the initial conditions by simulating first a cloud-like evolution of the chemistry for t_CE = 10⁶ yr does not improve the agreement. Table 2 contains a summary of the f_corr factors used in the different models. We highlight that it is not a contradiction to use different evolutionary stages for the distinct molecules, since the chemical model is run in post-processing with respect to the dynamical evolution. The fact that N₂D⁺ requires a shorter t with respectto the other species may be due to the fact that this deuterated species, requiring significant CO freeze-out for its efficient production, traces mainly the central parts of the pre-stellar core, where we expect the dynamical time to be the shortest due to the high volume densities.

In order to check if our observations are sensitive to the attenuation of ζ₂, or more to its average value across the core, we ran another chemical model using ζ₂ = 3 × 10⁻¹⁷ s⁻¹ (the average ζ₂ of model $\mathcal{L}$ at column densities of N = [0.2−5] × 10²² cm⁻²). The estimated abundance profiles are very similar to model $\mathcal{L}$, especially atradii r > 0.05 pc = 10⁴ AU. As a consequence, the synthetic spectra are also similar, and they do not allow us to distinguish between the cases of constant vs. attenuated ζ₂. We tested also the attenuation model from Neufeld & Wolfire (2017), which predicts a steeper decrease of ζ₂ with increasing density (ζ₂ ∝ N^−0.92). This model is in worse agreement with the observations with respect to model $\mathcal{L}$, since it predicts CR ionisation rate values lower than model $\mathcal{L}$ in the innermost parts of the core (<0.01 pc = 2000 AU), from where the targeted transitions mostly arise. In order to constrain the decrease of ζ₂ with increasing density, two complementary approaches appear necessary: (i) performing observations at higher spatial resolution (hundreds of astronomical units) of molecular lines with high critical densities, such as the one used for this work; (ii) targeting transitions with lower critical densities (10³ −10⁴ cm⁻³), which could trace the ζ₂ value in the low-density extended envelope.

Further improvements on this topic could come from expanding the sample by repeating the analysis on other pre- and proto-stellarcores. In this respect, we highlight that the two species that present the strongest sensitivity in their simulated spectra to changing ζ₂ are HC¹⁸O⁺ and N₂D⁺. The peak intensity of the HC¹⁸O⁺ (1-0) transition, in fact, changes by a factor of ≈10 going from model $\mathcal{F}$ to $\mathcal{H}$. Deuterated diazenylium is also sensitive to high ζ₂ values, which makes the abundance of this molecule drop (see for instance model $\mathcal{H}$). This, combined with the fact that N₂D⁺ has three transitions easily accessible at millimetre wavelengths and with high critical densities, makes this species a good probe of the CR ionisation rate in dense cores.

Our results appear to be in agreement with those of Bovino et al. (2020), who adopted an analytical approach and obtained ζ₂ = 2−3 × 10⁻¹⁷ s⁻¹ in L1544, and also with those of van der Tak & van Dishoeck (2000), obtained in a different molecular cloud. Moreover, they are consistent with Maret & Bergin (2007), who estimated ζ₂ = [1.5−9] × 10⁻¹⁷ s⁻¹ in the starless core Barnard 68, and within uncertainties with Maret et al. (2013), who found that ζ₂ ≈ 4.6 × 10⁻¹⁷ s⁻¹ reproduces reasonably the C¹⁸O and H¹³CO⁺ emission in two pre-stellar cores. It is worth noticing that in the recent work from Ivlev et al. (2019), the authors estimated an upper limit of ζ₂ ≈ 10⁻¹⁶ s⁻¹ using the gas temperature measurements in L1544. This upper limit, which is still consistent with our findings, was derived assuming unevolved dust with the MRN size distribution (Mathis et al. 1977); if small grains are eliminated due to a rapid coagulation, which is expected to occur in molecular clouds (Silsbee et al. 2020), the upper limit on ζ₂ is a factor of ≈3 lower.

One of the limitations of this work is that we have not considered the effects of magnetic fields. As charged particles, CR are affected by magnetic forces, and they gyrate around the magnetic field lines during their propagation. Hence, depending on the magnetic field morphology, and in particular on the ratio between the toroidal and the poloidal components of the magnetic field, CRs experience an effective column density which can be much larger than the observed one. In particular at the core centre, where the field lines are expected to be strongly entangled, the effective column density ‘seen’ by CRs can be significantly higher than N, causing a strong decrease of ζ₂ (see for instance Padovani et al. 2013). Future polarisation observations aimed at recovering the morphology of the magnetic field in L1544 will provide useful constraints in this sense.

Fig. 11 Top panels: abundance profiles obtained from the model $\mathcal{L}$ (solid lines), and from the model with constant ζ2 = 3 × 10−17 s−1 (dashed lines). The different colours refer to the four time steps investigated in this work (labelled in the first panel, in year). Bottom panels: ratios between the abundance profiles from model $\mathcal{L}$ and from themodel with constant ζ2 = 3 × 10−17 s−1. The colour code is the same as for the upper panels.

Appendix A CR ionisation rate and temperature of the models

Fig. A.1 Abundance profiles of N2H+, N2D+, HCO+, and DCO+ (from top-left to bottom-right panel) in the model $\mathcal{F}$, using the original kinetic temperature profile from Keto et al. (2015) (in black), and increasing or decreasing the temperature by 1 K (red and blue curves respectively). With respect to all other figures in this work, here we use the linear scale, since in the logarithmic scale the small variations would not be visible. The x-axis is limited to the central 0.2 pc.

Fig. A.2 Same as Fig. A.1, but for model $\mathcal{L}$.

As mentioned in the main text, increasing (or decreasing) ζ₂ affects the gas temperature(T_K). This, on theother hand, changes the excitation conditions for the different species. In principle one should self-consistently compute the temperature profiles (and the density profile) every time the CR ionisation rate is changed, which is beyond the scope of this work. Furthermore, we do not expect a significant temperature increase in the dense gas traced by the targeted molecular transitions, as the gas in this part of the core is mainly cooled by dust. This cooling mechanism is efficient enough that we do not expect important changes (see also Ivlev et al. 2019). A higher or lower ζ₂ will affect the temperature mainly in the outer part of the core, where the kinetic temperature can be constrained for instance using ammonia inversion transitions (Schmiedeke et al. in prep.).

Nevertheless, in order to check how small temperature variations impact our radiative transfer modelling and conclusions, for the three models discussed in Sects. 4.1 and 4.2 ($\mathcal{F}$, $\mathcal{L}$, and $\mathcal{H}$), we ran the chemical model in two more cases, obtained by increasing or decreasing the whole gas temperature profile uniformly by 1 K.

Figures A.1 to A.3 show the resulting abundance profiles of the four species of interest obtained after changing the gas temperature. In general, an increase of the temperature leads to an increase of the abundance for all the molecules. The changes are however small, especially for HCO⁺ (4 − 6% on average) and DCO⁺ (6 − 7%), whilst they are higher for N₂H⁺ (11 − 13%) and N₂D⁺ (12 − 18%). These small changes, together with the fact that varying the temperature in different models leads to similar abundance variations for a given molecule, leads us to conclude that taking into account these effects will not alter significantly our conclusions.

Fig. A.3 Same as Fig. A.1, but for model $\mathcal{H}$.

Appendix B Complete set of abundance profiles for models $\mathcal{L}$, $\mathcal{F}$, and $\mathcal{H}$

In Fig. B.1 we report the abundance profiles for the four targeted species obtained with model $\mathcal{L}$, $\mathcal{F}$, and $\mathcal{H}$ at the time steps not presented in Fig. 2.

Fig. B.1 Abundance profiles for the four analysed species (from left to right: N2H+, N2D+, HC18O+, DCO+) at t = 105 (top row), t = 5 × 105 yr (middle row), and t = 5 × 106 yr (bottom row). The colours refer to different models used for the CR ionisation rate: the fiducial model $\mathcal{F}$ with ζ2 = 1.3 × 10−17 s−1 is shown in blue, whilst the red and green curves show the models $\mathcal{L}$ and $\mathcal{H}$ of PI18, respectively.

Appendix C Uncertainties in the estimation of the f_corr values

As mentioned in Sect. 4.1, the non-LTE radiative transfer code requires long computational times to converge, which prevents us to perform a full exploration of the parameter space, in particular to identify the best-fit f_corr value for each species. Such a study is beyond the scope of this work, in which we focus on finding a general agreement between the distinct chemical models and the observed transitions. Finally, we highlight that this work is intrinsically characterised by large uncertainties, for instance in the physical model (L1544 is in reality not spherical), and in the chemical network (regarding e.g. rate coefficients and initial conditions). Nevertheless, it is worth to verify with a more quantitative method if the approach used in the main text to find the f_corr values is valid.

To this aim, we implemented a similar approach to that of Redaelli et al. (2018). Focussing on the results obtained for model $\mathcal{L}$, we test five f_corr factors around the value that provides the best-fit according to the analysis done in Sect. 4. For each molecular species we then compute the residuals R of the model with respect to the observations, following: $$R={\sum }_{i,j}{\left(T_{\text{MB},\text{mod}}^{i,j}(f_{\text{corr}})-T_{\text{MB},\text{obs}}^{i,j}\right)}^2,$$(C.1)

where $T_{\text{MB},\text{mod}}^{i,j}$ or $T_{\text{MB},\text{obs}}^{i,j}$ represents the modelled or observed main beam temperature in the spectral channel i for the rotational transition j. For each species we then fit a quadratic function to the R − f_corr relation. The $f_{\text{corr}}^{\text{min}}$ value that corresponds to the minimum value of the residuals is considered the best-fit value. Using standard error propagation from the uncertainties on the parameters of the quadratic relation, we estimate uncertainties on $f_{\text{corr}}^{\text{min}}$.

Figure C.1 shows the results of this analysis. The $f_{\text{corr}}^{\text{min}}$ values are shown as black vertical lines, with the shaded area indicating the corresponding ± 1σ intervals. The f_corr values found instead with the approach of Sect. 4 are shown as red vertical lines. For N₂D⁺, HC¹⁸O⁺, and DCO⁺ the differences between the two approaches are below 10% and well below the uncertainty level. Only for N₂H⁺ the two methods are not consistent. However, the discrepancy is only at a 1.2σ level. Furthermore, the value found ($f_{\text{corr}}^{\text{min}}=1.45$) is still equal or smaller to the corrective factors found in all the other tested models. We conclude that, given the uncertainties, our method is robust and equivalent to a more quantitative, but time consuming, analysis.

Fig. C.1 Sum of the residuals, over all available transitions, computed with models obtained by changing the fcorr value in the radiative transfer analysis. The synthetic spectra are produced using model $\mathcal{L}$ at t = 106 yr for all species but N2D+, for which we use the best-fit time step t = 5 × 105 yr. The dashed, blue line shows the quadratic function fitted to the R − fcorr relation. In each panel, the vertical black line shows the $f_{\text{corr}}^{\text{min}}$ value that minimises the residuals, and the shaded grey area shows the corresponding 1σ interval. The fcorr value evaluated in Sect. 4 is shown by the vertical red line. The panels refer to N2H+, N2D+, HC18O+, and DCO+ (from left toright).

Appendix D Cloud-evolution model: radiative transfer results

In this appendix we report the radiative transfer results obtained from the model with initial cloud evolution and constant ζ₂ = 1.3 × 10⁻¹⁷ s⁻¹. The modelled spectra are shown in Fig. D.1. It can be seen that, in general, this model requires equal or higher corrective factors with respect for instance to model $\mathcal{L}$ or model $\mathcal{F}$.

Fig. D.1 Radiative transfer results (in red) overlaid to the observations (in black) for all the targeted transition obtained with the chemical model with initial cloud evolution to set the initial conditions, obtained at t = 106 yr (see Sect. 4.4). The model implements constant ζ2 equal to the fiducial model (ζ2 = 1.3 × 10−17 s−1). The top row shows the DCO+ lines, while the bottom row presents the two N2H+ transitions and the HC18O+ (1-0) line, from left to right. The corresponding panels for N2D+ are presented in Fig. 10. The used fcorr value is shown on top of each panel.

Appendix E Model with ζ₂ = 3× 10⁻¹⁷ s⁻¹: resulting line profiles.

Figure E.1 shows the results of the radiative transfer applied to the chemical model with constant ζ₂ = 3× 10⁻¹⁷ s⁻¹.

Fig. E.1 Radiative transfer results (in red) overlaid to the observations (in black) for all the targeted transition obtained with the chemical model with constant ζ2 = 3 × 10−17 s−1, using the same evolutionary time steps and fcorr factors that provide the best fit for the fiducial model $\mathcal{F}$. The top row shown the N2D+ line, the middle row shows the DCO+ lines, and the bottom row presents the two N2H+ transitions and the HC18O+ (1-0) line, from left to right. The used time step and fcorr value are shown in the title of each panel.

References

All Tables

All Figures

Fig. 1 Left panel: H2 volume density profile (n) in thephysical model of L1544. Central panel: column density profile. Right panel: ζ2 profiles for the $\mathcal{F}$ (dashed bluecurve), $\mathcal{L}$ (solid red), and $\mathcal{H}$ (dashed-dotted green) models. The profiles for the models $\mathcal{L}$ and $\mathcal{H}$ PI18 have been derived integrating the volume density profile of the core model to obtain the column density profile, and then using the relation between ζ2 and N.

In the text

Fig. 2 Top panels: abundance profiles for the four analysed species (from left to right: N2H+, N2D+, HC18O+, DCO+) at t = 106 yr. The colours refer to different models used for the CR ionisation rate: model $\mathcal{F}$ (blue), model $\mathcal{L}$ (red), and model $\mathcal{H}$ (green). Bottom panels: the green and red solid curves show the ratio of model $\mathcal{H}$ and model $\mathcal{L}$ to model $\mathcal{F}$, respectively. The dashed, horizontal lines represent the inverse ratios of the corrective factors fcorr used to reproduce the observed lines, as reported in Sect. 4.1 and Table 2. To give an example, the radiative transfer of DCO+ lines requires fcorr = 0.9 in model $\mathcal{L}$, and fcorr = 1.5 in model $\mathcal{F}$; therefore thered solid curve $\mathcal{L}$ / $\mathcal{F}$ in the rightmost panel has to be compared with the red dashed line at ${(0.9/1.5)}^{-1}=1.67$.

In the text

Fig. 3 Spectra obtained with the radiative transfer code (model, in red) overlaid to the observations (black) for N2 H+ (1–0) and (3–2) lines. The columns refer to the three models: $\mathcal{F}$ in the left panels, $\mathcal{L}$ in the central panels, and $\mathcal{H}$ in the right ones. The corrective factors fcorr used to modify the whole abundance profile to match the observed fluxes is indicated above the top panel of each column. All models are at t = 106 yr.

In the text

	Fig. 4 As Fig. 3, but for N2D+.
In the text

	Fig. 5 As Fig. 3, but for HC18O+.
In the text

	Fig. 6 As Fig. 3, but for DCO+.
In the text

	Fig. 7 Abundance profiles for N2D+ in model $\mathcal{L}$ at four different time steps, from t = 105 yr to t = 5 × 106 yr. One can see that as the chemistry evolves, N2D+ is affected by stronger depletion, with a higher contrast between the central abundance and the peak value.
In the text

	Fig. 8 x(e) profiles in model $\mathcal{F}$ (dashed bluecurve), $\mathcal{L}$ (solid red), and $\mathcal{H}$ (dashed-dotted green). All the models are shown at t = 106 yr.
In the text

Fig. 9 Top panels: abundance profiles obtained with the model with constant ζ2 = 1.3 × 10−17 s−1, with no initial cloud evolution (model $\mathcal{F}$, solid curves) and with a cloud evolution of tCE = 106 yr to set the initial conditions (dashed curves). The different colours refer to the four time steps investigated in this work (labelled in the first panel, in year). Bottom panels: ratios between the abundance profiles of the model with and without cloud evolution. The colour code is the same as for the upper panels.

In the text

Fig. 10 Best-fit model for N2D+ lines obtained with the fiducial model with initial cloud evolution and ζ2 = 1.3 × 10−17 s−1, at t = 106 yr (red curves),overlaid to the observations (black). The used fcorr value is indicated at the top. The modelled lines produced using model $\mathcal{F}$ with no cloud evolution at t = 106 yr (and fcorr = 5) are shown in green. It is clear that the model with cloud evolution underestimates the (2–1) and (3–2) lines.

In the text

Fig. 11 Top panels: abundance profiles obtained from the model $\mathcal{L}$ (solid lines), and from the model with constant ζ2 = 3 × 10−17 s−1 (dashed lines). The different colours refer to the four time steps investigated in this work (labelled in the first panel, in year). Bottom panels: ratios between the abundance profiles from model $\mathcal{L}$ and from themodel with constant ζ2 = 3 × 10−17 s−1. The colour code is the same as for the upper panels.

In the text

Fig. A.1 Abundance profiles of N2H+, N2D+, HCO+, and DCO+ (from top-left to bottom-right panel) in the model $\mathcal{F}$, using the original kinetic temperature profile from Keto et al. (2015) (in black), and increasing or decreasing the temperature by 1 K (red and blue curves respectively). With respect to all other figures in this work, here we use the linear scale, since in the logarithmic scale the small variations would not be visible. The x-axis is limited to the central 0.2 pc.

In the text

	Fig. A.2 Same as Fig. A.1, but for model $\mathcal{L}$.
In the text

	Fig. A.3 Same as Fig. A.1, but for model $\mathcal{H}$.
In the text

Fig. B.1 Abundance profiles for the four analysed species (from left to right: N2H+, N2D+, HC18O+, DCO+) at t = 105 (top row), t = 5 × 105 yr (middle row), and t = 5 × 106 yr (bottom row). The colours refer to different models used for the CR ionisation rate: the fiducial model $\mathcal{F}$ with ζ2 = 1.3 × 10−17 s−1 is shown in blue, whilst the red and green curves show the models $\mathcal{L}$ and $\mathcal{H}$ of PI18, respectively.

In the text

Fig. C.1 Sum of the residuals, over all available transitions, computed with models obtained by changing the fcorr value in the radiative transfer analysis. The synthetic spectra are produced using model $\mathcal{L}$ at t = 106 yr for all species but N2D+, for which we use the best-fit time step t = 5 × 105 yr. The dashed, blue line shows the quadratic function fitted to the R − fcorr relation. In each panel, the vertical black line shows the $f_{\text{corr}}^{\text{min}}$ value that minimises the residuals, and the shaded grey area shows the corresponding 1σ interval. The fcorr value evaluated in Sect. 4 is shown by the vertical red line. The panels refer to N2H+, N2D+, HC18O+, and DCO+ (from left toright).

In the text

Fig. D.1 Radiative transfer results (in red) overlaid to the observations (in black) for all the targeted transition obtained with the chemical model with initial cloud evolution to set the initial conditions, obtained at t = 106 yr (see Sect. 4.4). The model implements constant ζ2 equal to the fiducial model (ζ2 = 1.3 × 10−17 s−1). The top row shows the DCO+ lines, while the bottom row presents the two N2H+ transitions and the HC18O+ (1-0) line, from left to right. The corresponding panels for N2D+ are presented in Fig. 10. The used fcorr value is shown on top of each panel.

In the text

Fig. E.1 Radiative transfer results (in red) overlaid to the observations (in black) for all the targeted transition obtained with the chemical model with constant ζ2 = 3 × 10−17 s−1, using the same evolutionary time steps and fcorr factors that provide the best fit for the fiducial model $\mathcal{F}$. The top row shown the N2D+ line, the middle row shows the DCO+ lines, and the bottom row presents the two N2H+ transitions and the HC18O+ (1-0) line, from left to right. The used time step and fcorr value are shown in the title of each panel.

In the text