Open Access
Issue
A&A
Volume 656, December 2021
Article Number A109
Number of page(s) 16
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202141776
Published online 10 December 2021

© E. Redaelli et al. 2021

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Open Access funding provided by Max Planck Society.

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 ≳ 1021 cm−2 (corresponding to visual extinction in magnitude AV ≳ 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: (1)

The ionised quickly reacts with another hydrogen molecule, producing the fundamental trihydrogen cation (). 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 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 H2D+, which is produced by isotopic exchange reactions between 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 N2, 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 . These, however, are possible only in diffuse gas via infrared spectroscopy. Lacking a permanent electric dipole, in fact, is not observable in emission in the interstellar gas at high densities. In dense clouds, van der Tak & van Dishoeck (2000) used absorption lines in combination with H13CO+ 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 towards 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 and other ionised tracers. Their data support a steep dependency of ζ2 on the gas column density (ζ2N−1).

In regions so dense that no background source is visible at the near infrared wavelengths of the lines, no 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-H2D+. 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 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 ζ2, 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: (2)

where k indicates the primary CR species (e.g. protons, electrons, positrons, and heavier nuclei), Eion = 15.44 eV is the energy threshold for H2 ionisation, is the ionisation cross section of H2 by the species k (Rudd et al. 1992; Kim et al. 2000), and jk(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 H2 ionisation. Their contribution is described by Φ(N), which is the ratio between the secondary to primary H2 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 ζ2 ≳ 10−16 s−1 up to N ≈ a few × 1023 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 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 N2H+ and carbon monoxide spectra (e.g. Keto et al. 2015). This model predicts low temperatures (T ≈6 K) and high volume densities (n ≳ 106 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 N2H+ (diazenylium) across the source. We now perform a similar analysis, focussing 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.

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 C18O, H2O, and N2H+ 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 N2H+ 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 N2H+ (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 N2D+), 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.

thumbnail 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 (dashed bluecurve), (solid red), and (dashed-dotted green) models. The profiles for the models and 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−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 was simulated by assuming an external visual extinction of AV = 2 mag1, which attenuates the external UV field. We tested four different time steps: t = 105, 5 × 105, 106, and 5 × 106 yr.

The chemical code takes as input the CR ionisation-rate ζ2, 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 (), has a constant value ζ2 = 1.3 × 10−17 s−1; the second one reproduces the high model () from PI18, whilst the third one corresponds to the low model () from PI18. The first model is similar to that used in RB19. In order to implement models and , 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 and , the ζ2 profile is computed assuming a column density value corresponding to AV = 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 , , and then .

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 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−2 pc = 2000 AU), and therefore we can safely employ the results of PI18.

We stress here that the models , , and 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 ≳ 105 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 HC18O+, which is the observed isotopologue, are derived diving those of HCO+ by the isotopic ratio 16O∕18O = 557 (Wilson 1999).

Table 1

Initial abundances with respect to the total number density of hydrogen nuclei used in the chemical model.

thumbnail 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 (blue), model (red), and model (green). Bottom panels: the green and red solid curves show the ratio of model and model to model , 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 , and fcorr = 1.5 in model ; therefore thered solid curve / in the rightmost panel has to be compared with the red dashed line at .

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: N2H+ (1–0) and (3–2); N2D+ (1–0), (2–1), and (3–2); DCO+ (1–0), (2–1), and (3–2); HC18O+ (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) = 05h04′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−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 N2H+ (1–0) and (3–2) transitions are nc = 6 × 104 cm−3 and nc = 1.4 × 106 cm−3, 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 = 106 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, N2H+ 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 HC18O+, but it is less evident because these molecules are already significantly depleted early on (t ≲105 yr) at radii smaller than0.04 pc.

The effect of varying ζ2 is species-dependent. HCO+ (and as a consequence HC18O+) 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 N2H+ isotopologues is more complex, and time-dependent. In general, we note that the higher the overall ζ2, the stronger is the depletion, and the earlier this phenomenon starts showing up. Model predicts a strong drop in the abundance of N2H+ and N2D+ starting at t = 106 yr. At t = 5 × 106 yr, also model (which has an intermediate ζ2 value between models and ) predicts an enhanced decrease of the abundances. For N2D+ 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 = 106 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 fcorr 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 N2H+ and N2D+), and a full parameter-space exploration is not feasible, in particular to determine fcorr. As in RB19, we therefore made several tests to determine by trial-and-error the fcorr 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 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 fcorr 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 = 106 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 fcorr 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 fcorr 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 = 106 yr are shown in Figs. 36, and they are discussed in detail in the following.

N2H+

None of the models at t = 106 yr is able to reproduce the intensity ratio of the two observed rotational transitions for N2H+. Models and 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 . 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 needs to be increased by fcorr = 1.5, whilst model 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.

N2D+

All models tend to underestimate the fluxes of the three N2D+ transitions, similarly to what was found in RB19. Model is characterised by the strongest underestimation, and the N2D+ 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 requires the smallest fcorr 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%.

HC18O+

Also for this molecule, model 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 fcorr factors among the targeted species, and is the one for which the intensity ratios of the different J-transitions are better reproduced. Model has to be reduced by a factor of ≈3 to reproduce the observed fluxes, whereas model has to be enhanced by fcorr = 1.5. For model , 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 , with the highest ζ2 profile, is excluded by our data. In particular, it severely underestimates the abundance of N2D+, and to a lesser extent also that of N2H+, whilst it overestimates the abundances of HCO+ isotopologues. Models and are instead in better agreement with the observations The latter appears better than model , since it requires the smallest fcorr 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 N2H+, HC18O+, and DCO+ (with a correction of only 10%), at the same time underestimating the fluxes of the N2D+ lines by the smallest factor.

The bottom panels of Fig. 2 show the ratio between the models ( / and / ), compared with the inverse ratio of the corresponding fcorr values used to obtain the best-fit solutions. Qualitatively, the range of radii where the model ratios are similar to the fcorr 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 ≈105 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 104to ≈106 cm−3.

thumbnail 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: in the left panels, in the central panels, and 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 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 at three further time steps: t = 105 yr, t = 5 × 105 yr, and t = 5 × 106 yr. These abundance profiles are shown in Appendix B (together with the ones for models and ). As previously described, we tested a few fcorr 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 nc) 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 N2H+, no time step significantly improves the agreement with the observations. In particular, the model at t = 105 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 × 106 yr underestimates the observed peak intensity by more than 50% for both transitions. The models at t = 5 × 105 yr and t = 106 yr are similar in terms of the ratio of the two rotational lines intensities, with the model at t = 106 yr providing a slightly better agreement.

