The cosmic-ray ionisation rate in the pre-stellar core L1544

Context. Cosmic rays (CRs) play an important role in the chemistry and dynamics of the interstellar medium. In dense environments, they represent the main ionising agent, driving the rich chemistry of molecular ions and determining the ionisation fraction, which regulates the degree of coupling between the gas and magnetic fields. Estimates of the CR ionisation rate ($\zeta_2$) span several orders of magnitude, depending on the targeted sources and on the used method. Aims. Recent theoretical models have characterised the CR attenuation with increasing density. We aim to test these models for the attenuation of CRs in the low-mass pre-stellar core L1544. Methods. We use a state-of-the-art gas-grain chemical model, which accepts the CR ionisation rate profile as input, to predict the abundance profiles of four ions: $\rm N_2H^+$, $\rm N_2D^+$, $\rm HC^{18}O^+$, and $\rm DCO^+$. Non-LTE radiative transfer is performed to produce synthetic spectra based on the derived abundances. These are compared with observations obtained with the Institut de Radioastronomie Millim\'etrique (IRAM) 30m telescope. Results. Our results indicate that a model with $\zeta_2>10^{-16} \rm \, s^{-1}$ is excluded by the observations. Also the model with the standard $\zeta_2 = 1.3 \times 10^{-17} \rm \, s^{-1}$ produces a worse agreement with respect to the attenuation model based on Voyager observations, which has an average $\zeta_2 = 3 \times 10^{-17} \rm \, s^{-1}$ at the column densities typical of L1544. The single-dish data, however, are not sensitive to the attenuation of the CR profile, which changes only by a factor of two in the range of column densities spanned by the core model. Interferometric observations at higher spatial resolution, combined with observations of transitions with lower critical density are needed to observe a decrease of $\zeta_2$ with density.


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 21 cm −2 (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: The ionised H + 2 quickly reacts with another hydrogen molecule, producing the fundamental trihydrogen cation (H + 3 ). H + 3 in turn is the starting point for the chain of reactions between charged and neutral species, producing several other key molecules, such as HCO + , via reaction with CO. Indirectly, CRs regulate also another important chemical process in the cold phases of the ISM: deuteration. Arguably the most important deuterated species in the gas phase is H 2 D + , which is produced by isotopic exchange This work is based on observations carried out with the IRAM 30m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). reactions between 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 2 , 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 magnetic fields (see e.g. the results of Planck Collaboration et al. 2016 for the magnetic properties of molecular clouds), which exert their influence directly on charged particles. 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, 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 it is represented by direct observations of H + 3 . These, however, are possible only in diffuse gas via infrared spectroscopy. Lacking a permanent electric dipole, in fact, H + 3 is not observable in emission in the interstellar gas at high densities. In dense clouds, van der Tak & van Dishoeck (2000) used H + 3 absorption lines in combination with H 13 CO + rotational lines towards background massive protostars, and they determined a CR ionisation rate of ≈ 3 × 10 −17 s −1 . More recently, Indriolo & McCall (2012) -expanding the sample from Indriolo et al. (2007)-reported 21 detections toward diffuse clouds, determining CR ionisation rate per hydrogen molecule (ζ 2 ) in the range [1.7 − 10.6] × 10 −16 s −1 . Neufeld & Wolfire (2017) derived a CR ionisation rate per hydrogen atom (ζ 1 ) of the order of a few 10 −16 s −1 in the Galactic disc, using H + 3 and other ionised tracers. Their data support a steep dependency of ζ 2 on the gas column density (ζ 2 ∝ N −1 ).
In regions so dense that no background source is visible at the near infrared wavelengths of the H + 3 lines, no 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 −18 − 10 −16 s −1 . Bovino et al. (2020) proposed a new analytical method to assess ζ 2 using observations of CO, of the deuterium fraction (for instance from HCO + and DCO + ), and of ortho-H 2 D + . The method was applied by Sabatini et al. (2020) to a sample of massive clumps, finding ζ 2 ≈ [0.7 − 6] × 10 −17 s −1 .
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 ζ 1 ≈ 6.8 × 10 −18 s −1 . This can be transformed into ζ 2 using the relation ζ 2 = (2.3/1.5) × ζ 1 (Glassgold & Langer 1974), obtaining the value ζ 2 ≈ 10 −17 s −1 which is still considered the standard value for molecular clouds. However, the assumption that ζ 2 is independent on the ISM density is incorrect, as also shown by the different observational results. CRs lose energy in their interactions with the ambient medium, which modifies their spectrum. These losses hence cause an attenuation of ζ 2 , which decreases with gas column density (N). This process has been studied in detail by Padovani et al. (2009), who computed the relation between the CR ionisation rate and the gas column density using: 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 2 ionisation, σ ion k (E) is the ionisation cross section of H 2 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 2 ionisation. Their contribution is described by Φ(N), which is the ratio between the secondary to primary H 2 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 focusing prac-tically 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 ζ 2 10 −16 s −1 up to N ≈ a few × 10 23 cm −2 , and reproduces the high ζ 2 measured in the diffuse medium. The "low" model, which instead implies ζ 2 ≈ a few 10 −17 s −1 , is the one that fits 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 from PI18 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 mm and submm 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 2 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 6 cm −3 ) 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 2 H + (diazenylium) across the source. We now perform a similar analysis, focusing this time on ζ 2 . 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.

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 18 O, H 2 O, and N 2 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. 2018Redaelli et al. , 2019. In the latest paper, we modelled the N 2 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 2 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 2 D + ), which are inconsistent with the observations. For this work we therefore use 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.

