© E. Redaelli et al. 2019

1 Introduction

Deuterium was formed during the primordial nucleosynthesis in the first minutes after the Big Bang with a fractional abundance of ≈1.5 × 10⁻⁵ (Linsky 2003). However, enhancements of this value of several orders of magnitude have been found both in different environments of the interstellar medium (ISM) and in several components of the solar system (Ceccarelli et al. 2014, and references therein). This has led to the idea of using the deuteration ratio D/H (the ratio between the abundance of a D-bearing molecule and its hydrogenated isotopologue) as a diagnostic tool to investigate the star formation process, and ultimately to understand how our own solar system was formed (Ceccarelli et al. 2014).

Prestellar cores represent the early phases of low-mass star formation (André et al. 2014). These are dense (central volume densities n ≳ 10⁵ cm⁻³) and cold (T ≲ 10 K) fragments of molecular clouds that are on the verge of gravitational collapse but still have not formed any central object (Keto & Caselli 2008). Prestellar cores are known to exhibit some of the highest levels of deuteration (D∕H > 0.1, Crapsi et al. 2005; Pagani et al. 2007), because their physical conditions greatly favour the deuteration processes. D-fractionation (i.e. the inclusion of a D-atom in hydrogenated species) is indeed driven by the reaction: $${\text{H}}_3^{+}+\text{HD}\leftrightharpoons {\text{H}}_2{\text{D}}^{+}+{\text{H}}_2,$$(1)

which proceeds more efficiently from left to right as the temperature decreases (assuming a low ortho-to-para H₂ ratio; e.g. Pagani et al. 1992) due to the lower zero-point energy of H₂D⁺, the progenitor of all the deuterated species in regions where the temperature is ≲ 30 K (Dalgarno & Lepp 1984). Thus, in cold environments, the increased abundance of H₂D⁺ and of the doubly and triply deuterated forms of ${\text{H}}_3^{+}$ rapidly leads to higher D/H ratios in all molecules produced by the reactions with the various H${}_3^{+}$ isotopologues. This process is also favoured by another mechanism. The ${\text{H}}_3^{+}$ main destruction route is its reaction with CO. Above n ≈a few 10⁴ cm⁻³, CO and other C- and O-bearing species catastrophically freeze onto the dust grain surfaces, leading to their depletion from the gas phase; in regions where the abundance of neutral species reacting with H${}_3^{+}$ is reduced, D-fractionation is enhanced (Dalgarno & Lepp 1984). In prestellar cores, where the depletion fraction of CO is > 90% in the central 6000–8000 AU (Caselli et al. 1999; Bacmann et al. 2002), the D-fraction in NH₃ reaches levels as high as 40% (Crapsi et al. 2005).

A quantity that plays a crucial role in driving the timescales of star formation is the ionisation degree, that is the ratio between the free electron density n(e) and the gas density n(H₂) [x(e) = n(e)∕n(H₂)] (McKee 1989). This factor determines the timescale for the ambipolar diffusion of neutrals across the magnetic field lines, which is the main process that allows the gravitational collapse of magnetised cores (Mouschovias & Spitzer 1976). In molecular clouds, the deuteration fraction is in general considered a good indicator of the ionisation degree (Guelin et al. 1977; Caselli et al. 1998, 2002a; Dalgarno 2006), as the D-fraction in species such as HCO⁺ is a function of x(e). As a consequence, accurate measurements of the D/H ratio allow the ionisation degree in prestellar cores to be determined. In particular, N₂H⁺ and HCO⁺ represent ideal tracers (Caselli et al. 1998; Caselli 2002); these molecules – and their deuterated counterparts – are abundant in low-mass cores. It has been shown that N₂H⁺ is usually highly concentrated in the central parts, whereas HCO⁺ traces the outer layers (e.g. Punanova et al. 2018). Combined, they can provide information on the ionisation degree in the whole source.

L1544 is one of the most studied low-mass prestellar cores. Its vicinity (d ≈ 135 pc, Schlafly et al. 2014)and bright continuum emission at mm and sub-mm wavelengths (Ward-Thompson et al. 1999; Doty et al. 2004; Spezzano et al. 2016)make it the ideal environment to study the early stages of low-mass star formation. The central part of the core is cold (T ≈6 K, Crapsi et al. 2007) and dense (n > 10⁶ cm⁻³, Tafalla et al. 2002; Keto & Caselli 2010). It is centrally peaked and shows signatures of contraction (Caselli et al. 2002b, 2012; Lee et al. 2001), suggesting that it is on the verge of gravitational collapse. In the last 20 yr, a vast observational campaign combined with theoretical efforts has led to a deep knowledge of its physical and chemical structure (see e.g. Crapsi et al. 2007; Caselli et al. 2012; Keto et al. 2015; Spezzano et al. 2017). Previous works (Caselli et al. 1998, 2002a) have already investigated the deuteration level and ionisation fraction in this source, using the telescope capabilities available at that time. We now have access to new high-sensitivity maps of N₂H⁺, N₂D⁺, HC¹⁸O⁺, and DCO⁺ across the whole core performed with the on-the-fly method. For the first time, we can analyse the first three rotational transitions of N₂D⁺ and DCO⁺ and the N₂H⁺ (1–0) and (3–2) lines¹ simultaneously. In this paper we report the deuteration fraction maps of N₂H⁺ and HCO⁺ in L1544, obtained with a detailed non-LTE modelling of the molecular spectra at the centre of the core. In a following work (Paper II) we will use the results presented here to accurately derive the ionisation fraction map of L1544.

In Sect. 2 we present the observations. The analysis is described in Sect. 3, divided into chemical modelling (Sect. 3.1), non-LTE modelling (Sect. 3.2), and column density derivation (Sect. 3.3). We discuss the results in Sect. 4, and draw a summary of the work in Sect. 5.

2 Observations

All the data presented in this paper were collected with the IRAM 30 m single dish telescope. We used the four receivers of the frontend EMIR (E0, E1, E2 and E3), which allow us to cover the required frequency range. The N₂H⁺ (1–0) data were acquired in summer 2015. Observations of N₂H⁺ (3–2), N₂D⁺ (1–0), (2–1), and (3–2), HC¹⁸O⁺ (1–0), and DCO⁺ (1–0), (2–1), and (3–2) were performed during winter 2015/2016 using four different spectral setups. The VESPA backend, with a spectral resolution of 20 kHz, was used for all data except for the DCO⁺ (1–0) transition, which was acquired with the FTS backend (resolution: 50 kHz). In all the sessions, the pointing was frequently checked on nearby planets or bright quasars, and found to be accurate within 2″ and 5″ (at the high and low frequencies, respectively). In order to encompass the whole core, all the transitions were acquired in the on-the-fly (OTF) map mode with a 2′ × 2′ (≈ 0.08 × 0.08 pc²) footprint centred at the L1544 dust peak (RA (J2000) = 05^h04′17.′′21, Dec (J2000) = 25°10′42.′′ 8, Ward-Thompson et al. 1999). The N₂H⁺ (1–0) data, observed in a different session, consist of a 4′ × 4′ (≈ 0.16 × 0.16 pc²) map across the dust peak.

