# Mapping the H 2 D + and N 2 H + emission toward prestellar cores. Testing dynamical models of the collapse using gas tracers

E. Koumpia, L. Evans, J. Di Francesco, F. van der Tak, R. Oudmaijer

► To cite this version:

E. Koumpia, L. Evans, J. Di Francesco, F. van der Tak, R. Oudmaijer. Mapping the H 2 D + and N 2 H + emission toward prestellar cores. Testing dynamical models of the collapse using gas tracers. Astronomy & Astrophysics - A&A, 2020, 643, pp.A61. 10.1051/0004-6361/202038457. hal-04895985

HAL Id: hal-04895985

https://hal.science/hal-04895985v1

Submitted on 20 Jan 2025

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
# Mapping the H₂D⁺ and N₂H⁺ emission toward prestellar cores. Testing dynamical models of the collapse using gas tracers*

E. Koumpia<fnref n="1" />, L. Evans<fnref n="1" />,<fnref n="2" />,<fnref n="3" />, J. Di Francesco<fnref n="4" />, F. F. S. van der Tak<fnref n="5" />,<fnref n="6" />, and R. D. Oudmaijer<fnref n="1" />

Received 20 May 2020 / Accepted 18 September 2020

## ABSTRACT

**Context.** The study of prestellar cores is critical as they set the initial conditions in star formation and determine the final mass of the stellar object. To date, several hypotheses have described their gravitational collapse. Deriving the dynamical model that fits both the observed dust and the gas emission from such cores is therefore of great importance.

**Aims.** We perform detailed line analysis and modeling of H₂D⁺ 1₀⁻<fnref n="1" />₁ and N₂H⁺ 4⁻3 emission at 372 GHz, using 2' × 2' maps (James Clerk Maxwell Telescope; JCMT). Our goal is to test the most prominent dynamical models by comparing the modeled gas kinematics and spatial distribution (H₂D⁺ and N₂H⁺) with observations toward four prestellar (L1544, L183, L694-2, L1517B) and one protostellar core (L1521f).

**Methods.** We fit the line profiles at all offsets showing emission using single Gaussian distributions. We investigate how the line parameters (V<sub>LSR</sub>, FWHM and T<sub>A</sub>) change with offset to examine the velocity field, the degree of nonthermal contributions to the line broadening, and the distribution of the material in these cores. To assess the thermal broadening, we derive the average gas kinetic temperature toward all cores using the non-LTE radiative transfer code RADEX. We perform a more detailed non-LTE radiative transfer modeling using RATRAN, where we compare the predicted spatial distribution and line profiles of H₂D⁺ and N₂H⁺ with observations toward all cores. To do so, we adopt the physical structure for each core predicted by three different dynamical models taken from literature: quasi-equilibrium Bonnor–Ebert sphere (QE-BES), singular isothermal sphere (SIS), and Larson–Penston (LP) flow. In addition, we compare these results to those of a static sphere, whose density and temperature profiles are based on the observed dust continuum. Lastly, we constrain the abundance profiles of H₂D⁺ and N₂H⁺ toward each core.

**Results.** We find that variable nonthermal contributions (variations by a factor of 2.5) are required to explain the observed line width of both H₂D⁺ and N₂H⁺, while the nonthermal contributions are found to be 50% higher for N₂H⁺. The RADEX modeling results in average core column densities of ~9 × 10¹<fnref n="2" /> cm⁻<fnref n="2" /> for H₂D⁺ and N₂H⁺. The LP flow seems to be the dynamical model that can reproduce the observed spatial distribution and line profiles of H₂D⁺ on a global scale of prestellar cores, while the SIS model systematically and significantly overestimates the width of the line profiles and underestimates the line peak intensity. We find similar abundance profiles for the prestellar cores and the protostellar core. The typical abundances of H₂D⁺ vary between 10⁻⁹ and 10⁻¹⁰ for the inner 5000 au and drop by about an order of magnitude for the outer regions of the core (2 × 10⁻¹⁰–6 × 10⁻¹¹). In addition, a higher N₂H⁺ abundance by about a factor of 4 compared to H₂D⁺ is found toward the two cores with detected emission. The presence of N₂H⁺ 4–3 toward the protostellar core and toward one of the prestellar cores reflects the increasing densities as the core evolves.

**Conclusions.** Our analysis provides an updated picture of the physical structure of prestellar cores. Although the dynamical models account for mass differences by up to a factor of 7, the velocity structure drives the shape of the line profiles, allowing for a robust comparison between the models. We find that the SIS model can be clearly excluded in explaining the gas emission toward the cores, but a larger sample is required to differentiate clearly between the LP flow, the QE-BES, and the static models. All models of collapse underestimate the intensity of the gas emission by up to several factors toward the only protostellar core in our sample, indicating that different dynamics take place in different evolutionary core stages. If the LP model is confirmed toward a larger sample of prestellar cores, it would indicate that they may form by compression or accretion of gas from larger scales. If the QE-BES model is confirmed, it means that quasi-hydrostatic cores can exist within turbulent ISM.

**Key words.** stars: formation – ISM: clouds – ISM: kinematics and dynamics – submillimeter: ISM

<footnote n="1">¹ School of Physics & Astronomy, University of Leeds, Woodhouse Lane, LS2 9JT Leeds, UK
e-mail: e.koumpia@leeds.ac.uk</footnote>

<footnote n="2">² IRAP, Université de Toulouse, CNRS, UPS, CNES, 31400 Toulouse, France</footnote>

<footnote n="3">³ INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Florence, Italy</footnote>

<footnote n="4">⁴ Herzberg Astronomy and Astrophysics Research Centre, National Research Council of Canada, Victoria, BC, Canada</footnote>

<footnote n="5">⁵ SRON Netherlands Institute for Space Research, Landleven 12, 9747 AD Groningen, The Netherlands</footnote>

<footnote n="6">⁶ Kapteyn Astronomical Institute, University of Groningen, Groningen, The Netherlands</footnote>

<footnote>* The reduced H₂D⁺ and N₂H⁺ JCMT data cubes of all cores are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/643/A61</footnote>
# 1. Introduction

Star formation begins within molecular clouds, where magnetic fields and turbulence dominate the dynamics, leading to the formation of filamentary structures of gas (Arzoumanian et al. 2011; Pudritz & Kevlahan 2013; André 2017). In particular, Herschel observations have revealed that prestellar cores and protostellar objects form in thermally critical and supercritical filaments (André et al. 2010; Tafalla & Hacar 2015). Understanding the physical and chemical processes that take place in these very early stages of star formation is of high importance. Not only do such cores set the initial conditions of star formation and determine the final mass of the star (Bergin & Tafalla 2007), but they also have a strong influence on the multiplicity (Pineda et al. 2015). Although low-mass star formation is better understood than high-mass star formation, the initial conditions of the collapse remain uncertain. Prestellar cores are known to contract under gravitational forces (Krumholz et al. 2005), while thermal and magnetic pressure and the presence of turbulence prevent the core from collapsing (Goodman et al. 1998).

The most broadly known dynamical models of a spherically symmetric collapsing core are: i) the Larson–Penston flow (LP flow; Larson 1969; Penston 1969), where at t~0 the inflow velocity at large radii reaches supersonic values (~3.3 × speed of sound), and therefore is far from equilibrium, and the density profile follows a power-law; ii) the singular isothermal sphere (SIS) evolving via an “inside-out” collapse (Shu 1977); and iii) the hydrostatic Bonnor–Ebert sphere (BES) supported by turbulence, which, when disturbed, can lead to a contraction in quasi-equilibrium (QE; Broderick & Keto 2010). Constraining the dynamical structure of prestellar cores observationally is very challenging. Hence, detailed studies focusing on the gas emission, as well as the age determination of observed prestellar cores, are crucial to distinguish between the proposed models.

Ward-Thompson et al. (1994) presented the first submillimeter continuum survey of a sample of cores without associated infrared emission shortward of 100 µm, therefore presumably starless cores, reporting their first detection using longer wavelengths (>450 µm). The lifetime of prestellar cores has been observationally estimated to be ~10⁴–10⁶ yr (e.g., L1544, L694-2, L183 Beichman et al. 1986; Caselli et al. 2008; Enoch et al. 2008). The prestellar cores are characterized by very low temperatures (T ≈ 10 K) that increase outward, and high densities (>10⁵ cm⁻³; Keto & Caselli 2008) that show a flattened density profile in the center, resulting in the unique gas chemistry presented in the current paper. Various studies have shown a strong correlation between CO depletion and the degree of deuteration of hydrogen-based species in prestellar cores (Caselli et al. 2003; Vastel et al. 2006). Below the critical temperature of T ≈ 25 K, gas-phase species, including almost all CNO species, are frozen out (≥98% Caselli et al. 2003) to dust grains and therefore are depleted (Caselli et al. 2008). The enrichment in H₂D⁺ and N₂H⁺ under those conditions can be understood if one takes a closer look at the production and destruction mechanisms of the relevant species (Eq. (1) to (4)):

$$H_3^+ + HD \rightleftharpoons H_2D^+ + H_2 + \Delta E, \quad (1)$$

$$H_3^+ + N_2 \rightarrow N_2H^+ + H_2, \quad (2)$$

$$H_3^+ + CO \rightarrow HCO^+ + H_2, \quad (3)$$

$$N_2H^+ + CO \rightarrow HCO^+ + N_2. \quad (4)$$

As soon as CO returns to the gas-phase (T > 20–30 K, i.e., after the formation of a protostar), the abundance of both H₂D⁺ and N₂H⁺ will subsequently decrease. We note, however, that

$N_2H^+$ survives with an abundance independent of the distance from the center (Pagani et al. 2012; Lique et al. 2015) and appears to remain longer (Tafalla et al. 2006).

Given the scarcity of available molecular tracers originating from prestellar cores, previous studies have naturally focused on H₂D⁺ (Caselli et al. 2002, 2008; van der Tak et al. 2005; Vastel et al. 2006; Harju et al. 2006) and N₂H⁺ (Tafalla et al. 2004; Pagani et al. 2007; Lique et al. 2015), using high resolution ground-based submm observatories (14''–22''; JCMT, CSO). Those studies provided important constraints on the temperature structure of such cores and the column densities and abundances of molecular species at individual cores. Since then, the collisional data for both H₂D⁺ (Hugo et al. 2009) and N₂H⁺ (Lique et al. 2015) have been revised and therefore the reported properties of those clouds also need to be revised.

In the past decade, there have been several attempts to trace substructure via the dust continuum and molecular line emission within the central 1000 au of prestellar cores using interferometers, but mostly without positive results. There have been no fruitful detections of dust continuum at mm wavelengths toward single objects or toward a sample of cores in Perseus or Ophiuchus in the past (CARMA, IRAM, SMA; Schnee et al. 2012, 2010; Crapsi et al. 2007), explained by the shallow density profiles of prestellar cores at less than few thousand astronomical units. Due to the spatial capabilities and high sensitivity of the Atacama Large Millimeter and submillimeter Array (ALMA), a few studies have been able to detect some compact emission but without resolving the substructure (Friesen et al. 2014, 2018; Kirk et al. 2017). Only very recently have the inner regions of a single prestellar core been resolved for the first time (L1544 with ALMA; Caselli et al. 2019). Given the very challenging nature of the continuum observations at those high angular resolutions (50–100 au), it is no surprise that H₂D⁺ detection with ALMA is limited to a single low-mass prestellar core (SM1N; Friesen et al. 2014). Therefore, the analysis and modeling of single-dish observations of dust and gas toward prestellar cores is still a powerful approach to probe these enigmatic stages of star formation.

In this paper, we present the first dedicated study to test the proposed dynamical models using advanced non-LTE radiative transfer modeling to simulate the gas emission, making use of the optically thin H₂D⁺ 1₀⁻¹₁ and N₂H⁺ 4–3 line emission mapped with the James Clerk Maxwell Telescope (JCMT) toward a sample of prestellar cores (4) and one protostellar core. In Sect. 2, we describe the target selection, observations and data reduction. In Sect. 3, we present the observed spatial distributions and line analyses of H₂D⁺ and N₂H⁺ toward all cores and we also investigated the thermal and nonthermal contributions to the observed line width. In Sect. 4, we present the average column densities of H₂D⁺ and N₂H⁺ as determined using a non-LTE radiative transfer code to fit the line emission at the central submm peak position toward all cores. In Sect. 5, we used a detailed radiative transfer model to simulate the 2' × 2' maps of the H₂D⁺ and N₂H⁺ line emission toward the cores, adopting three different dynamical models: quasi-equilibrium Bonnor–Ebert sphere (QE-BES), singular isothermal sphere (SIS), and Larston–Penston (LP) flow, and a static sphere. Our main results and conclusions are presented in Sect. 7.

# 2. Observations

## 2.1. Targets

In this study, we present the analysis of five cores. Their coordinates, distances and continuum brightnesses at 850 µm and
<table><caption>**Table 1.** Source sample: positions, continuum brightness, and distance.</caption><thead><tr><th>Name</th><th>RA (J2000)<br>(hms)</th><th>Dec (J2000)<br>(dms)</th><th>S<sub>850</sub><br>(Jy)</th><th>d<br>(pc)</th><th>Mass<br>(M<sub>☉</sub>)</th></tr></thead><tbody><tr><td>L1521f</td><td>04:28:39.3</td><td>+26:51:33</td><td>0.639</td><td>140</td><td>6<sup>(a)</sup></td></tr><tr><td>L694-2</td><td>19:41:04.5</td><td>+10:57:02</td><td>0.397</td><td>250</td><td>7<sup>(b)</sup></td></tr><tr><td>L183</td><td>15:54:08.6</td><td>-02:52:45</td><td>0.480</td><td>110</td><td>10<sup>(c)</sup></td></tr><tr><td>L1544</td><td>05:04:17.2</td><td>+25:10:44</td><td>0.470</td><td>140</td><td>8<sup>(d)</sup></td></tr><tr><td>L1517B</td><td>04:55:18.3</td><td>+30:37:48</td><td>0.261</td><td>140</td><td>4<sup>(e)</sup></td></tr></tbody></table>

