The Central 1000 au of a Prestellar Core Revealed with ALMA. II. Almost Complete Freeze-out

The initial conditions of the process of star and planet formation are to be found in dense cloud cores (Bergin & Tafalla 2007), especially in the so-called prestellar cores, which are gravitationally bound objects (André et al. 2014) with central number densities above 10⁵ cm⁻³ (Keto & Caselli 2008) and large deuterium fractions measured especially in N-bearing species (Crapsi et al. 2005), which preferentially trace the central regions (Hily-Blant et al. 2010; Spezzano et al. 2017). Detailed studies of prestellar cores in nearby clouds have provided clues on their physical and chemical structure: They are centrally concentrated, with low central temperatures (∼7 K, Crapsi et al. 2007; Pagani et al. 2007; Launhardt et al. 2013) and evidence of grain growth (Forbrich et al. 2015; Chacón-Tanarro et al. 2019). They display a clear chemical differentiation, where rare CO isotopologue emission maps present a valley within the central ∼8500 au, due to freeze-out onto dust grains (Caselli et al. 1999), while N-bearing species such as N₂H⁺, NH₃, and their deuterated forms appear to trace well the core center (Tafalla et al. 2002, 2004; Juárez et al. 2017; Spezzano et al. 2017), defined as the millimeter and submillimeter dust continuum emission peak. Although there is evidence of some central depletion of N₂H⁺ (Bergin et al. 2002), N-bearing deuterated species have been found to trace well the dust peak (Caselli et al. 2002). Only recently, high-sensitivity single-dish observations of multiple transitions of N₂D⁺, coupled with astrochemical and non-LTE radiative transfer modeling, have provided evidence of N₂D⁺ central depletion (Redaelli et al. 2019), but the size of the related depletion zone cannot be measured because of the limited angular resolution. The abundance of N₂H⁺ isotopologues is linked to N₂, which has a binding energy similar to that of CO (Bisschop et al. 2006), although recent studies have shown that the N₂ (but not the CO) binding energy is reduced in CO–N₂ ice mixtures (Nguyen et al. 2018). The smaller amount of N₂ freeze-out compared to CO can be explained if N₂ formation proceeds at a slower pace than CO formation, so that most of N₂ are still forming from (the more volatile) atomic nitrogen in the dense material when CO has already started to condense onto icy mantles (Flower et al. 2006; Hily-Blant et al. 2010; Pagani et al. 2012).

No evidence of ammonia freeze-out in starless and prestellar cores has been reported yet, despite the fact that the NH₃ binding energy is significantly larger than that of CO and N₂ and close to that of water molecules (Sipilä et al. 2019). This is puzzling, as the freeze-out timescale within the central 1000 au of prestellar cores, where the number density is about 10⁶ cm⁻³, is only ∼1000 yr, much shorter than dynamical timescales. Furthermore, it appears that the NH₃ abundance is slightly increasing toward the center when viewed with single-dish and interferometers (Tafalla et al. 2002; Crapsi et al. 2007; Caselli et al. 2017). Such an increasing abundance profile for NH₃ is at odds with state-of-the-art astrochemical modeling (Sipilä et al. 2019). Crapsi et al. (2007) also observed para-NH₂D(1₁₁−1₀₁) toward the Taurus prestellar core L1544 with the IRAM Plateau de Bure Interferometer, finding an integrated intensity map centered on the dust peak, thus ruling out any central depletion of NH₂D when viewed at an angular resolution of 58 × 45 (986 × 765 au, at the distance of 170 pc; ¹⁴ Galli et al. 2019). From this result, it was concluded that no complete depletion zone ¹⁵ (Walmsley & Flower 2004) was present in L1544.

With ALMA 1.3 mm dust continuum emission observations of L1544, we have recently revealed a compact central region, with mass ≃0.16 M_⊙, radius ≃1800 au, and average H₂ number density of ≃10⁶ cm⁻³, called the "kernel" (Caselli et al. 2019). This kernel is the stellar system cradle, which we are investigating in detail to shed light on the first steps toward star and protoplanetary disk formation. Here we present the ALMA (12 m+ACA) para-NH₂D(1₁₁−1₀₁) mosaic, which reveals for the first time the NH₂D depletion zone toward the L1544 kernel at an angular resolution of 25 (425 au; Sections 2 and 3). We also demonstrate that the current data are in agreement with our astrochemical model predictions (Section 4), implying that almost all (99.99%) species heavier than He reside on dust icy mantles within the kernel. Hereafter, pNH₂D indicates para-NH₂D. A detailed study of the kernel kinematics will be presented in a forthcoming paper.