The spectra, originally calibrated to antenna temperature ($T_{\text{A}}^{*}$), were converted to main beam temperature(T_MB) using the tabulated values for the 30m forward efficiency (F_eff) and main beam efficiency (B_eff). The main features of the observed lines, such as their frequency, angular resolution (θ_beam), and spectral resolution (ΔV_ch), are summarised in Table 1.

In addition to these data, we also make use of a single pointing observation of the HCO⁺ (1–0) line performed with the FCRAO 14 m telescope towards the dust peak of L1544. These data have a spectral resolution of 0.07 km s⁻¹ and an angular resolution of 55″, and were previously published in Tafalla et al. (1998). We refer to their paper for further details.

Results

Figures 1 and 2 show the observed spectra at the millimetre dust peak for each transition (in black). Considering only the emission-free channels, we computed the mean rms per channel for each transition, obtaining: 0.05 K for N₂D⁺ (1–0), 0.18 K for N₂D⁺ (2–1), 0.26 K for N₂D⁺ (3–2), 0.14 K for N₂H⁺ (1–0), 0.32 K for N₂H⁺ (3–2), 0.04 K for HC¹⁸O⁺ (1–0), 0.08 K for DCO⁺ (1–0), 0.07 K for DCO⁺ (2–1), and 0.23 K for DCO⁺ (3–2). Thus, at the dust peak the observed spectra have a signal-to-noise ratio (S/N) ranging from ≈ 4.0 to ≈ 40.0. Figure 3 shows the integrated intensity maps for each transition, computed integrating over all the hyperfine components when present.

Fig. 1 Observed spectra at the dust peak of N2D+ and DCO+, with the original angular resolution (in black histograms). The transition is labelled in the upper-left corner of each panel. In red, the best-fit solution found with MOLLIE is shown overlaying the observations (see Sect. 3.2). Underneath each line we present the residuals of the fit (difference between the model and the observation). The vertical blue bars show which hyperfine component has been used to derive the molecular column density, when the hyperfine structure has not been neglected.

3 Analysis

Our main goal is to derive the deuterium fraction across the source, and for this aim we need reliable measurements of the column density for each isotopologue. In principle, the knowledge of multiple transitions of a species can allow both its column density and excitation condition (i.e. the line excitation temperature, T_ex) to be determined. This approach, implemented in the population diagram method (Goldsmith & Langer 1999), has been extensively used in the literature, but it is valid only in the assumption that all the lines share the same T_ex. This assumption extends to all the hyperfine components, when present. However, we know from previous studies (Daniel et al. 2006, 2013) that this is often not the case for the molecules considered here. There are multiple possible reasons for this, such as selective photon trapping in very crowded hyperfine structure, or a significant difference incritical density for the different rotational transitions of the same molecule (the Einstein coefficient A scales with the third power of the frequency). In this case, the column density and the excitation temperatures of the different transitions are degenerate parameters, and cannot be constrained simultaneously. Therefore, an approach that does not assume LTE is needed to determine the T_ex of each transition and thus the molecular column density.

A full non-LTE analysis is on the other hand not feasible on maps, because it requires the knowledge of the full three-dimensional temperature and density structure of the core. For L1544, however, Keto et al. (2015) developed a one-dimensional physical model based on a collapsing Bonnor-Ebert sphere that provides a good fit to several molecules at the dust peak (Caselli et al. 2012, 2017; Bizzocchi et al. 2013; Redaelli et al. 2018). We thus decided to model the spectra at the dust peak of the core with a full non-LTE analysis. In this way, we are able to constrain the column density of each molecule at the dust peak. At the same time, we can use the equations of radiative transfer to derive N_col in the whole map. Since these two approaches must return the same result at the centre of the core, the T_ex (which is involved in the radiative transfer equations) is treated as a free parameter to be tuned in order to obtain the expected agreement of the two methods. The non-LTE modelling allows us to break the degeneracy between N_col and T_ex mentioned in the previous paragraph. In the following sections we provide a detailed description of the required steps.

Fig. 2 As in Fig. 1, but for N2H+, HC18O+ and HCO+ transitions.

3.1 Chemical models

The non-LTE modelling of line emission requires the knowledge of molecular abundance (X_mol) profiles. To derive these, we used the gas-grain chemical model discussed in Sipilä et al. (2015a,b), which has recently received an extensive update (Sipilä et al. 2019). Briefly, we employ the static physical core model from Keto et al. (2015), divided into concentric shells, and solve the chemical evolution separately in each shell. We then combine the results at different time steps to produce radius-dependent abundance profiles for the molecules of interest. The adopted initial conditions are summarised in Table 2. More information on typical model parameters can be found for instance in Sipilä et al. (2015a).

In this work we considered five different time steps: t = 10⁴ yr, t = 5 × 10⁴ yr, t = 10⁵ yr, t = 5 × 10⁵ yr, and t = 10⁶ yr. We tested the influence of external UV radiation by considering three different visual extinction values for the cloud embedding the core: A_V = 1, 2, or 5 mag. This amounts to a total of 15 modelled abundance profiles per molecule. We show the complete set of profiles in Appendix B. Simulated lines were produced based on these profiles as described in the following section.

3.2 Non-LTE modelling at the dust peak

The non-LTE approach was implemented using the radiative transfer code MOLLIE (Keto 1990; Keto et al. 2004), which allows us to derive synthetic spectra from a given physical model. MOLLIE is able to treat the line overlap that occurs in a crowded hyperfine structure, such as those of N₂H⁺ and N₂D⁺. As an output, it also provides the molecular column density, obtained by integrating the abundance profile multiplied by the gas density, eventually convolved to a desired beam size.

The physical model of the source, as already mentioned, is taken from Keto et al. (2015). It is an unstable quasi-equilibrium Bonnor-Ebert sphere, with a central density n₀ ≈ 10⁷ cm⁻³ and a central temperature of 6.5 and 6.3 K for the gas and the dust, respectively. Figure 4 shows the profiles of the main quantities for this model. Recently, Chacon-Tanarro et al. (2019) proposed a new density and (dust) temperature model for L1544, analysing new millimetre continuum observations of L1544. However, the authors did not develop a new model for the velocity field, a crucial quantity to correctly reproduce molecular spectra. Therefore, we prefer to use the previous model, which also includes the velocity profile.