**Notes.** The coordinates for all cores are the positions of the respective peak 850 µm emission found in the maps available in the SCUBA Legacy Catalogue data (see Di Francesco et al. 2008).

**References.** (a) Crapsi et al. (2004). (b) Di Francesco et al. (2008). (c) Crutcher et al. (2004). (d) Tafalla et al. (1998). (e) Fu et al. (2011).

the core masses (4–10 $M_{\odot}$) based on dust and gas mm observations are presented in Table 1. These cores have been previously observed by Crapsi et al. (2005) in lines of $N_2H^+$ and $N_2D^+$ to explore the use of deuterium enrichment to constrain the evolutionary status of starless cores. The five cores of the present sample were specifically selected from the 31 cores presented in Crapsi et al. (2005) to provide a sample whereby detection of $H_2D^+$ could be considered favorable, typically by a high ratio of $N(N_2D^+/N_2H^+)$ (see Table 7; Crapsi et al. 2005). The present sample consists of four cores that are considered to be more evolved based on a variety of chemical and kinematical probes (L1544, L1521f, L694-2, L183) and one less evolved core (L1517B).

Of the five cores in Table 1, four contain no known young stellar objects and, therefore, can be considered starless. Only toward one core, L1521f, a very low luminosity object (VELLO) was identified by Bourke et al. (2006), and therefore it is considered to be protostellar (see also, Tokuda et al. 2016, 2017).

## 2.2. Submillimeter single-dish observations

The five cores in this study were observed from January 2007 to January 2008 as part of three semesters of queue-mode observing at the JCMT<fnref n="1" /> on Mauna Kea, Hawaii, under observing programs M06BC11, M07AC15, and M07BC06. All observations were made under very dry weather conditions (i.e., $\tau_{225} < 0.05$), to ensure maximum sensitivity to the ortho-$H_2D^+$ $1_{10}-1_{11}$ line and the $N_2H^+$ 4–3 line (at 372.42 and 372.67 GHz, respectively; Pickett et al. 1998), which are adjacent to a broad atmospheric $H_2O$ feature. JCMT has a beamwidth of 14" at the relevant frequency band (350 GHz).

The observations were made using the 16-element Heterodyne Array Receiver Programme (HARP) focal-plane array that operates over 325–375 GHz and the AutoCorrelation Spectrometer Imaging System (ACSIS) back-end (Buckle et al. 2009). HARP was tuned to 372.5 GHz, and ACSIS configured to provide nominally 500 MHz total bandwidth with 61 kHz wide channels or 0.048 km s⁻<fnref n="1" /> spacing. This setup allowed both the $H_2D^+$ $1_{10}-1_{11}$ and $N_2H^+$ 4–3 lines to be observed simultaneously in the same spectral window. The velocity resolution of ACSIS

data is a factor of ~1.2 less than the configured channel spacing, or in this case, 0.058 km s⁻<fnref n="1" />.

The observations were defined in minimum schedulable blocks (MSBs) and were preceded by standard focus calibration observations at 345 GHz on nearby bright objects like R Aql. Additionally, pointing calibrations at 345 GHz were conducted on objects such as CRL618, both before the program MSBs and every ~60–90 min. Flux calibration was monitored by observing when possible spectral line standards, like W75N, at a variety of frequencies. The aperture efficiency, $\eta_{\alpha}$, of HARP at 372.5 GHz is ~0.5.

HARP consists of 16 receptors arranged in a 4 × 4 square pattern with an on-sky spacing between receptors of 30". Program MSBs were executed by pointing one of the inner four HARP receptors at the target positions of the cores (Table 1). Data were obtained by pointing HARP at the target positions, and integrating using a “stare-mode” fashion, through chopping between the target positions and positions 180" distant in azimuth at 7.8125 Hz. Each MSB consisted of five repeats of 300 s. Each target was observed for ~4 h in total. When $H_2D^+$ $1_{10}-1_{11}$ was detected at the target position by at levels $\ge \sim 5\sigma$, the telescope was shifted diagonally in position by ~22.5", namely, the target position was centered in the array between the center four receptors, and further staring observations were to be executed for up to another ~4 h. This strategy allowed the acquisition of high sensitivity data of faint lines in a “checkerboard” pattern across each core at a spatial sampling better than given by a single HARP footprint. Observing each core with HARP in a jiggle pattern would have resulted in better spatial sampling but at the cost of sensitivity in each source. Since HARP is aligned on the JCMT in azimuth and elevation, the final data per pointing sometimes include samples of more than 16 positions on the sky due to field rotation. Additionally, the total integration duration varies from core to core, given the nature of queue-mode observing and the scarcity of very dry weather conditions.

## 2.3. Data reduction

The data were reduced using standard routines and procedures in the STARLINK reduction package<fnref n="2" /> (Currie et al. 2014). Each integration was visually checked for baseline ripples, absent detectors, or large spikes, and specific integrations where these occurred were removed from the data ensemble. In particular, receptor H03 was unavailable during most of the observations; at various times, no more than two other receptors were also unavailable. Less strong spikes were identified and removed using a standard methodology created by developers at the Joint Astronomy Centre. Spectral baselines were subtracted, frequency axes converted to velocities and spectra trimmed using STARLINK scripts. To increase the signal-to-noise ratio (S/N), we re-binned the velocity axes of all the data by applying a factor of 2, resulting in a spectral resolution of ~0.12 km s⁻<fnref n="1" />.

# 3. Observational results

## 3.1. Spatial distribution of $H_2D^+$ and $N_2H^+$

The observed spatial distributions of $H_2D^+$ $1_{10}-1_{11}$ emission toward all cores are presented in Fig. 1. $H_2D^+$ in emission is clearly seen toward each core of the sample. In all cases, the peak line brightness is found at the position of the 850 µm peak

<footnote n="1">¹ At the time of the described observations, the JCMT was operated by the Joint Astronomy Centre on behalf of the Science and Technology Facilities Council of the United Kingdom, the National Research Council of Canada and the Netherlands Organisation for Scientific Research.</footnote>

<footnote n="2">² The Starlink software is currently supported by the East Asian Observatory.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 2. Spatial distribution of N₂H⁺ 4–3 as observed with JCMT toward L1544 and L1521f. These two are the only cores with N₂H⁺ 4–3 line emission detections. The area filled with light green indicates the clear (>3σ) detections and the possible detections (~2σ) based on the line shape and peak velocity position.</figcaption></figure>

continuum emission (i.e., within a ~14'' beam width). The line in emission is also clearly detected at several positions offset from the central position toward all cores. The observations of N₂H⁺ 4–3 from the cores where it was clearly detected (L1521f and L1544) are presented in Fig. 2. The detection of this emission is spatially limited to two peaks toward L1521f and possibly L1544, and in particular, the second peak is at a location ~15'' offset from the local peak of the submillimeter continuum emission.

## 3.2. Line profiles

We fit single Gaussians to the H₂D⁺ 1₁₀–1₁₁ and N₂H⁺ 4–3³ lines observed for each core and at each position where the emission is detected (T*A ≥ 3 × rms). Table 2 lists the peak brightness (T*A), central line velocity (VLSR), and the line width (FWHM) with the associated positions expressed in offsets from the local peak of the submillimeter 850 µm continuum emission.

Figures 3–7 show the derived line properties, VLSR, FHWM, and T*A of H₂D⁺ and N₂H⁺ emission with respect to their offsets from the 850 µm dust peak. Firstly, we plot the VLSR versus offsets and compare these to the known source velocities taken from the literature (e.g., Caselli et al. 2008). The upper plots in Figs. 3–7 show that the measured VLSR at the central position is in very good agreement and mostly within the errors when compared to the values reported by Caselli et al. (2008) toward most of the cores. Caselli et al. (2008) observed H₂D⁺ with CSO and got measurements of the emission at the central position of those cores at a similar spectral resolution to the current study but a lower angular resolution (22'' compared to our 14''). In contrast, the source core velocity of L1544 reported by Ho et al. (1978) using NH₃ measurements is systematically lower than those measured in our study and in Caselli et al. (2008) for the central positions. Instead, the measured velocities reported

in Ho et al. (1978) come to a better agreement with those we find at offsets further out in the cloud. We should note that the angular resolution in the 1978 study was larger by a factor of >4. Taking all the above into consideration, we conclude that our measurements are consistent with previous work.

The FWHMs of the lines are of interest since they can provide information on whether the broadening of the emission is a result of only a thermal contribution or if nonthermal motions are also present at those very early stages of star formation. The line profiles of all five cores are characterized by narrow emission (0.2–0.6 km s⁻¹). The thermal broadening (ΔV<sub>THERM</sub>) is estimated using:

$$ \Delta V_{\text{THERM}} = \sqrt{8\ln2 \frac{kT}{m}}, \qquad (5) $$

where T is the gas kinetic temperature we determine for each individual core (6–9 K; see RADEX analysis in Sect. 4, Table 3), and m is the mass of the ion.

ΔV<sub>THERM</sub> was thus found to be within 0.26–0.32 km s⁻¹ for H₂D⁺ and 0.09–0.11 km s⁻¹ for N₂H⁺. By determining the ΔV<sub>THERM</sub> and measuring the observed line width of the species at each position in the maps, we can quantify how much of the observed line broadening is due to nonthermal contributions (ΔV<sub>NT</sub> = √[ΔV<sub>OBS</sub>² - ΔV<sub>THERM</sub>²]; Myers et al. 1991) for each source (i.e., turbulence and magnetic fields). To address the significance of nonthermal contributions, we approached the temperature errors such that the temperature is within the range 5–12 K (typical values for dense cores).

Here, we present our findings regarding the entire sample, while more details on each individual core can be found in Appendix A. The middle plots in Figs. 3–7 show that thermal broadening alone is insufficient to explain the observed line widths toward most of the offset positions toward all cores. In particular, there is an indication that the nonthermal contributions become less significant for both H₂D⁺ and N₂H⁺ as one moves outward from the central position for three out of five

