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

Paola Caselli, Jaime E. Pineda, Olli Sipilä, Bo Zhao, Elena Redaelli, Silvia Spezzano, Maria José Maureira, Felipe Alves, Luca Bizzocchi, Tyler L. Bourke, et al.

► To cite this version:

Paola Caselli, Jaime E. Pineda, Olli Sipilä, Bo Zhao, Elena Redaelli, et al.. The Central 1000 au of a Prestellar Core Revealed with ALMA. II. Almost Complete Freeze-out. The Astrophysical Journal, 2022, 929, 10.3847/1538-4357/ac5913 . insu-03672028

HAL Id: insu-03672028

https://insu.hal.science/insu-03672028

Submitted on 19 May 2022

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The Central 1000 au of a Prestellar Core Revealed with ALMA. II. Almost Complete Freeze-out

Paola Caselli<fnref n="1" />, Jaime E. Pineda<fnref n="1" />, Olli Sipilä<fnref n="1" />, Bo Zhao<fnref n="1" />,<fnref n="2" />, Elena Redaelli<fnref n="1" />, Silvia Spezzano<fnref n="1" />, Maria José Maureira<fnref n="1" />, Felipe Alves<fnref n="1" />, Luca Bizzocchi<fnref n="1" />,<fnref n="3" />,<fnref n="4" />, Tyler L. Bourke<fnref n="5" />, Ana Chacón-Tanarro<fnref n="6" />, Rachel Friesen<fnref n="7" />, Daniele Galli<fnref n="8" />, Jorma Harju<fnref n="1" />,<fnref n="9" />, Izaskun Jiménez-Serra<fnref n="10" />, Eric Keto<fnref n="11" />, Zhi-Yun Li<fnref n="12" />, Marco Padovani<fnref n="8" />, Anika Schmiedeke<fnref n="1" />, Mario Tafalla<fnref n="6" />, and Charlotte Vastel<fnref n="13" />

Received 2021 October 2; revised 2022 February 25; accepted 2022 February 26; published 2022 April 8

## Abstract

Prestellar cores represent the initial conditions in the process of star and planet formation. Their low temperatures ($<10$ K) allow the formation of thick icy dust mantles, which will be partially preserved in future protoplanetary disks, ultimately affecting the chemical composition of planetary systems. Previous observations have shown that carbon- and oxygen-bearing species, in particular CO, are heavily depleted in prestellar cores due to the efficient molecular freeze-out onto the surface of cold dust grains. However, N-bearing species such as $NH_3$ and, in particular, its deuterated isotopologues appear to maintain high abundances where CO molecules are mainly in the solid phase. Thanks to ALMA, we present here the first clear observational evidence of $NH_2D$ freeze-out toward the L1544 prestellar core, suggestive of the presence of a “complete depletion zone” within a ≈1800 au radius, in agreement with astrochemical prestellar core model predictions. Our state-of-the-art chemical model coupled with a non-LTE radiative transfer code demonstrates that $NH_2D$ becomes mainly incorporated in icy mantles in the central 2000 au and starts freezing out already at ≈7000 au. Radiative transfer effects within the prestellar core cause the $NH_2D(l_{11} - l_{10})$ emission to appear centrally concentrated, with a flattened distribution within the central ≈3000 au, unlike the 1.3 mm dust continuum emission, which shows a clear peak within the central ≈1800 au. This prevented $NH_2D$ freeze-out from being detected in previous observations, where the central 1000 au cannot be spatially resolved.

*Unified Astronomy Thesaurus concepts: Interstellar medium (847); Interstellar molecules (849); Astrochemistry (75); Radiative transfer (1335)*