For the line radiative transfer, the collisional coefficients are of great importance, and we used the most recent ones. When they are available only for one isotopologue, we scaled them taking into account the different reduced mass of the collision system. The collisional coefficients for N₂H⁺ are taken from Lique et al. (2015), who computed them for the N₂H⁺/p-H₂ system. The rates for the DCO⁺/p-H₂ system were computed by Pagani et al. (2012), who take into account the hyperfine structure given by the deuterium nucleus. Finally, the collisions of HCO⁺ with p-H₂ were published by Flower (1999), from which we scaled the coefficients for HC¹⁸O⁺. In all cases, the o-H₂ is neglected,which is a reasonable assumption given that the ortho-to-para ratio (OPR) in L1544 is found to be as low as OPR = 10⁻³ (Kong et al. 2015).

We tested all the 15 models described in Sect. 3.1. The abundance profiles produced by the chemical model often systematically overestimate or underestimate the observed fluxes, and need to be multiplied by a numerical factor to find a good agreement with the observations. However, MOLLIE requires too much computational time to allow a full sampling of the parameter space. For each molecule, we therefore proceeded as follows: first, we test all the original models, varying the evolutionarystage and the thickness of the surrounding cloud, and we select the one that best reproduces the trend of the different transitions (e.g. the ratio between the peak temperatures of the brightest component of each line). We then multiply the whole profile by a numerical factor, testing approximately five values until the observed fluxes are matched. This is considered the best fit to the observations. We want to highlight that our goal is to obtain a simultaneous fit for all molecules that is roughly consistent with the observations, and a full parameter-space search for the best possible fit within the model uncertainties is beyond the scope of the present paper. In the following text, we focus on the abundance profiles that provide the best-fit solution for each molecule. As an example, wepresent in Sect. 3.2.2 the resulting spectra from several of the 15 abundance profiles for N₂D⁺.

Fig. 3 Integrated intensity maps of all the observed transitions, which are labelled in the bottom-right corners. The contours show 20, 40, 60, 80, and 90% of the peak values, which are in order from top left to bottom right panels: 2.0, 2.0, 0.7, 4.6, 0.7, 0.15, 1.8, 1.1, 0.7 K km s−1. The black cross represents the dust peak position, and the white circle is the beam size. All the maps have the same size (2′ × 2′), with the exception of N2H+ (1–0), which is 4′ × 4′: in this map the white rectangle shows the smaller coverage of the other transitions. Scalebars are shown in the top-left corners of the N2D+ (1–0) and N2H+ (1–0) maps.

Fig. 4 Profiles of dust temperature (green), gas temperature (red), H2 volume density (black, in logarithmic scale), and velocity (blue, in units of 0.1 km s−1) for the L1544 model developed in Keto et al. (2015). The velocity in the model is negative (it represents contraction motions), but it is shown here as positive for improved readability.

3.2.1 Best-fit solution for the observed molecules at the dust peak

The best-fit solutions for each molecule are shown in red in Figs. 1 and 2 overlaying the observations. We also include here the analysis of a single-pointing observation of HCO⁺ (1–0) performed at the dust peak (see Sect. 2 for details). Despite the fact that a map for this transition is not available, it is important to check the behaviour of the chemical model for HCO⁺, the main isotopologue of its family. The observed spectrum of HCO⁺, together with the best fit from MOLLIE, is shown in the bottom-right panel of Fig. 2.

N₂H⁺

The chemical network systematically underestimates the emission of N₂H⁺, for any combination of chemical time and external visual extinction. We find that in general the spectra are best reproduced at late evolutionary times (10⁶ yr), when the molecules at the centre of the core start to be affected by depletion. This allows the observed ratio between the (1–0) and the (3–2) line to be reproduced. Our simulations show that varying the thickness of the external layer does not have a significant impact on the synthetic spectra due to the high critical density of N₂H⁺ whose emission arises only from the central regions (n_crit ≳ 10⁵ cm⁻³; see Shirley 2015). Furthermore, N₂ is a late-type molecule that can take about ten times longer than CO to form (e.g. Hily-Blant et al. 2010). Therefore, the outer less dense envelope may not be rich enough in N₂ for N₂H⁺ to form efficiently. The abundance profile for t = 10⁶ yr and A_V = 1 had to be multiplied by 2.7 to match the observed fluxes.

N₂D⁺

Similar considerations to the ones made for the main isotopologue hold also for N₂D⁺. The best fit is found for t = 10⁶ yr and an external visual extinction of A_V = 1 mag, even thoughthis last parameter does not significantly affect the simulated spectra (see Sect. 3.2.2 for more details). The models again tend to underestimate the observed fluxes, and we had to increase the abundance profile by a factor of 3.0 to obtain a good agreement.

HCO⁺

In general, the original abundance profiles coming from the chemical network result in a better agreement with the observations of the HCO⁺ isotopologues with respect to the N₂H⁺ ones. The best agreement for the main isotopologue is found for the model with t = 10⁶ yr and an external visual extinction ofA_V = 5 mag, multiplied by a factor 0.5. This model presents a high molecular abundance in the external layers of the cores (X_mol ≳ 10⁻⁷), which is important to correctly reproduce the strong blue asymmetry of the (1–0) line due to self-absorption (Tafalla et al. 1998).

HC¹⁸O⁺

In the chemical model, the abundance profile for the HC¹⁸O⁺ is derived from that of HCO⁺ using the scaling factor valid for the local ISM ¹⁶O∕¹⁸O = 557 (Wilson 1999). This choice however leads to an underestimation of the observed flux. The best fit is found with the model at t = 10⁶ yr and A_V = 1 mag, multiplied by 6.0. This late evolutionary stage is important to be able to reproduce the double-peak feature seen in the (1–0) spectrum, but we had to decrease the external visual extinction with respect to the main isotopologue. If the abundance of HC¹⁸O⁺ is too high in the external layers of the core, as in the models with A_V = 5 mag, the resulting synthetic spectra present a blue asymmetry due to the self-absorption, which is not seen in the observations.

DCO⁺

The three transitions of DCO⁺ are best reproduced with the model at t = 10⁶ yr and A_V = 1 mag. In this case, we decrease the original abundance profile by 50% to match the observed peak temperatures.

The abundance profiles that provide the best-fit solution for each species are shown in Fig. 5, after the multiplication for the corrective factors. The adopted models and scaling factors, together with the molecular column density computed by MOLLIE at the centre of the core ($N_{\text{col}}^{\text{MOLLIE}}$), are summarised in Table 3. It is not surprising that the chemical model results do not provide a perfect match to the observations, and that the various lines cannot be fitted with a single model. There are many sources of uncertainty, such as the source model (density and gas/dust temperature structures), chemical modelling parameters, and rate coefficients. Our model does not consider the self-shielding of CO, which would affect the HCO⁺ abundance in the outer core, for example. A detailed parameter-space exploration of the abundance profiles is beyond the scope of the present paper.

3.2.2 N₂D⁺ depletion

In this section, we analyse in detail the case of N₂D⁺, showing the non-LTE radiative transfer results for 6 out of the 15 abundance profiles that we tested at the dust peak. The abundance profiles are all multiplied by a factor of 3.0 to allow an easier comparison with the observations. The obtained spectra are shown in Fig. 6 overlaying the observations.

