High sensitivity maps of molecular ions in L1544: I. Deuteration of N2H+ and HCO+ and first evidence of N2D+ depletion

Context. The deuterium fraction in low-mass prestellar cores is a good diagnostic indicator of the initial phases of star formations, and it is also a fundamental quantity to infer the ionisation degree in these objects. Aims. With the analysis of multiple transitions of $\rm N_2H^+$, $\rm N_2D^+$, $\rm HC^{18}O^+$ and $\rm DCO^+$ we are able to determine the molecular column density maps and the deuterium fraction in $\rm N_2H^+$ and $\rm HCO^+$ toward the prototypical prestellar core L1544. This is the preliminary step to derive the ionisation degree in the source. Methods. We use a non-local thermodynamic equilibrium (non-LTE) radiative transfer code, combined with the molecular abundances derived from a chemical model, to infer the excitation conditions of all the observed transitions, which allows us to derive reliable maps of each molecule's column density. The ratio between the column density of a deuterated species and its non-deuterated counterpart gives the searched deuteration level. Results. The non-LTE analysis confirms that, for the analysed molecules, higher-J transitions are characterised by excitation temperatures $\approx 1-2\,$K lower than the lower-J ones. The chemical model that provides the best fit to the observational data predict the depletion of $\rm N_2H^+$ and to a lesser extent of $\rm N_2D^+$ in the innermost region. The peak values for the deuterium fraction that we find are $\mathrm{D/H_{N_2H^+}} = 0.26^{+0.15}_{-0.14}$ and $\mathrm{D/H_{HCO^+}} = 0.035^{+0.015}_{-0.012}$, in good agreement with previous estimates in the source.