The lines of N2D+, on the other hand, are better reproduced when adopting an earlier evolutionary time (t =5 × 105 yr), and fcorr = 3. In order to better understand the evolution of the N2D+ profile, we show in Fig. 7 the abundance profiles at four time steps for model . 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 = 106 yr the line intensity ratios are well reproduced, but the abundance profile has to be increased by a factor fcorr = 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 = 105 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 × 106 yr). The intermediate time step t = 5 × 105 yr presents instead the best combination of total molecular density and depletion level.

For DCO+ and HC18O+, the model at t = 106 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 N2D+, 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 N2D+ 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 are the severe underestimation of N2D+ abundance (by more than one order of magnitude) and a large overestimation of HC18O+ abundance (by a factor of 4). At any time step, model presents the highest abundance of HCO+ (and thus of HC18O+), which will always be overestimated. Concerning N2D+, model always predicts the lowest abundance, with the exception of the earliest evolutionary time (t =105 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 at t = 106 yr are the underestimation of HC18O+ (by a factor of 2.6) and of N2D+ (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 N2D+, we tested the same time-step that provides the best fit for the model : t = 5 × 105 yr and fcorr = 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 and , with respect to model .

thumbnail Fig. 4

As Fig. 3, but for N2D+.

thumbnail Fig. 5

As Fig. 3, but for HC18O+.

thumbnail Fig. 6

As Fig. 3, but for DCO+.

thumbnail Fig. 7

Abundance profiles for N2D+ in model 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 , , and . We focus on the model results at t = 106 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 is x(e)= 10−9. 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.

thumbnail Fig. 8

x(e) profiles in model (dashed bluecurve), (solid red), and (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 ζ2 = 1.3 × 10−17 s−1, we changed the initial conditions of the chemical model. We simulated first a cloud-like evolution for tCE = 106 yr. For this initial simulation, we used the following physical conditions: n = 103 cm−3, T =15 K, and external extinction AV = 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 tCE).

In Fig. 9 we show the final abundance profiles for the species of interest from the two models, with and without cloud evolution. N2H+ and N2D+ 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 =105 yr).

As a result of these differences, the cloud-evolution model requires in general higher corrective factors fcorr to reproduce the observed fluxes. We report in Appendix D the complete set of comparisons between the synthetic and observed spectra at t = 106 yr. We discuss here the case of N2D+, 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 fcorr = 5, in Fig. 10, compared with the best-fit solution found for model . 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 N2D+, together with the fact that the corrective factors fcorr 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 tCE. Furthermore, neglecting the dynamical evolution from cloud to core is also a source of uncertainty.

4.5 Constant ζ2 = 3 × 10−17 s−1 model

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 , ⟨ ζ2 ⟩ = 3 × 10−17 s−1 in model , and ⟨ ζ2 ⟩ = 2 × 10−16 s−1 in model . The last two values refer to gas column densities in the range N = [0.2−5] × 1022 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 ran 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 ). The resulting abundance profiles are shown in Fig. 11, in comparison to those of model . The two models produce very similar profiles, in particular at the two time steps that provide the best-fit (t =5 × 105 yr and t = 106 yr). The median relative differences at these time steps are 6−9% for N2H+, 8− 10% for N2D+, 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 ) 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 ζ2N−0.92 (see their Fig. 6). The reader may note that in the range of column densities spanned by the physical model, the model predicts a shallow dependency on the column density, approximately following ζ2N−0.2. We tested the steeper relation derived from Neufeld & Wolfire (2017) in our chemical model3. The resulting abundances underestimate the observed transitions with respect to model . 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 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 HC18O+ and DCO+, for which the dependence on ζ2 is monotonic, as discussed in Sect. 3.