Chemical models
The physical model of L1544 is coupled to the gas-grain chemical model presented in Sipilä et al. (2015aSipilä et al. ( ,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 −3 pc in the central parts to 2 × 10 −2 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 −3 . The effect of the embedding cloud is simulated by assuming an external visual extinction of A V = 2 mag 1 , which attenuates the external UV field. We test four different time steps: t = 10 5 , 5 × 10 5 , 10 6 , and 5 × 10 6 yr. Table 1. Initial abundances with respect to the total number density of hydrogen nuclei used in the chemical model.
The first model is similar to that used in RB19. In order to implement models L and H , we have used the parametrisation presented in Appendix F of PI18 to compute the ζ 2 profile from the column density profile. The radial profiles of the three ζ 2 models are presented in the right panel of Fig. 1. In models L and H , the ζ 2 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 ζ 2 values that they span, and can be ordered by increasing ζ 2 as F , L , and then 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 −3 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 focusing 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 −2 pc = 2000 AU), and therefore we can safely employ the results of PI18.
We stress here that the models F , L , and H differ in the shape of the ζ 2 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 5 cm −3 ) 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 18 O + , which is the observed isotopologue, are derived from those of HCO + , by dividing the latter by the isotopic ratio 16 O/ 18 O = 557 (Wilson 1999).

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 use single-dish data obtained with the Institut de Radioastronomie Millimétrique (IRAM) 30m telescope of: N 2 H + (1-0) and (3-2); N 2 D + (1-0), (2-1), and (3-2); DCO + (1-0), (2-1), and (3-2); HC 18 O + (1-0). The latest is preferred to the optically thick main isotopologue. All the spectra are obtained as single pointings taken at the millimetre dust peak of L1544 (R.A.(J2000) = 05 h 04 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 these data, which makes 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 −1 ), 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 2 H + (1-0) and (3-2) transitions are n c = 6 × 10 4 cm −3 and n c = 1.4 × 10 6 cm −3 , respectively (Shirley 2015). This fact, combined with the improved angular resolution at higher frequencies, allows to test reliably the radial profiles of these molecules, even though the spatial resolution of the single-dish data is limited for the low-J transitions. Figure 2 shows the simulated abundance profiles for the different species of interest at t = 10 6 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 2 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 18 O + , but it is less evident because these molecules are already significantly depleted early on (t 10 5 yr) at radii smaller than 0.04 pc. The effect of varying ζ 2 is species-dependent. HCO + (and as a consequence HC 18 O + ) presents a clear, monotonic effect at every time step: the higher the value of ζ 2 , 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 2 H + isotopologues is more complex, and timedependent. In general, we note that the higher the overall ζ 2 , the stronger is the depletion, and the earlier this phenomenon starts showing up. Model H predicts a strong drop in the abundance of N 2 H + and N 2 D + starting at t = 10 6 yr. At t = 5 × 10 6 yr, also model L (which has an intermediate ζ 2 value between models F and H ) predicts an enhanced decrease of the abundances. For N 2 D + in particular, the high model tends to always predict the lowest abundances with respect to the other two models.