The first four rows from the top in Fig. 6 show increasing evolutionary times, keeping the external extinction fixed to A_V = 1 mag. The difference between the observations and the model can be quantified by comparing the relative ratio of the peak temperatures of the three lines. In the observations, this ratio is T_peak(1–0) : T_peak(2–1) : T_peak(3–2) = 1:1.5:1.1. In the models with the earlier times (t = 5 × 10⁴ and t = 10⁵ yr), the obtained ratios are 1:3.3:2.5 and 1:2.6:1.9, respectively. This means that these abundance profiles overestimate the (2–1) and (3–2) with respect to the (1–0) line. In the later evolutionary stages (t =5 × 10⁵ and t = 10⁶ yr), the ratios are 1:1.2:0.7 and 1:1.3:0.8, closer to theobserved one. This trend can be understood by looking at the behaviour of the different abundance profiles in the inner 10 000 AU, shown in Fig. 7. At earlier times, the N₂D⁺ abundance peaks within the central 2000–3000 AU, where the molecular hydrogen density is n(H₂) ≳ 10⁵ cm⁻³. Since the critical densities of the (1–0), (2–1), and (3-2) transitions are n_crit = 7.6 × 10⁴, 8.7 × 10⁵, and 3.8 × 10⁶ cm⁻³, respectively (Shirley 2015), this is also the region where the (2–1) and the (3–2) lines are emitting the most. At later evolutionary times, N₂D⁺ is depleted and the abundance decreases by more than a factor of five for r ≲ 2000 AU. This is due to the fact that its precursor species, N₂, starts to freeze-out onto dust grains, and as a consequence the formation rate of N₂D⁺ is reduced. The N₂D⁺ abundance peak moves outwards to r ≈ 5000 AU, where n(H₂) ≈ 5 × 10⁴ cm⁻³. Therefore, the (2–1) and (3–2) lines are less excited, while the (1–0) flux keeps rising. This is the first observational confirmation of N₂D⁺ partial depletion in cold and dense environments.

We now focus on the influence of the external embedding layer, keeping the evolutionary stage fixed (t = 10⁶ yr) and varying the external visual extinction (see the lower three rows of Fig. 6). The synthetic spectra are not very sensitive to this change, both in the flux level and in the spectral profiles, especially concerning the (2–1) and (3–2) lines. This is again due to a combination of molecular abundance and excitation conditions. In fact, increasing A_V significantly changes the abundance only for r ≳ 10 000 AU, where the molecular hydrogen density is n(H₂) ≲ 10⁴ cm⁻³. Therefore, the high-frequency transitions are not significantly affected by the change of A_V.

Fig. 5 Left panel: molecular abundances that provide the best fit to the spectra at the dust peak for the different species as a function of the core radius. Right panel: zoom-in of the molecular abundances in the inner 0.14 pc.

3.3 Column density and deuteration maps

In the following, we describe how we derive the column density map of each species, knowing the N_col value at the centre of the core from the non-LTE modelling, via the radiative transfer equations. The spectra of the analysed molecules greatly differ due to their opacity properties and hyperfine structure. For instance, we have optically thin transitions (HC¹⁸O⁺ (1–0)), or optically thick ones (HCO⁺ (1–0)). Some species present crowded hyperfine structures (such as N₂H⁺ and N₂D⁺), while others show a much simpler profile (e.g. DCO⁺, where the hyperfine splitting due to the deuterium spin is not detectable at the given spectral resolution; Caselli & Dore 2005). As a result, we tackle the problem of deriving the column density maps by developing a specific approach for each species individually. In all cases, we make use of the following well-known relations. $$T_{\text{MB}}={\eta }_{\text{bf}}\left[J_{\nu }(T_{\text{ex}})-J_{\nu }(T_{\text{bg}})\right]\left(1-{\text{e}}^{-{\tau }_{\nu }}\right),$$(2)

and $${\tau }_{\nu }=\sqrt{\frac{\ln 2}{16{\pi }^3}}\frac{c^3{A}_{\text{ul}}{g}_{\text{u}}}{{\nu }^{3}Q(T_{\text{ex}})\Delta V}{\text{e}}^{-\frac{E_{\text{u}}}{k_{\text{B}}{T}_{\text{ex}}}}\left({\text{e}}^{\frac{h\nu }{k_{\text{B}}{T}_{\text{ex}}}}-1\right)N_{\text{col}}.$$(3)

Equation (2) describes the main beam temperature observed for a line with a given optical depth τ_ν and a given excitation temperature T_ex. The function J_ν is the equivalent Rayleigh-Jeans temperature, T_bg = 2.73 K is the cosmic background temperature, and η_bf is the beam filling factor, considered to be equal to 1.0 due to the large extension of the core emission compared to the beam size. The second equation relates the optical depth of a given transition to the total column density of the molecule, N_col. The line is considered to be Gaussian, with a full width at half maximum equal to ΔV. Here, h is the Planck constant and k_B the Boltzmann constant. Equation (3) depends on several spectroscopic constants: the line frequency ν, the upper state energy E_u and degenerancyg_u, the Einstein coefficient A_ul, and the partition function Q. This latter equation holds both for an entire rotational transition and for an individual hyperfine component, depending on the splitting considered for the levels in the calculation of the partition function and the corresponding level degeneracy used. We emphasise that Eqs. (2) and (3) represent the integrated form of the (differential) radiative transfer equations. This is the appropriate approach since observations yield essentially integrated quantities (such as column density, fluxes, etc.).

Fig. 6 Results of the line radiative transfer for six chemical models on the three transitions of N2D+, respectively (1–0), (2–1), and (3–2) from left to right in each panel. The synthetic spectra are shown with red histograms overlaying the observations (in black). Each model is labelled with its evolutionary stage and external visual extinction.

Fig. 7 Zoom-in of the central 10 000 AU of the N2D+ abundances predicted by the chemical code for AV = 1 mag at four evolutionary stages: t = 5 × 104 yr (green), t = 105 yr (blue), t = 5 × 105 yr (black), and t = 106 yr (red). The abundance profiles are multiplied by a factor of 3.0.

3.3.1 Column density maps of individual molecules

To derive the column density of each molecule, we first smooth the angular resolution of the high-frequency transitions to the one with the largest beam. We then fit N_col pixel by pixel, minimising the residuals between the values of T_MB obtained through Eqs. (3) and (2) and the observed peak temperatures for all the available transitions.The observed peak temperatures are obtained with a Gaussian fit to the observed spectra, deriving in this way both T_MB and ΔV at the same time. When the line presents hyperfine structure, we select a single component. The choice of the T_ex value to beused in Eqs. (2) and (3) is made so that this approach at the centre of the core agrees with the N_col value obtained with MOLLIE, within the uncertainties. The T_ex values obtained in this way are used to model the whole map. The chosen method to select T_ex for each line and the assumption that this quantity is spatially constant are discussed in detail in Appendix C. We now describe in detail on which spectral features we focus to fit N_col for each species.