thumbnail 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 , 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.

thumbnail 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 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: N2H+, N2D+, HC18O+, 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 ζ2 (corresponding to the high model of 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 N2D+ abundance, while it overestimates by a factor of 3−4 the intensities of HC18O+ and DCO+ lines.

The model from PI18 is the one that provides the overall best agreement with the observations. It is remarkable that, with the exception of N2D+ – 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 , and requires corrections in the range fcorr = 1.5−5 at t = 106 yr to reproduce the observed rotational lines. Changing the initial conditions by simulating first a cloud-like evolution of the chemistry for tCE = 106 yr does not improve the agreement. Table 2 contains a summary of the fcorr 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 N2D+ 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 ζ2, or more to its average value across the core, we ran another chemical model using ζ2 = 3 × 10−17 s−1 (the average ζ2 of model at column densities of N = [0.2−5] × 1022 cm−2). The estimated abundance profiles are very similar to model , especially atradii r > 0.05 pc = 104 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 ζ2. We tested also the attenuation model from Neufeld & Wolfire (2017), which predicts a steeper decrease of ζ2 with increasing density (ζ2N−0.92). This model is in worse agreement with the observations with respect to model , since it predicts CR ionisation rate values lower than model 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 (103 −104 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 pre- and proto-stellarcores. In this respect, we highlight that the two species that present the strongest sensitivity in their simulated spectra to changing ζ2 are HC18O+ and N2D+. The peak intensity of the HC18O+ (1-0) transition, in fact, changes by a factor of ≈10 going from model to . Deuterated diazenylium is also sensitive to high ζ2 values, which makes the abundance of this molecule drop (see for instance model ). This, combined with the fact that N2D+ 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 C18O and H13CO+ 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 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 ζ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.

thumbnail Fig. 11

Top panels: abundance profiles obtained from the model (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 and from themodel with constant ζ2 = 3 × 10−17 s−1. The colour code is the same as for the upper panels.

Table 2

Summary of the fcorr values used in the different chemical models for each molecular species.

Appendix A CR ionisation rate and temperature of the models

thumbnail Fig. A.1

Abundance profiles of N2H+, N2D+, HCO+, and DCO+ (from top-left to bottom-right panel) in the model , 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.

thumbnail Fig. A.2

Same as Fig. A.1, but for model .

As mentioned in the main text, increasing (or decreasing) ζ2 affects the gas temperature(TK). 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 ζ2 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 (, , and ), 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 N2H+ (11 − 13%) and N2D+ (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.

thumbnail Fig. A.3

Same as Fig. A.1, but for model .

Appendix B Complete set of abundance profiles for models , , and

In Fig. B.1 we report the abundance profiles for the four targeted species obtained with model , , and at the time steps not presented in Fig. 2.

thumbnail 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 with ζ2 = 1.3 × 10−17 s−1 is shown in blue, whilst the red and green curves show the models and of PI18, respectively.

Appendix C Uncertainties in the estimation of the fcorr 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 fcorr 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 fcorr 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 , we test five fcorr 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: (C.1)

where or 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 Rfcorr relation. The 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 .

Figure C.1 shows the results of this analysis. The values are shown as black vertical lines, with the shaded area indicating the corresponding ± 1σ intervals. The fcorr values found instead with the approach of Sect. 4 are shown as red vertical lines. For N2D+, HC18O+, and DCO+ the differences between the two approaches are below 10% and well below the uncertainty level. Only for N2H+ the two methods are not consistent. However, the discrepancy is only at a 1.2σ level. Furthermore, the value found () 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.

thumbnail 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 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 Rfcorr relation. In each panel, the vertical black line shows the 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 ζ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 or model .

thumbnail 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 ζ2 = 3× 10−17 s−1: resulting line profiles.

Figure E.1 shows the results of the radiative transfer applied to the chemical model with constant ζ2 = 3× 10−17 s−1.

thumbnail 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 . 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

  1. Bizzocchi, L., Caselli, P., Leonardo, E., & Dore, L. 2013, A&A, 555, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Black, J. H., Hartquist, T. W., & Dalgarno, A. 1978, ApJ, 224, 448 [Google Scholar]
  3. Bovino, S., Ferrada-Chamorro, S., Lupi, A., Schleicher, D. R. G., & Caselli, P. 2020, MNRAS, 495, L7 [Google Scholar]
  4. Caselli, P., Walmsley, C. M., Terzieva, R., & Herbst, E. 1998, ApJ, 499, 234 [Google Scholar]
  5. Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002, ApJ, 565, 344 [Google Scholar]
  6. Caselli, P., Pineda, J. E., Zhao, B., et al. 2019, ApJ, 874, 89 [NASA ADS] [CrossRef] [Google Scholar]
  7. Chacón-Tanarro, A., Caselli, P., Bizzocchi, L., et al. 2017, A&A, 606, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Crapsi, A., Caselli, P., Walmsley, M. C., & Tafalla, M. 2007, A&A, 470, 221 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Cummings, A. C., Stone, E. C., Heikkila, B. C., et al. 2016, ApJ, 831, 18 [CrossRef] [Google Scholar]
  10. Dalgarno, A., & Lepp, S. 1984, ApJ, 287, L47 [NASA ADS] [CrossRef] [Google Scholar]
  11. Doty, S. D., Schöier, F. L., & van Dishoeck, E. F. 2004, A&A, 418, 1021 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Glassgold, A. E., & Langer, W. D. 1974, ApJ, 193, 73 [NASA ADS] [CrossRef] [Google Scholar]
  13. Indriolo, N., & McCall, B. J. 2012, ApJ, 745, 91 [NASA ADS] [CrossRef] [Google Scholar]
  14. Indriolo, N., Geballe, T. R., Oka, T., & McCall, B. J. 2007, ApJ, 671, 1736 [NASA ADS] [CrossRef] [Google Scholar]
  15. Ivlev, A. V., Padovani, M., Galli, D., & Caselli, P. 2015, ApJ, 812, 135 [Google Scholar]
  16. Ivlev, A. V., Silsbee, K., Sipilä, O., & Caselli, P. 2019, ApJ, 884, 176 [Google Scholar]
  17. Ivlev, A. V., Silsbee, K., Padovani, M., & Galli, D. 2021, ApJ, 909, 107 [NASA ADS] [CrossRef] [Google Scholar]
  18. Keto, E. R. 1990, ApJ, 355, 190 [CrossRef] [Google Scholar]
  19. Keto, E., Rybicki, G. B., Bergin, E. A., & Plume, R. 2004, ApJ, 613, 355 [NASA ADS] [CrossRef] [Google Scholar]
  20. Keto, E., Caselli, P., & Rawlings, J. 2015, MNRAS, 446, 3731 [NASA ADS] [CrossRef] [Google Scholar]
  21. Kim, Y.-K., Santos, J. P., & Parente, F. 2000, Phys. Rev. A, 62, 052710 [NASA ADS] [CrossRef] [Google Scholar]
  22. Maret, S., & Bergin, E. A. 2007, ApJ, 664, 956 [Google Scholar]
  23. Maret, S., Bergin, E. A., & Tafalla, M. 2013, A&A, 559, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [Google Scholar]
  25. Mouschovias, T. C., & Spitzer, L., J. 1976, ApJ, 210, 326 [NASA ADS] [CrossRef] [Google Scholar]
  26. Neufeld, D. A., & Wolfire, M. G. 2017, ApJ, 845, 163 [Google Scholar]
  27. Padovani, M., Galli, D., & Glassgold, A. E. 2009, A&A, 501, 619 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Padovani, M., Hennebelle, P., & Galli, D. 2013, A&A, 560, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Padovani, M., Ivlev, A. V., Galli, D., & Caselli, P. 2018, A&A, 614, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Phan, V. H. M., Schulze, F., Mertsch, P., Recchia, S., & Gabici, S. 2021, Phys. Rev. Lett., 127, 141101 [NASA ADS] [CrossRef] [Google Scholar]
  31. Planck Collaboration Int. XXXV. 2016, A&A, 586, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Redaelli, E., Bizzocchi, L., Caselli, P., et al. 2018, A&A, 617, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Redaelli, E., Bizzocchi, L., Caselli, P., et al. 2019, A&A, 629, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Rudd, M. E., Kim, Y. K., Madison, D. H., & Gay, T. J. 1992, Rev. Mod. Phys., 64, 441 [Google Scholar]
  35. Ruffle, D. P., Bettens, R. P. A., Terzieva, R., & Herbst, E. 1999, ApJ, 523, 678 [NASA ADS] [CrossRef] [Google Scholar]
  36. Sabatini, G., Bovino, S., Giannetti, A., et al. 2020, A&A, 644, A34 [EDP Sciences] [Google Scholar]
  37. Scherer, K., Fichtner, H., Ferreira, S. E. S., Büsching, I., & Potgieter, M. S. 2008, ApJ, 680, L105 [NASA ADS] [CrossRef] [Google Scholar]
  38. Schlafly, E. F., Green, G., Finkbeiner, D. P., et al. 2014, ApJ, 786, 29 [NASA ADS] [CrossRef] [Google Scholar]
  39. Shirley, Y. L. 2015, PASP, 127, 299 [Google Scholar]
  40. Silsbee, K., Ivlev, A. V., Padovani, M., & Caselli, P. 2018, ApJ, 863, 188 [Google Scholar]
  41. Silsbee, K., Ivlev, A. V., Sipilä, O., Caselli, P., & Zhao, B. 2020, A&A, 641, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Sipilä, O., Caselli, P., & Harju, J. 2015a, A&A, 578, A55 [Google Scholar]
  43. Sipilä, O., Harju, J., Caselli, P., & Schlemmer, S. 2015b, A&A, 581, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Sipilä, O., Caselli, P., & Harju, J. 2019, A&A, 631, A63 [Google Scholar]
  45. Spezzano, S., Bizzocchi, L., Caselli, P., Harju, J., & Brünken, S. 2016, A&A, 592, L11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Spitzer, L. J., & Tomasko, M. G. 1968, ApJ, 152, 971 [NASA ADS] [CrossRef] [Google Scholar]
  47. Stone, E. C., Cummings, A. C., Heikkila, B. C., & Lal, N. 2019, Nat. Astron., 3, 1013 [NASA ADS] [CrossRef] [Google Scholar]
  48. van der Tak, F. F. S., & van Dishoeck, E. F. 2000, A&A, 358, L79 [NASA ADS] [Google Scholar]
  49. Ward-Thompson, D., Motte, F., & Andre, P. 1999, MNRAS, 305, 143 [Google Scholar]
  50. Wilson, T. L. 1999, Rep. Prog. Phys., 62, 143 [Google Scholar]

1

RB19 tested AV values of 1, 2, and 5 mag, and showed that this parameter does not affect significantly the synthetic spectra of the targeted species, in particular of N2H+ isotopologues.

2

We highlight that multiplying the abundance profile by fcorr does not correspond to multiplying the fluxes modelled by radiative transfer by the same factor, due to opacity and non-LTE effects.

3

The relation derived from Neufeld & Wolfire (2017) reads log10 (ζ2∕s-1) = −0.92 × log10(N∕cm−2) + 3.9.

All Tables

Table 1

Initial abundances with respect to the total number density of hydrogen nuclei used in the chemical model.

Table 2

Summary of the fcorr values used in the different chemical models for each molecular species.

All Figures

thumbnail 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 (dashed bluecurve), (solid red), and (dashed-dotted green) models. The profiles for the models and 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
thumbnail 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 (blue), model (red), and model (green). Bottom panels: the green and red solid curves show the ratio of model and model to model , 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 , and fcorr = 1.5 in model ; therefore thered solid curve / in the rightmost panel has to be compared with the red dashed line at .

In the text
thumbnail 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: in the left panels, in the central panels, and 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
thumbnail Fig. 4

As Fig. 3, but for N2D+.

In the text
thumbnail Fig. 5

As Fig. 3, but for HC18O+.

In the text
thumbnail Fig. 6

As Fig. 3, but for DCO+.

In the text
thumbnail Fig. 7

Abundance profiles for N2D+ in model 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
thumbnail Fig. 8

x(e) profiles in model (dashed bluecurve), (solid red), and (dashed-dotted green). All the models are shown at t = 106 yr.

In the text
thumbnail 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 , 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
thumbnail 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 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
thumbnail Fig. 11

Top panels: abundance profiles obtained from the model (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 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
thumbnail Fig. A.1

Abundance profiles of N2H+, N2D+, HCO+, and DCO+ (from top-left to bottom-right panel) in the model , 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
thumbnail Fig. A.2

Same as Fig. A.1, but for model .

In the text
thumbnail Fig. A.3

Same as Fig. A.1, but for model .

In the text
thumbnail 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 with ζ2 = 1.3 × 10−17 s−1 is shown in blue, whilst the red and green curves show the models and of PI18, respectively.

In the text
thumbnail 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 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 Rfcorr relation. In each panel, the vertical black line shows the 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
thumbnail 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
thumbnail 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 . 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

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.