2.1. Deuterated Ammonia

We observed L1544 with ALMA during Cycle 4 (ESO Project ID 2016.1.00240.S; PI: Caselli) using the 12 m (Main array) and the 7 m (ACA) arrays. The 12 m array three-pointing observations were carried out on 2016 December 22, while the 7 m array single-pointing observations were carried out on 2016 October 2. The Main array observations included the quasar J0510+1800 for bandpass and phase calibration, while quasar J0423–0120 is observed for flux calibration. The ACA observations included the quasar J0510+1800 for flux, bandpass, and phase calibration. The three-point mosaic observations used a correlator with the continuum window centered at 100 GHz and a spectral window for pNH₂D (1₁₁–1₀₁) centered at 110.15355 GHz (these data will be presented in a future paper dedicated to the study of dust properties within the kernel). The bandwidth and spectral resolution for the pNH₂D observations are 58.59 MHz and 30.518 kHz (corresponding to a velocity resolution of 83 m s⁻¹), respectively. The 12 m and 7 m array data were calibrated using the Common Astronomy Software Applications package (CASA; McMullin et al. 2007) version 4.7.0. We assigned weights to the data sets with the statwt command before combining them (concat command). We perform a joint deconvolution of both data sets (ACA and 12 m array) in CASA (version 5.6.2) using tclean for mosaics, natural weighting, and restored with a circular beam of 25. We use a rest frequency of 110.153562 GHz for the imaging (De Lucia & Helminger 1975; Cohen & Pickett 1982). We use the multiscale clean technique, with the multiscale parameter of [0, 2.5, 7.5, 22.5]''. The final noise level is 7 mJy beam⁻¹ per 0.08 km s⁻¹ channel, which corresponds to a noise level of 0.13 K in the central 45'' diameter region of the primary-beam-corrected data. The integrated intensity map of pNH₂D(1₁₁−1₀₁) (ACA + 12 m array) is shown in the top panel of Figure 1, and the full mosaic is shown in Figure A1. We note that the combined ACA + 12 m data perfectly recover the single-dish flux within the IRAM 30 m beam (see Figure B1), demonstrating that mosaic pointing observations with ALMA are a powerful technique to carry out detailed studies of prestellar cores.

Zoom In Zoom Out Reset image size

2.2. 1.3 mm Continuum

We perform the line fit of pNH₂D(1₁₁−1₀₁) using the nh2d model implemented in pyspeckit (Ginsburg & Mirocha 2011). The pNH₂D (1₁₁–0₀₁) line is modeled using the hyperfine structure and relative intensity from the recent laboratory works of Melosso et al. (2020, 2021). The model assumes a Gaussian velocity distribution (with centroid velocity and velocity dispersion V_lsr and σ_v ), equal excitation temperatures (T_ex) for all hyperfine components, and total optical depth (τ₀). We perform a fit for the whole cube with all four free parameters. Following the procedure in Pineda et al. (2021), we discard all centroid velocity determinations if the uncertainty in V_lsr is larger than 0.02 km s⁻¹, all velocity dispersion determinations if the uncertainty σ_v is larger than 0.015 km s⁻¹, and the excitation temperature if the uncertainty is larger than 1 K.

The column density has been derived following Mangum & Shirley (2015):

$$\begin{eqnarray}\begin{array}{rcl}N({\mathrm{pNH}}_{2}{\rm{D}}) & = & \displaystyle \frac{8\pi {\nu }^{3}}{{c}^{3}{A}_{\mathrm{ul}}}\displaystyle \frac{Q({T}_{\mathrm{ex}})}{{g}_{u}{e}^{-{E}_{u}/{{kT}}_{\mathrm{ex}}}}\\ & & \times {\left[{e}^{h\nu /{{kT}}_{\mathrm{ex}}}-1\right]}^{-1}\displaystyle \int \tau {dv},\end{array}\end{eqnarray} \tag{ 1 }$$

where $\int \tau {dv}=\sqrt{2\pi }\ {\tau }_{0}{\sigma }_{v}$, ν is the frequency of the transition, g_u is the upper-level degeneracy, E_u is the upper energy level of the transition, A_ul is the Einstein coefficient, and Q(T_ex) is the partition function. The values for E_u , g_u , and A_ul are obtained from the Leiden Atomic and Molecular Database (LAMDA; Schöier et al. 2005; van der Tak et al. 2020), ¹⁶ and the partition function is calculated up to the first 30 energy levels. This formula has been used previously (Daniel et al. 2016; Harju et al. 2017). T_ex has been estimated at each position using the hfs fitting within pyspeckit, with average uncertainties of 2% (and a maximum of 5% at the edge of the map). The error associated with the column density is about 10%, close to the calibration error.