N₂H⁺

The (1–0) transition of N₂H⁺ presents a single, isolated component at the high-frequency end of its spectrum, which is a perfect candidate for our analysis, also providing reliable values for the line width. On the contrary, the hyperfine structure of the (3–2) line collapses in a very complicated structure, with a strong shoulder on the red side. For this transition we thus use the peak temperature of the observed spectra combined with the ΔV value obtained from the analysis of the (1–0) transition². For the spectroscopic constants needed in Eq. (3), we selected the brightest hyperfine component according to the spectroscopic catalogues (see Appendix A).

N₂D⁺

The N₂D⁺ (1–0) transition also presents a single isolated component suitable to our purposes. This is not the case for the other two transitions. For the (2–1) line, we find that the best solution is to consider the brightest feature of the spectra, which is composed by four components separated in total by only ≈ 90 kHz (see Appendix A for details). We fit a single Gaussian to this spectral feature in the observations, and we use the additive property of Eq. (3), according to which ${\tau }_{\nu }^{\text{tot}}={\sum }_i{\tau }_{\nu }^i$, where i labels the hyperfine components. Finally, in the (3–2) spectra the hyperfine structure is closely blended, and therefore we consider it a single line³.

HC¹⁸O⁺

This is the only molecule for which we have a single transition, and thus we cannot fit N_col using the procedure described above. However, this transition is optically thin, as confirmed by the radiative transfer performed by MOLLIE, which predicts a maximum value of τ_ν ≈ 0.20 at the dust peak. Therefore, we use the optically thin approximation of Eqs. (2) and (3), following Caselli et al. (2002a): $$N_{\text{col}}=\frac{8\pi W{\nu }^3}{c^3{A}_{\text{ul}}}\frac{Q(T_{\text{ex}})}{J_{\nu }(T_{\text{ex}})-J_{\nu }(T_{\text{bg}})}\frac{{\text{e}}^{\frac{E_{\text{u}}}{k_{\text{B}}{T}_{\text{ex}}}}}{g_{\text{u}}\left({\text{e}}^{\frac{h\nu }{k_{\text{B}}{T}_{\text{ex}}}}-1\right)},$$(4)

where W is the integrated intensity of the line. We therefore fit a Gaussian to the spectrum in each pixel; we subsequently compute its integrated intensity and then derive N_col using Eq. (4).

DCO⁺

The spectra of DCO⁺ are the simplest to analyse since the hyperfine structure due to the deuterium atom is compact if compared to our spectral resolution. We consider the three transitions as single lines, fit a Gaussian profile to each of them, and fit the molecular column density minimising the residuals between the synthetic spectra and the Gaussian fits to the observations.

Figure 8 shows the column density maps obtained for each molecule as described immediately above, while their peak values are summarised in Table 4. Our results are consistent with those of Caselli et al. (2002a), which are quoted to be accurate within a factor of two. We however highlight that a direct comparison with other works is difficult (see Sect. 4). We report in Appendix A the spectroscopic values used for each transition. The values of T_ex that ensure that the above method provides the same peak column density as computed by MOLLIE are reported in the last column of Table 3. These values confirm that the assumption of a single T_ex value for all transitions is incorrect, since the (3–2) lines present excitation temperatures that are ≈ 2 K lower than the respective (1–0) ones. For most transitions, T_ex is lower than the gas kinetic temperature, suggesting that many of the detected lines are subthermally excited.

One strong assumption of our method is the line Gaussianity. This neglects any spectral asymmetry due for example to the contraction motions, which are clearly visible in some spectra. We want to emphasise that any improvement of our analysis is strongly related to the development of a better physical model for L1544. In fact, the one-dimensional model of Keto et al. (2015) allows us to solve the radiative transfer problem only at the dust peak, since the fact that L1544 is not spherically symmetric becomes critical moving away from the central point. Since our goal is to investigate the deuteration properties in the whole extended core, we are forced to make further assumptions and simplifications.

Once the column density maps are ready, we can compute the deuterium fraction, dividing the column density of the D-bearing species and that of the main isotopologue pixel by pixel. First, the beam sizes of the maps are matched to allow a fair comparison. The column density of HCO⁺ is derived from that of HC¹⁸O⁺ using the standard isotopic ratio: N_col(HCO⁺) = N_col(HC¹⁸O⁺) × 557. Figures 9 and 10 show the resulting images for the D/H ratio.

3.3.2 Evaluation of the uncertainties

In the approach adopted to compute the molecular column densities, the main source of uncertainty is found in the excitation temperature values. In fact, they depend on several other parameters, such as the assumed physical and chemical models. We choose a conservative approach, and decide that the T_ex values in Table 3 are accurate within 0.5 K. This choice is supported by the following considerations. First of all, such a variation in T_ex translates into an average variation of 20–30% in the computed column densities, which corresponds to an equal variation of the molecular abundance. We ran MOLLIE again, modifying the abundance profiles by 30%, and the resulting modelled spectra deviate significantly from the observed ones. This means that our data are indeed sensitive to a variation of 0.5 K in T_ex. Furthermore, this variation in T_ex is the limit suggested by our analysis on the possible spatial variation of this quantity (see Appendix C).

The column density maps, computed with the T_ex values modified by 0.5 K, therefore provide the sought-after uncertainties. The uncertainties on the D/H ratios are derived through standard error propagation from the ones on the molecular column density. The median relative errors are 41% for ${\text{D/H}}_{{\text{N}}_2{\text{H}}^{+}}$ and 22% for ${\text{D/H}}_{{\text{HCO}}^{+}}$. The lower mean uncertainty on the HCO⁺ results depends on how sensitive N_col is to T_ex variations, according to Eq. (3). The column density of N₂D⁺, for example, is more heavily affected than that of HC¹⁸O⁺ by a T_ex change of 0.5 K, and therefore its relative error is higher.

Fig. 8 Column density maps obtained for each molecule (labelled in the bottom-right corner of each panel). The beam sizes are shown in the bottom-left corners. The black cross represents the position of the millimetre dust peak.

4 Discussion

The column density maps presented in Fig. 8 show that different molecules exhibit different morphologies. In particular, one can see differences in the position of the column density maxima with respect to the dust peak. Both N_col(N₂H⁺) and N_col(N₂D⁺) peak close to the millimetre dust peak, and while the former has a more extended distribution, the latter is highly concentrated in the densest part of the core (as already noted in Caselli et al. 2002a). The abundance profiles in Fig. 5, despite referring to a one-dimensional core model, confirm these trends. The peak of N₂H⁺ abundance is found at a larger radius with respect to N₂D⁺.