Introduction
Deuterium was formed during the primordial nucleosynthesis in the first minutes after the Big Bang with a fractional abundance of ≈ 1.5×10 −5 (Linsky 2003). However, enhancements of several orders of magnitude of this value 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 5 cm −3 ) and cold (T 10 K) fragments of molecular clouds, which 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 Datom in hydrogenated species), in fact, is driven by the reaction: This work is based on observations carried out with the IRAM 30 m Telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). which proceeds more efficiently from left to right as the temperature decreases (assuming a low ortho-to-para H 2 ratio; e.g. Pagani et al. 1992) due to the lower zero-point energy of H 2 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 2 D + and of the doubly and triply-deuterated forms of H + 3 rapidly brings to higher D/H ratios in all molecules produced by the reactions with the various H + 3 isotopologues. This process is favoured also by another mechanism. The H + 3 main destruction route is its reaction with CO. Above n ≈ a few 10 4 cm −3 , CO and other C-and Obearing 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 3 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, i.e. the ratio between the free electron density n(e) and the gas density n(H 2 ) [x(e) = n(e)/n(H 2 )] (McKee 1989). In fact, this factor determines the time scale 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. 1998Caselli et al. , 2002bDalgarno 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 to determine the ionisation degree in prestellar cores. In particular, N 2 H + and HCO + represent ideal tracers (Caselli et al. 1998;Caselli 2002). These molecules -and their deuterated counterparts-are abundant in low-mass cores. N 2 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 lowmass star formation. The central part of the core is cold (T ≈ 6 K, Crapsi et al. 2007) and dense (n > 10 6 cm −3 , Tafalla et al. 2002;Keto & Caselli 2010). It is centrally peaked and shows signatures of contraction (Caselli et al. 2002a;Lee et al. 2001;Caselli et al. 2012), suggesting that it is on the verge of gravitational collapse. In the last 20 years, a vast observational campaign combined with theoretical efforts has led to a deep knowledge of its physical and chemical structure (see for instance Crapsi et al. 2007;Caselli et al. 2012;Keto et al. 2015;Spezzano et al. 2017). Previous works (Caselli et al. 1998(Caselli et al. , 2002b) have already investigated the deuteration level and ionisation fraction in this source, using the telescope capabilities available at that time. We have now access to new, high sensitivity maps of N 2 H + , N 2 D + , HC 18 O + , and DCO + across the whole core performed with the On-The-Fly method. For the first time, we can analyse at the same time the first three rotational transitions of N 2 D + and DCO + and the N 2 H + (1-0) and (3-2) lines 1 . In this paper we report the deuteration fraction maps of N 2 H + and HCO + in L1544, obtained with a detailed non-LTE modelling of the molecular spectra at the core's centre. In a following work (Paper II) we will use the results presented here to accurately derive the ionisation fraction map of L1544.
In Section 2 we present the observations. The analysis is described in Sec. 3, divided in chemical modelling (Sec. 3.1), non-LTE modelling (Sec. 3.2) and column density derivation (Sec. 3.4). We discuss the results in Sec. 4, and draw a summary of the work in Sec. 5.
The spectra, originally calibrated in antenna temperature (T * 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 14m telescope towards the dust peak of L1544. These data have a spectral resolution of 0.07 km s −1 and an angular resolution of 55 , and were previously published in Tafalla et al. (1998). We refer to their paper for further details.

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 could allow to determine both its column density and excitation condition (i.e. the line excitation temperature, T ex ). 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 assumptions extends to all the hyperfine components, when present. However, we know from previous works (Daniel et al. 2006(Daniel et al. , 2013 that this is often not the case for the molecules considered here. This can be due to multiple reasons, such as selective photon trapping in very crowded hyperfine structure, or a significant difference in critical 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 they cannot be constrained simultaneously. Therefore, an approach that does not assume LTE is needed to determine each transition's T ex 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 threedimensional 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;Bizzocchi et al. 2013;Caselli et al. 2017;Redaelli et al. 2018). We thus decided to model the spectra at the core's dust peak 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 core's centre, 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 Subsections we provide a detailed description of the needed steps.

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 4 yr, t = 5 × 10 4 yr, t = 10 5 yr, t = 5 × 10 5 yr, and t = 10 6 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 makes it 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 Subsection.

Non-LTE modeling 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 Na + 2.00 × 10 −9 Mg + 7.00 × 10 −9 Fe + 3.00 × 10 −9 P + 2.00 × 10 −10 Cl + 1.00 × 10 −9 hyperfine structure, such as those of N 2 H + and N 2 D + . As an output, it also provides the molecular column density, obtained 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 0 ≈ 10 7 cm −3 and a central temperature of 6.5 K 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. Fig. 1. Observed spectra at the dust peak of N 2 D + and DCO + , with the original angular resolution (in black histograms). The transition is labeled in the upper-left corner of each panel. In red, the best-fit solution found with MOLLIE are overlaid to the observations (see Sec. 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.
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 2 H + are taken from Lique et al. (2015), who computed them for the N 2 H + /p-H 2 system. The rates for the DCO + /p-H 2 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 2 were published by Flower (1999), from whom we scaled the coefficients for HC 18 O + . In all cases, the o-H 2 is neglected, which is a reasonable assumption given that the ortho-to-para ra-tio (OPR) in L1544 is found to be as low as OPR = 10 −3 (Kong et al. 2015).
We tested all the 15 models described in Sec. 3. The abundance profiles produced by the chemical model often overestimate or underestimate systematically the observed fluxes, and need to be multiplied by a numerical factor to find a good agreement with the observations. However, MOLLIE requires too long computational times to allow a full sampling of the parameter space. For each molecule, we thus proceeded as follows: first, we test all the original models, varying the evolutionary stage and the thickness of the surrounding cloud, and we select the one that best reproduces the trend of the different transitions (for instance, the ratio between the peak temperatures of the brightest component of each line). Then, we multiply the whole profile by a numerical factor, testing ≈ 5 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, we present in Sec. 3.3 the resulting spectra from several of the 15 abundance profiles for N 2 D + .

Best fit solution for the observed molecules at the dust peak
The best fit solutions for each molecule are shown in red in Figures 1 and 2, overlaid to the observations. We also include here the analysis of a single-pointing observation of HCO + (1-0) performed at the dust peak (see Sec. 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. Its observed spectrum, together with the best fit from MOLLIE, is shown in the bottom-right panel of Figure 2.
N 2 H + The chemical network systematically underestimates the emission of N 2 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 6 yr), when the molecules at the core's centre start to be affected by depletion. This allows to reproduce the observed ratio between the (1-0) and the (3-2) line. Our simulations show that varying the thickness of the external layer does not have a great impact on the synthetic spectra, due to the high critical density of N 2 H + , whose emission arises thus only from the central regions (n crit 10 5 cm −3 , see Shirley 2015). Furthermore, N 2 is a latetype molecule that can take ∼ 10 times longer than CO to form (e.g. Hily-Blant et al. 2010). Thus, the outer less dense envelope may not be rich in N 2 enough for N 2 H + to form efficiently. The abundance profile for t = 10 6 yr and A V = 1 had to be multiplied by 2.7 to match the observed fluxes.
A&A proofs: manuscript no. PaperDeuteration_2resub 5h04m14.0s 16.0s 18.0s 20.0s RA (J2000) N 2 D + Similar considerations to the ones made for the main isotopologue hold also for N 2 D + . The best fit is found for t = 10 6 yr and an external visual extinction of A V = 1 mag, even though this last parameter does not affect significantly the simulated spectra (see Sec. 3.3 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 2 H + ones. The best agreement for the main isotopologue is found for the model with t = 10 6 yr and an external visual extinction of A 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 −7 ), which is important to correctly reproduce the strong blue-shifted asymmetry of the (1-0) line, due to self-absorption (Tafalla et al. 1998).
HC 18 O + In the chemical model, the abundance profile for the HC 18 O + is derived from the one of HCO + , using the scaling factor valid for the local ISM 16 O/ 18 O = 557 (Wilson 1999). This choice, however, leads to the underestimation of the observed flux. The best fit is found with the model at t = 10 6 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. In fact, if the abundance of HC 18 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 6 yr and A V = 1 mag. In this case, we decrease by 50% the original abundance profile to match the observed peak temperatures. The abundance profiles that provide the best fit solution for each species are shown in Figure 5, after the multiplication for the corrective factors. The adopted models and scaling factors, together with the molecular column density computed by MOL-LIE at the core's centre (N MOLLIE col ), 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 fit with a single model. There are many sources Table 3. Summary of the models that provide the best fit of each molecule, with the adopted scaling factor for the abundance profile.
The velocity in the model is negative (it represents contraction motions), but it is shown here as positive for better readability.
of uncertainty, such as the source model (density and gas/dust temperature structures), chemical modeling parameters and, of course, rate coefficients. Our model does not consider the selfshielding 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.

N 2 D + depletion
In this Subsection, we analyse in detail the case of N 2 D + , showing the non-LTE radiative transfer results for 6 out of the 15 abundances 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, overlaid to the observations, are shown in Figure 6.
The first four rows from the top in Figure 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 4 and t = 10 5 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 5 and t = 10 6 yr), the ratios are 1 : 1.2 : 0.7 and 1 : 1.3 : 0.8, closer to the observed one. This trend can be understood by looking at the behaviour of the different abundance profiles in the inner 10000 AU, shown in Figure  7. At earlier times, the N 2 D + abundance peaks within the central 2000 − 3000 AU, where the molecular hydrogen density is n(H 2 ) 10 5 cm −3 . Since the critical densities of the (1-0), (2-1) and (3-2) transitions are n crit = 7.6 × 10 4 cm −3 , 8.7 × 10 5 cm −3 , and 3.8 × 10 6 cm −3 , 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 2 D + is depleted and the abundance to decreases by more than a factor of five for r 2000 AU. This is due to the fact that its precursor species, N 2 , starts to freezeout onto dust grains, and as a consequence the formation rate of N 2 D + is reduced. The N 2 D + abundance peak moves outwards to r ≈ 5000 AU, where n(H 2 ) ≈ 5 × 10 4 cm −3 . 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 2 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 6 yr) and varying the external visual extinction (see the lower three rows of Figure 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 changes significantly the abundance only for r 10000 AU, where the molecular hydrogen density is n(H 2 ) 10 4 cm −3 . Therefore, the high frequency transitions are not significantly affected by the change of A V .
Article number, page 7 of 19 A&A proofs: manuscript no. PaperDeuteration_2resub The 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.

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 core's centre 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 18 O + (1-0)), or optically thick ones (HCO + (1-0)). Some species present crowded hyperfine structures (such as N 2 H + and N 2 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 well known relations: and τ ν = ln 2 16π 3 Eq.
(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, and T bg = 2.73 K is the cosmic background temperature. η 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-half-maximum equal to ∆V. h is the Planck constant and k B the Boltzmann constant. Eq. (3) depends on several spectroscopic constants: the line frequency ν , the upper state energy E u and degenerancy g u , the Einstein coefficient A ul and the partition function Q. Eq. (3) 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 emphasize that Eq. (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.).

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. Then, we fit N col pixel by pixel, minimising the residuals between the values of T MB obtained through Eq. (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 be used in Eq. (2) and (3) is made so that this approach at the core's centre 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 2 H + The (1-0) transition of N 2 H + presents a single, isolated component at the high-frequency end of its spectrum, which is a perfect candidate for our analysis, providing also reliable values for the linewidth. On the contrary, the hyperfine structure of the (3-2) line collapses in a much 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) 2 . For the spectroscopic constants needed in Eq. 3, we selected the brightest hyperfine component according to the spectroscopic catalogues (see Appendix A). A&A proofs: manuscript no. PaperDeuteration_2resub Fig. 7. Zoom-in of the central 10000 AU of the N 2 D + abundances predicted by the chemical code for A V = 1 mag, at four evolutionary stages: t = 5 10 4 yr (green), t = 10 5 yr (blue), t = 5 10 5 yr (black), and t = 10 6 yr (red). The abundance profiles are multiplied by a factor of 3.0.
N 2 D + The N 2 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 4 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 τ tot ν = i τ i ν , where i labels the hyperfine components. Finally, in the (3-2) spectra the hyperfine structure is closely blended, so that we consider it a single line 3 .
HC 18 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 Eq. (2) and (3), following Caselli et al. (2002b): where W is the integrated intensity of the line. We therefore fit a gaussian to the spectrum in each pixel, we 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 just described, while their peak values are summarised in Table 4. Our results are consistent with the ones from Caselli et al. (2002b), which are quoted to be accurate within a factor of 2. We however highlight that a direct comparison with other works is difficult (see Sec. 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. They 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 in fact 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 pixel by pixel the column density of the D-bearing species and the one of the main isotopologue. First, the maps' beam sizes are matched to allow a fair comparison. The column density of HCO + is derived from the one of HC 18 O + using the standard isotopic ratio: N col (HCO + ) = N col (HC 18 O + ) · 557. Figures 9 and 10 show the resulting images for the D/H ratio.

Evaluation of the uncertainties
In the approach adopted to compute the molecular column densities, the main source of uncertainty are 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 in an average variation of 20 ∼ 30% in the computed column densities, which corresponds to an equal variation of the molecular abundance. We ran again MOLLIE modifying the abundance profiles by 30%, and the resulting modeled 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, provide thus the searched uncertainties. The uncertainties on the D/H ratios are derived through standard er-ror propagation from the ones on the molecular column density. The median relative error are 41% for D/H N 2 H + and 22% for D/H HCO + . The lower mean uncertainty on HCO + results depends on how sensitive N col is to T ex variations, according to Eq. (3). The column density of N 2 D + , for instance, is affected by a T ex change of 0.5 K more heavily than HC 18 O + , and thus its relative error is higher.

Discussion
The column density maps presented in Figure 8 show that different molecules exhibit different morphologies. In particular, one can notice differences in the position of the column density maxima with respect to the dust peak. N col (N 2 H + ) and N col (N 2 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. 2002b). The abundance profiles in Figure 5, despite referring to a one-dimensional core model, confirms these trends. The peak of N 2 H + abundance is found at a larger radius with respect to N 2 D + .
Despite the fact that N 2 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 Sec. 3.3. This is to our knowledge the first time that N 2 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 too poor to investigate this point exhaustively. However, the fact that we need a model with some central depletion to reproduce the observed spectra is a clear confirmation. As a further test, we tried to model all the N 2 H + and N 2 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 reproduce decently 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 2 D + analysing only the (3-2) transition in the optically thin approximation, using Eq. (4). This allows us to obtain a higher resolution map, since the beam size at 1 mm is ≈ 12 (corresponding to ≈ 0.01 pc at L1544 distance). The results are shown in Figure  11. Understandably, the absolute values of N col are different from the ones obtain via our complete analysis. However, in this case we are more interested in the morphology of the molecular distribution. Fig. 11 shows that N 2 D + column density is not so concentrated around the core's centre and 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 resolutions will help us to better understand the extent of N 2 H + and N 2 D + depletion, which is of particular interest since both molecules (and N 2 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 2 , which is not observable 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). Ncol(N2D + ) 1e12 Fig. 11. N 2 D + column density obtained analysing only the (3-2) transition, thus obtaining a spatial resolution almost three times better 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 ).
DCO + morphology is also quite concentrated, similarly to N 2 D + , but its peak is shifted in the north-west direction. This is confirmed by the integrated intensity maps shown in the bottom panels of Figure 3. In particular the one 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: (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); (ii) deuterated molecules prefer exactly these physical conditions, as described in the Introduction. This explanation gives also a straightforward interpretation of the N(HC 18 O + ) map. This molecule presents the highest depletion level of the ones here analysed, 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 North-West direction.
Another interesting point is that, as described in Sec. 3.2, we could not reproduce the HCO + and HC 18 O + (1-0) spectra using the same abundance profile, simply scaled by the oxygen isotopic ratio. In fact, 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 −7 ). 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 18 O + abundance in the envelope. We believe that this is a consequence of the selective-photodissociation of HC 18 O + . While the main isotopologue is so abundant that it is able to self-shield effectively at large radii, the HC 18 O + abundance is not high enough to block the photodissociating photons that can thus penetrate deeper in the core.
The morphology 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 D/H N 2 H + = 0.26 +0.15 −0.14 and D/H HCO + = 0.035 +0.015 −0.012 , which are in good agreement with previous estimations (see for instance Caselli et al. 2002b;Crapsi et al. 2005). It is however difficult to make a direct comparison with previous works. For example, Caselli et al. (2002b) used a constant T ex approach with different excitation temperatures with respect to ours (for instance, T ex = 5.0 K for N 2 H + and 4.9 K for N 2 D + ). They also lacked some important transitions such as N 2 D + (1-0) and DCO + (1-0), being their frequencies too low for the receivers of the time.
We find that N 2 H + is more deuterated than HCO + , as expected being N 2 (and N 2 H + in turn) a late-type molecule and less affected by depletion (Nguyen et al. 2018). Therefore, it appears in physical condition where CO is already depleted and H 2 D + is more abundant. The deuteration level of N 2 H + decreases from 26% to 10% in 45 , while D/H HCO + varies from 3.5% to 2.0%. Figure 12 show the D/H values for the two tracers as a function of radius for a cut crossing the dust peak and the D/H HCO + maximum, shown with a dashed line in Figure 10. The deuteration level thus decreases faster for N 2 H + than for HCO + , mainly as a consequence of the fact that DCO + is more extended than N 2 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 2 H + and N 2 D + .  Figure 10. The data points of HCO + have been shifted upwards by 0.10 to allow an easier comparison.

Conclusions
In this work we performed a detailed analysis of multiple transitions of N 2 H + , N 2 D + , HC 18 O + , DCO + with recent, high sensitivity observations, which allow us to investigate the molecular properties with high SNR 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 2 H + and N 2 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 2 H + and to a lesser extent for N 2 D + . This is of particular interest because it is the first time that N 2 D + depletion, which is often predicted by chemical models, finds confirmation in observational data. Further, higherresolution observations will help us investigate this point. The D/H ratios are peaked very close to dust peak of L1544, and shows different values, being the deuteration of N 2 H + higher than that of HCO + (26% and 3.5%, respectively). In an upcoming paper, we will use the results derived here to obtain spatiallyresolved information on the ionisation degree, a key parameter in the early stages of star formation.

Appendix A: The spectroscopic constants
In Table A.1 we report the values of the spectroscopic parameters used to derive the column densities, according to the method explained in Sec. 3.4. 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 2 H + and HCO + (and isotopologues), respectively. In the considered range of temperatures, the correlation is almost linear.

Rot. Transition
Hyperfine comp. a ν (GHz) Notes. (a) When no hyperfine component is indicated, the transition is intended as single transition. (b) The notation J, F 1 , F ← J , F 1 , * indicates that there are multiple but degenerate lower states, which have the same energy but different F quantum number. (c) The spectroscopic data are taken from: 1) Our calculation based on data from (Cazzoli et al. 2012). 2) Our calculation based on data from (Dore et al. 2004;Amano et al. 2005;Yu et al. 2015). 3) Our calculation based on data from (Caselli & Dore 2005;Lattanzi et al. 2007). 4) From Bizzocchi et al. (in prep.).

Appendix B: Molecular abundance profiles
In this Appendix, we plot the abundance profiles computed with our chemical network and tested with the non-LTE approach. Appendix C: Discussion on the excitation temperature As described in Sec. 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 thus necessary to constrain these parameters (N col and T ex ) independently.
One could argue that MOLLIE can be used to derive directly the T ex values. This is only partially true. MOLLIE, as 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 the statistical 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 4 . The results is thus 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 2 column density as a weighting function. Since H 2 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 2 ) 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 Eq. (2) and (3), must be consistent and provide the same result at the core's centre. With this constrain, 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 remark that this approach is possible only at the dust peak, because the physical and chemical models are spherically symmetric and cannot reproduce the core's asymmetries that the observations unveils. Furthermore, the model of Keto et al. (2015) has been specifically developed to reproduce the dust peak. Appendix C.1: T ex spatial variations A strong assumption that we make in Sec. 3.4 is that the T ex values found as just described are constant across all our map coverage. This can be objected, since one can expect that molecules are 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 2 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 maps' coverage) is enough to make our N col fitting routine to break down. With these profiles, the resulting N col maps present unphysical features, such as values at the core's border several times higher that at the integrated intensity peak, and/or unrealistic morphologies.
In order to further test this assumption, we focus on N 2 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 2 D + is sufficiently centrally concentrated and symmetric around the core's centre to allow a decent 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 core's centre. 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 Sec. 3.4.1), respectively at the dust peak and at 40 offset. The shaded area indicate the spectra rms 5 . The solid curves indicate the T MB of each line as a function of T ex , obtained using Eq. (2) and (3) using the column density values predicted by MOLLIE at the two offsets (3.5 and 1.6 × 10 12 cm −2 , 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. As one can notice, 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 small coverage of our maps, which are 2 × 2 in size, meaning that in a core radius we have at most 2 − 3 independent beams. All this said, we assume that our T ex values reported in Table 3 are accurate within 0.5 K. In each panel, the horizontal dashed lines are the observed peak temperatures of the selected N 2 D + transition at the dust peak (red), and at 40" of offset (green). The shadowed areas represent observational uncertainties. The solid curves indicate T MB as a function of T ex , obtained via the radiative transfer equations and using the column density values predicted by MOLLIE at the two offsets. The vertical dashed line is the T ex value used in the analysis of the maps, and the grey shade is its uncertainty (0.5 K). The panels refer to the isolated component of N 2 D + (1-0) (left), the central component of N 2 D + (2-1) (centre), and the N 2 D + (3-2) line (right).