The column density map in Figure 1 shows a complex structure, partly due to noise, without a clear peak toward the 1.3 mm dust continuum peak, which identifies the kernel (Caselli et al. 2019). This implies that the pNH₂D column density does not follow the H₂ column density, thus a drop in pNH₂D fractional abundance (with respect to H₂) within the kernel radius of ≃1800 au (≃10'') must be present. The pNH₂D column density flattening within the central ≃3000 au, which includes the L1544 kernel, is well illustrated in Figure 2. The circles represent individual measurements, while the black and gray horizontal lines are the mean and its associated uncertainty in a 25 bin. The projected radius pr is calculated as ${pr}=\sqrt{{r}_{\mathrm{maj}}^{2}+{({r}_{\min }{b}_{\mathrm{axes}})}^{2}},$ where r_maj and r_min are the distances along the semimajor and -minor axes, respectively; b_axes = 22/12 is the semimajor and -minor axis ratio for the pNH₂D emission, where the semimajor axis is along the position angle of 160° (measured east from north; see the bottom panel of Figure 1 for an illustration of the elliptical shape and orientation used to define pr).

Zoom In Zoom Out Reset image size

We carried out chemical and radiative transfer modeling to reproduce the observed pNH₂D column density and line emission distribution. First, we produced simulated pNH₂D abundance distributions using the most recent version of our gas-grain chemical code (Sipilä et al. 2020), which includes an extensive description of deuterium and spin-state chemistry, of critical importance for the present work. Very briefly, we used the one-dimensional (1D) source model for L1544 presented by Keto et al. (2015), which was divided into concentric shells; chemical simulations were carried out separately in each shell. This process yields time-dependent abundance profiles. The initial abundances used in the chemical modeling are displayed in Table 1. We assume a constant grain radius of 0.1 μm, and a grain material density of 2.5 g cm⁻³. Other model details can be found in Sipilä et al. (2020) and are omitted here for brevity.

We searched for the best fit to the observed line map by extracting the pNH₂D abundance profile from the chemical model at various time steps (between 10⁴ and 10⁶ yr) and using the abundance data as input for the non-LTE line radiative transfer code LOC (Juvela 2020) to simulate the line emission. We adopted the pNH₂D–pH₂ collisional coefficients from Daniel et al. (2014) and assumed that the hyperfine components in each rotational level have the same excitation temperature. The hyperfine component frequencies and relative intensities have been adopted from Melosso et al. (2021). The velocity profile of the contracting cloud (including a uniform level of turbulence of 0.1 km s⁻¹) was taken from Keto et al. (2015). With this setup, the best fit between the model and the observations is reached after a chemical evolution of t ∼ 2.5 × 10⁵ yr, as determined by a χ² analysis applied to modeled and observed spectra toward six positions at and near the center of the core. ¹⁷ We stress that the chemical best-fit time of 2.5 × 10⁵ yr cannot be used as a chemical clock or to deduce the dynamical age of the core. The reason for this is that the chemistry is run using a static physical structure as input (a BE sphere with a density, temperature, and velocity profile deduced by Keto et al. 2015), so chemical processes are not followed simultaneously with the dynamical evolution. Future work will focus on self-consistent chemical-dynamical evolution models, using hydro and magnetohydrodynamics simulations of the types mentioned in Appendix D.