Despite the fact that N₂D⁺ is concentrated in the central region, it is interesting to notice that the chemical model that provides the best fit to its transitions predicts some extent of molecular depletion due to freeze-out, as described in Sect. 3.2.2. This is to our knowledge the first time that N₂D⁺ depletion, which is often predicted by chemical models, is confirmed by observational data. We are aware that our spatial resolution (35″ corresponds to ≈0.03 pc at the distance of L1544) is not sufficiently high to investigate this point exhaustively. However, the fact that we need a model with some central depletion to reproduce the observed spectra is clear confirmation. As a further test, we tried to model all the N₂H⁺ and N₂D⁺ spectra with MOLLIE using a flat abundance profile with no depletion, which only decreases in the outskirts of the source due to photodissociation. This model is able to accurately reproduce the (1–0) lines, but overestimates the fluxes of the higher-J transitions, which have a higher critical density and are thus sensitive to the abundance in the central part of the core. Furthermore, we derived the molecular column density of N₂D⁺ analysing only the (3–2) transition in the optically thin approximation, using Eq. (4). This allows us to obtain a map of higher resolution, since the beam size at 1 mm is ≈ 12″ (corresponding to ≈0.01 pc at L1544 distance). The results are shown in Fig. 11. Understandably, the absolute values of N_col are different from those obtained via our complete analysis. However, in this case we are more interested in the morphology of the molecular distribution. Figure 11 shows that N₂D⁺ column density is not particularly concentrated around the centre of the core and that it does not peak exactly at this position. Instead, it presents multiple secondary peaks, which are significant with respect to our noise level. These features support our hypothesis of an abundance profile which is not increasing or constant in the innermost part of the core, but that on the contrary is decreasing due to the depletion of the molecule from the gas phase.

Future observational follow-ups with higher resolution will help us to better understand the extent of N₂H⁺ and N₂D⁺ depletion, which is of particular interest since both molecules (and N₂D⁺ in particular) are expected to preferentially trace the central regions of prestellar cores. Moreover, this provides clues on the freeze-out of the parent molecule N₂, which is notobservable in cold and dense environments. It is also interesting to mention that ammonia, in comparison, does not exhibit any sign of depletion at the centre of L1544 (Crapsi et al. 2007).

The DCO⁺ morphologyis also relatively concentrated, similarly to N₂D⁺, but its peak is shifted in the northwest direction. This is confirmed by the integrated intensity maps shown in the bottom panels of Fig. 3. In particular the integrated intensity map for DCO⁺ (3–2), which has the best angular resolution (12″), shows that the molecular emission peaks ≈15″ away from the dust peak. This morphology is due to the counter-balance of two chemical mechanisms: namely (i) the freeze-out onto dust grains of the C- and O-bearing species, which makes them depleted at very high densities and low temperatures (Caselli et al. 1999; Tafalla et al. 2002); and (ii) the fact that deuterated molecules prefer exactly these physical conditions, as described inSect. 1. This explanation additionally gives a straightforward interpretation of the N(HC¹⁸O⁺) map. This molecule presents the highest depletion level of those analysed here, as also confirmed by the abundance profile predicted by our chemical network. Its column density peak is found ≈ 35″ away from the dust peak, in the northwest direction.

Another interesting point is that, as described in Sect. 3.2, we could not reproduce the HCO⁺ and HC¹⁸O⁺ (1–0) spectra using the same abundance profile, simply scaled by the oxygen isotopic ratio. To reproduce the strong blue asymmetry shown by the main isotopologue, we needed to adopt a model with a thick envelope, where the molecule is very abundant (X_col(HCO+) ≳ 10⁻⁷). On the contrary, if we use this dense outer layer also for the rare isotopologue, we are not able to correctly fit the spectral profile, since we overestimate the HC¹⁸O⁺ abundance in the envelope. We believe that this is a consequence of the selective photodissociation of HC¹⁸O⁺. While the main isotopologue is so abundant that it is able to self-shield effectively at large radii, the HC¹⁸O⁺ abundance is not high enough to block the photodissociating photons that can therefore penetrate deeper into the core.

The morphological features previously described for the different column density maps reflect directly on the deuterium fraction maps. Since the D-bearing isotopologues peak closer to the dust peak than the corresponding main species, the D/H maps exhibit a compact morphology around the dust peak. The maximum values that we found are ${\text{D/H}}_{{\text{N}}_2{\text{H}}^{+}}={0.26}_{-0.14}^{+0.15}$ and ${\text{D/H}}_{{\text{HCO}}^{+}}={0.035}_{-0.012}^{+0.015}$, which are in good agreement with previous estimations (see e.g. Caselli et al. 2002a; Crapsi et al. 2005). It is however difficult to make a direct comparison with previous works. For example, Caselli et al. (2002a) used a constant T_ex approach with different excitation temperatures with respect to ours (e.g. T_ex = 5.0 K for N₂H⁺ and 4.9 K for N₂D⁺). These latter study also lacked some important transitions such as N₂D⁺ (1–0) and DCO⁺ (1–0), their frequencies being too low for the receivers used at the time.

We find that N₂H⁺ is more deuterated than HCO⁺, as expected since N₂ (and N₂H⁺ in turn) is a late-type molecule and is less affected by depletion (Nguyen et al. 2018). Therefore, N₂H⁺ appears in conditions where CO is already depleted and H₂D⁺ is more abundant. The deuteration level of N₂H⁺ decreases from 26 to 10% in 45″, while ${\text{D/H}}_{{\text{HCO}}^{+}}$ varies from 3.5 to 2.0%. Figure 12 shows the D/H values for the two tracers as a function of radius for a cut crossing the dust peak and the ${\text{D/H}}_{{\text{HCO}}^{+}}$ maximum, shown with a dashed line in Fig. 10. The deuteration level thus decreases faster for N₂H⁺ than for HCO⁺, mainly as a consequence of the fact thatDCO⁺ is more extended than N₂D⁺ with respect to their main isotopologues and that HCO⁺ and DCO⁺ are more affected by freeze-out in the core centre compared to N₂H⁺ and N₂D⁺.

Fig. 9 D/H ratio obtained for N2H+ in L1544. The beam is shown in red, and the black cross represents the dust peak position.

Fig. 10 Deuterium fraction of HCO+ obtained in L1544. The beam is shown in red, and the black cross represents the dust peak position. The dashed line is the cut usedto produce Fig. 12.

Fig. 11 N2D+ column density obtained analysing only the (3–2) transition, thus obtaining a spatial resolution that is almost three times higher with respect to the corresponding panel in Fig. 8. The white contours represent the integrated intensity of the line, at levels of [5, 9, 11, 12]σ (1σ = 0.05 K km s−1).

Fig. 12 Comparison of the trends of ${\text{D/H}}_{{\text{N}}_2{\text{H}}^{+}}$ (red) and ${\text{D/H}}_{{\text{HCO}}^{+}}$ (blue) along the cut shown in Fig. 10. The data points of HCO+ have been shifted upwards by 0.10 to allow an easier comparison.