Non-LTE modelling at t = 10 6 yr
In order to produce the synthetic spectra to be compared with the observations we have used the non-local-thermodynamicequilibrium (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 2 ; 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 −1 , 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 2 H + and N 2 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 −3 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 convolve the results from MOLLIE to the corresponding beam size for each transition, and we extract 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 10 2 10 1 10 13 10 12 10 11 10 10 10 9 X mol N 2 H + 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 focus on the evolutionary time of t = 10 6 yr, which in RB19 provided the best fit for all the molecules of interest. We hence run 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 particular on the parameters listed above at the 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 6 yr are shown in Figs. 3 to 6, and they are discussed in detail in the following.
N 2 H + None of the models at t = 10 6 yr is able to reproduce the intensity ratio of the two observed rotational transitions for N 2 H + . Models F and 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 H . This latter, furthermore, requires an enhancement of 400% of the overall abundance profile in order to match the observed line intensities. The abundance obtained with model F needs to be increased by f corr = 1.5, whilst model 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 2 D + All models tend to underestimate the fluxes of the three N 2 D + transitions, similarly to what was found in RB19. Model H is characterised by the strongest underestimation, and the N 2 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 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 18 O + Also for this molecule, model L provides the best agreement in terms of flux. None of the models however is able to reproduce the symmetric double-peak features 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 H has to be reduced by a factor of ≈ 3 to reproduce the observed fluxes, whereas model F has to be enhanced by f corr = 1.5. For model L , a modification of 10% is enough to obtain a good agreement with the observations. Based on these results, we conclude first of all that the model H , with the highest ζ 2 profile, is excluded by our data. In particular, it underestimates severely the abundance of N 2 D + , and to a lesser extent also that of N 2 H + , whilst it overestimates the abundances of HCO + isotopologues. Models F and L are instead in better agreement with the observations The latter appears better than model 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 2 H + , HC 18 O + , and DCO + (with a correction of only 10%), at the same time underestimating the fluxes of the N 2 D + lines by the smallest factor. The bottom panels of Fig. 2 show the ratio between the models (L /F and H /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 −2 pc, corresponding to volume densities higher than ≈ 10 5 cm −3 (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 4 cm −3 to ≈ 10 6 cm −3 .

Testing different evolutionary stages
The previous subsection has shown that model 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 L at three further time steps: t = 10 5 yr, t = 5 × 10 5 yr, and t = 5 × 10 6 yr. These abundance profiles are shown in Appendix B (together with the ones for models F and 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 modify 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 2 H + , no time step improves significantly the agreement with the observations. In particular, the model at t = 10 5 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 6 yr underestimates the observed peak intensity by more than 50% for both transitions. The models at t = 5 × 10 5 yr and t = 10 6 yr are similar in terms of the ratio of the two rotational lines intensities, with the model at t = 10 6 yr providing a slightly better agreement.
The lines of N 2 D + , on the other hand, are better reproduced when adopting an earlier evolutionary time (t = 5 × 10 5 yr), and f corr = 3. In order to better understand the evolution of the N 2 D + profile, we show in Fig. 7 the abundance profiles at four time steps for model L . As the chemistry evolves, the peak increases in value and it shifts outwards, whilst at the centre the abundance decreases. In Sect. 4.1 we showed how at t = 10 6 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 5 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 6 yr). The intermediate time step t = 5 × 10 5 yr presents instead the best combination of total molecular density and depletion level.
For DCO + and HC 18 O + , the model at t = 10 6 yr requires the smallest corrective factors to reproduce the observed lines, and 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 2 D + , which is formed only when CO is depleted from the gas phase  Abundance profiles for N 2 D + in model L at four different time steps, from t = 10 5 yr to t = 5×10 6 yr. One can see that as the chemistry evolves, N 2 D + is affected by stronger depletion, with a higher contrast between the central abundance and the peak value. 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 2 D + is better reproduced at an earlier chemical timescale with respect to the other species. A better performance of the two remaining ζ 2 models can be excluded based on the following arguments. The main problems of model H are the severe underestimation of N 2 D + abundance (by more than one order of magnitude) and a large overestimation of HC 18 O + abundance (by a factor of 4). At any time step, model H presents the highest abundance of HCO + (and thus of HC 18 O + ), which will always be overestimated. Concerning N 2 D + , model H always predicts the lowest abundance, with the exception of the earliest evolutionary time (t = 10 5 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 are not able to reproduce the intensity ratios of the distinct J transitions.
The main problems of model F at t = 10 6 yr are the underestimation of HC 18 O + (by a factor of 2.6) and of N 2 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 2 D + , we tested the same time-step that provides the best fit for the model L : t = 5 × 10 5 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 H and F , with respect to model L . parameter, which is shown in Fig. 8 for the three models F , L , and H . We focus on the model results at t = 10 6 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 ζ 2 , as expected due to the nature of the ionisations. At the core centre, the value obtained in model L is x(e) = 10 −9 . This result is in agreement with the one from Caselli et al. (2002), who used more simplified chemical models on a subset of the molecular transitions presented in this work to estimate this quantity.

Simulating the cloud evolution
As a further test on the fiducial model with constant ζ 2 = 1.3 × 10 −17 s −1 , we changed the initial conditions of the chemical model. We simulate first a cloud-like evolution for t CE = 10 6 yr. For this initial simulation, we use the following physical conditions: n = 10 3 cm −3 , T = 15 K, and external extinction A V = 2 mag. The final abundances of this model are then taken as initial condition 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 2 H + and N 2 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 5 yr).
As a result of these differences, the cloud-evolution model requires in general higher corrective factor 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 6 yr. We discuss here the case of N 2 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 L . 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 2 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.

Constant
The models so far explored are characterised by average ζ 2 values that are significantly different: ζ 2 = const = 1.3 × 10 −17 s −1 in model F , ζ 2 = 3 × 10 −17 s −1 in model L , and ζ 2 = 2 × 10 −16 s −1 in model H . The last two values refer to gas column densities in the range N = [0.2 − 5] × 10 22 cm −2 (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 ζ 2 .
In order to test this, we have run another chemical model, identical to the fiducial one but adopting a constant CR ionisation rate equal to ζ 2 = 3 × 10 −17 s −1 (the average value of model L ). The resulting abundance profiles are shown in Fig. 11, in comparison to those of the model L . The two models produce very similar profiles, in particular at the two time steps that provide the best-fit (t = 5 × 10 5 yr and t = 10 6 yr). The median relative differences at these time steps are 6 − 9% for N 2 H + , 8 − 10% for N 2 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 resolution 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 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 ζ 2 ∝ N −0.92 (see their Figure 6). The reader may note that in the range of column densities spanned by the physical model, the model L predicts a shallow dependency on the column density, approximately following ζ 2 ∝ N −0.2 . We tested the steeper relation derived from Neufeld & Wolfire (2017) in our chemical model 3 . The resulting abundances underestimate the observed transitions with respect to model 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 ζ 2 < 3 × 10 −17 s −1 . As a Model with Cloud Evolution, f corr = 5.0 Fig. 10. Red curves show the best-fit model for N 2 D + lines obtained with the fiducial model with initial cloud evolution and ζ 2 = 1.3 × 10 −17 s −1 , at t = 10 6 yr, overlaid to the observations (black). The used f corr value is indicated at the top. The modelled lines produced using model F with no cloud evolution at t = 10 6 yr (and f corr = 5) are shown in green. It is clear that the model with cloud evolution underestimates the (2-1) and (3-2) lines.
consequence, in the central regions -where the targeted transitions are excited, see Sect. 4.1-the lower CR ionisation rate leads in general to a lower predicted abundance, in particular for HC 18 O + and DCO + , for which the dependence on ζ 2 is monotonic, as discussed in Sect. 3.

Summary and Conclusions
In this work we have studied the effect of the CR ionisation rate on the results of gas-grain chemical simulations, focusing on four species: N 2 H + , N 2 D + , HC 18 O + , and DCO + . The predicted abundance profiles have been used to perform a full non-LTE radiative transfer analysis, and the synthetic spectra have been 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 ζ 2 (H from PI18), characterised by an average value of ζ 2 = 2 × 10 −16 s −1 , is excluded by our observations. In particular, this model strongly underestimates the N 2 D + abundance, while it overestimates by a factor of 3 − 4 the intensities of HC 18 O + and DCO + lines.
The model L from PI18 is the one that provides the overall best agreement with the observations. It is remarkable that, with the exception of N 2 D + -which is underestimated by all the models here considered-the other targeted molecules are reproduced within 10%. The fiducial model, with standard ζ 2 = 1.3 × 10 −17 s −1 , performs worse than model L , and requires corrections in the range f corr = 1.5 − 5 at t = 10 6 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 6 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 E. Redaelli et al.: CR ionisation rate in L1544 evolution. The fact that N 2 D + requires a shorter t with respect to 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 ζ 2 , or more to its average value across the core, we have run another chemical model using ζ 2 = 3×10 −17 s −1 (the average ζ 2 of model L at column densities of N = [0.2−5]×10 22 cm −2 ). The estimated abundance profiles are very similar to model L , especially at radii r > 0.05 pc = 10 4 AU. As a consequence, the synthetic spectra are also similar, and do not allow us to distinguish between the cases of constant vs. attenuated ζ 2 . We tested also the attenuation model from Neufeld & Wolfire (2017), which predicts a steeper decrease of ζ 2 with increasing density (ζ 2 ∝ N −0.92 ). The latter is in worse agreement with the observations with respect to model L , since it predicts CR ionisation rate values lower than model 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 ζ 2 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 3 − 10 4 cm −3 ), which could trace the ζ 2 value in the low-density extended envelope.
Further improvements on this topic could come from expanding the sample by repeating the analysis on other prestellar/protostellar cores. In this respect, we highlight that the two species that present the strongest sensitivity in their simulated spectra to changing ζ 2 are HC 18 O + and N 2 D + . The peak intensity of the HC 18 O + (1-0) transition, in fact, changes by a factor of ≈ 10 going from model F to H . N 2 D + is also sensitive to high ζ 2 values, which makes the abundance of this molecule drop (see for instance model H ). This, combined with the fact that N 2 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 = 2−3×10 −17 s −1 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 ζ 2 = [1.5−9]×10 −17 s −1 in the starless core Barnard 68, and within uncertainties with Maret et al. (2013), who found that ζ 2 ≈ 4.6×10 −17 s −1 reproduces reasonably the C 18 O and H 13 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 ζ 2 ≈ 10 −16 s −1 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 ζ 2 is a factor of ≈ 3 lower.
One of the limitations of this work is that we are not considering 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 ζ 2 (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.
A&A proofs: manuscript no. L1544_Ionisation_IIIsub  As mentioned in the main text, increasing (or decreasing) ζ 2 affects the gas temperature (T K ). This, on the other 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 ζ 2 will affect mainly the temperature 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 (F , L , and H ), we have run the chemical model in two more cases, obtained by increasing/decreasing the whole gas temperature profile uniformly by 1 K. Figs. 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 2 H + (11 − 13%) and N 2 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. 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 Article number, page 13 of 17 A&A proofs: manuscript no. L1544_Ionisation_IIIsub   Fig. B.1. Abundance profiles for the four analysed species (from left to right: N 2 H + , N 2 D + , HC 18 O + , DCO + ) at t = 10 5 (top row), t = 5 × 10 5 yr (middle row), and t = 5 × 10 6 yr (bottom row). The colours refer to different models used for the CR ionisation rate: the fiducial model F with ζ 2 = 1.3 × 10 −17 s −1 is shown in blue, whilst the red and green curves show the models L and H of PI18, respectively. 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). Focusing on the results obtained for model 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: where T i, j MB,mod/obs represents the modelled/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 min corr 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 min corr . Figure C.1 shows the results of this analysis. The f min corr 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 2 D + , HC 18 O + , and DCO + the differences between the two approaches are below 10% and well below the uncertainty level. The synthetic spectra are produced using model L at t = 10 6 yr for all species but N 2 D + , for which we use the best-fit time step t = 5 × 10 5 yr. The dashed, blue line shows the quadratic function fitted to the R − f corr relation. In each panel, the vertical black line shows the f min corr value that minimises the residuals, and the shaded grey area shows the corresponding 1σ interval. The f corr value evaluated in Sect. 4 is shown by the vertical red line. The panels refer to N 2 H + , N 2 D + , HC 18 O + , and DCO + (from left to right).
Only for N 2 H + the two methods are not consistent. However, the discrepancy is only at a 1.2σ level. Furthermore, the value found ( f min corr = 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.

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 ζ 2 = 1.3 × 10 −17 s −1 . 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 L or model F .