Figure 3 shows the simulated pNH₂D column density map and the (1₁₁−1₀₁) line emission profiles at the best-fit time. The modeled column densities and line emission maps agree well with the observations, showing a relatively flat column density map within the central ≃13'' (or about 2200 au). In fact, although the pNH₂D column density map in Figure 3, left panel, shows a central valley within a bright ring of ≃5'' (≃850 au), while our observations show a flattened structure (see Figure 2), the model column density at the dust peak is only 4% lower than at the bright ring (1.22 × 10¹⁴ cm⁻² instead of 1.27×10¹⁴ cm⁻²), well within our calibration uncertainties. As we will show in Section 5, this implies that our observations are consistent with NH₂D freeze-out in the central core. We have verified that the (1₁₁−1₀₁) line emission profiles cannot be reproduced by the models without NH₂D depletion toward the center (see Appendix C). We note that the observed average pNH₂D column density in Figure 2 is slightly (less than a factor of 2) larger than that in Figure 3 (≃2 × 10¹⁴ cm⁻² versus ≃1.2 × 10¹⁴ cm⁻²), which is once again underlying the good agreement between model and observations, despite the simplicity of the model. The measured and modeled excitation temperatures T_ex are also similar, as can be seen in Figure 4, which compares the measured T_ex (as deduced from the pyspeckit hfs fit of all the spectra within the map in Figure 1) as a function of projected radius (pr; see Section 3) with the excitation temperature profile within the core (T_ex(r), with r the core radius) predicted by the radiative transfer LOC. Although a one-to-one comparison between the two excitation temperatures cannot be made, as one should instead compare the measured T_ex with the integral along the line of sight of the model T_ex(b) for each impact parameter b, the two values are comparable, with the model T_ex(r) within 20%–30% of the measured T_ex.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A final note concerns the adopted physical structure. The Keto et al. (2015) 1D model used here provides self-consistently volume density and velocity profiles across the core, and it can be easily input in comprehensive gas-grain chemical models such as the present one and the previous work by, e.g., Vasyunin et al. (2017) on complex organic molecules in L1544. Naturally, this 1D model cannot reproduce the elongated structure observed at 1.3 mm (see Figure 3 in Caselli et al. 2019 and Appendix D). In fact, ALMA observations are better reproduced by 3D nonideal magnetohydrodynamic (niMHD) simulations of a contracting cloud with a peak volume density of 10⁷ cm⁻³ (see their Figure 4). The 3D niMHD simulation gives rise to a prestellar core with an elongated flattened central structure, in harmony with the observed 1.3 mm ALMA dust continuum emission. The adoption of this 3D model for chemistry is, however, beyond the scope of this paper, and it will be used in a future publication, where only simple species will be followed in a reduced chemical network to be included in the dynamical simulation. Despite the difference between the two (1D and 3D) models, it is interesting to note that the density and velocity profiles along the major and minor axes of the simulated 3D core do not differ substantially from those predicted by the 1D model of Keto et al. (2015), as shown in Figure D1. In particular, the density and velocity profiles along the major axis of the 3D model are especially close to our adopted 1D model (see Figure D1), suggesting similarities that, together with the radiative transfer results shown in Figure 3, make our 1D model a very good approximation of the L1544 structure.

In the previous sections, we showed that ALMA has allowed us to resolve for the first time the distribution of deuterated ammonia in the central 1000 au of the prestellar core L1544. Despite the complex morphology of the L1544 core as viewed at the frequency of the pNH₂D (1₁₁−1₀₁) line (see Figure 1), it is clear that the pNH₂D column density has a flattened distribution in correspondence with the dust continuum peak, where the kernel is located (see also Figure 2). The radius of the kernel is about ≃1800 au, which is ∼5 times smaller than the CO depletion zone observed toward the same prestellar core by Caselli et al. (1999). This is the reason why previous interferometric observations from Crapsi et al. (2007), done with poorer angular resolution (58 × 45), could not reveal the abundance drop of deuterated ammonia, as demonstrated in Appendix E.

In the following, we will present the fractional abundance profiles predicted by our chemical model, at 2.5 × 10⁵ yr, of all ammonia isotopologues as well as the profile of the total depletion factor, ${f}_{{{\rm{D}}}_{\mathrm{tot}}}$, defined as

$$\begin{eqnarray}&&{f}_{{{\rm{D}}}_{\mathrm{tot}}}=\displaystyle \frac{{\sum }_{i}[x{(i)}_{\mathrm{gas}}+x{(i)}_{\mathrm{solid}}]}{{\sum }_{i}x{(i)}_{\mathrm{gas}}},\end{eqnarray} \tag{ 2 }$$

where i represents any species in the chemical model containing an element heavier than He. Thus, ${f}_{{{\rm{D}}}_{\mathrm{tot}}}$ quantifies the amount of elements heavier than He subsisting in the gas phase, with the higher values corresponding to higher levels of depletion.