<footnote n="3">³ There is no detectable hyperfine structure of N₂H⁺ 4–3, therefore it can be fit with a single Gaussian (see Splatalogue: https://www.cvs.nrao.edu/php/splat/advanced.php).</footnote>
<table><caption>**Table 2.** Line measurements at all offsets with clear detections for each core: JCMT data.</caption><thead><tr><th>Core</th><th>RA offset</th><th>Dec offset</th><th>Species</th><th>T<sup>*</sup><sub>A</sub></th><th>Error</th><th>V<sub>LSR</sub></th><th>Error</th><th>FWHM</th><th>Error</th><th>rms</th></tr><tr><th></th><th>('')</th><th>('')</th><th></th><th>(K)</th><th></th><th>(km s<sup>-1</sup>)</th><th></th><th>(km s<sup>-1</sup>)</th><th></th><th>(K)</th></tr></thead><tbody><tr><td rowspan="15">L1544</td><td>0</td><td>0</td><td rowspan="5">H<sub>2</sub>D<sup>+</sup></td><td>0.55</td><td>0.02</td><td>7.14</td><td>0.01</td><td>0.57</td><td>0.03</td><td>0.05</td></tr><tr><td>0</td><td>30</td><td>0.26</td><td>0.03</td><td>7.07</td><td>0.02</td><td>0.48</td><td>0.06</td><td>0.06</td></tr><tr><td>-30</td><td>30</td><td>0.17</td><td>0.02</td><td>7.10</td><td>0.03</td><td>0.61</td><td>0.08</td><td>0.05</td></tr><tr><td>15</td><td>15</td><td>0.20</td><td>0.03</td><td>7.12</td><td>0.04</td><td>0.45</td><td>0.08</td><td>0.06</td></tr><tr><td>-15</td><td>15</td><td>0.37</td><td>0.03</td><td>7.11</td><td>0.02</td><td>0.49</td><td>0.05</td><td>0.06</td></tr><tr><td>-45</td><td>15</td><td>0.16</td><td>0.04</td><td>7.21</td><td>0.04</td><td>0.4</td><td>0.1</td><td>0.07</td></tr><tr><td>-30</td><td>0</td><td>0.18</td><td>0.02</td><td>7.17</td><td>0.02</td><td>0.46</td><td>0.05</td><td>0.04</td></tr><tr><td>15</td><td>-15</td><td>0.27</td><td>0.02</td><td>7.21</td><td>0.02</td><td>0.53</td><td>0.05</td><td>0.05</td></tr><tr><td>-15</td><td>-15</td><td>0.28</td><td>0.05</td><td>7.21</td><td>0.03</td><td>0.38</td><td>0.07</td><td>0.05</td></tr><tr><td>0</td><td>-30</td><td>0.10</td><td>0.02</td><td>7.20</td><td>0.06</td><td>0.5</td><td>0.1</td><td>0.05</td></tr><tr><td>0</td><td>0</td><td>N<sub>2</sub>H<sup>+</sup></td><td>0.11</td><td>0.02</td><td>7.20</td><td>0.04</td><td>0.4</td><td>0.1</td><td>0.03</td></tr><tr><td>-15</td><td>15</td><td></td><td>0.12</td><td>0.03</td><td>7.26</td><td>0.05</td><td>0.4</td><td>0.1</td><td>0.07</td></tr><tr><td rowspan="15">L183</td><td>15</td><td>45</td><td rowspan="5">H<sub>2</sub>D<sup>+</sup></td><td>0.31</td><td>0.04</td><td>2.32</td><td>0.03</td><td>0.56</td><td>0.08</td><td>0.09</td></tr><tr><td>-15</td><td>45</td><td>0.30</td><td>0.04</td><td>2.48</td><td>0.04</td><td>0.58</td><td>0.09</td><td>0.09</td></tr><tr><td>30</td><td>30</td><td>0.17</td><td>0.04</td><td>2.23</td><td>0.03</td><td>0.32</td><td>0.08</td><td>0.06</td></tr><tr><td>0</td><td>30</td><td>0.45</td><td>0.03</td><td>2.36</td><td>0.01</td><td>0.3</td><td>0.2</td><td>0.06</td></tr><tr><td>-30</td><td>30</td><td>0.21</td><td>0.03</td><td>2.41</td><td>0.04</td><td>0.1</td><td>0.1</td><td>0.07</td></tr><tr><td>15</td><td>15</td><td>0.37</td><td>0.04</td><td>2.38</td><td>0.03</td><td>0.58</td><td>0.07</td><td>0.09</td></tr><tr><td>-15</td><td>15</td><td>0.43</td><td>0.03</td><td>2.39</td><td>0.02</td><td>0.57</td><td>0.05</td><td>0.07</td></tr><tr><td>0</td><td>0</td><td>0.48</td><td>0.03</td><td>2.40</td><td>0.01</td><td>0.45</td><td>0.03</td><td>0.06</td></tr><tr><td>15</td><td>-15</td><td>0.34</td><td>0.05</td><td>2.38</td><td>0.03</td><td>0.50</td><td>0.08</td><td>0.1</td></tr><tr><td>-15</td><td>-15</td><td>0.43</td><td>0.04</td><td>2.38</td><td>0.02</td><td>0.42</td><td>0.05</td><td>0.05</td></tr><tr><td>0</td><td>-30</td><td>0.29</td><td>0.04</td><td>2.36</td><td>0.03</td><td>0.42</td><td>0.07</td><td>0.09</td></tr><tr><td>15</td><td>-45</td><td>0.19</td><td>0.05</td><td>2.25</td><td>0.05</td><td>0.4</td><td>0.1</td><td>0.09</td></tr><tr><td>0</td><td>-60</td><td>0.12</td><td>0.03</td><td>2.32</td><td>0.05</td><td>0.4</td><td>0.1</td><td>0.06</td></tr><tr><td>-15</td><td>-45</td><td>0.19</td><td>0.06</td><td>2.40</td><td>0.04</td><td>0.3</td><td>0.1</td><td>0.07</td></tr><tr><td rowspan="15">L694-2</td><td>0</td><td>30</td><td rowspan="5">H<sub>2</sub>D<sup>+</sup></td><td>0.18</td><td>0.02</td><td>9.65</td><td>0.03</td><td>0.45</td><td>0.07</td><td>0.05</td></tr><tr><td>15</td><td>15</td><td>0.18</td><td>0.02</td><td>9.60</td><td>0.03</td><td>0.55</td><td>0.07</td><td>0.05</td></tr><tr><td>-15</td><td>15</td><td>0.35</td><td>0.02</td><td>9.64</td><td>0.02</td><td>0.48</td><td>0.04</td><td>0.05</td></tr><tr><td>0</td><td>0</td><td>0.41</td><td>0.02</td><td>9.56</td><td>0.02</td><td>0.51</td><td>0.04</td><td>0.06</td></tr><tr><td>-30</td><td>0</td><td>0.21</td><td>0.03</td><td>9.63</td><td>0.03</td><td>0.50</td><td>0.08</td><td>0.06</td></tr><tr><td>15</td><td>-15</td><td>0.27</td><td>0.02</td><td>9.58</td><td>0.02</td><td>0.48</td><td>0.04</td><td>0.04</td></tr><tr><td>-15</td><td>-15</td><td>0.20</td><td>0.02</td><td>9.52</td><td>0.02</td><td>0.41</td><td>0.04</td><td>0.03</td></tr><tr><td>30</td><td>-30</td><td>0.12</td><td>0.03</td><td>9.57</td><td>0.05</td><td>0.4</td><td>0.1</td><td>0.06</td></tr><tr><td>0</td><td>-30</td><td>0.15</td><td>0.02</td><td>9.47</td><td>0.04</td><td>0.63</td><td>0.09</td><td>0.05</td></tr><tr><td rowspan="5">L1517B</td><td>-15</td><td>15</td><td rowspan="5">H<sub>2</sub>D<sup>+</sup></td><td>0.10</td><td>0.02</td><td>5.79</td><td>0.03</td><td>0.38</td><td>0.07</td><td>0.03</td></tr><tr><td>0</td><td>0</td><td>0.17</td><td>0.02</td><td>5.82</td><td>0.03</td><td>0.07</td><td>0.08</td><td>0.05</td></tr><tr><td>15</td><td>-15</td><td>0.12</td><td>0.01</td><td>5.73</td><td>0.03</td><td>0.49</td><td>0.06</td><td>0.03</td></tr><tr><td>-15</td><td>-15</td><td>0.10</td><td>0.01</td><td>5.77</td><td>0.02</td><td>0.36</td><td>0.06</td><td>0.03</td></tr><tr><td>0</td><td>-30</td><td>0.10</td><td>0.03</td><td>5.88</td><td>0.04</td><td>0.4</td><td>0.1</td><td>0.04</td></tr><tr><td rowspan="15">L1521f</td><td>-15</td><td>45</td><td rowspan="5">H<sub>2</sub>D<sup>+</sup></td><td>0.17</td><td>0.02</td><td>6.41</td><td>0.02</td><td>0.34</td><td>0.05</td><td>0.07</td></tr><tr><td>0</td><td>30</td><td>0.30</td><td>0.05</td><td>6.45</td><td>0.03</td><td>0.35</td><td>0.07</td><td>0.09</td></tr><tr><td>0</td><td>-30</td><td>0.21</td><td>0.06</td><td>6.43</td><td>0.03</td><td>0.24</td><td>0.07</td><td>0.06</td></tr><tr><td>15</td><td>15</td><td>0.18</td><td>0.02</td><td>6.53</td><td>0.02</td><td>0.53</td><td>0.05</td><td>0.04</td></tr><tr><td>15</td><td>45</td><td>0.10</td><td>0.02</td><td>6.35</td><td>0.06</td><td>0.6</td><td>0.1</td><td>0.03</td></tr><tr><td>-15</td><td>15</td><td>0.16</td><td>0.02</td><td>6.54</td><td>0.03</td><td>0.55</td><td>0.07</td><td>0.03</td></tr><tr><td>0</td><td>0</td><td>0.27</td><td>0.05</td><td>6.54</td><td>0.03</td><td>0.35</td><td>0.07</td><td>0.09</td></tr><tr><td>15</td><td>-15</td><td>0.18</td><td>0.02</td><td>6.51</td><td>0.03</td><td>0.46</td><td>0.07</td><td>0.05</td></tr><tr><td>-15</td><td>-15</td><td>0.10</td><td>0.02</td><td>6.50</td><td>0.05</td><td>0.6</td><td>0.1</td><td>0.04</td></tr><tr><td>-45</td><td>15</td><td rowspan="5">N<sub>2</sub>H<sup>+</sup></td><td>0.07</td><td>0.02</td><td>6.62</td><td>0.08</td><td>0.6</td><td>0.2</td><td>0.03</td></tr><tr><td>-15</td><td>15</td><td>0.11</td><td>0.02</td><td>6.47</td><td>0.03</td><td>0.36</td><td>0.08</td><td>0.04</td></tr><tr><td>0</td><td>0</td><td>0.30</td><td>0.04</td><td>6.56</td><td>0.04</td><td>0.6</td><td>0.1</td><td>0.06</td></tr><tr><td>0</td><td>30</td><td>0.22</td><td>0.06</td><td>6.53</td><td>0.04</td><td>0.29</td><td>0.09</td><td>0.07</td></tr></tbody></table>

**Notes.** To convert to main beam temperature, $T_{mb}$, one should divide the antenna temperature, $T_A^*$, with the beam efficiency $\eta = 0.6 (T_{mb} = \frac{T_A^*}{\eta})$.
<figure><img src="image_3.png" /></figure>

**Fig. 3.** $V_{LSR}$, FWHM, and $T_A^*$ measurements of H₂D⁺ plotted over offset from the 850 µm dust peak of L183. The shaded region in the middle plot corresponds to the thermal broadening within the 5–12 K temperature range. The shaded regions in purple correspond to the errors reported in Caselli et al. (2008). The square, triangle, and circle refer to north, south, and 0'' declination offsets, respectively. The blue and red symbols refer to offsets toward the east and west, respectively, while the black symbols correspond to 0'' Right Ascension offsets. For example, a 40'' offset can be represented by a red square (40'' north-west), a black triangle (40'' south), a blue circle (40'' east), and so on.

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

**Fig. 4.** Same as Fig. 3, but for L694-2.

cores. This finding contradicts the theory and what was previously observed toward molecular clouds, where the velocity dispersion increases with scale as the density decreases (e.g., Maloney 1988).

For the two cores for which both H₂D⁺ and N₂H⁺ measurements are available, we compared the nonthermal contribution at the central position. In particular, we find nonthermal contributions for H₂D⁺ to be 0.47<sup>+0.05</sup><sub>-0.06</sub> km s⁻¹ and 0.23<sup>+0.15</sup><sub>-0.10</sub> km s⁻¹ for L1544 and L1521f, respectively, while for N₂H⁺ we find these

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

**Fig. 5.** Same as Fig. 3, but for L1517B.

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

**Fig. 6.** Same as Fig. 3, but for L1544.

to be 0.42<sup>+0.10</sup><sub>-0.08</sub> km s⁻¹ and 0.6 ± 0.1 km s⁻¹, respectively. For L1544, the regions traced by the two species are characterized by similar nonthermal contributions (within the uncertainties) and, therefore, originate from regions with similar internal motions. For the more evolved, L1521f core, the nonthermal contributions are higher by a factor of 3 for the N₂H⁺ line compared to the H₂D⁺, which is beyond the uncertainties. This latter difference indicates that even though both lines are observed at a central position of the protostellar core, the H₂D⁺ emission originates from more quiescent gas than the N₂H⁺ emission does. This difference can be explained by the more complex structure of protostellar cores (i.e., outflow, disk) compared with prestellar cores.

Lastly, the $T_A^*$ at each offset can provide valuable information regarding the gas density and temperature properties toward the cores. The lower plots in Figs. 3–7 suggest that $T_A^*$ generally decreases with increasing offset. This behavior can be explained
<figure><img src="image_7.png" /><figcaption>Fig. 7. Same as Fig. 3, but for L1521f.</figcaption></figure>

by the decreasing density or decreasing excitation temperature, $T_{ex}$ (i.e., sub-thermal excitation), as one moves outward from the center of the core, and it is suggestive of the low optical depth of the line (see also Sect. 5.3). We note that the kinetic temperature in prestellar cores generally increases outward (up to 5 K difference) and therefore it cannot be behind the observed decrease in $T_A^*$. The observed small change in temperature does not seem to be the one driving the brightness of the emitting line. Fluctuations in the $T_A^*$ values of the lines either at the same offsets (see, e.g., L1544) or at increasing offsets (see, e.g., L183) indicate the presence of a nonhomogeneous medium (e.g., clumps, cavities).

# 4. Column densities and gas temperatures

## 4.1. Radiative transfer modeling

The determination of the temperature and column densities of the targeted species is a first step to help us study the chemistry, and, therefore, evolutionary status of the cores. We estimated the column density of $H_2D^+$ and $N_2H^+$ and kinetic temperature toward all cores by modeling the velocity-integrated main beam temperature at the peak/central position using the non-LTE radiative transfer code RADEX (van der Tak et al. 2007). We assumed an isothermal homogeneous sphere. RADEX predicts the line intensities of the molecular transitions of interest for a given set of physical parameters: kinetic temperature, $H_2$ density, molecular column density, Rayleigh-Jeans equivalent background temperature of a black body shaped radiation field (in this case, the cosmic microwave background (CMB) temperature of 2.73 K was adopted), and the observed line width. To convert to the main beam temperature, we adopted the main beam efficiency, $\eta_{mb}$, of 0.6 ($\eta_\alpha/\eta_{mb} \sim 0.8$) measured by Buckle et al. (2009).

### 4.1.1. $H_2D^+$

The observed spatial distribution of the $H_2D^+$ emission is extended for all cores, and therefore the implicit assumption of unit beam filling factor is proper. The most recent available collisional properties of the $H_3^+-H_2$ system and all its isotopic

<table><caption>**Table 3.** Column density ($N_{mol}$) and gas kinetic temperature ($T_{KIN}$) estimates for the adopted volume density ($n_{H_2}$) at the central position of core using RADEX.</caption><thead><tr><th>Species</th><th>Core</th><th><i>T</i><sub>KIN</sub><sup>(a)</sup><br>(K)</th><th><i>n</i><sub>H<sub>2</sub></sub> (10<sup>6</sup>)<br>(cm<sup>-3</sup>)</th><th>N(mol)<sup>(a)</sup><br>(10<sup>13</sup> cm<sup>-2</sup>)</th><th>Δ<i>V</i><sub>th</sub><sup>(b)</sup><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="5">H<sub>2</sub>D<sup>+</sup></td><td>L1544</td><td>9/7</td><td>2.0</td><td>2/3.2</td><td>0.32/0.28</td></tr><tr><td>L183</td><td>8/7</td><td>2.0</td><td>1.4/2.5</td><td>0.30/0.28</td></tr><tr><td>L694-2</td><td>9/7</td><td>0.9</td><td>0.8/3.2</td><td>0.32/0.28</td></tr><tr><td>L1517B</td><td>6/9.5</td><td>0.2</td><td>0.4/0.4</td><td>0.26/0.33</td></tr><tr><td>L1521f</td><td>6/9.3</td><td>1.1</td><td>0.5/0.7</td><td>0.26/0.33</td></tr><tr><td rowspan="2">N<sub>2</sub>H<sup>+</sup></td><td>L1544</td><td>5</td><td>2.0</td><td>0.8</td><td>0.09</td></tr><tr><td>L1521f</td><td>8</td><td>1.1</td><td>0.9</td><td>0.11</td></tr></tbody></table>

**Notes.** The volume density for each core is adopted by Caselli et al. (2008). <sup>(a)</sup>The first values of $T_{KIN}$ and $N_{mol}$ refer to this work followed by the values reported in Caselli et al. (2008) for direct comparison. <sup>(b)</sup>The thermal broadening is computed based on $T_{KIN}$. Accounting for kinetic temperatures within a range of 5–12 K, the reported column densities can vary by up to ~60–80%.

variations were studied and presented in Hugo et al. (2009). To determine the excitation of $H_2D^+$, we adopted the rates presented in Hugo et al. (2009) with two simplifications. First, our calculations only consider the lowest two energy levels of $o-H_2D^+$, the $l_{11}$ ground state, and the $l_{10}$ excited state which lies 17 K above the ground. This two-level approximation is valid at the low temperatures of prestellar cores ($T \sim 10$ K), and therefore for the present sample, since the next highest level, $l_{21}$, lies at 109 K above the ground and its excitation is therefore negligible. The second simplification is that our calculations ignore reactive collisions. We note that in general the collisions between $H_3^+$ and $H_2$ may lead to reaction (i.e., the fully elastic case), excitation (i.e., the fully inelastic case), or both. Our assumption is justified by the work of Hugo et al. (2009), who show that at the low temperatures of interest here, the reactive collision rates are factors >100 lower than the inelastic ones.

Our calculations include inelastic collisions of $o-H_2D^+$ with $o-H_2$ and $p-H_2$. Since we are dealing with very low temperatures, the ortho-to-para ratio of $H_2$ is assumed to be equal to its thermal value (~$10^{-4}$ in our case) instead of the high-temperature limit of 3. In several previous studies (van der Tak et al. 2005; Vastel et al. 2006; Caselli et al. 2008), the $H_2D^+$ abundances toward the cores of our sample (when studied) were calculated using scaled radiative rates adopted by Black et al. (1990), and therefore the column densities and abundances derived in those studies should be updated with the most recent rates, as we did here.

### 4.1.2. $N_2H^+$

In this work, we used the newest available collisional data of $N_2H^+$ (Lique et al. 2015). As opposed to former collisional data, where the collisional partner was Helium (He), the new ones are calculated based on collisions with the most abundant partner in the interstellar medium, which is $H_2$, and are found to be significantly different (Lique et al. 2008).

## 4.2. RADEX results

Table 3 presents the adopted gas densities $n(H_2)$ of the cores along with their derived column densities ($N(H_2D^+)$ and
N(N₂H⁺)), and derived temperatures (T_KIN) at the central posi-
tion of each core. The average volume density, at least at the
core center where we are focusing in this part of the analysis, was
adopted as the values presented in Caselli et al. (2008). In partic-
ular for this sample, we report low kinetic temperatures between
6 and 9 K, and column densities for both species lying within the
range of 0.4–2 × 10¹³ cm⁻². The derived kinetic temperatures are
used to assess the line thermal broadening per core presented in
Sect. 3.2.

Comparing our observables with those presented in Caselli
et al. (2008), we see that the reported T<sub>mb</sub> values are similar in
both studies, which indicates that the difference in beam size
between JCMT and CSO plays a minor role, supporting our con-
clusion that the emission is extended and the beam is filled. The
line widths are also similar except for L1521f and L1544. The
observed difference is likely due to the more limited spectral
resolution of the CSO backends. The column densities among
the studies are consistent (mostly within a factor of 2), given
the differences in adopted kinetic temperatures (up to 3 K),
observed line fluxes, and molecular data. Increasing the num-
ber of observed transitions per species, and at multiple offsets,
would allow the accurate determination of the kinetic gas tem-
perature, column, and volume density, and the production of the
column and gas kinetic temperature maps of the cores.

5. Modeling the spatial distribution of H₂D⁺ and N₂H⁺

A primary goal of this paper is to distinguish between dynamical
models that turn a prestellar core into a protostar that have been
proposed in the literature. To do so, we modeled the spatial dis-
tribution and line profiles of H₂D⁺ and N₂H⁺ (detected toward
only two cores) for each core, as predicted from each dynamical
model using advanced radiative transfer modeling, and compare
the predicted emission with our observed maps.

**5.1. RATRAN: initializing the model**

To model the spatial distribution and line profiles of H₂D⁺ and N₂H⁺ (where present) emission from the cores, we ran the Monte Carlo radiative transfer code RATRAN (Hogerheijde & van der Tak 2000). RATRAN calculates synthetic line emission for the species of interest, taking into account the temperature and density gradients of the source, the dust emission and absorption properties, and the small-scale (e.g., thermal motions, turbulence) and large-scale (e.g., infall motions) source kinematics. We fix the physical structure of the sources by adopting the corresponding predicted temperature, density, and velocity structure from proposed dynamical models in the literature. The dynamical models examined in this study are described in Sect. 5.2. For our models, we assumed that the dust temperature is equal to the gas temperature, which is valid for environments of high volume density (>10⁴ cm⁻³), and a gas-to-dust ratio of 100.

For our calculations, and to simulate the 2' × 2' observa-
tions, we defined a grid of 13 logarithmically spaced shells for
each core, extending to an outer radius of 10 000 au for all
cores except L694-2, which was extended to 15 000 au due to
its larger distance (see Table 1). The density and velocity struc-
tures are shown in Fig. 8 and were adopted by Keto et al. (2015)
for the different models which are described in more detail in
Sects. 5.2.1–5.2.3.

The temperature profiles were adopted from the literature
(L183, L694-2, L1517B, L1544, L1521f Pagani et al. 2007;

<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Density (solid line) and velocity (dashed line) profiles of a prestellar core predicted by the three dynamical models: LP flow (green), SIS (red), and QE-BES (blue), as presented in Keto et al. (2015). Each of the modeled profiles was adopted to model the H₂D⁺ emission using RATRAN.</figcaption></figure>

Redaelli et al. 2018; Maret et al. 2013; Chacón-Tanarro et al. 2019; Crapsi et al. 2004). These profiles were derived using radiative transfer models to fit the observed gas (e.g., N₂H⁺) or dust emission (e.g., mm IRAM observations), and are specific for each core. The LP and SIS models are by definition isothermal, while the QE-BES model predicts a range of temperatures, and is more consistent with the observed temperature variations (~5 K). Our approach to adopting the same temperature profile in all dynamical models may affect the brightness of the line, but not the dynamics (Keto & Caselli 2010). In this study, we examined the line morphology, which is highly dependent on the dynamics, meaning, the distribution of the velocities of the gas, rather than the small variations in the gas kinetic temperature; therefore, our general conclusions are robust.

The Doppler parameter, *b*, was calculated for each shell
as equal to 0.6 × *FWHM*, where *FWHM* is the result of the
Gaussian fitting process for each core and position presented in
Table 2. In Table 4, we present the initial abundance profiles
adopted during the modeling process for H₂D⁺ and N₂H⁺ for
each core with the associated references, and the relevant modifi-
cations per model to achieve a good fit. In particular, we proceed
to a more accurate determination of the abundance profiles by fit-
ting the emission at the central position and expand the solution
to the entire map (see Sect. 5.3). Lastly, to compare the mod-
eled emission with observations, the resulting modeled spatial
distribution is convolved to the beam size of the JCMT (14" at
372 GHz).

**5.2. Physical structure: dynamical models**

Multiple hypotheses exist about the dynamical processes that
take place for a prestellar core to collapse. To date, the determi-
nation of the most accurate model has been challenging. In this
study, we examined three different dynamical models based on
spherical symmetry that result in similar density profiles but very
different velocity profiles across the cores. These models are the
QE contraction of a BES, the LP flow (LP), and the inside-out
collapse of the SIS.
<table><caption>**Table 4.** H₂D⁺ and N₂H⁺ abundances estimated with RATRAN for each core based on the JCMT observations.</caption><thead><tr><th rowspan="2">Core</th><th rowspan="2">Species</th><th colspan="3">Literature</th><th rowspan="2">Authors</th><th rowspan="2">QE-BES</th><th colspan="3">This work</th></tr><tr><th>X<sub>IN</sub></th><th>X<sub>INTERMEDIATE</sub></th><th>X<sub>OUT</sub></th><th>LP</th><th>SIS</th><th>Static</th></tr></thead><tbody><tr><td rowspan="2">L183</td><td>H₂D⁺</td><td>(r &lt; 300 au)<br>1 × 10<sup>-10</sup></td><td>(300 au &lt; r &lt; 5000 au)<br>2 × 10<sup>-10</sup></td><td>(5000 au &lt; r &lt; 10 000 au)<br>4 × 10<sup>-11</sup></td><td>α</td><td>2</td><td>Multiplication factor<br>0.5</td><td>100</td><td>1</td></tr><tr><td>L1544</td><td>(r &lt; 2500 au)<br>1 × 10<sup>-9</sup><br>(r &lt; 2500 au)<br>0</td><td>n/a<br>n/a<br>n/a<br>n/a</td><td>(2500 au &lt; r &lt; 10 000 au)<br>5 × 10<sup>-10</sup><br>(2500 au &lt; r &lt; 10 000 au)<br>2 × 10<sup>-10</sup></td><td>β<br>γ</td><td>1<br>20</td><td>Multiplication factor<br>0.3<br>Multiplication factor<br>0.8</td><td>60<br>500</td><td>0.8<br>4</td></tr><tr><td>L694-2</td><td>H₂D⁺</td><td>n/a<br>n/a</td><td>n/a<br>n/a</td><td>n/a<br>n/a</td><td>n/a<br>n/a</td><td>4.4 × 10<sup>-10</sup></td><td>Abundance (const.)<br>5.2 × 10<sup>-11</sup></td><td>2 × 10<sup>-8</sup></td><td>2.5 × 10<sup>-10</sup></td></tr><tr><td>L1517B</td><td>H₂D⁺</td><td>n/a<br>n/a</td><td>n/a<br>n/a</td><td>n/a<br>n/a</td><td>n/a<br>n/a</td><td>5 × 10<sup>-11</sup></td><td>Abundance (const.)<br>1.2 × 10<sup>-11</sup></td><td>5 × 10<sup>-9</sup></td><td>1.5 × 10<sup>-10</sup></td></tr><tr><td>L1521f</td><td>H₂D⁺</td><td>(r &lt; 3500 au)<br>3 × 10<sup>-10</sup><br>(r &lt; 3500 au)<br>1 × 10<sup>-10</sup></td><td>n/a<br>n/a<br>n/a<br>n/a</td><td>(3500 au &lt; r &lt; 10 000 au)<br>6 × 10<sup>-11</sup><br>(3500 au &lt; r &lt; 10 000 au)<br>2 × 10<sup>-10</sup></td><td>δ</td><td>1.3<br>1</td><td>Multiplication factor<br>0.75<br>Multiplication factor<br>0.3</td><td>100<br>15</td><td>1<br>1</td></tr></tbody></table>

**Notes.** The table contains also the initial abundance profiles adopted from literature when available for direct comparison. The first two columns of the table report the core name and the observed species. The abundance profiles of the observed species per core as retrieved by the literature are presented in Cols. 3–5. In particular, we report the inner, intermediate and outer abundances (X<sub>IN</sub>, X<sub>INTERMEDIATE</sub> and X<sub>OUT</sub>) following their corresponding shell sizes. Lastly, the resulted abundance profiles of this work and per model (QE-BES, LP, SIS, and Static) are reported in Cols. 7–10. The multiplication factor is the one applied to the abundance profiles taken from literature (when available) resulting in the best fit between the modeled and observed line intensity at the central position of the maps. For L1517B and L694-2 we report the resulted constant abundance. Authors: α Pagani et al. (2009); for 300 au < r < 5000 au, the average value of the reported abundances was adopted, β van der Tak et al. (2005), γ Caselli et al. (2002), δ Crapsi et al. (2004).

We note that despite the similarity of the slopes in the density distribution, the different density profiles account for a total cloud mass of ~10 M<sub>☉</sub> for the QE-BES model, ~27 M<sub>☉</sub> for the LP model, and ~4 M<sub>☉</sub> for the SIS model, for a radius up to 30 000 au. Based on mass estimations alone, the SIS and QE-BES models appear to be the most suitable in reproducing the estimated mass for the cores in our sample (based on dust continuum observations; see Table 1). Here, we investigate the behavior of the gas. The quite distinct velocity profiles resulting from the different dynamical models are shown in Fig. 8. Gas emission observations in the optically thin regime that can trace those kinematics are therefore crucial. The models of interest are described below.

### 5.2.1. LP flow

Larson (1969) attempted to solve numerically the dynamical equations concerning the collapse of a spherically symmetric isothermal core, known as Larson core, which is gravitationally unstable, yet all rotational, magnetic, and turbulent motions are neglected (free-fall collapse). The collapse, in this case, is non-homologous, with the central region collapsing to a high density and reaching hydrostatic equilibrium and the outer regions collapsing to a much lesser extent, creating a first hydrostatic core object. In the meantime, Penston (1969) presented the dynamics of self-gravitating gaseous spheres, focusing on the isothermal case, which collapses in free-fall. The analytical solutions to the equations of free-fall collapse are used to produce the various symmetries (including spherical) at the instant when a dense hydrostatic object forms in the innermost dense regions of the core. For the hydrodynamic simulations, the diffusion approximation for radiative transfer was used.

The density and velocity profiles of the LP model (Fig. 8) were adopted by Keto et al. (2015), and are products of

integrating Eqs. (11) and (12) in Shu (1977). To produce a similar density profile (slope ~r<sup>-2</sup>) for all models (LP, SIS, QE-BES), the sound speed was 0.2 km s<sup>-1</sup> at ~11 K, and the evolutionary time was 2.8 × 10<sup>11</sup> s. The density follows a flat distribution with a central density n = 6.3 × 10<sup>6</sup> cm<sup>-3</sup> out to ~850 au and follows an r<sup>-2</sup> slope outward, with n = 3.5 × 10<sup>4</sup> cm<sup>-3</sup> at ~10 000 au. The velocity gradually increases from 5 × 10<sup>-2</sup> km s<sup>-1</sup> at 100 au, reaching a maximum of 10<sup>-1</sup> km s<sup>-1</sup> at ~850 au, where it reaches and remains at a constant velocity through to the outermost regions of the core.

### 5.2.2. SIS

Shu (1977) proposed the SIS model to describe the gravitational collapse of isothermal gas clouds. In this case, the collapse occurs from the inside out, with the inner and outer envelopes modeled separately. The inner envelope is assumed to collapse under free-fall conditions while the outer envelope is modeled as essentially static. The proposed density distribution in this case follows the r<sup>-3/2</sup> for the inner envelope, and r<sup>-2</sup> for the outer envelope, adopting a sound speed, α = 0.2 km s<sup>-1</sup> at ~11 K, and an evolutionary time t = ~2.5 × 10<sup>5</sup> yr. In particular, the inner density is 2.5 × 10<sup>6</sup> cm<sup>-3</sup> at 100 au, decreasing down to 1.4 × 10<sup>4</sup> cm<sup>-3</sup> at ~10 000 au. The velocity profile in this dynamical model is characterized by a higher velocity at the center of the core (compared to the other dynamical models) of ~4 km s<sup>-1</sup> at ~100 au and decreases down to 0.5 km s<sup>-1</sup> at ~10 000 au, before starting to behave as a static sphere without infall in the outer parts of the core.

### 5.2.3. QE-BES

This model treats the prestellar core as a BES undergoing quasi-static contraction, in which the pressure in the system remains
<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Example of the abundance profile adjustment during the fitting process. The observed line profiles of H₂D⁺ are overplotted with the modeled ones using RATRAN toward the central position of the L183 core. Three different dynamical models were adopted, QE-BES, SIS and LP flow. The original abundance profiles (left panel) were adopted by Pagani et al. (2009) and cannot reproduce the central peak of the observations. Adjusting the H₂D⁺ abundance profile by applying a multiplication factor, resulted in a better reproduction of the observations for the LP and QE-BES models, while the SIS model failed to match the shape of the observed line profile.</figcaption></figure>

uniform and constant, due to a low rate of volume change. The pressure is not affected by small oscillations connected with small-scale processes within or near the cores. The notion of a BES was originally put forward in 1956 when it was shown that an isothermal sphere of gas in an external pressurized medium has a critical size above which it is gravitationally unstable (Bonnor 1956; Ebert 1957). This model predicts a density profile following $r^{-2}$ in the outer region of the core with an approximately constant density inward (e.g., Keto & Caselli 2010). The central density in this case is $n = 10^7$ cm⁻³ (Keto et al. 2015). The velocity profile predicted by this contraction model shows a so-called Λ-shape, starting from 0 km s⁻¹ in the center (<100 au) and increasing outward, while reaching a maximum value (~1 km s⁻¹) near the characteristic radius, $r_f$, (~1000 au in this case), where we observe the transition from a constant/flat density distribution to that of $r^{-2}$. After the velocity reaches its maximum value, it decreases again when moving toward the edge of the core.

## 5.3. Abundances

Before determining the most representative dynamical model per source, we proceed with fine-tuning the abundance profiles of H₂D⁺ and N₂H⁺.

### 5.3.1. Abundance derivation: dynamical approach

We started the fitting process by adopting an abundance profile for the observed H₂D⁺ per core from the literature and proceed to a more accurate determination of the abundance profiles by fitting the line emission at the central position for each core. To do so, we let the abundance vary by up to 2 orders of magnitude higher and lower from the initial abundance and using a step of a factor of 2. Table 4 lists the initial H₂D⁺ abundance profiles adopted from the literature with their associated references and the multiplication factor we applied in this study to fit the peak line intensity of the central position. A constant abundance profile was also tested toward these cores but without resulting in a qualitative change on the results. Therefore, we adopted the

more detailed abundance profiles when available throughout the paper. We note that for two cores, L694-2 and L1517B, initial abundance profiles could not be retrieved from the literature; therefore, we assumed a constant initial abundance of $1 \times 10^{-10}$ which we varied as described above. Figure 9 shows an example of the fitting process toward L183.

Our observations concern only single transitions of both H₂D⁺ and N₂H⁺, while the physical structure is fixed. Therefore, although our fitting process results in an uncertainty of ~14%-18% (~3 × $\chi^2_{min}$) for the derived abundances, this process somewhat underestimates the real uncertainties. Once we find the abundance profile that can reproduce the observed line profile in the central position of the core, we expand the solution to the entire map. In Table 5, we present the resulting optical depth for each line per model toward the central positions of each core. We find that the emission from both H₂D⁺ l₁₀–l₁₁ and N₂H⁺ 4–3 is mostly in the optically thin regime (τ < 1) justifying their good performance in tracing the dynamics of the cores.

### 5.3.2. Static sphere

To provide an independent determination of the abundances, we proceed with applying the same radiative transfer technique, but this time we adopted the physical structure of each source based on their previous dust and gas observations. In this case we assumed a static envelope, without infall or expansion motions. The temperature and density profiles were adopted from the literature (L183, L694-2, L1517B, L1544, L1521f Pagani et al. 2007; Redaelli et al. 2018; Maret et al. 2013; Chacón-Tanarro et al. 2019; Crapsi et al. 2004), and are based on continuum and gas observations, and they are, therefore, unique to each core. Again, the basic assumptions for all cores are that the gas and dust are thermally coupled, their temperatures are therefore set to be equal, and the gas to dust ratio is 100. As before, we defined a grid of 13 logarithmically spaced shells for each core, while the density and temperature profiles per core are described as follows.

L183. For this core, we adopted the density and temperature profiles presented in Pagani et al. (2007), and are based on
<table><caption>**Table 5.** Optical depth, $\tau$, of H$_2$D$^+$ and N$_2$H$^+$ estimated with RATRAN at the central position of each core.</caption><thead><tr><th>Core</th><th>Species</th><th>τ (QE-BES)</th><th>τ (LP)</th><th>τ (Static)</th><th>τ (SIS)</th></tr></thead><tbody><tr><td>L183</td><td>H<sub>2</sub>D<sup>+</sup> 1<sub>10</sub>-1<sub>11</sub></td><td>0.4</td><td>0.4</td><td>0.4</td><td>3.6</td></tr><tr><td>L694-2</td><td>H<sub>2</sub>D<sup>+</sup> 1<sub>10</sub>-1<sub>11</sub></td><td>0.7</td><td>0.2</td><td>0.3</td><td>0.4</td></tr><tr><td>L1517B</td><td>H<sub>2</sub>D<sup>+</sup> 1<sub>10</sub>-1<sub>11</sub></td><td>0.1</td><td>0.05</td><td>0.08</td><td>0.8</td></tr><tr><td>L1544</td><td>H<sub>2</sub>D<sup>+</sup> 1<sub>10</sub>-1<sub>11</sub></td><td>1.5</td><td>0.4</td><td>0.4</td><td>0.2</td></tr><tr><td></td><td>N<sub>2</sub>H<sup>+</sup> 4-3</td><td>1.6</td><td>0.5</td><td>0.5</td><td>0.2</td></tr><tr><td>L1521f</td><td>H<sub>2</sub>D<sup>+</sup> 1<sub>10</sub>-1<sub>11</sub></td><td>0.08</td><td>0.07</td><td>0.08</td><td>0.02</td></tr><tr><td></td><td>N<sub>2</sub>H<sup>+</sup> 4-3</td><td>0.7</td><td>0.3</td><td>0.6</td><td>0.5</td></tr></tbody></table>

**Notes.** The optical depth, $\tau$, at the central position of each core is found to be <1 for both species for most models.

a detailed 1D radiative transfer modeling to fit both N$_2$H$^+$ and N$_2$D$^+$ emission (IRAM-30m cut observations), consistently with the dust emission (MAMBO; Pagani et al. 2004). In particular, the density profile follows a broken power law with an inner density $\rho_0 = 2.34 \times 10^6$ cm$^{-3}$, with the density proportional to $r^{-1}$ for a radius up to ~4000 au, and to $r^{-2}$ for the rest of the core. The temperature profile has a constant value of 7 K up to 5600 au and increases outward up to 12 K.

**L694-2.** For this core, we adopted the density and temperature profiles presented in Redaelli et al. (2018). According to Arzoumanian et al. (2011), the density profile is given by:

$$ \rho(r) = \left[ \frac{\rho_0}{1 + \left(\frac{r}{r_0}\right)^2} \right]^{\frac{p}{2}} \quad (6) $$

which Redaelli et al. (2018) adjusted for L694-2 to fit its observed Herschel/SPIRE dust continuum at 250, 350, and 500 µm, using radiative transfer modeling. In particular, the inner density $\rho_0 = 1.8 \times 10^5$ cm$^{-3}$, the half maximum density radius $r_0 = 2500$ au, and $p = 3.3$. The temperature profile is 8 K in the center of the core and increases outward to 12 K.

**L1517B.** Here, we adopted the density and temperature profiles presented in Maret et al. (2013) and Tafalla et al. (2004), respectively, based on radiative transfer modeling of the 1.2 mm continuum emission. In particular, the density profile of L1517B follows the form:

$$ \rho(r) = \frac{\rho_0}{1 + \left(\frac{r}{r_0}\right)^\alpha} \quad (7) $$

with a $\rho_0 = 2.2 \times 10^5$ cm$^{-3}$, a half maximum density radius $r_0 = 4.9 \times 10^3$ au, and $\alpha = 2.5$. The temperature profile is also taken from those studies and is assumed to be isothermal and equal to 10 K throughout the core.

**L1544.** For this core, we adopted the most recent density and temperature profiles presented in Chacón-Tanarro et al. (2019), based on fitting the 1.1 and 3.3 mm continuum emission. The temperature profile follows the form:

$$ T(r) = T_{\text{OUT}} - \frac{T_{\text{OUT}} - T_{\text{IN}}}{1 + \left(\frac{r}{r_0}\right)^{1.7}} \quad (8) $$

with outer temperature $T_{\text{OUT}} = 12$ K, inner temperature $T_{\text{IN}} = 6.9$ K, and critical radius $r_0 = 28.07''$. The density profile follows the same form as in Eq. (7), but for $\rho_0 = 1.6 \times 10^6$ cm$^{-3}$, a half maximum density radius $r_0 = 17.3''$, and $\alpha = 2.6$.

**L1521f.** For L1521f, we adopted the density and temperature profiles presented in Crapsi et al. (2004). In particular, the density profile follows the same form as in Eq. (7), but for a central density $\rho_0 = 1 \times 10^6$ cm$^{-3}$, a half maximum density radius $r_0 = 2800$ au, and $\alpha = 2$, while the temperature is assumed to be isothermal and equal to 10 K for the entire core.

### 5.3.3. Abundance derivation: static approach

We find that the initial abundance profiles reproduce the line strength of the central peak position, requiring no modification for L183 and L1521f, and a modification of a factor between 0.8 and 4 for the H$_2$D$^+$ and N$_2$H$^+$ emission, respectively, toward L1544. The static model performs well in reproducing the central peak for all cores with minimum modification of the initial abundance profile. Therefore, we consider the H$_2$D$^+$ abundance determination to be more accurate compared to the one derived from the dynamical models toward the two cores for which we did not have prior knowledge of the H$_2$D$^+$ abundance. Our best fit results in [H$_2$D$^+$] of $2.5 \times 10^{-10}$ for L694-2 and $1.5 \times 10^{-10}$ for L1517B, and similar (within a factor of 2) to what was previously found for the rest of the cores. These results are summarized in Table 4.

## 5.4. RATRAN: results

In Figs. 10–16, we present the azimuthally averaged observed and modeled emission of H$_2$D$^+$ toward the cores for the different physical and kinematic structures predicted by the three dynamical models, LP, SIS, and QE-BES, and the static model. The corresponding plots using the original, before azimuthal averaging, maps of each core can be found in Appendix B. The synthetic spectrum is compared to the observed one at each map position using a simple $\chi^2$ approximation given by:

$$ \chi^2 = \sum_i \frac{(T_{\text{MB,obs}}^i - T_{\text{MB,mod}}^i)^2}{\sigma_{\text{OBS}}^2}, \quad (9) $$

where $T_{\text{MB,obs}}^i$ is the observed main beam temperature per velocity channel (i), $T_{\text{MB,mod}}^i$ is the modeled one, and $\sigma_{\text{OBS}}$ is the observed rms.

### 5.4.1. Individual cores

**L183.** For L183, as seen in Figs. 10 and B.1, the static, QE-BES and LP models all can reproduce the central (0'', 0'') observed peak intensity and line shape with <5% deviation. This similarity is easily justified since the abundance profiles are adjusted accordingly in order to fit the central position. The SIS model, however, predicts ~2 times larger line widths than observed, resulting in an order of magnitude larger $\chi^2$ (~900) compared to the rest of the models and demonstrating the poor quality of the fit produced by this model.

Moving outward, we see at 15'' offsets, that the LP model performs better in reproducing the observed peak intensities with absolute deviations of 10–30%, followed by both the static and QE-BES models with absolute deviations of 30–50%. The SIS model is again significantly worse (>3 times higher $\chi^2$ compared to the other three models) in reproducing the shape and intensity of the observed line profiles. At 30'' offsets, the LP model is still superior to the other models with typical deviations of ~15% from the observed emission, which is significantly
<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Observed azimuthally averaged line profiles for each position of H₂D⁺ maps toward L183, overplotted with the modeled ones for Static, LP, SIS and QE-BES models from RATRAN.</figcaption></figure>

closer to the observations compared to QE-BES and static models that show >40% deviations. The exception to this rule is the 30'' north, where all models underestimate the emission by 50–65%. Similarly, all models underestimate the observed emission by 50–90% at 45'' offsets northward (both east, west, and center). On the other hand, at 45'' offsets southward, the

core is characterized by weaker emission and the LP model can reproduce reasonably well the observations (<15% deviations), which is indicative of a diversion from the spherical symmetry assumption (see Fig. B.1). We note that Pagani et al. (2007) did not report such diversions for other species (e.g., N₂D⁺), but in that paper only the west–east cut is presented. Noticeably,
<figure><img src="image_11.png" /><figcaption>Fig. 11. Same as Fig. 10, but for L694-2.</figcaption></figure>

the QE-BES model performs better than the static sphere at the larger offsets.

We conclude that the LP model performs better in reproducing the observed line profiles and intensities globally for L183, followed by the QE-BES model, while the static model significantly underestimates the emission toward the outer parts of the

core. The SIS model cannot reproduce the intensity and morphology of the line profiles globally, and therefore it is the least suitable among the models.

L694-2. For L694-2, as seen in Figs. 11 and B.2, the static model reproduces the central (0", 0") observed peak intensity
<figure><img src="image_12.png" /><figcaption>**Fig. 12.** Same as Fig. 10, but for L1517B.</figcaption></figure>

and line shape with only 1% deviation, followed by the LP model with <10% deviation. While the QE-BES model can also fairly reproduce the central emission, the SIS model fails to match the observations by producing a line profile with much stronger line wings than the observed, resulting in a χ² >1000, which is two orders of magnitude higher than the LP and static models.

Interestingly, when we examine the modeled versus the observed emission globally, we see that the LP, QE-BES and static models can all reproduce the observed peak intensities equally well with absolute deviations of 2–12% at 15" north-east and south-west from the central peak, while all three underestimate the observed line intensity at 15" south-east and
<figure><img src="image_13.png" /><figcaption>Fig. 13. Same as Fig. 10, but for L1544.</figcaption></figure>

north-west by 30–50%. The SIS model reproduces neither the strength nor the shape of the lines at those offsets. At offsets >30'', the QE-BES model predicts emission that is 30–60% higher compared to the other models. This weak emission, however, is mostly within the rms limits, and therefore not possible to validate with the dataset in hand.

We conclude that for L694-2 we cannot clearly distinguish between LP, QE-BES, and the static models, as their predicted emission is very similar (mostly similar resulting $\chi^2$). Future deeper observations toward the outer parts of the core (e.g., >30'') will reveal if there is a weak emission similar to what the QE-BES model predicts, which will allow us to tell those models
<figure><img src="image_14.png" /><figcaption>**Fig. 14.** Same as Fig. 10, but for N₂H⁺ toward L1544.</figcaption></figure>

apart. The SIS model appears to be the least suitable among the models, failing to reproduce both the intensity and the shape of the line profiles globally.

**L1517b.** For L1517B, as seen in Figs. 12 and B.3, all three QE-BES, LP, and static models can reproduce the central (0'',

0'') observed peak intensity and line profile with <5% deviation, while the SIS model predicts ~15% lower intensity at the source velocity.

When we examine the modeled versus the observed emission globally, we see that both the LP and the static model can reproduce the observed peak intensities with absolute deviations
<figure><img src="image_15.png" /><figcaption>Fig. 15. Same as Fig. 10, but for the protostellar core L1521f.</figcaption></figure>

of 0.5–8% at offsets within 15'', where the detections are more
apparent, while the deviations of the models are >30% at larger
offsets where the observed emission gets weaker. For the same
grid of offsets, the QE-BES model does not perform as well,
showing larger deviations of 25–75%, with the only exception
toward the central position (0'', 0''). Lastly, the SIS model

appears to reproduce the peak intensities with deviations of
5–30% for within 15'', while it gets significantly worse at larger
offsets (up to 80% for >30''). When we look at the line pro-
file of the modeled emission though, it becomes apparent that
the SIS model fails to reproduce the observed emission glob-
ally. The unsuitability is reflected in high χ² at most positions
<figure><img src="image_16.png" /><figcaption>**Fig. 16.** Same as Fig. 10, but for N₂H⁺ toward L1521f.</figcaption></figure>

that are about an order of magnitude higher compared to those of the other three models at the same offsets (Fig. 12). At off-sets with no detected emission (most offsets >30′′), all four models predict some weak emission within the rms limits, which cannot be directly addressed in the current paper. Future observations of higher sensitivity will be able to reveal such emission.

We conclude that for L1517B, both the static and LP models are able to reproduce the observed peak intensities, line profiles, and spatial distribution of the emission equally well, and therefore we cannot clearly discriminate them. The QE-BES model underestimates the emission at most offsets for >50%, while the SIS appears to be the least suitable among the models, and therefore can be safely excluded.
L1544. As seen in Figs. 13 and B.4 for L1544, both the QE-BES and LP models can reproduce the central (0'', 0'') observed peak intensity of $H_2D^+$ and line profile with <3% deviation, followed by the static model (~8%), while the SIS model predicts ~55% lower intensity at the source velocity.

When we examine the modeled versus the observed emission globally, we see that the LP model performs better in reproducing the observed peak intensities of detected emission with absolute deviations of 3–23% at offsets within 30''. For the same grid of offsets, the QE-BES model is the next most suitable model showing deviations of 0.1–52%, with the 0.1% only toward the central position (0'', 0''), followed by the static model with deviations of 8–70%. Lastly, the SIS model does a fairly good job in reproducing the peak intensities with deviations of 2–50%. When we examine the line shape (profile and spatial distribution) of the modeled emission though, it becomes apparent that the SIS model fails to reproduce the observations globally. This failure is clearly depicted by the resulting $\chi^2$, which accounts for all velocity channels of the emission (not just the central peak), and is found to be systematically and significantly higher than the rest of the models at most positions where emission is present (e.g., up to 1000 at 0'' offset, compared to <100 for the rest of the models; Fig. 13). The situation is different at offsets with no detected emission (i.e., most offsets >45''), where the static model is the only one to predict no emission, while all dynamical models predict weak emission but mostly within the rms limits. This discrepancy can be addressed in the future by obtaining higher sensitivity observations, which can reveal the predicted weak emission of the dynamical models. If the static model predictions are correct and there is indeed no emission toward the outer parts of the core, this result can be explained by the more realistic approach on density and temperature profiles, for this core, which are based on dust observations and not simulations.

The situation is less clear when simulating the $N_2H^+$ 4–3 emission. As we see in Figs. 14 and B.5, all models apart from the SIS model can reproduce the observed peak intensity and line profile at the (0'', 0'') and (15'', 15'') offsets (where emission is present). The QE-BES model is also tested by Redaelli et al. (2019) toward L1544, and is found to reproduce the $N_2H^+$ 1–0 and 3–2 transitions at the central dust peak position reasonably well. All models in our study systematically and significantly overestimate the emission by an order of magnitude at all other inner positions and at offsets of 30'' regardless of the direction. Therefore, no model is able to explain the observed $N_2H^+$ 4–3 spatial distribution toward L1544 globally. We conclude that the LP model performs better in reproducing the observed line profiles and intensities of $H_2D^+$ globally for L1544, while the SIS model appears to be the least suitable among the models.

L1521f. For this core, as seen in Figs. 15 and B.6, both QE-BES and LP models can reproduce the central (0'', 0'') observed peak intensity and line profile of $H_2D^+$ with <12% deviation, followed by the static model (~17%), while the SIS model predicts ~75% lower intensity at the source velocity.

When we examine the modeled versus the observed emission globally, we see that apart from two offsets at (0'', 0'') and (-15'', -15''), where the LP, QE-BES, and static models can reproduce the observed peak intensity similarly and with <35% deviations, they systematically and significantly underestimate the peak intensity between 52% and 96% at the rest of the offsets with detections. The SIS is the only model providing a good representation of the observed peak intensity for the majority of offsets within 45'' from the center, with deviations of 8–40%.

Nevertheless, similarly to the other cores, the SIS model fails to reproduce the line profiles resulting in very high $\chi^2$ (>100), when all velocity channels are taken into account (Fig. 15). The exceptions to this trend are the two observed offsets at 45'' north, where the SIS model appears to be the only one to predict the observed emission.

When we examine the $N_2H^+$ emission in Figs. 16 and B.7 we see that all models apart from the SIS model can reproduce the central peak emission within 10%, but overestimate the emission at 15'' offsets by >60%, while they are all consistent with showing no emission at larger offsets (>30''). The amount of available information connected to $N_2H^+$ is not sufficient to draw clear conclusions.

We find that none of the models perform well in reproducing the observed emission globally, both line profiles, and the intensity and distribution of $H_2D^+$, toward L1521f. This outcome may be due to the fact that L1521f is a more evolved (protostellar) object compared to the rest of the sample of prestellar cores, and therefore it is characterized by more complex dynamics and physical structure. Protostellar cores, for example, are known to drive powerful outflows, and therefore the tested dynamical models that account only for infall motions may not be expected to reproduce the observations consistently with what we see. The situation is similar for the model of the static envelope. The good quality fit of the SIS model at the reported offsets may indicate that the dynamical evolution from a prestellar to a protostellar core is better represented by the SIS model rather than the LP model in opposition to what we see for prestellar cores. Given the poor ability of the model to reproduce the emission at a more global scale, however, we consider this alignment to be rather coincidental. The less problematic fit we get toward $N_2H^+$ indicates that $H_2D^+$ is dynamically decoupled from $N_2H^+$, with the former requiring a more complex structure.

### 5.4.2. Summary of RATRAN results

We find that the LP flow can reproduce the observed line profile and intensity of $H_2D^+$ fairly well on a global scale (most offsets), and appears to be superior to the static, SIS, and QE-BES models for two out of the four prestellar cores, L183 and L1544. The influence of the abundance profiles on the simulated line intensity makes clear differentiation among those models challenging. For L1517B, the situation is less clear, with the static sphere reproducing the observations equally as well as the LP model. For L694-2, we cannot clearly distinguish between the LP, QE-BES, and static sphere models. On the other hand, all models fail to reproduce the observed line profiles and distribution of $H_2D^+$ toward the protostellar core L1521f, indicating that the predicted model dynamics taking place during the prestellar phase are not suitable for a more evolved object, where outflows and a central heating source are present. These findings are summarized in Fig. 17 at offsets along a cut in the south-west/north-east direction.

In addition, our modeling clearly demonstrates that the SIS dynamical model fails to reproduce the observed line profile for all cores, independent of the initial abundances. While the density profile of the SIS model shows a similar slope to the other two models (all models follow ~1/r² profile), the exact density profile results in a smaller mass in the same volume of gas compared to the other two dynamical models. As a result, the overall density distribution of the SIS model is not sufficient to excite the $H_2D^+$ transition in the amount required to match the observations. This behavior is clearly seen in the abundance calculations, where we find the $H_2D^+$ abundance needs to be 2–4 orders of
<figure><img src="image_17.png" /><figcaption>**Fig. 17.** Observed line profiles of H₂D⁺, overplotted with the modeled ones for LP, SIS, and QE-BES models from RATRAN at selected offsets: (0", 0"), (15", 15"), and (-15", -15").</figcaption></figure>

magnitude higher for the SIS model compared to the other two models for all sources. The observed differences on the abundances can be depicted in the resulting column densities as well. For a given volume density an increase of a factor of 100 in the abundance corresponds to an increase in the molecular density by the same factor. For example, in the case of the H₂D⁺ emission toward L1544 the LP model results in a column density of 1.6 × 10¹³ cm⁻², while the SIS model results in a column density of 3.2 × 10¹⁵ cm⁻². Referring back to the RADEX analysis (Sect. 4.2), we find N(H₂D⁺) = 2.0 × 10¹³ cm⁻², which is in very good agreement with the LP model but 2 orders of magnitude lower than the column density from the SIS model. This result demonstrates once more the poor performance of the SIS model compared to the LP and QE-BES models.

Moreover, the QE-BES model gives a core mass of ~10 $M_\odot$, the LP model of ~27 $M_\odot$ and the SIS model of 4 $M_\odot$, while the masses of the cores derived from dust emission are 4–10 $M_\odot$. We conclude that the differences in the modeled versus observed masses affect our abundance determinations, resulting in higher or lower H₂D⁺ abundance depending on whether the modeled mass for a defined volume is higher or lower, respectively, than observed. To examine the effects of the assumed mass in the reproduction of the observations we notice that L1517B is as massive as 4 $M_\odot$ (i.e., the same mass as the SIS model has), but the SIS model still fails to reproduce the observations; therefore, our conclusion on excluding the SIS model as a possible dynamical model to explain the observed H₂D⁺ emission is independent from the mass assumptions. Similarly, the LP model has

a mass ~3 times higher than the core masses in our sample, yet it is able to reproduce the observations at all prestellar cores. We therefore conclude that although the mass differences introduce higher uncertainties in the abundance determinations, they do not affect the line profiles. The line profiles depend on the velocity structure of the models, making our analysis and conclusions in differentiating between the dynamical models robust. Referring to the results presented in Sect. 3.2, we found that there are significant nonthermal contributions toward all cores, which are traditionally attributed to turbulence and magnetic pressure. The ability of the LP flow to reproduce the observations without taking those phenomena into account, unlike the adopted QE-BES model (i.e., turbulence), indicates that they are not significant in driving the core dynamics.

## 6. Discussion

In this section, we discuss the main findings of this paper. In particular, we attempt to explain the observed line broadening, the presence of N₂H⁺ toward only two of the five cores in our sample, and, lastly, the core collapse dynamics, which is a debatable topic both in theory and observations.

### 6.1. Line broadening

We first discuss our results from the line profile analysis presented in Sect. 3.2. We find that the line width of both N₂H⁺ and H₂D⁺ cannot be explained by thermal broadening alone, and the
nonthermal contributions are variable with offsets and overall higher for N₂H⁺ than H₂D⁺ (Sect. 3.2). The observed varia- tions with offsets of the nonthermal contributions are as high as ~130% and beyond the uncertainties. Traditionally, nonther- mal motions are explained by presence of magnetic pressure or turbulence. In this study, we show that the LP flow can repro- duce the line width of the H₂D⁺ observations in multiple offset locations (e.g., see Fig. 10), which by definition (free-fall) does not take into account turbulence or magnetic fields; therefore, the nonthermal broadening can be merely due to infall/compression motions, characteristic of the dynamical structure of prestellar cores.

## 6.2. N₂H⁺ 4–3 emission as an evolution tracer

Here, we discuss the origin of the N₂H⁺ 4–3 emission. We detect the N₂H⁺ 4–3 transition toward two out of the five cores in our sample, in contrast to the ground ortho-H₂D⁺ 1₁₀–1₁₁ emission which is present toward all cores. In addition, H₂D⁺ is more spa- tially extended than N₂H⁺, which is traced only within 15" of the core centers.

To understand the spatial differences, and the presence or absence of these species toward these cores, we considered the critical densities of the observed transitions, calculated based on the available molecular data. In particular, the observed *ortho*-H₂D⁺ 1₁₀–1₁₁ at 372.4 GHz is characterized by a critical density $n_{cr} = 1.2 \times 10^5$ cm⁻³ at 10 K, while the observed N₂H⁺ 4–3 transition at 372.67 GHz is characterized by a critical density $n_{cr} = 7.7 \times 10^6$ cm⁻³ at 10 K, which is more than an order of magnitude higher. We also note that the N₂H⁺ 1–0 transition at 93.17 GHz, which has a comparable critical density to the ground H₂D⁺ transition ($n_{cr} \sim 1.4 \times 10^5$ cm⁻³ at 10 K), has been previ- ously detected toward all cores in our sample (Crapsi et al. 2005). We conclude that the observed spatial differences between the two species are due to the higher densities traced by N₂H⁺ 4–3, which are unsurprisingly located toward the more central regions of the cores (see, e.g., density profiles in Sect. 5.3.2). Indeed, this interpretation is supported by the observed lower J transitions 3-2 and 1–0 of N₂H⁺ toward L1544 (Redaelli et al. 2019), which are found to be up to a factor of four more extended compared to the 4–3 transition.

The absence of N₂H⁺ 4–3 emission for L694-2 and L1517B can be explained by their lower central densities on the order of ~10⁵ cm⁻³, namely, an order of magnitude lower than the critical density of N₂H⁺ 4–3. Similarly, the presence of this tran- sition can be explained for L1521f and L1544, given their higher central densities of ~10⁶ cm⁻³. Along this line of reasoning, the non-detection of N₂H⁺ 4–3 toward L183 is more puzzling, given that its central density is found to be similar to those of L1521f and L1544, indicating that its central density may be overestimated.

We conclude that N₂H⁺ 4–3 could be used as a chemical clock, in the sense that it traces a fairly high density environment achievable only in more evolved prestellar/early protostellar core phases. Our conclusion is supported by the presence of this emis- sion toward the only protostellar (more evolved) core in our sample (L1521f). This approach would mean that L694-2 and L1517B are less-evolved cores.

## 6.3. Differentiating between the theories of collapse

Here, we discuss our findings from the RATRAN analysis pre- sented in Sect. 5, where we compare the performance of the different dynamical models of collapse. We find that the LP

model can mostly reproduce the observed line profiles and spa- tial distribution of H₂D⁺, followed in suitability by the QE-BES model. Meanwhile, the SIS model fails to reproduce the line width globally, predicting a very prominent double peak line profile toward the central position of all cores, in contrast to the observed single peak. Although the LP model seems to be superior to the QE-BES model for at least two of the prestel- lar cores in the sample, we cannot distinguish between the two models with high confidence. The widths of the line profiles can be mostly reproduced by both models, while in addition the emission resulting from adopting LP conditions is closer to the observed line intensity on a more global scale compared to the QE-BES model. All models fail to reproduce the emission from the one protostellar core in our sample, L1521f, which is indicative of the more complex structure of a protostellar object with respect to a prestellar core, namely, presence of winds and outflows, disk, internal radiation, different grain sizes, and chemistry.

The abundance profiles have a significant effect on the simu- lated line intensity. We consider the abundance derivation from the dynamical models to be less accurate compared to those of the static model, due to the differences in the total core mass that each model implies. As we explained in Sect. 5.4 the exact density distribution over a specific volume (i.e., different mass) affects the degree of excitation of H₂D⁺ and therefore its line intensity. For example, to match the H₂D⁺ line intensity of L1544 (~8 M☉), the H₂D⁺ abundance had to be lower in the case of the LP model (>10 M☉) compared to those of the QE-BES (~10 M☉) or SIS (<10 M☉) models. The static model, on the other hand, is based on continuum observations and therefore it accounts for the observed total mass of each individual core. Breaking the observed discrepancy between masses derived from continuum and gas observations will allow proper differentiation between the LP and QE-BES models. To achieve this goal, the dynam- ical models need to be “tailored” to match the observed mass per individual core. Such a procedure, however, requires prior knowledge of the exact evolutionary status of each core, and is beyond the scope of this work.

Unlike the line intensity, which in the optically thin regime ($\tau < 1$; see Table 5) highly depends on the adopted abund- ances, the shape of the line profile for a given abundance profile depends mostly on the velocity structure, which is unique per dynamical model. Therefore, we can safely exclude the SIS model as a possible dynamic model to explain H₂D⁺ emis- sion for all cores. The situation is similar for N₂H⁺, where present. Our qualitative results remain even when a constant abundance profile is adopted for each core.

To understand the differences in the shape of the line profile between the dynamical models, we explore the physical loca- tion of the H₂D⁺ emission and we compare it to other tracers from previous work. The observed 1₁₀–1₁₁ ground-state transi- tion of ortho-H₂D⁺ at 372.4 GHz is characterized by a critical density $n_{cr} = 1.2 \times 10^5$ cm⁻³ at 10 K, serving as a rather dense gas tracer. Looking at Fig. 8, we see that in the case of the SIS model, the observed H₂D⁺ transition would better probe the innermost regions (<1000 au). In the case of the QE-BES and LP models, however, it traces a greater volume of the gas (<10 000 au), but not the outer layers up to 100 000 au. The obvi- ous mismatch between the predicted and observed H₂D⁺ line profiles produced assuming an SIS collapse suggests that high- velocity inflow in the inner parts of the cores fails to explain the H₂D⁺ emission. The line profiles produced by the LP and QE-BES models can mostly reproduce the observed line profiles, suggesting that more moderate contraction motions in the inner
regions are more applicable (e.g., by three orders of magnitude
difference in velocity). Keto et al. (2015) performed a similar
analysis for H₂O and C¹⁸O toward L1544 and found that the
QE-BES model reproduces the shape and intensity of the line
profiles for both species, while the LP model produced a more
complex profile for the C¹⁸O similar to what we see for H₂D⁺ in
the SIS model in this study. To understand better the observed
differences in the line profiles produced from different models,
we need to consider the origin of the emission. In particular, the
C¹⁸O transition presented in Keto et al. (2015) is characterized
by a critical density about two orders of magnitude lower than
that of the H₂D⁺ line examined here. Therefore, it is better at
tracing less dense gas compared to the ground H₂D⁺ transition.
Looking at Fig. 8, we see that the outer regions of the core in
the LP hypothesis are characterized by much higher velocities,
explaining the broad wings modeled in C¹⁸O. In conclusion,
the QE-BES model can explain the observed line profiles of
C¹⁸O, H₂D⁺ and H₂O, both in shape and intensity, for the cen-
tral position of L1544, while the LP and SIS models cannot
simultaneously explain the observed emission.

In addition, we examine the performance of the static sphere.
We find that the static model, which is based on the continuum
observations of each core, can also reproduce the observa-
tions fairly well and similarly to the LP model for at least one
core (L 1517b), demonstrating the importance of the adopted
abundance profiles. This model, however, does not take into
account infall motions. With only a single available transition
of H₂D⁺, a more accurate determination of the abundance pro-
file and its proper coupling to the specific core dynamics is very
challenging.

We conclude that coupling gas core dynamics and chem-
istry with dust properties is crucial to differentiating between the
QE-BES and LP models, and therefore, should be considered in
future studies. In particular, chemical modeling is essential to
constrain the chemical age of each core, which will subsequently
lead to the coupling of the dynamical models to the observed
mass of each individual core.

To be able to draw stronger conclusions it is necessary to
expand this analysis to a larger sample of known prestellar
and protostellar cores. If the LP model is confirmed toward a
larger sample of cores, its specific initial conditions indicate that
the prestellar cores may form by compression or accretion of
gas from larger scales (>10 000 au). On the other hand, if the
QE-BES model is confirmed, it means that quasi-hydrostatic
cores can exist within turbulent ISM, which would be against
predictions from numerical simulations (Tilley & Pudritz 2004).

**7. Summary–conclusions**

This paper presents an analysis of H₂D⁺ 1₀⁻₁₁ and N₂H⁺ 4⁻3 toward one protostellar and four prestellar cores using 2' × 2' maps from JCMT. We present a basic line analysis for all positions in the maps where emission is detected. We estimated the average physical conditions (T<sub>KIN</sub> and N<sub>COL</sub>) toward the center of the cores using the non-LTE radiative transfer code, RADEX, incorporating the latest collisional data available for both species. The derived kinetic temperatures per core are used to assess the thermal contributions on the observed line widths. Finally, we tested three different dynamical scenarios, using the advanced radiative transfer code RATRAN to model the spatial distribution of the gas traced by the optically thin H₂D⁺ 1₀⁻₁₁ and N₂H⁺ 4⁻3 line emission and reproduce the sizes of the observed maps. The latter allowed us to not only test which of the available dynamical scenarios is the most favorable for the

prestellar core phase but also to derive abundance profiles of
the observed species. Our main conclusions are summarized as
follows.

- Thermal broadening alone is not sufficient to explain the line widths of H₂D⁺ and N₂H⁺ measured toward the majority of the positions in the maps with clear detections. This finding suggests that nonthermal effects play a significant role at the very early stages of star formation. The nonthermal contributions can be merely attributed to the infall/compression motions predicted by the LP flow. Therefore, in that scenario, turbulent and magnetic pressure contributions are not required to explain the observed line width.

- The nonthermal contributions are generally found to be stronger in the innermost parts of L694-2 and L1544 (<15"; with the exceptions of few offsets which are similar) and are found to be higher for N₂H⁺ than H₂D⁺. We report variations to nonthermal contributions up to a factor of 2.5 at all cores for different offsets and beyond the uncertainties. The lack of sufficient N₂H⁺ detections at offset positions does not allow a meaningful investigation on nonthermal variations with radius for that emission.

- N₂H⁺ 4–3 was only detected toward two of the cores, L1521f and L1544. Given that L1521f is known to be a VeLLO, and therefore more evolved than a prestellar core, and given that both cores are located at the same distance, these detections indicate that L1544 may also be more evolved compared to the rest of the sample. In that case, N₂H⁺ 4–3 can be used as a chemical clock, tracing the higher densities achieved during core evolution.

- The RADEX analysis results in similar mean column densities of H₂D⁺ and N₂H⁺ (~10¹³ cm⁻²) for L1521f and L1544.

- The RATRAN analysis reveals that the LP flow is superior in reproducing the H₂D⁺ spatial distribution and line profiles toward two of the prestellar cores, while for the other two it is followed by the QE-BES and static scenarios, without being able to differentiate clearly between those.

- The SIS model fails to reproduce the shape and strength of the line profiles for the majority of the positions toward all cores and, therefore, does not appropriately describe the dynamical nature of these cores.

- The RATRAN analysis reveals that none of the scenarios can reproduce the H₂D⁺ emission toward the protostellar core, L1521f, indicating that more complex dynamical and chemical structures need to be considered, including outflow activity and radiation or wind pressure from the central heating source.

We explored the most prominent dynamical scenarios based on
spherical contraction to explain the H₂D⁺ and N₂H⁺ (where
present) emission toward four prestellar and one protostellar
core. With the observations and sample in hand we cannot clearly
differentiate between the LP flow and QE-BES models, but can
clearly exclude the SIS model. Our conclusions remain consis-
tent after adopting a jump-drop or a constant abundance profile.
To differentiate between the other collapse scenarios is neces-
sary to construct dynamical models adapted to the individual
observed core mass. In addition, testing the models will benefit
from transitions of more molecular species in both the optically
thin and optically thick regimes (e.g., van der Tak et al. 2005)
and a larger sample of cores. Such data will inform us if the gas
collapses as a free-fall (LP) or if turbulence is indeed crucial for
the gas dynamics of cores.

In future, non-spherical models, such as oblate ones resulting
from ambipolar diffusion, should also be considered. Applying
chemical modeling tailored to the physical structure of each
source to fit the H₂D⁺ and N₂H⁺ observed abundances, will pro-
vide us with the age (e.g., Brünken et al. 2014) and therefore the
evolutionary status of each core, allowing us to link the chemi-
cal to the dynamical evolution of prestellar cores overall. Lastly,
high angular resolutions to reach the emission in the inner 100 au
(e.g., ALMA) of the prestellar cores should prove to be a power-
ful tool to explore the gas dynamics, but the observed chemical
inner holes toward those objects and the heavy freeze-out, make
such observations very challenging.

*Acknowledgements.* We would like to thank the anonymous referee for providing helpful comments and suggestions that improved the paper. EK is funded by the STFC (ST/P00041X/1). We thank Tim van Kempen and Werner Salomons for useful contributions during data reduction and analysis.

References

André, P. 2017, *Comptes Rendus Geosci.*, **349**, 187

André, P., Men'shchikov, A., Bontemps, S., et al. 2010, A&A, 518, A102

Arzoumanian, D., André, P., Didelon, P., et al. 2011, A&A, 529, A6

Beichman, C. A., Myers, P. C., Emerson, J. P., et al. 1986, ApJ, 307, 337

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

Black, J. H., van Dishoeck, E. F., Willner, S. P., & Woods, R. C. 1990, ApJ, 358, 459

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

Bourke, T. L., Myers, P. C., Evans, Neal J., I., et al. 2006, ApJ, 649, L37

Broderick, A. E., & Keto, E. 2010, ApJ, 721, 493

Brünken, S., Sipilä, O., Chambers, E. T., et al. 2014, Nature, 516, 219

Buckle, J. V., Hills, R. E., Smith, H., et al. 2009, MNRAS, 399, 1026

Caselli, P., Benson, P. J., Myers, P. C., & Tafalla, M. 2002, ApJ, 572, 238

Caselli, P., van der Tak, F. F. S., Ceccarelli, C., & Bacmann, A. 2003, A&A, 403, A37

Caselli, P., Vastel, C., Ceccarelli, C., et al. 2008, A&A, 492, 703

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

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

Crapsi, A., Caselli, P., Walmsley, C. M., et al. 2004, A&A, 420, 957

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

Crutcher, R. M., Nutter, D. J., Ward-Thompson, D., & Kirk, J. M. 2004, ApJ, 600, 279

Currie, M. J., Berry, D. S., Jenness, T., et al. 2014, ASP Conf. Ser., 485, 391

Di Francesco, J., Johnstone, D., Kirk, H., MacKenzie, T., & Ledwosinska, E. 2008, ApJS, 175, 277

Ebert, R. 1957, ZAp, 42, 263

Enoch, M. L., Evans, Neal J., I., Sargent, A. I., et al. 2008, ApJ, 684, 1240

Friesen, R. K., Di Francesco, J., Bourke, T. L., et al. 2014, ApJ, 797, 27

Friesen, R. K., Pon, A., Bourke, T. L., et al. 2018, ApJ, 869, 158

Fu, T.-M., Gao, Y., & Lou, Y.-Q. 2011, ApJ, 741, 113

Goodman, A. A., Barranco, J. A., Wilner, D. J., & Heyer, M. H. 1998, ApJ, 504, 223

Harju, J., Haikala, L. K., Lehtinen, K., et al. 2006, A&A, 454, L55

Ho, P. T. P., Martin, R. N., & Barrett, A. H. 1978, ApJ, 221, L17

Hogerheijde, M. R., & van der Tak, F. F. S. 2000, A&A, 362, 697

Hugo, E., Asvany, O., & Schlemmer, S. 2009, J. Chem. Phys., 130, 164302

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

Kirk, H., Dunham, M. M., Di Francesco, J., et al. 2017, ApJ, 838, 114

Krumholz, M. R., McKee, C. F., & Klein, R. I. 2005, Nature, 438, 332

Larson, R. B. 1969, MNRAS, 145, 271

Lique, F., Tobola, R., Klos, J., et al. 2008, A&A, 478, 567

Lique, F., Daniel, F., Pagani, L., & Feautrier, N. 2015, MNRAS, 446, 1245

Maloney, P. 1988, ApJ, 334, 761

Maret, S., Bergin, E. A., & Tafalla, M. 2013, A&A, 559, A53

Myers, P. C., Ladd, E. F., & Fuller, G. A. 1991, ApJ, 372, L95

Pagani, L., Bacmann, A., Motte, F., et al. 2004, A&A, 417, 605

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

Pagani, L., Vastel, C., Hugo, E., et al. 2009, A&A, 494, 623

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

Penston, M. V. 1969, MNRAS, 144, 425

Pickett, H. M., Poynter, R. L., Cohen, E. A., et al. 1998, J. Quant. Spec. Rad. Transf., 60, 883

Pineda, J. E., Offner, S. S. R., Parker, R. J., et al. 2015, Nature, 518, 213

Pudritz, R. E., & Kevlahan, N. K. R. 2013, Phil. Trans. Roy. Soc. London A, 371, 20120248

Redaelli, E., Bizzocchi, L., Caselli, P., et al. 2018, A&A, 617, A7

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

Schnee, S., Enoch, M., Johnstone, D., et al. 2010, ApJ, 718, 306

Schnee, S., Sadavoy, S., Di Francesco, J., Johnstone, D., & Wei, L. 2012, ApJ, 755, 178

Shu, F. H. 1977, ApJ, 214, 488

Tafalla, M., & Hacar, A. 2015, A&A, 574, A104

Tafalla, M., Mardones, D., Myers, P. C., et al. 1998, ApJ, 504, 900

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

Tafalla, M., Santiago-García, J., Myers, P. C., et al. 2006, A&A, 455, 577

Takahashi, S., Ohashi, N., & Bourke, T. L. 2013, ApJ, 774, 20

Tilley, D. A., & Pudritz, R. E. 2004, MNRAS, 353, 769

Tokuda, K., Onishi, T., Matsumoto, T., et al. 2016, ApJ, 826, 26

Tokuda, K., Onishi, T., Saigo, K., et al. 2017, ApJ, 849, 101

Tokuda, K., Onishi, T., Saigo, K., et al. 2018, ApJ, 862, 8

van der Tak, F. F. S., Caselli, P., & Ceccarelli, C. 2005, A&A, 439, 195

van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627

Vastel, C., Caselli, P., Ceccarelli, C., et al. 2006, ApJ, 645, 1198

Ward-Thompson, D., Scott, P. F., Hills, R. E., & Andre, P. 1994, MNRAS, 268, 276
# Appendix A: Line analysis per core

Our line profile findings per source are described as follows.

**L183.** For L183 (Fig. 3), the $V_{LSR}$ ranges between $2.25 \pm 0.05$ and $2.48 \pm 0.04$ km s⁻¹, and within the upper and lower velocities reported by Ho et al. (1978) and Caselli et al. (2008). The fact that the measured $V_{LSR}$ is not constant but varies with offset indicates local changes in the global kinematics of the core. The FWHM ranges from $0.32 \pm 0.08$ km s⁻¹ to $0.58 \pm 0.07$ km s⁻¹. With the exception of only few offsets (>40'' southwest and northwest), the observed FWHM cannot be explained by thermal broadening alone (beyond the errors). The FWHM is overall up to ~90% larger than that expected from thermal broadening alone. Lastly, the $T_A^*$ varies between $0.12 \pm 0.03$ K and $0.48 \pm 0.03$ K, with the highest measurements found toward the center and a gradual decrease outward. The observed variations at similar offsets, as well as the fact that there is an increase in $T_A^*$ at the most outward position, can be explained by inhomogeneities present in the medium.

**L694-2.** The $V_{LSR}$ for L694-2 (see Fig. 4) ranges from $9.47 \pm 0.04$ km s⁻¹ to $9.65 \pm 0.03$ km s⁻¹. There is not a specific trend on the observed $V_{LSR}$ variations, but the difference of up to ~0.2 km s⁻¹ indicates local kinematic differences throughout the core. In this case, the FWHM ranges from $0.39 \pm 0.1$ to $0.63 \pm 0.09$ km s⁻¹. With the exception of a single offset at ~42'', where thermal broadening alone can explain the observed FWHM, the nonthermal contributions are significant at all other offsets inward (<40''). The line is generally up to ~96% broader than expected from thermal broadening alone. The decreased FWHM at a larger offset, may indicate a decrease in nonthermal contributions at the outermost parts of the core. Finally, the $T_A^*$ values gradually decrease with offset, with local variations in $T_A^*$ pointing toward medium inhomogeneities.

**L1517B.** If one compares Fig. 5 to Figs. 3–4, 6–7, it is noticeable that the H₂D⁺ emission of L1517B is less spatially extended than seen toward the rest of the cores; Indeed, the line emission in this core is rather concentrated around the central position of the 850 µm dust peak (<30''). The $V_{LSR}$ of L1517B ranges from $5.73 \pm 0.03$ to $5.82 \pm 0.03$ km s⁻¹, and decreases with increasing offset. The FWHM ranges from $0.36 \pm 0.06$ to $0.49 \pm 0.06$ km s⁻¹. Thermal broadening alone can explain the observed FWHM at most offsets for core temperatures >10 K. This temperature is higher by a factor of 2 compared to the one we determined in Sect. 4.2 for this source, but still within the acceptable temperature range of prestellar cores. The line is ~85% broader than expected from thermal broadening alone. Lastly, the $T_A^*$ is once again the highest in the central position and decreases with increasing offset without signs of inhomogeneities.

**L1544.** The H₂D⁺ analysis for this source (Fig. 6) shows that the $V_{LSR}$ ranges between $7.07 \pm 0.03$ and $7.21 \pm 0.02$ km s⁻¹, with a mean $V_{LSR}$ value of $7.15 \pm 0.03$ km s⁻¹. All measurements are within the range previously reported by Ho et al. (1978) and Caselli et al. (2008). The FWHM ranges from $0.38 \pm 0.07$ to $0.61 \pm 0.08$ km s⁻¹, which is systematically higher than the thermal broadening contributions (e.g., a mean value of $0.48 \pm 0.07$ km s⁻¹ is ~40% higher). The observed widths are statistically larger than the expected thermal width in five different offsets, including the central position. As mentioned earlier, we attribute this difference to nonthermal contributions. With

the exception of one offset at ~40'' north, the nonthermal contributions are higher in the center compared to the rest of the offsets. The line is generally up to ~90% broader than expected from thermal broadening alone. Finally, the highest value of $T_A^*$ (0.55 ± 0.02 K) is found in the central position, while it generally decreases with increasing offset. The fact that at certain offsets, the $T_A^*$ varies by a factor >2 and beyond the errors is indicative of local inhomogeneities in the medium (in density and/or temperature).

Given that N₂H⁺ 4–3 emission was detected toward only two positions, we do not show a plot associated with its analysis but summarize our results as follows. The $V_{LSR}$ is found to be ~7.20 ± 0.04 km s⁻¹, that is, within the uncertainties of the corresponding H₂D⁺ value. As in the case of H₂D⁺, thermal broadening does not fully account for the FWHM of N₂H⁺ (0.4 ± 0.1 km s⁻¹), as the observed FWHM is broader by a factor of ~4. In the case of N₂H⁺, nonthermal widths are higher by at least a factor of ~2.5 compared to those of H₂D⁺. Therefore, N₂H⁺ emission appears to stem from a medium of higher nonthermal motions compared to H₂D⁺ or that these mechanisms somehow affect N₂H⁺ more. Finally, both detections of N₂H⁺ have very low $T_A^*$ values, with a maximum value of $0.12 \pm 0.03$ K, a factor of 5 lower than the H₂D⁺ peak intensities.

**L1521f.** This core (see Fig. 7) is the only known protostellar core in the sample (Bourke et al. 2006) and the one out of two for which N₂H⁺ 4–3 is detected. The line analysis of H₂D⁺ is presented in Fig. 7. The $V_{LSR}$ ranges between 6.41 ± 0.02 and 6.54 ± 0.03 km s⁻¹. The FWHM values range from 0.34 ± 0.05 to 0.55 ± 0.07 km s⁻¹. The FWHM can be explained by thermal broadening at the central position and higher offsets (>30''). At intermediate offsets (~20''), however, the FWHM is up to a factor of 2 larger than that expected from thermal broadening. While magnetic and turbulent contributions are also likely to play a role toward this core, we recall that this core is at a later evolutionary state compared to the rest of the sample. Indeed, its known compact outflow (Takahashi et al. 2013) may also contribute to the observed line broadening, although H₂D⁺ is not known to trace outflow activity. $T_A^*$ fluctuates between 0.10 ± 0.02 and 0.30 ± 0.05 K without showing a specific trend. We attribute the observed variations to local changes in density and/or temperature (e.g., a clumpy medium; Tokuda et al. 2018).

As in L1544, N₂H⁺ emission was detected toward only two positions, and therefore we do not show an associated offset plot. The $V_{LSR}$ took values of 6.47 ± 0.03 and 6.56 ± 0.04 km s⁻¹, are again within the uncertainties of the corresponding H₂D⁺ values. At the central position, the FWHM is 0.6 ± 0.1 km s⁻¹, therefore broader by a factor of 5 compared to the width expected from thermal broadening. The nonthermal contributions decrease by 65% at the offset position. Therefore, the nonthermal effects are higher in the center of the core. This behavior is in contrast to what we saw in the H₂D⁺ case, where the central position did not require nonthermal effects to explain the observed FWHM. Therefore, as in the case of L1544, the nonthermal contributions seem to have a more significant effect on N₂H⁺ emission compared to H₂D⁺ emission. In the case of L1521f, this behavior could indicate that N₂H⁺ emission originates from outflow cavities, as opposed to H₂D⁺ emission, or that the magnetic pressure and turbulence affect the two species differently. Finally, at the offset position, the measured $T_A^*$ of the N₂H⁺ line (0.30 ± 0.04 K) is very similar to that of the H₂D⁺ line.
**Appendix B: RATRAN results: original dataset**

In Figs. 10–16, we present the original, before azimuthal averaging, maps of H₂D⁺ 1₁₀–1₁₁ and N₂H⁺ 4–3 toward the cores, overplotted with the modeled ones for Static, LP, SIS and QE-BES models from RATRAN.

<figure><img src="image_18.png" /><figcaption>**Fig. B.1.** Observed line profiles before azimuthal averaging for each position of H₂D⁺ maps toward L183, overplotted with the modeled ones for Static, LP, SIS and QE-BES models from RATRAN.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>L694-2 (H₂D⁺)</figcaption></figure>

Fig. B.2. Same as Fig. B.1, but for L694-2.
<figure><img src="image_20.png" /><figcaption>L1517b ($H_2D^+$)</figcaption></figure>

Fig. B.3. Same as Fig. B.1, but for L1517B.
<figure><img src="image_21.png" /><figcaption>L1544 (H₂D⁺)</figcaption></figure>

**Fig. B.4.** Same as Fig. B.1, but for L1544.
<figure><img src="image_22.png" /><figcaption>L1544 (N₂H⁺)</figcaption></figure>

Fig. B.5. Same as Fig. B.1, but for N₂H⁺ toward L1544.
<figure><img src="image_23.png" /><figcaption>L1521f (H₂D⁺)</figcaption></figure>

Fig. B.6. Same as Fig. B.1, but for L1521f.
<figure><img src="image_24.png" /><figcaption>L1521f ($N_2H^+$)</figcaption></figure>

Fig. B.7. Same as Fig. B.1, but for $N_2H^+$ toward L1521f.