5 Conclusions

In this work we performed a detailed analysis of multiple transitions of N₂H⁺, N₂D⁺, HC¹⁸O⁺, DCO⁺ with recent high-sensitivity observations, which allow us to investigate the molecular properties with high-S/N data across the whole L1544 core. Using a non-LTE approach combined with the molecular abundances computed with chemical models, we derived the excitation conditions of the molecules at the dust peak. Our resulting T_ex values confirm that the assumption of a common excitation temperature for different rotational transitions or even for different hyperfine components does not hold for these molecular ions. The molecular abundances derived with our chemical models are in general not able to reproduce the observations, especially for N₂H⁺ and N₂D⁺ for which we had to increase the predicted abundance by up to 300% to obtain a good fit to the observed spectra. Despite these discrepancies, we are able to draw some important conclusions.

All the analysed molecules show some level of depletion in the very central parts of the core due to freeze-out onto the dust grains. This phenomenon is most prominent for HCO⁺ and its isotopologues, but it is seen also for N₂H⁺ and to a lesser extent for N₂D⁺. This is of particular interest because it is the first time that N₂D⁺ depletion, which is often predicted by chemical models, finds confirmation in observational data. Further higher-resolution observations will help us investigate this point. The D/H ratios are peaked very close to the dust peak of L1544, and show different values, the deuteration of N₂H⁺ being higherthan that of HCO⁺ (26 and 3.5%, respectively). In an upcoming paper, we will use the results derived here to obtain spatially resolved information on the ionisation degree, a key parameter in the early stages of star formation.

Acknowledgements

The authors thank the anonymous referee, whose comments helped to increase the quality of the manuscript. E.R. thanks Dr. Jacob Laas for his help in checking the overall language quality.

Appendix A The spectroscopic constants

In Table A.1 we report the values of the spectroscopic parameters used to derive the column densities, according tothe method explained in Sect. 3.3. The references for the shown values are reported in the last column of the table. Figures A.1 and A.2 show the partition function for N₂H⁺ and HCO⁺ (and isotopologues), respectively. In the considered range of temperatures, the correlation is almost linear.

Fig. A.1 Partition function of N2H+ (red circles) and N2D+ (blue crosses) as a function of excitation temperature. Q for rotational transition, neglecting the hyperfine structure, are indicated with solid lines, while the dashed curves take into accountall the hyperfine levels.

Fig. A.2 Rotational partition function of DCO+ (blue crosses) and HC18O+ (red dots), as a function of Tex.

Appendix B Molecular abundance profiles

Here, we plot the abundance profiles computed with our chemical network and tested with the non-LTE approach. Figures B.1–B.3 show the models at the five different time-steps, respectively, for A_V = 1, 2, and 5.

Fig. B.1 Molecular abundances in the model with AV = 1 mag for N2H+ (top left), N2D+ (top right), HCO+ (bottom left), and DCO+ (bottom right). The colours represent the different time-steps: 104 yr (red), 5 × 104 yr (green), 105 yr (blue), 5 × 105 yr (purple), and 106 yr (black).

Fig. B.2 Same as Fig. B.1, but for AV = 2 mag.

Fig. B.3 Same as Fig. B.1, but for AV = 5 mag.

Appendix C Discussion on the excitation temperature

As described in Sect. 3, a main issue to derive the molecular column density when its rotational lines are observed consists in determining the excitation temperatures of the transitions. When T_ex is not the same for all the lines (or when just one line is available), the problem becomes degenerate. A non-LTE analysis is therefore necessary to constrain these parameters (N_col and T_ex) independently.

One could argue that MOLLIE can be used to derive the T_ex values directly. This is only partially true. MOLLIE, as in any non-LTE radiative transfer code that performs ray tracing, does not operate directly with the excitation temperature quantity; it propagates and counts photons, and it solves thestatistical equilibrium equations for the energy level populations. The latter are indeed related to T_ex through the Boltzmann equation. However, the T_ex obtained in this way are only local values, since they can be computed cell by cell⁴. The result is therefore a one-dimensional radial profile of T_ex for each transition.

On the other hand, the T_ex that appears in the equations of radiative transfer is an average quantity, integrated along the line of sight and over the beam size, which relates the observed fluxes with the molecular column density for the beam of that specific observation. It is not straightforward to link the local T_ex profiles obtained in MOLLIE with the T_ex values needed to compute N_col, especially because T_ex is an intensive quantity. An approach that we have tested is to integrate the T_ex profiles along the line of sight, using the H₂ column density as a weighting function. Since H₂ is the main collisional partner, this choice seems reasonable. However, it neglects other important parameters that affect T_ex, such as the gas temperature and the density of the molecule itself. Another point that must be taken into account is that, contrary to T_ex, the molecular abundance is an extensive quantity that can be multiplied by n(H₂) and integrated along the line-of-sight to derive the molecular column density without any ambiguity. It is a natural requirement that the two methods to derive N_col, i.e. the physical and chemical modelling done with MOLLIE and the solution of Eqs. (2) and (3), must be consistent and provide the same result at the centre of the core. With this constraint, T_ex becomes a free parameter that can be tuned to ensure the required consistency. This is the selected method to derive the T_ex values reported in Table 3.

We reiterate that this approach is possible only at the dust peak, because the physical and chemical models are spherically symmetric and cannot reproduce the asymmetries of the core that the observations unveil. Furthermore, the model of Keto et al. (2015) has been specifically developed to reproduce the dust peak.

A strong assumption that we make in Sect. 3.3 is that the T_ex values found, as described immediately above, are constant throughout our map coverage. This is questionable, since one might expect molecules to be less excited in the outskirts of the core, where the gas is less dense. In order to investigate this point, we tried to derive the DCO⁺ column density using a decreasing T_ex profile for each line. DCO⁺ is a representative case, because it has a more extended distribution than for instance N₂D⁺ and three transitions are available at the same time. We realised that a decrease of more than 0.5– 1.0 K within ≈ 60″ (the limit of our map coverage) is enough to cause our N_col fitting routine to break down. With these profiles, the resulting N_col maps present unphysical features, such as values at the border of the core that are several times higher than at the integrated intensity peak, and/or unrealistic morphologies.

Fig. C.1 Isolated component of N2D+ (1–0) (left panel), the central component of N2D+ (2–1) (centre), and the N2D+ (3–2) line (right panel). In each panel, the horizontal dashed lines are the observed peak temperatures of the selected N2D+ transition at the dust peak (red), and at 40″ of offset (green). The shadowed areas represent observational uncertainties. The solid curves indicate TMB as a functionof Tex, obtainedvia the radiative transfer equations and using the column density values predicted by MOLLIE at the two offsets. The vertical dashed line is the Tex value used in the analysis of the maps, and the grey shaded area is its uncertainty (0.5 K).