The abundance profile of pNH₂D at the best-fit time is displayed in Figure 5 along with ${f}_{{{\rm{D}}}_{\mathrm{tot}}}$. It is interesting to note that the pNH₂D abundance starts to drop already within about 7000 au, whereas the integrated intensity map, as well as the simulated observations, shows a centrally concentrated structure up to the central ≃2000 au. The abundance drop is in fact partially hidden to observations due to the radiative transfer of the pNH₂D(1₁₁−1₀₁) line. In particular, despite a "catastrophic" pNH₂D freeze-out present within the central 2000 au (see Figure 5), the corresponding decrease in column density is only ∼4%; this is shown in the left panel of Figure 3, where a relatively bright ring of about 5'' in radius, within the pNH₂D depletion zone, is clearly visible. This is due to the fact that the column density is dominated by larger scales, where pNH₂D is abundant. The critical density of the (1₁₁−1₀₁) line is ∼1.4 × 10⁵ cm⁻³ at T = 10 K, and so collisional excitation dominates in a region roughly 13'' in radius, enclosing the bright ring in the left panel of Figure 3. The inferred column density ring and shallow inner valley (Figure 3) are due to the drop in the pNH₂D abundance toward the center of the core. The simulated main hyperfine component peak has an optical depth close to 1, in agreement with our observations, and the strength of the whole line emission follows the morphology of the column density distribution.

Zoom In Zoom Out Reset image size

The abundance distributions of the modeled normal and deuterated ammonia (summed over the spin states) as well as that of CO are also shown in Figure 5 (right panel); it is evident that the depletion of singly deuterated ammonia within the central 2000 au is not due to conversion into multiply deuterated ammonia, but is instead a result of freeze-out onto grain surfaces. These results, including the ortho-to-para ratios of NH₂D, are consistent with those obtained by Hily-Blant et al. (2018), who modeled the molecular composition and the nuclear spin chemistry of a collapsing Bonnor–Ebert (BE; Ebert 1955; Bonnor 1956) sphere. According to our model, the total depletion factor ${f}_{{{\rm{D}}}_{\mathrm{tot}}}$ increases from 1000 to 10,000 in the central 2000 au, implying that between 99.9% and 99.99% of all the species heavier than He are predicted to reside on the surface of dust grains within the kernel. This has the consequence that thick icy mantles are expected to form on top of dust grains within the central regions of prestellar cores, just before the formation of a protostar and protostellar disk. Thick icy mantles could be responsible for the change in dust opacity measured toward the center of L1544 by Chacón-Tanarro et al. (2019).

We note that, prior to this work, coarser interferometric observations of pNH₂D showed a clear centrally peaked distribution toward L1544 (Crapsi et al. 2007). For this reason, pNH₂D was considered the species most resilient to freeze-out (compared with all other species observed so far toward this core, which all show direct or indirect evidence of depletion, including N₂D⁺; e.g., Spezzano et al. 2017; Redaelli et al. 2019). Now that even pNH₂D freeze-out has been unveiled and our modeling reproduces the observed column density map and line profiles of pNH₂D(1₁₁−1₀₁), we can conclude that the almost complete freeze-out predicted by our chemical model is then consistent with the data available for this source. Moreover, as already mentioned, previous single-dish and interferometric observations suggested an NH₃ abundance profile increasing toward the central regions (e.g., Crapsi et al. 2007). To investigate this apparently contradicting result, we also simulate the NH₃(1,1) line profiles using the pNH₃ abundance profile predicted by our model, after 2.5 × 10⁵ yr of chemical evolution, and compare them with the VLA-only observations in Figure 2 of Crapsi et al. (2007). We also consider the case without NH₃ freeze-out in the central 6000 au (which differs from our predicted profile shown in Figure 5, right panel) and see if this better reproduces the observations, in line with the conclusions from Crapsi et al. (2007). The assumed beam for the simulated observations of NH₃(1,1) is 4'', close to the VLA resolution of 435 × 345. The spectra are reported in Figure F1, showing that our model prediction of NH₃ freeze-out toward the center (see Figure 5) is consistent with the Crapsi et al. (2007) observations, including the relative intensity of the groups of hyperfines and the negligible difference between the spectrum centered at the dust peak position and the one at the 223 offset. The main difference is that our simulated lines are about two times brighter than the observed ones, but this is expected due to the fact that we are not simulating interferometric observations, and the core envelope is significantly contributing to the simulated line flux (unlike in the case of interferometric observations). In fact, higher-sensitivity VLA data of NH₃(1,1) toward L1544, combined with GBT single-dish data, show a peak temperature in better agreement with our modeling; these data will be presented in an upcoming paper focusing on a new gas temperature map of L1544 (A. Schmiedeke et al. 2022, in preparation). Another interesting result of this comparison is that variations in the abundance profile of ammonia in the central 5000 au do not change the spectra (see dashed histograms in Figure F1, which refer to a model with a flat fractional abundance of NH₃ of 2 × 10⁻⁸ within the central 5000 au), demonstrating that NH₃(1,1) cannot be used to measure its abundance profile within the central few thousand astronomical units because of its relatively low critical density (see Appendix F). Interestingly, new high-resolution and high-sensitivity NH₃ observations of the prestellar core H-MM1 in Ophiuchus show a clear depletion zone toward its dust peak (J. E. Pineda et al., 2022, in preparation), thus confirming our predictions.