## 1. Introduction

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^5$ cm$^{-3}$ (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 ($\sim$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_2H^+$, $NH_3$, 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_2H^+$ (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_2D^+$, coupled with astrochemical and non-LTE radiative transfer modeling, have provided evidence of $N_2D^+$ 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_2H^+$ isotopologues is linked to $N_2$, which has a binding energy similar to that of CO (Bisschop et al. 2006), although recent studies have shown that

<figure><img src="image_2.png" /></figure>

<footnote n="1">¹Max-Planck-Institut für extraterrestrische Physik, Gießenbachstrasse 1, D-85748 Garching bei München, Germany; caselli@mpe.mpg.de</footnote>

<footnote n="2">²Department of Physics and Astronomy, McMaster University, 1280 Main St. W, Hamilton, ON, L8S 4M1, Canada</footnote>

<footnote n="3">³Dipartimento di Chimica “Giacomo Ciamician”, Università di Bologna, Via F. Selmi 2, I-40126 Bologna, Italy</footnote>

<footnote n="4">⁴Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy</footnote>

<footnote n="5">⁵SKA Observatory, Jodrell Bank, Lower Withington, Macclesfield, SK11 9FT, UK</footnote>

<footnote n="6">⁶Observatorio Astronómico Nacional (OAG-IGN), Alfonso XII 3, E-28014, Madrid, Spain</footnote>

<footnote n="7">⁷Department of Astronomy & Astrophysics, University of Toronto, 50 St. George St., Toronto, ON, M5S 3H4, Canada</footnote>

<footnote n="8">⁸INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy</footnote>

<footnote n="9">⁹Department of Physics, P.O. BOX 64, FI-00014 University of Helsinki, Finland</footnote>

<footnote n="10">¹⁰Centro de Astrobiología (CSIC-INTA), Ctra. de Torrejón a Ajalvir km 4, E-28850, Torrejón de Ardoz, Spain</footnote>

<footnote n="11">¹¹Harvard-Smithsonian Center for Astrophysics, 160 Garden Street, Cambridge, MA 02420, USA</footnote>

<footnote n="12">¹²Department of Astronomy, University of Virginia, Charlottesville, VA 22904, USA</footnote>

<footnote n="13">¹³IRAP, Université de Toulouse, CNRS, CNES, UPS, Toulouse, France</footnote>
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 5′′8 × 4′′5 (986 × 765 au, at the distance of 170 pc;<fnref n="14" /> Galli et al. 2019). From this result, it was concluded that no complete depletion zone<fnref n="15" /> (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 2′′5 (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. Observations

### 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, band-pass, 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 2′′5. 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.

### 2.2. 1.3 mm Continuum

We observed L1544 with ALMA during Cycle 2 (ESO Project ID 2013.1.01195.S; PI: Caselli) using both the 12 m (Main array) and the 7 m (ACA) arrays. The data were already presented in Caselli et al. (2019), but here we briefly describe the data reduction process. The 12 m array observations were carried out on 2014 December 27 and 29, while the 7 m array observations were carried out on 2014 June 15, July 20 and 29, and August 6 and 11. The single-pointing observations used a correlator with the continuum centered at 228.973 GHz and a 2 GHz bandwidth. The 12 m array data were calibrated using CASA version 4.2.2, while the 7 m array data were calibrated using CASA version 4.2.1. We run the statwt command on both data sets before the combination. We perform a joint deconvolution of both data sets (ACA and 12 m array) in CASA (version 5.6.2) using tclean multifrequency synthesis (mfs) for mosaics, a natural weighting, a taper of 2′′1 × 1′′45 and a PA of 147°2, and restored with a circular beam of 2′′5. Because the continuum emission is extended, we use the multiscale clean technique, with the multiscale parameter of [0, 1.25, 3.75, 11.25]′. The 1.3 mm dust continuum emission map can be seen in the bottom panel of Figure 1, overlaid onto the pNH₂D column density map, which will be presented in Section 3.

## 3. Column Density Calculation

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

<footnote n="14">¹⁴ We note that previous papers on L1544 adopted a distance to the source of 140 pc or 135 pc, based on previous distance measurements (Schlafly et al. 2014; Galli et al. 2018). Here we adjust previous results to the new distance estimate.</footnote>

<footnote n="15">¹⁵ A “complete depletion zone” is defined as a region within a prestellar core where, due to freeze-out, species heavier than He have left the gas phase and reside in the icy mantles of dust grains.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Top: integrated intensity map of the pNH₂D(1₁₁-1₀₁) transition toward L1544. Bottom: pNH₂D column density map derived from fitting the hyperfine structure. The contours show the 14, 17, and 20 × σ levels ($\sigma = 36 \mu$Jy beam⁻¹) of the 1.3 mm dust continuum emission map presented in Caselli et al. (2019), which locates the kernel. The scale bar and beam size are shown at the bottom-right and -left corners, respectively. In the top left of the bottom panel, we show the orientation and shape of the ellipse used to define the projected radius in Section 3 and Figure 2. The typical uncertainty on the pNH₂D column density is about 10%.</figcaption></figure>

laboratory works of Melosso et al. (2020, 2021). The model assumes a Gaussian velocity distribution (with centroid velocity and velocity dispersion $V_{lsr}$ and $\sigma_v$), equal excitation temperatures ($T_{ex}$) for all hyperfine components, and total optical depth ($\tau_0$). 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 $\sigma_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):

$$N(\text{pNH}_2\text{D}) = \frac{8\pi\nu^3 Q(T_{\text{ex}})}{c^3 A_{\text{ul}} g_u e^{-E_u/kT_{\text{ex}}}} \times [e^{h\nu/kT_{\text{ex}}} - 1]^{-1} \int \tau dv, \quad (1)$$

where $\int \tau dv = \sqrt{2\pi} \tau_0 \sigma_v$, $\nu$ is the frequency of the transition, $g_u$ is the upper-level degeneracy, $E_u$ is the upper energy level of

<figure><img src="image_4.png" /><figcaption>**Figure 2.** pNH₂D column density as a function of projected radius (see the main text for its definition). Blue circles correspond to individual measurements across the column density map in the bottom panel of Figure 1. The black and gray horizontal lines are the mean and corresponding uncertainty within a 2"/5 bin. Note the flat column density profile within the central 10', where the 1.3 mm dust continuum map shows a clear peak (see Figure 1, bottom panel).</figcaption></figure>

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),<fnref n="16" /> 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 2"/5 bin. The projected radius *pr* is calculated as $pr = \sqrt{r_{\text{maj}}^2 + (r_{\text{min}} b_{\text{axes}})^2}$, where $r_{\text{maj}}$ and $r_{\text{min}}$ are the distances along the semimajor and -minor axes, respectively; $b_{\text{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*).

<footnote n="16">¹⁶ https://home.strw.leidenuniv.nl/~moldata/datafiles/p-nh2d.dat</footnote>
<table><caption>**Table 1**
Initial Fractional Abundances (with Respect to the Total H Nuclei, $n_H$)</caption><thead><tr><th>Species</th><th>Abundance</th></tr></thead><tbody><tr><td>H<sub>2</sub></td><td>5.00 × 10<sup>-1a</sup></td></tr><tr><td>He</td><td>9.00 × 10<sup>-2</sup></td></tr><tr><td>C<sup>+</sup></td><td>1.20 × 10<sup>-4</sup></td></tr><tr><td>N</td><td>7.60 × 10<sup>-5</sup></td></tr><tr><td>O</td><td>2.56 × 10<sup>-4</sup></td></tr><tr><td>S<sup>+</sup></td><td>8.00 × 10<sup>-8</sup></td></tr><tr><td>Si<sup>+</sup></td><td>8.00 × 10<sup>-9</sup></td></tr><tr><td>Na<sup>+</sup></td><td>2.00 × 10<sup>-9</sup></td></tr><tr><td>Mg<sup>+</sup></td><td>7.00 × 10<sup>-9</sup></td></tr><tr><td>Fe<sup>+</sup></td><td>3.00 × 10<sup>-9</sup></td></tr><tr><td>P<sup>+</sup></td><td>2.00 × 10<sup>-10</sup></td></tr><tr><td>Cl<sup>+</sup></td><td>1.00 × 10<sup>-9</sup></td></tr></tbody></table>

**Note.**

## 4. Chemical and Radiative Transfer Modeling

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.<fnref n="17" /> 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}$.

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

<footnote>a The initial H₂ ortho/para ratio is 1 × 10⁻³.</footnote>

<footnote n="17">¹⁷ The six observed spectra have been extracted from the six 2''5 circular areas, along the major axis and toward the southeast direction, displayed in Figure A1. We also compared with the six spectra extracted from the perpendicular direction along the minor axis (see Figure A1) and found similar results within the central 1800 au of the kernel (see Appendix A and Figure A2), while the outer positions show some deviation, suggesting that the (spherically symmetric) model is in better agreement with the physical structure along the major axis (see also Appendix D).</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Left: map of the pNH₂D column density distribution predicted by the chemical model at the best-fit time. The column densities have been convolved to a beam of 2''5 (425 au). The black circles with numbers are the areas within which the red spectra in the right panel, with the corresponding numbers, have been extracted. Right: simulated pNH₂D (1₁₁–1₀₁) line emission profiles (red) at ~2''5 intervals from the core center. The observed lines toward equivalent positions in the southeast direction in Figure A1, parallel to the major axis, are shown in black.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** pNH₂D(1₁₁–1₀₁) excitation temperature ($T_{ex}$) as a function of projected radius (see Section 3 for its definition). Blue circles correspond to individual measurements. The black and gray horizontal lines are the mean and corresponding uncertainty within a 2''5 bin. The red curve is the model $T_{ex}$ profile, the predicted excitation temperature as a function of core radius, $T_{ex}(r)$.</figcaption></figure>

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.

## 5. Discussion: Almost Complete Freeze-out in the Kernel

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 (5''8 × 4''5), 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_{D_{tot}}$, defined as

$$f_{D_{tot}} = \frac{\sum_{i}[x(i)_{gas} + x(i)_{solid}]}{\sum_{i}x(i)_{gas}}, \qquad (2)$$

where i represents any species in the chemical model containing an element heavier than He. Thus, $f_{D_{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_{D_{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
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Left: model fractional abundance of pNH₂D (blue) and overall gas-phase depletion factor (red) at the best-fit time as a function of radius. Right: fractional abundances w.r.t. H₂ of normal and deuterated ammonia (colors, indicated in the figure), as well as of CO (multiplied by 0.01, black), at the best-fit time as a function of distance from the core center. The ammonia fractional abundances represent sums over the nuclear spin forms (ortho and para; also meta for ND₃).</figcaption></figure>

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.

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_{D_{\text{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 4.′′35 × 3.′′45. 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 22.′′3 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).

## 6. Conclusions

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.

## Appendix A: Full Mosaic Map

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.
<figure><img src="image_8.png" /><figcaption>**Figure A1.** Integrated intensity map of the pNH₂D (1₁₁–1₀₁) transition on L1544, showing the three-point mosaic. The white circles overlaid on the integrated intensity map mark the positions where the spectra of Figure 3 (southeast direction, along the major axis) and Figure A2 (southwest direction, along the minor axis) have been extracted.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure A2.** Same as the right panel of Figure 3 but now the black histograms are the spectra extracted from the six circles in the southwest direction in Figure A1 (along the L1544 minor axis). Note that the observed spectra in the outer positions (numbers 5 and 6) are now less intense than those predicted by our model. This suggests that the major axis of L1544 has a physical and chemical structure in better agreement with our 1D model.</figcaption></figure>

# Appendix B

## Comparison with the Single-dish NH₂D(1₁₁–1₀₁) Spectrum

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
<figure><img src="image_10.png" /><figcaption>**Figure B1.** Comparison between the pNH₂D(1₁₁–1₀₁) IRAM 30 m spectrum (red histogram) and the ALMA 12 m + ACA spectrum within the same 22" beam of the 30 m observations (black histogram), obtained toward the L1544 kernel. This clearly shows that the ALMA 12 m + ACA combined mosaic fully recovers the line flux.</figcaption></figure>

with the combination of 12 m and ACA data, provides a powerful technique to study prestellar cores in detail.

## Appendix C
Modeled Spectra in the Case of No Central Freeze-out

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.

<figure><img src="image_11.png" /><figcaption>**Figure C1.** Same as Figure 3, right panels, but assuming that the abundance of pNH₂D is constant within the central 2000 au instead of dropping as predicted by our chemical models (see Figure 5). The modeled lines clearly overestimate the observed spectra in the central four positions, demonstrating that freeze-out is needed.</figcaption></figure>
## Appendix D
### Density and Velocity Profiles

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^7$ cm$^{-3}$, in agreement with the peak density at the center of the BE sphere, while the average density within the kernel is about $10^6$ cm$^{-3}$, 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.

<figure><img src="image_12.png" /></figure>

**Figure D1.** Comparison between the H₂ number density (left panel) and velocity (right panel) radial profiles of the BE sphere in quasi-static contraction from the HD simulations of Keto et al. (2015), after $6 \times 10^5$ yr of evolution of an initial BE sphere with central density $1 \times 10^4$ cm$^{-3}$, radius $2 \times 10^5$ au, and total mass $10 M_{\odot}$ (dashed red curves), and those along the major (black curves) and minor (blue curves) axes of the niMHD simulation of a magnetized contracting core, after $1.4 \times 10^6$ yr of evolution starting from a uniform $8.1 M_{\odot}$ cloud with volume density $2 \times 10^4$ cm$^{-3}$, radius $5 \times 10^4$ au, and initial mass-to-flux ratio of 1, which best fits the ALMA 1.3 mm dust continuum emission (see Caselli et al. 2019 for details). Note the overall similarities despite the very different simulations (HD vs. niMHD) and different initial conditions.
Appendix E

Model Predictions Assuming a 5" × 5" Beam

Crapsi et al. (2007) observed the same transition of
deuterated ammonia using the IRAM Plateau de Bure
interferometer with an angular resolution of 5′′8 × 4′′5 . To
demonstrate that this is not sufficient to detect the pNH$_{2}$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 2′′5). 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 2′′5
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.

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

$t = 2.5 \times 10^5 \text{yr}$

<figure><img src="image_13.png" /><figcaption>**Figure E1.** Simulated spectra of the same L1544 core model presented in Section 4 and Figure 3 but now assuming a telescope beam of 5" × 5", similar to the one used by Crapsi et al. (2007). The number in the top left of each panel is the integrated intensity in K km s⁻¹, while the number in the top right shows the corresponding position across the strip labeled in Figure 3, left panel, but with a spacing of 5" instead of 2".5. Although the simulated spectra look similar to those in Figure 3, the difference in integrated intensity between the spectrum extracted toward the core center (position 1) and the one 5" away (position 2, the integrated intensity peak) is very small. Thus, the flattening zone cannot be detected with the ~5" beam of the Crapsi et al. (2007) observations.</figcaption></figure>
$p\text{NH}_2\text{D}(1_{11}-1_{01})$ 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 $p\text{NH}_2\text{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.

## Appendix F
### $\text{NH}_3(1,1)$ Simulated Spectra

Given that previous work has claimed that the fractional abundance of $\text{NH}_3$ is actually increasing toward the L1544 center (e.g., Crapsi et al. 2007), we would like to test if our model predictions about the $\text{NH}_3$ fractional abundance across the core (see Figure 5) is consistent with previous interferometric work. For this, we simulate VLA $\text{NH}_3(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, 4"$34 \times 3"$45), 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 $\text{NH}_3(1,1)$ line obtained by combining VLA with GBT data will be presented. Therefore, we conclude that our model, where $\text{NH}_3$ presents a sharp drop in abundance toward the prestellar core center, is consistent with VLA observations. We also considered a case where the $\text{NH}_3$ 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 \times 10^{-8}$. The resultant spectra are also in Figure F1 (see dashed histograms); they are almost identical to the solid histograms, underlining the fact that $\text{NH}_3(1,1)$ data cannot be used to constrain the abundance profile of $\text{NH}_3$ within the central region of the core. This is due to the relatively low critical density of $\text{NH}_3(1,1)$ ($\approx 4 \times 10^3 \text{ cm}^{-3}$; 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 $p\text{NH}_2\text{D}(1_{11}-1_{01})$ line, with the important difference that $p\text{NH}_2\text{D}(1_{11}-1_{01})$ has a significantly higher critical density than $\text{NH}_3(1,1)$, and then it is more sensitive to variations in the $\text{NH}_2\text{D}$ fractional abundance within the central regions (as demonstrated by Figure C1).

<figure><img src="image_14.png" /></figure>

**Figure F1.** Simulated $\text{NH}_3(1,1)$ spectra toward the dust peak (top panel) and an offset 22"/3 (bottom panel) away from the center, to be compared with the two spectra in Figure 2 of Crapsi et al. (2007). A beam of 4" (close to the angular resolution of the VLA observations of Crapsi et al. 2007) has been assumed. The solid blue histogram is the spectrum obtained with the $\text{NH}_3$ abundance profile predicted by our model and shown in Figure 5, while the dashed blue histogram (only visible in the top panel, as in the bottom panel it coincides with the solid histogram) is the spectrum predicted assuming $\text{NH}_3$ constant abundance within the central 5000 au. Both spectra reproduce well the observed ones, within the observational uncertainties. It is then clear that $\text{NH}_3(1,1)$ interferometric spectra of L1544 cannot provide constraints on the fractional abundance distribution of ammonia, unlike ALMA observations of $\text{NH}_2\text{D}(1_{11}-1_{01})$ (see Figure C1).

## ORCID iDs

Paola Caselli https://orcid.org/0000-0003-1481-7911
Jaime E. Pineda https://orcid.org/0000-0002-3972-1978
Olli Sipilä https://orcid.org/0000-0002-9148-1625
Elena Redaelli https://orcid.org/0000-0002-0528-8125
Silvia Spezzano https://orcid.org/0000-0002-6787-5245
Maria José Maureira https://orcid.org/0000-0002-7026-8163
Felipe Alves https://orcid.org/0000-0002-7945-064X
Luca Bizzocchi https://orcid.org/0000-0002-9953-8593
Tyler L. Bourke https://orcid.org/0000-0001-7491-0048
Rachel Friesen https://orcid.org/0000-0001-7594-8128
Daniele Galli https://orcid.org/0000-0001-7706-6049
Jorma Harju https://orcid.org/0000-0002-1189-9790
Izaskun Jiménez-Serra https://orcid.org/0000-0003-4493-8714
Zhi-Yun Li https://orcid.org/0000-0002-7402-6487
Marco Padovani https://orcid.org/0000-0003-2303-0096
Anika Schmiedeke https://orcid.org/0000-0002-1730-8832
Mario Tafalla https://orcid.org/0000-0002-2569-1253

## References

* Altwegg, K., Balsiger, H., & Fuselier, S. A. 2019, ARA&A, 57, 113

* André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (Tucson, AZ: Univ. Arizona Press), 27

* Bergin, E. A., Alves, J., Huard, T., & Lada, C. J. 2002, ApJL, 570, L101

* Bergin, E. A., & Tafalla, M. 2007, ARA&A, 45, 339

* Bisschop, S. E., Fraser, H. J., Öberg, K. I., van Dishoeck, E. F., & Schlemmer, S. 2006, A&A, 449, 1297

* Bizzocchi, L., Caselli, P., Leonardo, E., & Dore, L. 2013, A&A, 555, A109

* Bonnor, W. B. 1956, MNRAS, 116, 351

* Caselli, P., Bizzocchi, L., Keto, E., et al. 2017, A&A, 603, L1
Caselli, P., & Ceccarelli, C. 2012, A&ARv, 20, 56

Caselli, P., Keto, E., Bergin, E. A., et al. 2012, ApJL, 759, L37

Caselli, P., Pineda, J. E., Zhao, B., et al. 2019, ApJ, 874, 89

Caselli, P., Walmsley, C. M., Tafalla, M., Dore, L., & Myers, P. C. 1999, ApJL, 523, L165

Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002, ApJ, 565, 344

Ceccarelli, C., Caselli, P., Bockelée-Morvan, D., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 859

Chacón-Tanarro, A., Pineda, J. E., Caselli, P., et al. 2019, A&A, 623, A118

Cleeves, L. I., Bergin, E. A., Alexander, C. M. O. D., et al. 2014, Sci, 345, 1590

Cohen, E. A., & Pickett, H. M. 1982, JMoSp, 93, 83

Crapsi, A., Caselli, P., Walmsley, C. M., et al. 2005, ApJ, 619, 379

Crapsi, A., Caselli, P., Walmsley, M. C., & Tafalla, M. 2007, A&A, 470, 221

Daniel, F., Coudert, L. H., Punanova, A., et al. 2016, A&A, 586, L4

Daniel, F., Faure, A., Wiesenfeld, L., et al. 2014, MNRAS, 444, 2544

De Lucia, F. C., & Helminger, P. 1975, JMoSp, 54, 200

Draine, B. T. 1985, in Protostars and Planets II, ed. D. C. Black & M. S. Matthews (Tucson, AZ: Univ. Arizona Press), 621

Drozdovskaya, M. N., Schroeder, I. R. H. G., Rubin, M., et al. 2021, MNRAS, 500, 4901

Ebert, R. 1955, ZAp, 37, 217

Flower, D. R., Pineau Des Forêts, G., & Walmsley, C. M. 2006, A&A, 456, 215

Forbrich, J., Lada, C. J., Lombardi, M., Román-Zúñiga, C., & Alves, J. 2015, A&A, 580, A114

Galli, P. A. B., Loinard, L., Bouy, H., et al. 2019, A&A, 630, A137

Galli, P. A. B., Loinard, L., Ortiz-Léon, G. N., et al. 2018, ApJ, 859, 33

Ginsburg, A., & Mirocha, J. 2011, PySpecKit: Python Spectroscopic Toolkit, ascl:1109.001

Harju, J., Daniel, F., Sipilä, O., et al. 2017, A&A, 600, A61

Hily-Blant, P., Faure, A., Rist, C., Pineau des Forêts, G., & Flower, D. R. 2018, MNRAS, 477, 4454

Hily-Blant, P., Walmsley, M., Pineau Des Forêts, G., & Flower, D. 2010, A&A, 513, A41

Juárez, C., Girart, J. M., Frau, P., et al. 2017, A&A, 597, A74

Juvela, M. 2020, A&A, 644, A151

Keto, E., & Caselli, P. 2008, ApJ, 683, 238

Keto, E., & Caselli, P. 2010, MNRAS, 402, 1625

Keto, E., Caselli, P., & Rawlings, J. 2015, MNRAS, 446, 3731

Launhardt, R., Stutz, A. M., Schmiedeke, A., et al. 2013, A&A, 551, A98

Mangum, J. G., & Shirley, Y. L. 2015, PASP, 127, 266

Maret, S., Faure, A., Scifoni, E., & Wiesenfeld, L. 2009, MNRAS, 399, 425

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Melosso, M., Bizzocchi, L., Dore, L., et al. 2021, JMoSp, 377, 111431

Melosso, M., Dore, L., Gauss, J., & Puzzarini, C. 2020, JMoSp, 370, 111291

Mumma, M. J., & Charnley, S. B. 2011, ARA&A, 49, 471

Nguyen, T., Baouche, S., Congiu, E., et al. 2018, A&A, 619, A111

Öberg, K. I., Guzmán, V. V., Furuya, K., et al. 2015, Natur, 520, 198

Ormel, C. W., Paszun, D., Dominik, C., & Tielens, A. G. G. M. 2009, A&A, 502, 845

Pagani, L., Bacmann, A., Cabrit, S., & Vastel, C. 2007, A&A, 467, 179

Pagani, L., Bourgoin, A., & Lique, F. 2012, A&A, 548, L4

Pineda, J. E., Schmiedeke, A., Caselli, P., et al. 2021, ApJ, 912, 7

Redaelli, E., Bizzocchi, L., Caselli, P., et al. 2019, A&A, 629, A15

Redaelli, E., Sipilä, O., Padovani, M., et al. 2021, A&A, 656, A109

Schlafly, E. F., Green, G., Finkbeiner, D. P., et al. 2014, ApJ, 786, 29

Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369

Sipilä, O., Caselli, P., Redaelli, E., Juvela, M., & Bizzocchi, L. 2019, MNRAS, 487, 1269

Sipilä, O., Zhao, B., & Caselli, P. 2020, A&A, 640, A94

Spezzano, S., Caselli, P., Bizzocchi, L., Giuliano, B. M., & Lattanzi, V. 2017, A&A, 606, A82

Tafalla, M., Myers, P. C., Caselli, P., & Walmsley, C. M. 2004, A&A, 416, 191

Tafalla, M., Myers, P. C., Caselli, P., Walmsley, C. M., & Comito, C. 2002, ApJ, 569, 815

van der Tak, F. F. S., Lique, F., Faure, A., Black, J. H., & van Dishoeck, E. F. 2020, Atoms, 8, 15

van Dishoeck, E. F., Bergin, E. A., Lis, D. C., & Lunine, J. I. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 835

Vasyunin, A. I., Caselli, P., Dulieu, F., & Jiménez-Serra, I. 2017, ApJ, 842, 33

Walmsley, C. M., Flower, D. R., & Pineau des Forêts, G. 2004, A&A, 418, 1035