In order to further test this assumption, we focus on N₂D⁺. As already mentioned, MOLLIE and our physical and chemical models are not able to reproduce the observed spectra outside the dust peak. However, N₂D⁺ is sufficiently centrally concentrated and symmetric around the centre of the core to allow a reasonable fit to the observations also outside dust peak. In particular, we were able to compare the synthetic spectra extracted from the model at 40″ offset with the observed signals obtained by averaging the data in a ring at a distance of 40″ from the centre of the core. We can therefore constrain N_col at two different offsets, and verify that the T_ex values that reproduce these column densities are equal within the uncertainties. Figure C.1 shows the results of this analysis. The red and green dashed lines represent the T_MB observed values of the line components used to derive N_col (see Sect. 3.3.1), at the dust peak and at 40″ offset, respectively. The shaded area indicates the spectra rms⁵. The solid curves indicate the T_MB of each line as a function of T_ex, obtained using Eqs. (2) and (3) using the column density values predicted by MOLLIE at the two offsets (3.5 and 1.6 × 10¹² cm⁻², respectively). Finally, the black dashed lines show the T_ex values adopted in this work, and the grey shadows indicate a 0.5 K variation of this quantity. All the observations are smoothed at the same beam size (34″), which is the resolution of the final N_col map. For all three transitions the solid lines cross the interceptions of the shaded areas, suggesting that within our uncertainties we are not able to detect a significant change of T_ex.

We want to highlight that our analysis does not imply that T_ex is indeed constant across the whole core, but only that our observations are not sensitive to its variation and that they can be analysed assuming a constant value. Most likely, this is also due to the limited coverage of our maps, which are 2′ × 2′ in size, meaning that in a core radius we have at most two or three independent beams. All this said, we assume that our T_ex values reported in Table 3 are accurate to within 0.5 K.

References

All Tables

All Figures

Fig. 1 Observed spectra at the dust peak of N2D+ and DCO+, with the original angular resolution (in black histograms). The transition is labelled in the upper-left corner of each panel. In red, the best-fit solution found with MOLLIE is shown overlaying the observations (see Sect. 3.2). Underneath each line we present the residuals of the fit (difference between the model and the observation). The vertical blue bars show which hyperfine component has been used to derive the molecular column density, when the hyperfine structure has not been neglected.

In the text

	Fig. 2 As in Fig. 1, but for N2H+, HC18O+ and HCO+ transitions.
In the text

Fig. 3 Integrated intensity maps of all the observed transitions, which are labelled in the bottom-right corners. The contours show 20, 40, 60, 80, and 90% of the peak values, which are in order from top left to bottom right panels: 2.0, 2.0, 0.7, 4.6, 0.7, 0.15, 1.8, 1.1, 0.7 K km s−1. The black cross represents the dust peak position, and the white circle is the beam size. All the maps have the same size (2′ × 2′), with the exception of N2H+ (1–0), which is 4′ × 4′: in this map the white rectangle shows the smaller coverage of the other transitions. Scalebars are shown in the top-left corners of the N2D+ (1–0) and N2H+ (1–0) maps.

In the text

	Fig. 4 Profiles of dust temperature (green), gas temperature (red), H2 volume density (black, in logarithmic scale), and velocity (blue, in units of 0.1 km s−1) for the L1544 model developed in Keto et al. (2015). The velocity in the model is negative (it represents contraction motions), but it is shown here as positive for improved readability.
In the text

	Fig. 5 Left panel: molecular abundances that provide the best fit to the spectra at the dust peak for the different species as a function of the core radius. Right panel: zoom-in of the molecular abundances in the inner 0.14 pc.
In the text

	Fig. 6 Results of the line radiative transfer for six chemical models on the three transitions of N2D+, respectively (1–0), (2–1), and (3–2) from left to right in each panel. The synthetic spectra are shown with red histograms overlaying the observations (in black). Each model is labelled with its evolutionary stage and external visual extinction.
In the text

	Fig. 7 Zoom-in of the central 10 000 AU of the N2D+ abundances predicted by the chemical code for AV = 1 mag at four evolutionary stages: t = 5 × 104 yr (green), t = 105 yr (blue), t = 5 × 105 yr (black), and t = 106 yr (red). The abundance profiles are multiplied by a factor of 3.0.
In the text

	Fig. 8 Column density maps obtained for each molecule (labelled in the bottom-right corner of each panel). The beam sizes are shown in the bottom-left corners. The black cross represents the position of the millimetre dust peak.
In the text

	Fig. 9 D/H ratio obtained for N2H+ in L1544. The beam is shown in red, and the black cross represents the dust peak position.
In the text

	Fig. 10 Deuterium fraction of HCO+ obtained in L1544. The beam is shown in red, and the black cross represents the dust peak position. The dashed line is the cut usedto produce Fig. 12.
In the text

	Fig. 11 N2D+ column density obtained analysing only the (3–2) transition, thus obtaining a spatial resolution that is almost three times higher with respect to the corresponding panel in Fig. 8. The white contours represent the integrated intensity of the line, at levels of [5, 9, 11, 12]σ (1σ = 0.05 K km s−1).
In the text

	Fig. 12 Comparison of the trends of ${\text{D/H}}_{{\text{N}}_2{\text{H}}^{+}}$ (red) and ${\text{D/H}}_{{\text{HCO}}^{+}}$ (blue) along the cut shown in Fig. 10. The data points of HCO+ have been shifted upwards by 0.10 to allow an easier comparison.
In the text

	Fig. A.1 Partition function of N2H+ (red circles) and N2D+ (blue crosses) as a function of excitation temperature. Q for rotational transition, neglecting the hyperfine structure, are indicated with solid lines, while the dashed curves take into accountall the hyperfine levels.
In the text

	Fig. A.2 Rotational partition function of DCO+ (blue crosses) and HC18O+ (red dots), as a function of Tex.
In the text

	Fig. B.1 Molecular abundances in the model with AV = 1 mag for N2H+ (top left), N2D+ (top right), HCO+ (bottom left), and DCO+ (bottom right). The colours represent the different time-steps: 104 yr (red), 5 × 104 yr (green), 105 yr (blue), 5 × 105 yr (purple), and 106 yr (black).
In the text

	Fig. B.2 Same as Fig. B.1, but for AV = 2 mag.
In the text

	Fig. B.3 Same as Fig. B.1, but for AV = 5 mag.
In the text

Fig. C.1 Isolated component of N2D+ (1–0) (left panel), the central component of N2D+ (2–1) (centre), and the N2D+ (3–2) line (right panel). In each panel, the horizontal dashed lines are the observed peak temperatures of the selected N2D+ transition at the dust peak (red), and at 40″ of offset (green). The shadowed areas represent observational uncertainties. The solid curves indicate TMB as a functionof Tex, obtainedvia the radiative transfer equations and using the column density values predicted by MOLLIE at the two offsets. The vertical dashed line is the Tex value used in the analysis of the maps, and the grey shaded area is its uncertainty (0.5 K).

In the text