A final comment should be added about the velocity structure used as input to the chemical model. In some previous studies of L1544, the Keto & Caselli (2010) infall velocity profile had to be scaled up by a factor of ∼2 to reproduce observations obtained with single-dish telescopes (Bizzocchi et al. 2013; Redaelli et al. 2019). In the present case, a similar scaling leads to double-peaked line profiles, even for the satellite components that are completely excluded by the observations. This discrepancy with the earlier single-dish studies is likely related to the 1D spherically symmetric structure of the source model (while the source itself is certainly not spherical on the large scales traced by single-dish observations) and to the possible presence of more complex velocity fields at the core scale (see Appendix D). The present interferometric observations probe a relatively small volume, where the number density is close to the critical density of the observed transition (1.4 × 10⁵ cm⁻³), i.e., a region with a radius of about 4000 au, considering the adopted density profile shown in Figure D1 (red dashed curve); hence, the 1D symmetry assumption is adequate, as already discussed. Therefore, the Keto & Caselli (2010) physical model, which already successfully reproduced line profiles of several species observed with single-dish telescopes toward L1544 (e.g., Caselli et al. 2012; Keto et al. 2015), still remains valid for these new ALMA interferometric observations, as well as recent work from Redaelli et al. (2021).

Our ALMA spectroscopic observations, coupled with chemical and radiative transfer modeling, of the prototypical prestellar core L1544 have unveiled the depletion zone of deuterated ammonia, with size very close to the kernel previously found with ALMA continuum observations (radius ≃1800 au). This has not been possible before because of the poorer angular resolution (factor of ≥2) and quality of previous observations of the same line. Our chemical model, applied to the physical structure of L1544, reproduces remarkably well the pNH₂D(1₁₁−1₀₁) emission. Simulated observations of the modeled prestellar core show that the high pNH₂D fractional abundance in regions outside the kernel, with volume densities close to the critical density of the (1₁₁–1₀₁) transition, prevent a clear view of the actual catastrophic NH₂D freeze-out region without high-enough angular resolution. We find that more than 99.9% of all species heavier than He reside in thick icy mantles within the central kernel, the future stellar cradle. In prestellar cores like L1544, with central volume densities ≥10 ⁶ cm⁻³ and corresponding freeze-out timescales ≤10³ yr, an almost complete freeze-out is therefore inevitable, just before star formation. The thick icy mantles (Draine 1985) may promote coagulation of dust grains (e.g., Ormel et al. 2009; Chacón-Tanarro et al. 2019) and allow preservation and delivery of the frozen prestellar chemistry into the next stages of evolution, when the protostar and the protoplanetary disk will form. Some of this prestellar ice, especially those trapped within icy pebbles in the outer part of the disk, may survive later stages of planet formation and evolution. Indeed, similarities in molecular composition and isotopic fractionation between star- and planet-forming regions and primitive material in our solar system, such as comets and carbonaceous chondrites (e.g., Mumma & Charnley 2011; Caselli & Ceccarelli 2012; Ceccarelli et al. 2014; Cleeves et al. 2014; van Dishoeck et al. 2014; Öberg et al. 2015; Altwegg et al. 2019; Drozdovskaya et al. 2021), strongly hint at a crucial role of prestellar chemistry in shaping the composition and further chemical/mineralogy evolution of the building blocks of planets.

We are grateful to Malcolm Walmsley, for the many inspiring discussions on complete depletion that lead to this work. We thank the anonymous referee and Doug Johnston for a final review, which improved the clarity of the paper. We gratefully acknowledge the support of the Max Planck Society. I.J.-S. has received partial support from the Spanish State Research Agency (PID2019-105552RB-C41). Z.Y.L. is supported in part by NASA 80NSSC18K1095 and NSF AST-1910106. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.01195.S and ADS/JAO.ALMA#2016.1.00240.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

In Figure 1, a zoom-in view of the pNH₂D(1₁₁−1₀₁) integrated intensity map is shown in order to highlight the structure and better compare with the column density map. In this appendix, we present the full three-point mosaic of the same map to show the extent of the mapped area (see Figure A1). The six white circles along the major axis in the southeast direction in Figure A1 enclose the areas where the spectra have been extracted for comparison with the model spectra (see Figure 3). The spectra from the other six white circles, those along the minor axis, are reported in Figure A2 and compared with the same model spectra as in Figure 3. It is apparent that the outer positions in Figure A2 (numbers 5 and 6) show model spectra brighter than the observed ones, suggesting that the major axis has a physical and chemical structure closer to our spherically symmetric model. This is supported also by Figure D1, where the density and velocity profiles of the 3D model are compared with our 1D model profiles.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure B1 shows the comparison between the pNH₂D (1₁₁−1₀₁)spectrum obtained with IRAM 30 m single-dish observations (with a 22'' beam) and the one obtained by ALMA (12 m + ACA) within the same area. The IRAM 30 m spectrum of pNH₂D(1₁₁−1₀₁) has been extracted at the dust peak of an On-The-Fly map covering the inner 22'' × 22'' of L1544 around the dust peak. The emission map was obtained in the framework of project 013–13 (PI: S. Spezzano) in 2013 October . Position switching was used, with the reference position set at (−180'', 180'') offset with respect to the map center. The EMIR E090 receiver was used with the Fourier Transform Spectrometer backend with a spectral resolution of 50 kHz. The mapping was carried out in good weather conditions (τ ∼ 0.03) and a typical system temperature of ∼90–100 K. It is clear from the figure that all the flux is fully recovered, demonstrating that ALMA mosaicking, together with the combination of 12 m and ACA data, provides a powerful technique to study prestellar cores in detail.

Zoom In Zoom Out Reset image size

As shown in Figure 3, our model reproduces very well the pNH₂D(1₁₁ − 1₀₁) spectra toward the L1544 kernel. This model predicts that NH₂D (as well as other species containing elements heavier than He) quickly disappears from the gas phase within the central 2000 au because of freeze-out onto dust grains (see Figure 5). To test that central freeze-out is indeed needed to reproduce observations, Figure C1 shows the predicted pNH₂D(1₁₁−1₀₁) spectra when the NH₂D fractional abundance is kept at the constant value of 6 × 10⁻⁹ within the central 2000 au, instead of dropping as in Figure 5. The predicted spectra clearly overestimate the observed intensities in the central four positions, thus demonstrating that observations can only be reproduced when freeze-out is taken into account, as predicted by our chemical models.

Zoom In Zoom Out Reset image size

The physical structure of the L1544 core used for the chemistry code has been described in Section 4. As discussed there, the density profile from the self-consistent hydrodynamic simulation of a contracting BE (Ebert 1955; Bonnor 1956) sphere in quasi-static equilibrium (Keto & Caselli 2010) reproduces well line profiles previously observed toward the L1544 dust peak (e.g., Keto & Caselli 2010; Caselli et al. 2012; Keto et al. 2015). However, because of the spherical symmetry of the BE model, the elongated 1.3 mm dust continuum emission observed at high angular resolution with ALMA cannot be reproduced (Caselli et al. 2019). Indeed, Caselli et al. (2019) showed that the central kernel is consistent with flattening produced by nonideal magnetohydrodynamic (niMHD) simulations of the contraction of a rotating magnetized cloud core. The peak density within the flattened structure in the 3D niMHD simulation reaches 10⁷ cm⁻³, in agreement with the peak density at the center of the BE sphere, while the average density within the kernel is about 10⁶ cm⁻³, consistent with previous single-dish observations (Crapsi et al. 2007; Chacón-Tanarro et al. 2019), which did not have high-enough angular resolution to detect the central density enhancement.

Despite the structural difference between the central 1000 au of the Keto et al. (2015) model (used as input for our chemical model) compared with the 3D niMHD simulated core (which better reproduces the ALMA 1.3 mm dust continuum emission), the agreement between model and observed spectra in Figure 3 is excellent. This is partly due to the choice of the cut, taken along the (projected) major axis of the kernel, as this is the direction where the density profile of the 3D structure more closely resembles the BE sphere (see Figure D1, left panel). The agreement with observed spectra along the minor axis becomes poorer (than that with observed spectra along the major axis) in the outer two positions, as seen in Figure A2. To reproduce in detail the spectra along the minor axis, we will have to include the chemistry in the 3D niMHD simulation, and this will be done in a future paper. For now, it is interesting to point out that the BE sphere in quasi-static contraction of Keto et al. (2015) has a density profile just in between the minor and major axes of the 3D niMHD simulated core (Figure D1, left panel). Also, the radial velocity profile of the contracting BE sphere is very similar to that along the minor and major axes within a radius of 500 au (Figure D1, right panel); above this radius, the velocity profile of the contracting BE sphere follows more closely that along the major axis of the 3D niMHD core, and it is only a factor of 2 away from that along the minor axis. Considering the simplicity of the contracting BE sphere HD simulation compared with the niMHD simulation of the contracting magnetized core, it is striking to see how closely the two resemble each other.

Zoom In Zoom Out Reset image size

Crapsi et al. (2007) observed the same transition of deuterated ammonia using the IRAM Plateau de Bure interferometer with an angular resolution of 58 × 45 . To demonstrate that this is not sufficient to detect the pNH₂D depletion zone, we perform simulated observations of the same modeled core described in Section 4, assuming a spherical half-power beamwidth of 5''. Figure E1 presents six simulated spectra extracted from six positions similar to those shown in Figure 3, but now separated by 5'' (instead 25). A 5'' beam does not have a substantial effect on the simulated lines toward the center of the core (when compared with the red spectra in Figure 3). The most noticeable difference is that the satellite components become slightly less bright (and less optically thick) compared to the simulations with a 25 beam; the line brightness now peaks toward the central two positions (number 1 and 2 in Figure E1). An integrated intensity map would then show a peak at the center and a monotonic drop toward the outer parts, consistent with the Crapsi et al. (2007) observations. Therefore, with 5'' resolution, as for the integrated intensity map, the central flattening in the column density profile could not be detected because it is only present within the central resolution elements of the Crapsi et al. (2007) observations.

Zoom In Zoom Out Reset image size

We would also like to point out another important difference between our ALMA observations and those carried out by Crapsi et al. (2007): Our new ALMA image of pNH₂D(1₁₁−1₀₁) is the result of a three-point mosaic (see Figure A1), while the IRAM-PdBI observations consist of a single-point image without including a primary-beam correction. This may cause an artificially steeper flux drop away from the central region, mimicking a centrally concentrated pNH₂D flux profile consistent with the distribution of the dust continuum emission. Finally, this leads to the misleading conclusion that central freeze-out is not needed for deuterated ammonia, pointing out the importance of high-angular resolution and high-sensitivity mosaic observations for a correct interpretation of prestellar core observations.

Given that previous work has claimed that the fractional abundance of NH₃ is actually increasing toward the L1544 center (e.g., Crapsi et al. 2007), we would like to test if our model predictions about the NH₃ fractional abundance across the core (see Figure 5) is consistent with previous interferometric work. For this, we simulate VLA NH₃(1,1) observations of our model cloud using the LOC radiative transfer code and assuming a telescope beam of 4'' (close to the VLA synthesized beam in the observations of Crapsi et al. 2007, 434 × 345), at two different positions: the dust peak and the offset (10'', −20''), as in Figure 2 of Crapsi et al. (2007). The resultant spectra, focusing on the observed central three groups of hyperfine components, are displayed in Figure F1 (solid histogram). Both spectra are in good agreement with the VLA observations, taking into account their 15% flux uncertainty. Also, we are not simulating interferometric-only observations, as reported in Crapsi et al. (2007), as we take into account all the flux from the various scales, including the extended envelope, filtered out by the VLA observations. This point will be made clear in A. Schmiedeke et al. (2022, in preparation), where a higher-sensitivity map of the NH₃(1,1) line obtained by combining VLA with GBT data will be presented. Therefore, we conclude that our model, where NH₃ presents a sharp drop in abundance toward the prestellar core center, is consistent with VLA observations. We also considered a case where the NH₃ fractional abundance does not drop in the central 6000 au (unlike in our models, shown in Figure 5), and instead, it maintains a constant value of 2 × 10⁻⁸. The resultant spectra are also in Figure F1 (see dashed histograms); they are almost identical to the solid histograms, underlining the fact that NH₃(1,1) data cannot be used to constrain the abundance profile of NH₃ within the central region of the core. This is due to the relatively low critical density of NH₃(1,1) (≃4 × 10³ cm⁻³; Maret et al. 2009), which then makes the (1,1) inversion transition mainly sensitive to the outer parts of the prestellar core, hiding abundance variations in the central few thousand astronomical units. In a sense, this result is similar to that found for the pNH₂D(1₁₁−1₀₁) line, with the important difference that pNH₂D(1₁₁−1₀₁) has a significantly higher critical density than NH₃(1,1), and then it is more sensitive to variations in the NH₂D fractional abundance within the central regions (as demonstrated by Figure C1).

Zoom In Zoom Out Reset image size