# Hunting pre-stellar cores with APEX IRAS16293E (Oph464)

S. Spezzano, E. Redaelli, P. Caselli, O. Sipilä, J. Harju, F Lique, D

Arzoumanian, J E Pineda, F. Wyrowski, A. Belloche

► To cite this version:

S. Spezzano, E. Redaelli, P. Caselli, O. Sipilä, J. Harju, et al. Hunting pre-stellar cores with APEX IRAS16293E (Oph464). Astronomy & Astrophysics - A&A, 2025, 694, pp.A27. 10.1051/0004-6361/202452625. hal-04984623

HAL Id: hal-04984623

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

Submitted on 10 Mar 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
# Hunting pre-stellar cores with APEX
IRAS16293E (Oph464)

S. Spezzano<fnref n="1" />,*, E. Redaelli<fnref n="1" />,<fnref n="2" />, P. Caselli<fnref n="1" />, O. Sipilä<fnref n="1" />, J. Harju<fnref n="1" />,<fnref n="3" />, F. Lique<fnref n="4" />, D. Arzoumanian<fnref n="5" />, J. E. Pineda<fnref n="1" />, F. Wyrowski<fnref n="6" />, and A. Belloche<fnref n="6" />

Received 15 October 2024 / Accepted 16 December 2024

## ABSTRACT

**Context.** Pre-stellar cores are the first steps in the process of star and planet formation. However, the dynamical and chemical evolution of pre-stellar cores is still not well understood. Our partial knowledge of the chemical and physical structure of pre-stellar cores, as well as how they are fed and influenced by the surrounding environment, limits the level of knowledge that we can achieve at later stages in the star and planet formation process, from protostellar cores to exoplanets.

**Aims.** Our aims are to estimate the central density of the pre-stellar core IRAS16293E and to carry out an inventory of molecular species towards the density peak of the core.

**Methods.** We observed high-J rotational transitions of $N_2H^+$ and $N_2D^+$, and several other molecular lines towards the dust emission peak using the Atacama Pathfinder EXperiment (APEX) telescope, and derived the density and temperature profiles of the core using far-infrared surface brightness maps from Herschel. The $N_2H^+$ and $N_2D^+$ lines were analysed by non-local thermodynamic equilibrium (non-LTE) radiative transfer modelling.

**Results.** Our best-fit core model consists of a static inner region, embedded in an infalling envelope with an inner radius of approximately 3000 au (21" at 141 pc). The observed high-J lines of $N_2H^+$ and $N_2D^+$ (with critical densities greater than $10^6 \text{ cm}^{-3}$) turn out to be very sensitive to depletion; the present single-dish observations are best explained with no depletion of $N_2H^+$ and $N_2D^+$ in the inner core. The $N_2D^+/N_2H^+$ ratio that best reproduces our observations is 0.44, one of the highest observed to date in pre-stellar cores. Additionally, half of the molecules that we observed are deuterated isotopologues, confirming the high level of deuteration towards this source.

**Conclusions.** Non-LTE radiative transfer modelling of $N_2H^+$ and $N_2D^+$ lines proved to be an excellent diagnostic of the chemical structure (i.e. molecular freeze-out) and dynamics (infall velocity profile) of a pre-stellar core. Probing the physical conditions immediately before the protostellar collapse is a necessary reference for theoretical studies and simulations with the aim of understanding the earliest stages of star and planet formation and the timescale of this process.

**Key words.** ISM: clouds – ISM: lines and bands – ISM: molecules

## 1. Introduction

Molecular emission is of paramount importance when we try to understand the process of star formation. The wealth of information contained in molecular spectra from star-forming regions encompasses the physical processes driving the phenomena and the chemical budget available for the forming planets. Pre-stellar cores are gravitationally bound objects on the verge of forming a protostar. They represent the initial conditions in stars and planets formation (e.g. Crapsi et al. 2007). Most of the ingredients that build up a planet are already present at this stage, from the organic molecules that form pre-biotic material to the solids that form the refractory part of the planet (Caselli & Ceccarelli 2012). Despite their importance, the dynamical and chemical evolution of pre-stellar cores is still not well understood. Our partial knowledge of the chemical and physical structure of pre-stellar cores, and how they are fed and influenced by the surrounding environment, limits the level of knowledge that we can achieve at later stages in the star and planet formation process, from protostellar cores to exoplanets.

This paper is part of a systematic effort aiming to provide a representative sample of dynamically evolved cores that can also be used as targets for higher-resolution studies of the inner core structure, dynamics, and chemistry. The whole sample is presented in Caselli et al. (in prep.). We performed a survey of nearby pre-stellar cores with the Atacama Pathfinder EXperiment (APEX) telescope in the rotational lines of $N_2H^+$ and $N_2D^+$ around 300 GHz and 460 GHz. The pre-stellar cores targeted in this survey represent the densest starless cores within 200 pc of the Sun identified in the Herschel Gould Belt Survey (HGBS, Andre et al. 2010). Through these observations, we want to confirm the high central density derived from Herschel data, and to use the $N_2D^+/N_2H^+$ ratio to estimate their evolutionary stage. With this sample we are probing the physical conditions immediately before the protostellar birth, necessary to guide theoretical studies and simulations with the overall aim of understanding the earliest stages of star and planet formation and the timescale of this process. For chemical reasons, $N_2H^+$ and $N_2D^+$ are indicators of an advanced evolutionary stage of a dense core. These molecules form relatively slowly in dense gas from molecular nitrogen, $N_2$, and benefit from the depletion of CO. They too

<footnote n="1">¹ Max-Planck-Institut für Extraterrestrische Physik, Giessenbachstrasse 1, 85748 Garching, Germany</footnote>

<footnote n="2">² European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching, Germany</footnote>

<footnote n="3">³ Department of Physics, PO Box 64, 00014 University of Helsinki, Finland</footnote>

<footnote n="4">⁴ Univ. Rennes, CNRS, IPR (Institut de Physique de Rennes) – UMR 6251, 35000 Rennes, France</footnote>

<footnote n="5">⁵ National Astronomical Observatory of Japan, Osawa 2-21-1, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="6">⁶ Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany</footnote>

<footnote>* Corresponding author; spezzano@mpe.mpg.de</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** APEX observations of several transitions of N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup> towards the dust peak of the pre-stellar core I16293E (black histograms) and results of the LOC non-LTE radiative transfer model (red curve). The physical model of I16293E computed from *Herschel* SPIRE and PACS data is shown in the upper right panel. The black curve (left y-axis) shows the total volume density (n<sub>H</sub><sub>2</sub>), the red curve (right y-axis) shows the dust temperature, and the vertical dashed line indicates the 3000 au (21′′) inner radius of the infalling envelope (see Sect. 3).</figcaption></figure>

deplete in the dense kernel, but the effect is not as severe as for
CO and other molecules. Deuterium fractionation is enhanced
when CO disappears from the gas phase (Roberts et al. 2003).
As shown by Crapsi et al. (2005), a high N<sub>2</sub>D<sup>+</sup>/N<sub>2</sub>H<sup>+</sup> abundance
ratio (≥ 0.1) is a clear signature of a centrally concentrated core
with a very high degree of CO depletion.

This paper focuses on the pre-stellar core IRAS 16293E
(hereafter I16293E), located within the L1689N molecular cloud
in the Rho Ophiuchi cloud complex at a distance of ~141 pc
(Dzib et al. 2018). I16293E shows some of the highest deutera-
tion levels in the interstellar medium; several multiply deuterated
molecules have been observed there, including triply deuter-
ated ammonia, ND<sub>3</sub> (Lis et al. 2016). The environment where
I16293E is embedded is quite active. The low-mass protostellar
system IRAS16923-2422 A and B is located just 90" (12700 au)
to the west of the pre-stellar core (Stark et al. 2004). Fur-
thermore, there is observational evidence that the large-scale
outflows that are driven by the protostellar system have inter-
acted with the pre-stellar core. The area around (and including)
I16293E has been extensively studied with both interferometer
and single-dish facilities (e.g. Lis et al. 2016 and Kahle et al.
2023, and references therein), showing a pattern of velocity that
can be explained by the interaction between the pre-stellar core
and the molecular outflow emanating from the nearby solar-
type protostars. Through Herschel observations of water and
singly deuterated ammonia, Lis et al. (2016) showed that the
outflow wraps around and largely avoids the pre-stellar core.
These authors suggest that the shock associated with the out-
flow has already propagated through the north-eastern part of
the core and what can be observed is quiescent, cold, shock-
compressed, dense gas, that is blue-shifted with respect to the
systemic velocity of the cloud. The broader line profiles in

the western–south-western part of the core suggest instead an
ongoing interaction with the outflow towards the south.

For this paper, we used a non-local thermodynamic equilib-
rium (non-LTE) radiative transfer code to model several rota-
tional lines of N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup>, exhibiting critical densities up to
5×10<sup>6</sup> cm<sup>-3</sup>, thus probing the dense central part of the pre-stellar
core. Additionally, we investigated the rich chemistry towards
I16293E by comparing the line shapes and velocity components
in several molecular emission lines of about 20 different species.

This paper is structured as follows. In Sect. 2 we present the
observations; the non-LTE modelling of the N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup>
lines is described in Sect. 3. The line emission observed in sev-
eral molecules is presented in Sect. 4. We discuss the results in
Sect. 5, and summarise our conclusions in Sect. 6.

2. Observations

2.1. APEX

The spectra observed with APEX towards the dust emission peak
of I16293E are shown in Figures 1 and 2. The molecular emis-
sion data presented in this work were collected with the APEX
telescope as part of the large core sample described in Caselli
et al. (in prep.). Two tunings of the SEPIA345 receiver were used
to cover the N<sub>2</sub>H<sup>+</sup> 3–2 and N<sub>2</sub>D<sup>+</sup> 4–3 transitions at 279.5 GHz
and 308.4 GHz, respectively. These observations were performed
in July and October 2022 (Proposal ID: O-0110.F-9310A-2022).
The N<sub>2</sub>D<sup>+</sup> 3–2 line was covered with the nFLASH230 receiver,
whilst the N<sub>2</sub>D<sup>+</sup> 6–5 and N<sub>2</sub>H<sup>+</sup> 5–4 transitions were covered
simultaneously with one tuning of nFLASH460. These obser-
vations, collected under project M-0110.F-9501C-2022, were
performed in September and November 2022.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Molecular lines observed towards the dust peak of I16293E with APEX. The vertical red dotted line shows the $v_{LSR}$ of the core, 3.7 km s$^{-1}$. Given the complexity of the environment around I16293E, the rest velocity slightly changes depending on the molecule, as already pointed out in Lis et al. (2016) and Kahle et al. (2023). Both lines present in the panel of fully deuterated ammonia (ND$_3$), top right corner, are hyperfine components of the same rotational transition of ND$_3$. The stick spectrum at the bottom of the panel of ND$_3$ shows the hyperfine structure of the transition. Their spectroscopic parameters are reported in Table 2.</figcaption></figure>

All observations were performed in ON-OFF position switching (with the OFF position being +212" in RA and Dec) towards the core central position (RA(J2000) = 16h32m29s.07 and Dec(J2000) = -24°29'09".0), consistent with the continuum peak emission at 970 µm observed with the ALMA compact array (Lis et al. 2016). We used the FFTS backend with a spectral resolution of 64 kHz. This translates into velocity resolution from $\Delta V_{ch} \sim 0.08$ km s$^{-1}$ at 230 GHz to $\sim 0.04$ km s$^{-1}$ at 460 GHz. The calibration accuracy for the APEX observations is $\sim 10\%$.

We reduced the data using the CLASS package of the GILDAS software, and subtracted the baseline using a first-order polynomial. The intensity scale was converted into main beam

temperature $T_{MB}$ assuming the forward efficiency $F_{eff} = 0.95$ and computing the main-beam efficiency at each frequency according to $\eta_{MB} = 1.00797 - 0.000857 * v(\text{GHz})^1$.

The spectra in Figures 1 and 2 show the presence of a very rich chemistry as well as a very complex physical structure. In particular, we can see different velocity components in all bright lines, with the exception of the spectra of $N_2H^+$ and $N_2D^+$ likely because they only trace the pre-stellar core. The relevant information on the targeted transitions of $N_2H^+$ and $N_2D^+$ is given in Table 1.

<footnote n="1">¹ https://www.apex-telescope.org/telescope/efficiency/index.php</footnote>
<table><caption>**Table 1.** Parameters of the targeted N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup> lines.</caption><thead><tr><th>Transition</th><th>Frequency<br>MHz</th><th>A<sub>ul</sub><br>× 10<sup>-3</sup>/s</th><th>g<sub>up</sub></th><th>n<sub>crit</sub><br>cm<sup>-3</sup></th><th>E<sub>up</sub><br>K</th><th>θ<sub>MB</sub><br>''</th><th>η<sub>MB</sub></th><th>Δ V<sub>ch</sub><br>km s<sup>-1</sup></th><th>rms<br>mK</th><th>V<sub>lsr</sub><br>km s<sup>-1</sup></th><th>FWHM<br>km s<sup>-1</sup></th></tr></thead><tbody><tr><td>N<sub>2</sub>H<sup>+</sup> 3–2</td><td>279511.7491</td><td>1.26</td><td>63</td><td>1.4×10<sup>6</sup></td><td>27</td><td>22.3</td><td>0.76</td><td>0.07</td><td>43</td><td>3.610(2)</td><td>0.544(3)</td></tr><tr><td>N<sub>2</sub>H<sup>+</sup> 5–4</td><td>465824.7770</td><td>6.18</td><td>99</td><td>5.5×10<sup>6</sup></td><td>67</td><td>13.3</td><td>0.60</td><td>0.04</td><td>36</td><td>3.600(3)</td><td>0.41(1)</td></tr><tr><td>N<sub>2</sub>D<sup>+</sup> 3–2</td><td>231321.8283</td><td>0.71</td><td>63</td><td>7.4×10<sup>5</sup></td><td>22</td><td>26.0</td><td>0.80</td><td>0.08</td><td>10</td><td>3.630(1)</td><td>0.410(1)</td></tr><tr><td>N<sub>2</sub>D<sup>+</sup> 4–3</td><td>308422.2672</td><td>1.75</td><td>81</td><td>1.5×10<sup>6</sup></td><td>37</td><td>20.2</td><td>0.74</td><td>0.06</td><td>38</td><td>3.600(1)</td><td>0.381(7)</td></tr><tr><td>N<sub>2</sub>D<sup>+</sup> 6–5</td><td>462603.8520</td><td>6.15</td><td>117</td><td>4.4×10<sup>6</sup></td><td>78</td><td>13.5</td><td>0.60</td><td>0.04</td><td>35</td><td>3.620(1)</td><td>0.42(2)</td></tr></tbody></table>

**Notes.** The spectroscopic information refers to the unsplit rotational transitions and is taken from the CDMS catalogue (Müller et al. 2005). The last two columns report the kinematic parameters obtained with the HFS fit. The line critical density is computed with Eq. (4) of Shirley (2015), using the collisional rates of the unsplit lines.

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

**Fig. 3.** Far-infrared maps and averaged intensity profiles of I16293E. Top: Herschel SPIRE HiRes and PACS surface brightness maps at 500 µm, 350 µm, 250 µm, and 160 µm. The latter three are smoothed to the SPIRE HiRes resolution 15". The hatched areas are masked out in the calculation of intensity profiles to avoid a contribution from IRAS16293-2422. Bottom: circularly averaged surface brightness profiles. The error bars show standard deviations within annular regions. The solid curves show Plummer-type functions fitted to the averages.

## 2.2. Herschel

The top row of Fig. 3 shows the thermal dust emission maps of the core at four Herschel wavelengths. The bottom row shows the circularly averaged intensity profiles as a function of the impact parameter. The calibration accuracy for SPIRE and PACS has been estimated to be ~15% (Ladjelate et al. 2020). To perform the non-LTE radiative transfer analysis, we computed a spherically symmetric model of the core based on these intensity profiles. The radial density and dust temperature profiles, $\rho(r)$ and $T_d(r)$, of the core model were determined using the method described in Roy et al. (2014). According to Eq. (3) of Roy et al. (2014), the product $\rho(r) B_\lambda(T_d) \kappa_\lambda$, where $B_\lambda$ is the Planck function at wavelength $\lambda$ and $\kappa_\lambda$ is the dust opacity at this wavelength, can be obtained from the inverse Abel transformation of the surface brightness gradient $dI_\lambda/dr$ as a function of the impact parameter. The density and temperature profiles were fitted by determining the surface brightness gradient at four wavelengths. For the dust opacity, we assumed $\kappa_{250\mu m} = 0.1 \text{ cm}^2 \text{ g}^{-1}$ and a dust opacity index of $\beta = 2.0$ (Hildebrand 1983). We used deconvolved Herschel/SPIRE HiRes images at $\lambda = 250 \mu m$, 350 $\mu m$, and 500 $\mu m$ (1.2, 0.86, and 0.60 THz, respectively), and the Herschel/PACS image of the region at 160 $\mu m$ (1.9 THz). The SPIRE HiRes images at 500 $\mu m$ have a resolution of

approximately 15"<fnref n="2" />. The images at shorter wavelengths were convolved to this resolution. The surface brightness profiles were determined by fitting Plummer-type functions to the concentric circular averages. The purpose of the fitting was to provide a smooth intensity gradient $dI_\lambda/dr$. Figure 3 shows the surface brightness maps of a $10' \times 10'$ region around I16293E. The star-forming core IRAS 16293-2422 is masked out in the images, and this region was not used in the calculation of the intensity profiles. The error bars in the bottom panel represent the standard deviations of the averages, and the solid curves show Plummer-type fits to the average values. The resulting volume density and dust temperature distributions are shown in the top right panel of Fig. 1. The logarithmic fits of the temperature and density profiles beyond the apparent flat radius at 3000 au (21") until 50 000 au (350") are $T \sim r^{0.12}$ and $n \sim r^{-2.45}$. The dust temperature ranges from 10.5 to 15 K, slightly warmer compared to typical pre-stellar cores, where the temperatures range from ~6 to ~13 K. The volume density values go from a few times $10^3 \text{ cm}^{-3}$ in the outer regions to a few times $10^6 \text{ cm}^{-3}$ towards the centre.

<footnote n="2">² http://herschel.esac.esa.int/Docs/SPIRE/spire_handbook.pdf</footnote>
### 3. Non-LTE modelling of N₂H⁺ and N₂D⁺

We used the radiative transfer code LOC (Juvela 2020) to predict the observed N₂H⁺ and N₂D⁺ lines taking into account the overlap of the hyperfine components. As physical models of the source, we used the ρ(r) and T(r) profiles derived from Herschel data, as described in Section 2.2. We assume T<sub>gas</sub> = T<sub>dust</sub> (Sokolov et al. 2017), where T<sub>dust</sub> is derived from the modelling of the Herschel data shown in Figure 3. The output spectral cube is convolved with the corresponding observational beam and compared to the observed spectrum. To model our observations we used the hyperfine-resolved collisional coefficients taken from the EMAA catalogue<fnref n="3" />. The collisional coefficients of N₂H⁺ were computed by Lique et al. (2015), and those of N₂D⁺ by Lin et al. (2020). We tested different abundance profiles for N₂H⁺ and N₂D⁺, and different distributions of the radial velocity and velocity dispersion.

Simulated spectra from the best-fit model are shown in red in Figure 1. Here, the turbulent velocity dispersion is set to σ<sub>turb</sub> = 0.195 km s<sup>-1</sup>, and the radial velocity profile is a step function where the velocity is set to zero in the inner 3000 au (21'') and -0.5 km s<sup>-1</sup> at larger radii (i.e. a static inner core surrounded by infalling envelope). The abundance profiles are also step functions with values of 4.5×10<sup>-10</sup> for N₂H⁺ and 2×10<sup>-10</sup> for N₂D⁺ in the inner 3000 au (21'') and a drop of 70% of the central values at radii larger than 3000 au.

**Velocity profile.** To reproduce the asymmetric line profile of the N₂H⁺ 3–2 line, it is necessary to use a velocity radial profile with a velocity of -0.5 km s⁻¹ at radii larger than 3000 au (21'') and equal to zero in the central part of the core. These values were determined by testing a range of typical values for similar objects. A negative velocity means that the gas is moving towards the centre of the core. Such infalling motions are responsible for red-shifted absorption. While it is possible to reproduce well the high-J transitions with a static core (i.e. v=0 in the whole core), the asymmetry in the N₂H⁺ 3–2 central hyperfine transition can only be reproduced with an infall velocity profile. The results of our test with a static core are shown as a green curve in Figure A.1. We also tried to invert the step function of our best-model velocity profile, assuming our core to have an infalling central region and a static envelope. In this case the velocity was set to -0.5 km s⁻¹ in the central part of the core and equal to zero at radii larger than 3000 au (21''). The results of this test are shown as a blue curve in Figure A.1. In this case, neither the low-J nor the high-J transitions are well reproduced. In particular, the high-J transitions of both N₂H⁺ and N₂D⁺ appear to have separated velocity components because of the infall velocity, which is not the case for the observed spectra.

A radius of 3000 au (21'') for the static core proves to be a good value to reproduce the APEX observations. When decreasing the radius to 1000 au (7''), the intensity of the high-J transitions decreases significantly as the lines become broader because of the effect of the infall. If we instead increase the radius of the static core to 4000 au (28''), the asymmetry in the N₂H⁺ 3–2 line decreases and disappears with the further increase in this value. We note that a step function is obviously a very crude description of the velocity profile of I16293E. Nevertheless, our results show how powerful multi-line non-LTE modelling of N₂H⁺ and N₂D⁺ is when constraining the physical and dynamical structure of a pre-stellar core.

**Abundance profile.** Our N₂H⁺ and N₂D⁺ observations can be well reproduced using a step-function abundance profile with

a drop of 70% at radii larger than 3000 au (21''). The central abundances that best reproduce our observations are 4.5×10<sup>-10</sup> and 2×10<sup>-10</sup> for N₂H⁺ and N₂D⁺, respectively, one of the highest D/H ratios in N₂H⁺ observed to date (0.44). The central regions of the pre-stellar cores are characterised by very high densities (of the order of a few 10<sup>6</sup> cm<sup>-3</sup>) and very low temperatures (down to ~6 K). One of the consequences of such low temperatures and high densities is the freeze-out of molecules on top of the icy mantles of dust grains in the central regions of the cores, and hence a depletion of molecules from the gas phase. Recent ALMA observations of deuterated ammonia in the low-mass pre-stellar core L1544 in the Taurus molecular cloud showed an almost complete freeze-out in the inner 1800 au of the core, where ≥99% of the gas is frozen on the dust grains (Caselli et al. 2022; Pineda et al. 2022). To examine whether the same effect is observable in I16293E, we ran models with an abundance drop in the central part of the core. It turned out that such a drop would dramatically weaken the high-J lines of N₂H⁺ and N₂D⁺, and that models with a constant abundance in the core nucleus reproduce best the observations. Scaling up the fractional abundances, but keeping the centrally decreasing profiles would not help, because the high-J lines would still be too weak and the low-J lines would become optically thick. The results of our tests using different depletion radii are reported in the Appendix. A complete depletion in the inner 2000 au (14'') of the core significantly affects the high-J lines of both molecules (see the blue curve in Figure A.2). The high-J lines are much better reproduced if we consider a complete depletion in the inner 1000 au (7'') of the core (green curve in Figure A.2). The spatial resolution of the present observations is approximately 2000 au (14''), meaning that observations with higher angular resolution are necessary to make a conclusive statement about the freeze-out of N₂H⁺ and N₂D⁺ in the centre of I16293E. We note that the weaker hyperfine components of the 3–2 and 4–3 lines of N₂D⁺ are overestimated by our best model shown in Figure 1. While scaling up the volume density profile and simultaneously decreasing the molecular abundance would help in this regard, we prefer to keep the physical structure of the core derived from the Herschel data fixed. Daniel et al. (2016) derived an abundance profile for N₂H⁺ by modelling the maps of the 1–0 and 3–2 transitions. The abundance they derive at radii larger than 15'' (~2000 au) is 3–4×10<sup>-10</sup>, and it agrees quite well with the abundance derived in this work (3–4.5×10<sup>-10</sup>). Lis et al. (2002) showed that the emission of some simple hydrogenated and deuterated molecules in L1689N (the cloud where I16293E is embedded) comes from distinct physical regions. Bacmann et al. (2016) showed that the ND emission and the NH absorption in I16293E originate from different layers in the cloud. However, the authors point out that this might only be a consequence of the difference in excitation for the two species and that the molecules are both present in the inner dense region of the core, while only NH is present in the outer region. Our work here confirms this, as the emission of N₂H⁺ and N₂D⁺ is very well modelled assuming that the two molecules are co-spatial. As dissociative recombination of N₂H⁺ and N₂D⁺ with electrons is the main pathway to form NH and ND (Roueff et al. 2015), we can safely assume that NH and ND will also be present in the same gas in I16293E in the inner dense region of the core. This is also confirmed by the abundance profiles from the chemical model described in Sipilä et al. (2019), computed at 1×10<sup>5</sup> yr, as shown in Figure A.4.

As discussed for the velocity profile, using a step function for the abundance profile is also a crude approximation of the true distribution of N₂H⁺ and N₂D⁺ across I16293E. Thus, we have also used a different approach and coupled the physical structure

<footnote n="3">³ Available at https://emaa.osug.fr/species-list</footnote>
<table><caption>**Table 2.** Spectroscopic parameters of the observed lines.</caption><thead><tr><th>Molecule</th><th>Transition</th><th>Rest frequency (MHz)</th><th>E<sub>up</sub> (K)</th><th>T<sup>*</sup><sub>A</sub> K</th><th>v<sub>LSR</sub> km s<sup>-1</sup></th><th>FWHM km s<sup>-1</sup></th><th>W K km s<sup>-1</sup></th></tr></thead><tbody><tr><td>C<sup>17</sup>O</td><td>N=2–1 F=7/2–5/2 and 5/2–3/2</td><td>224714.19(8)<sup>(a)</sup></td><td>16</td><td>1.43(1)</td><td>3.787(2)</td><td>1.040(4)</td><td>1.586(6)</td></tr><tr><td>C<sup>17</sup>O</td><td>N=2–1 F=3/2–3/2 and 5/2–3/2</td><td>224715.24(8)<sup>(a)</sup></td><td>16</td><td>0.48(1)</td><td>3.824(1)</td><td>0.85(1)</td><td>0.434(5)</td></tr><tr><td>CN</td><td>N=2–1, J=3/2–1/2, F=1/2–1/2</td><td>226663.693(2)</td><td>16</td><td>0.37(8)</td><td>3.58(3)</td><td>0.58(9)</td><td>0.23(3)</td></tr><tr><td>CS</td><td>N=5–4</td><td>244935.556(3)</td><td>35</td><td>0.81(4)</td><td>3.830(9)</td><td>1.13(2)</td><td>0.97(2)</td></tr><tr><td>CS</td><td>N=6–5</td><td>293912.086(3)</td><td>49</td><td>0.16(2)</td><td>3.89(2)</td><td>1.12(5)</td><td>0.193(8)</td></tr><tr><td>SO</td><td>N J=2 2–1 2</td><td>309502.44(6)</td><td>19</td><td>0.219(3)</td><td>3.68(2)</td><td>0.57(5)</td><td>0.134(8)</td></tr><tr><td>DCO<sup>+</sup></td><td>N=4–3</td><td>288143.858(3)</td><td>35</td><td>3.16(2)</td><td>3.589(1)</td><td>0.573(2)</td><td>1.929(6)</td></tr><tr><td>HDS</td><td>J<sub>K<sub>u</sub>K<sub>c</sub></sub> = 1<sub>01</sub>–0<sub>00</sub></td><td>244555.6(1)</td><td>12</td><td>0.20(1)</td><td>3.723(7)</td><td>0.44(2)</td><td>0.110(3)</td></tr><tr><td>HDO</td><td>J<sub>K<sub>u</sub>K<sub>c</sub></sub> = 1<sub>01</sub>–0<sub>00</sub></td><td>464924.5200</td><td>22</td><td>0.343(3)</td><td>3.530(8)</td><td>0.40(2)</td><td>0.145(6)</td></tr><tr><td>HNC</td><td>J=3–2</td><td>271981.1420<sup>(b)</sup></td><td>26</td><td>2.34(2)</td><td>3.45(1)</td><td>0.51(3)</td><td>1.27(6)</td></tr><tr><td>DNC</td><td>J=3–2</td><td>228910.4810<sup>(b)</sup></td><td>22</td><td>2.103(1)</td><td>3.600(1)</td><td>0.640(2)</td><td>1.433(2)</td></tr><tr><td>DNC</td><td>J=4–3</td><td>305206.22(3)<sup>(b)</sup></td><td>37</td><td>0.69(2)</td><td>3.534(4)</td><td>0.518(9)</td><td>0.380(6)</td></tr><tr><td>o-NH<sub>2</sub>D</td><td>J<sub>K<sub>u</sub>K<sub>c</sub></sub>=2<sub>21</sub>–2<sub>11</sub> v<sub>f</sub> = 1–0, F=1–1, 2–2, and 3–3</td><td>305723.78(3)<sup>(c)</sup></td><td>72</td><td>0.14(2)</td><td>3.63(2)</td><td>0.55(6)</td><td>0.083(6)</td></tr><tr><td>ND<sub>3</sub></td><td>J K F=1 0 1–0 0 1</td><td>309908.46(10)</td><td>14</td><td>0.215(3)</td><td>3.69(1)</td><td>0.5(1)</td><td>0.11(2)</td></tr><tr><td>ND<sub>3</sub></td><td>J K F=1 0 2–0 0 1</td><td>309909.69(10)</td><td>14</td><td>0.439(3)</td><td>3.67(1)</td><td>0.40(3)</td><td>0.19(1)</td></tr><tr><td>ND<sub>3</sub></td><td>J K F=1 0 0–0 0 1</td><td>309911.53(10)</td><td>14</td><td>0.05(2)</td><td>3.68<sup>(d)</sup></td><td>0.45<sup>(d)</sup></td><td>0.02(1)</td></tr><tr><td>H<sub>2</sub>CO</td><td>J<sub>K<sub>u</sub>K<sub>c</sub></sub>=4<sub>04</sub>–3<sub>03</sub></td><td>290623.40(1)</td><td>35</td><td>0.72(2)</td><td>3.623(4)</td><td>0.62(1)</td><td>0.469(7)</td></tr><tr><td>HDCO</td><td>J<sub>K<sub>u</sub>K<sub>c</sub></sub>=4<sub>13</sub>–3<sub>13</sub></td><td>246924.6(1)</td><td>38</td><td>0.43(1)</td><td>3.563(3)</td><td>0.528(8)</td><td>0.241(3)</td></tr><tr><td>D<sub>2</sub>CO</td><td>J<sub>K<sub>u</sub>K<sub>c</sub></sub>=4<sub>04</sub>–3<sub>03</sub></td><td>231410.23(5)</td><td>28</td><td>0.493(1)</td><td>3.703(3)</td><td>0.544(8)</td><td>0.134(8)</td></tr><tr><td>D<sub>2</sub>CO</td><td>J<sub>K<sub>u</sub>K<sub>c</sub></sub>=4<sub>13</sub>–3<sub>12</sub></td><td>245532.75(5)</td><td>35</td><td>0.117(9)</td><td>3.723(7)</td><td>0.44(2)</td><td>0.110(3)</td></tr><tr><td>CH<sub>3</sub>OH</td><td>J<sub>K<sub>u</sub>K<sub>c</sub></sub>=5<sub>15</sub>–4<sub>14</sub></td><td>241767.234(4)</td><td>40</td><td>0.35(4)</td><td>3.69(2)</td><td>0.72(4)</td><td>0.27(1)</td></tr><tr><td>CH<sub>3</sub>OH</td><td>J<sub>K<sub>u</sub>K<sub>c</sub></sub>=5<sub>05</sub>–4<sub>04</sub></td><td>241791.352(4)</td><td>35</td><td>0.44(4)</td><td>3.68(1)</td><td>0.76(3)</td><td>0.36(1)</td></tr><tr><td>CH<sub>3</sub>OH</td><td>J<sub>K<sub>u</sub>K<sub>c</sub></sub>=6<sub>16</sub>–5<sub>15</sub></td><td>290069.747(5)</td><td>54</td><td>0.10(2)</td><td>3.62(3)</td><td>0.70(7)</td><td>0.074(6)</td></tr><tr><td>CH<sub>3</sub>OH</td><td>J<sub>K<sub>u</sub>K<sub>c</sub></sub>=6<sub>06</sub>–5<sub>05</sub></td><td>290110.637(4)</td><td>49</td><td>0.12(2)</td><td>3.68(2)</td><td>0.78(6)</td><td>0.104(7)</td></tr><tr><td>c-C<sub>3</sub>H<sub>2</sub></td><td>J<sub>K<sub>u</sub>K<sub>c</sub></sub>=3<sub>21</sub>–2<sub>12</sub></td><td>244222.151(4)</td><td>18</td><td>0.133(9)</td><td>3.66(1)</td><td>0.77(3)</td><td>0.110(3)</td></tr></tbody></table>

**Notes.** The numbers in brackets denote $1\sigma$ uncertainties in units of the last quoted digit. The frequencies in the table are reported in the CDMS and JPL catalogues (Müller et al. 2005); Pickett et al. (1998) and are derived from the laboratory work in the following: Bocquet et al. (1999); Bechtel et al. (2006); Clark & De Lucia (1976); Helminger & Gordy (1969); Messer et al. (1984); Cazzoli et al. (2002); Klapper et al. (2003); Dixon & Woods (1977); Caselli & Dore (2005); Helminger et al. (1971); Brünken et al. (2006); Saykally et al. (1976); Melosso et al. (2021); Zakharenko et al. (2015); Müller & Lewen (2017); Xu et al. (2008); Bogey et al. (1986)). $^{(a)}$The frequency is averaged between two components, not resolvable in this work, reported in Cazzoli et al. (2002). $^{(b)}$This frequency does not take into account the hyperfine splitting. $^{(c)}$The frequency is averaged between three hyperfine components. $^{(d)}$For the fit of this line, the FWHM and $v_{LSR}$ were fixed to the average values of the other two hyperfine components.

of I16293E derived from the Herschel observations (cf. Section 2) with the chemical model presented in Sipilä et al. (2019). The results of the LOC modelling using the abundance profiles from the chemical models are not satisfactory and in particular cannot reproduce the high-$J$ lines of either $N_2H^+$ or $N_2D^+$ (see Figure A.3). The models in fact predict a severe depletion towards the inner part of the core starting already at radii larger than 5000 au (45"). Because the high-$J$ transitions have high critical densities, the depletion in the central part of the core has a severe effect on the amount of gas that is traced by the $N_2H^+$ 5–4 and $N_2D^+$ 6–5 lines. The results of the modelling using the abundance profiles computed at $1\times10^5$ yr are shown in Figure A.3 in blue. The abundance profiles extracted at later times are more affected by the depletion in the central part of the core, and consequently reproduce the observed high-$J$ lines even less accurately. To decrease the effect of the depletion in the inner part of the core, we tested the effects of a different prescription for the cosmic-ray ionisation rate ($\zeta$), described in Padovani et al. (2018). A higher $\zeta$ could be justified by the proximity to an active site of star formation, as energetic particles can be accelerated along the outflows driven by young stellar objects (Padovani et al. 2016). Additionally, a higher $\zeta$ is expected to keep more $N_2$ in the gas phase in the central part of the core. We used two models with a radial dependence of $\zeta$, one with the

high model ($\mathcal{H}$) and one with the low model ($\mathcal{L}$) from Padovani et al. (2018), as done in Redaelli et al. (2021) to model the emission of $N_2H^+$ in the pre-stellar core L1544. The results of the modelling using these abundances are shown in Figure A.3, in green ($\mathcal{H}$) and red ($\mathcal{L}$). These models did not translate into a big improvement when reproducing the high-$J$ lines because the depletion is still significant at radii larger than 3000 au. Additionally, these models predict $N_2H^+$ and $N_2D^+$ to be abundant in the outer parts of the core, resulting in overly bright (and broad) low-$J$ lines of these molecules. Further studies on the effect of an enhanced cosmic-ray ionisation rate on the abundance of $N_2H^+$ and $N_2D^+$ are necessary to assess the effect of the environment on the depletion radius of these two molecules. More details on our tests with the chemical models are presented in the Appendix.

## 4. Molecular emission towards the dust peak of I16293E

While the main objective of the present observations was to use $N_2H^+$ and $N_2D^+$ lines to confirm the high central density of the core, the large bandwidths of the APEX receivers allowed us to simultaneously detect several additional molecular lines. They
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Zoomed-in image of the blue-shifted and red-shifted velocity components present in some of the observed lines. The vertical red dotted line shows the $v_{LSR}$ of the core, 3.7 km s$^{-1}$.</figcaption></figure>

are listed in Table 2, and are shown in Figure 2. The results of the Gaussian fits (performed with the CLASS software) on all observed lines are also presented in Table 2. We observed simple species, such as diatomics and triatomics ($C^{17}$O, CS, CN, SO, HDO, HDS, DCO$^+$, HNC, and DNC), as well as larger molecules ($H_2CO$, HDCO, $D_2CO$, $CH_3OH$, $c-C_3H_2$, $NH_2D$, and ND$_3$). Among the species that we observed, deuterated molecules are strikingly dominant. Half of the molecules reported in Table 2 are in fact deuterated, and we also observe relatively bright lines of a triply deuterated species, ND$_3$.

Some of the lines shown in Figure 2 have single-Gaussian profiles and do not show hints of multiple velocity components. This is the case for HDO, HDS, $NH_2D$, and ND$_3$. Other lines instead show multiple velocity components. The brightest methanol and CS lines have a quite broad (~1.5 km s$^{-1}$) velocity component blue-shifted from the $v_{LSR}$, between 0.5 and 2 km s$^{-1}$. The lines of CS, CN, $c-C_3H_2$, and $D_2CO$ instead have an additional component that is red-shifted with respect to the core $v_{LSR}$, at ~4.5 km s$^{-1}$. The blue-shifted and red-shifted velocity components are both shown in Figure 4. In order to understand where these velocity components might originate from, it is necessary to inspect the crowded area surrounding the pre-stellar core. Maps of over 30 molecules in the area around I16293E have recently been reported by Kahle et al. (2023), with an exhaustive study of the large-scale molecular emission in the pre-stellar core, the envelope surrounding the IRAS 16293-2422 A and B protostars, and the shocked dense gas in the area. The CO J=6–5 emission map shown in Figure 2 of Kahle et al. (2023) is dominated by the presence of two outflows that originate from the IRAS 16293-2422 A and B protostars and that interact with the pre-stellar core. In particular, the blue side of the east–west outflow broadens around the I16293E pre-stellar core, appearing to wrap around the dense core structure. The red side of the

north-east–south-west outflow instead contributes to two emission peaks of shocked dense gas in the north of I16293E, E1 and E2. Given their velocity and vicinity to the core, we suggest that the red- and blue-shifted components observed in some of the emission lines towards the dust peak of I16293E trace gas from outflows that generated the E1 and HE2 emission peaks. This gas is not part of the pre-stellar core, but it is in the line of sight.

**5. Discussion: A dense pre-stellar core in an active environment**

The non-LTE radiative transfer modelling of the $N_2H^+$ and $N_2D^+$ lines showed how sensitive these lines are to small changes in the velocity and abundance profiles, and have proved to be excellent diagnostics for the gas they trace. To reproduce the asymmetries observed in the low-$J$ transitions, the gas needed to be infalling. The best velocity profile to reproduce all the observed lines of $N_2H^+$ and $N_2D^+$ is a static core in the inner 3000 au (21") surrounded by an infalling gas layer. The velocity difference between the inner static core and the envelope in our model (0.5 km s$^{-1}$) is consistent with the ND and NH observations in Bacmann et al. (2016), where NH absorption is red-shifted with respect to ND emission by ~0.4 km s$^{-1}$.

The $N_2H^+$ and $N_2D^+$ lines proved to be very sensitive to depletion, in particular the high-$J$ transitions. As discussed in Section 3, an abundance profile with no depletion in the centre reproduces our observations best. We also get good results with a depletion radius of 1000 au (7"), which is however a spatial scale that we cannot resolve with our Herschel observations, and hence in the physical structure that we use for the modelling. Our results using a depletion radius of 1000 au (7") are in agreement with the very low abundance of $N_2D^+$ reported
in Lis et al. (2016) with ALMA observations within a 1300 au (9′′) beam. Testing the abundance profiles predicted by chemical models using the physical structure of I16293E gave us interesting insights into how to further develop the chemical models. The chemical model that we used (Sipilä et al. 2019) predicts too much depletion, and consequently cannot reproduce the high-J lines of N₂H⁺ and N₂D⁺. We also used two models with a radial dependence of the cosmic-ray ionisation rate, following the work of Padovani et al. (2018), and the one with the high cosmic-ray ionisation rate (H) does a good job reproducing the high-J lines, but not the low-J lines. A higher cosmic-ray ionisation rate keeps more molecular nitrogen, N₂, in the gas phase in the central part of the core, but it also predicts high abundances of N₂H⁺ and N₂D⁺ in the outer parts of the core (see Figure A.4). The presence of a considerable amount of N₂H⁺ and N₂D⁺ in the outer gas layers, where the gas temperature is 13–15 K, translates into an additional broader line component in the N₂H⁺ and N₂D⁺ 3–2 lines (quite noticeable in the green and red curves in the left panels in Figure A.3), which is not consistent with the observed spectra.

To better reproduce the observed lines in a pre-stellar core such as I16293E, we plan to include the effects of the shock in the chemical models. A realistic model of I16293E should also combine chemistry with the dynamic evolution of the core. In particular, I16293E has likely evolved quickly into its present high-density state (e.g. due to shock compression), which is not captured by the present model, and hence the depletion is overestimated.

The N₂D⁺/N₂H⁺ ratio that best reproduces our observations is 0.44, one of the highest observed to date in pre-stellar cores. Using the LTE approach, Kahle et al. (2023) derived a ratio of ~10% (Table F.2 in their paper). The large discrepancy between the two values strongly suggests that non-LTE effects have quite an impact on these spectra. I16293E is well known to have some of the largest deuteration fractions among pre-stellar cores (Lis et al. 2016). How much of the observed deuteration is a consequence of the environment where the core is embedded still needs to be explored in a quantitative manner. It is important to note that half of the molecules that we observed in our spectra are deuterated isotopologues, confirming the high level of overall deuteration towards this source. A quantitative study of the deuteration of different classes of molecules will be relevant, but it is beyond the scope of this paper. Our results on N₂H⁺ show that the complex morphology of the source warrants synergetic approaches that include observations, radiative transfer modelling, and chemical modelling, when possible. The presence of line components that are red-shifted and/or blue-shifted from the rest velocity of the core is a consequence of the complex morphology of the source as discussed in Section 4. We note that the N₂D⁺ and N₂H⁺ lines are the only bright lines in our sample that show no additional velocity components, thus proving to be excellent tracers of the dense core, also in the case of a crowded environment similar to the one where I16293E is embedded.

In a recent paper, Pagani et al. (2024) present para-D₂H⁺ and ortho-H₂D⁺ maps of the IRAS16293E core. The emission of both molecular ions peak approximately 25'' north-east from the dust emission peak (our centre position). Because H₂D⁺ and D₂H⁺ are predicted to be most abundant in cold, dense, and heavily depleted gas, Pagani et al. suggest that the true centre of the core lies at the D₂H⁺ peak, and that the dust emission peak, which roughly coincides with the N₂H⁺(4–3) and ND₃(1₀ – 0₀) emission peaks as observed with ALMA (Lis et al. 2016) indicates a hotspot resulting from compression by outflows from IRAS 16293-2422. In Fig. B.1 we show the dust colour

temperature (T~C~) and H~2~ column density (N(H~2~)) maps of the region, derived from Herschel PACS and SPIRE/HiRes images at the wavelengths ~λ = 70, 100, 160, 250, 350, and 500 µm. The maps are obtained by fitting a grey-body function pixel-by-pixel to the aligned surface brightness maps smoothed to a common resolution, with the T~C~ and N(H~2~) as free parameters. The T~C~ represents the line-of-sight average of the dust temperature. Both the dust emission peak and the D₂H⁺ peak, indicated with the diamond in the maps, are found in a valley of relatively low colour temperatures (T~C~ < 13 K), where the H₂ column densities exceed 10²³ cm⁻². The N(H~2~) at the D₂H⁺ peak is 1.5 × 10²³ cm⁻², which agrees with the core model of Pagani et al. (2024). The H₂ column density reaches its maximum at the centre position adopted here. This is no surprise, because the position was chosen from Herschel HiRes column density maps provided by the HGBS archive. The N(H~2~) and T~C~ distributions derived from Herschel reproduce very well the 850 µm surface brightness map of the core observed with SCUBA (Pattle et al. 2015), which also peaks at our (0,0). This position probably represents a true maximum in the total column density. Usually the highest densities are found in the direction of N(H~2~) maxima, and here it probably also contains a local density maximum. The present line observations, in particular the strong detection of the N₂H⁺(5–4) and N₂D⁺(6–5) transitions confirm that the density in this position exceeds 10⁶ cm⁻³, which agrees with the estimate from the Abel transformation method described in Sect. 2.2. We note that according to the 1D physical model used by Pagani et al. (2024), the volume density at the offset 25'' from the D₂H⁺ peak is ~4 × 10⁵ cm⁻³. This density would be too low to produce the observed intensities of N₂H⁺(5–4) and N₂D⁺(6–5) in this position. However, according to our tests, another core resembling the inner portions (~r < 10¹⁷ cm) of the model of Pagani et al. (2024) (see their Fig. 4) with a central temperature of 8 K can be placed in the direction of the D₂H⁺ peak without substantial changes in the observable 850 µm map.

## 6. Summary and conclusions

Several rotational transitions of N₂H⁺ and N₂D⁺ have been observed towards the dust peak of the pre-stellar core I16293E using APEX, namely the N₂H⁺ 3–2 and 5–4 and the N₂D⁺ 3–2, 4–3, and 6–5 transitions. I16293E is a dense pre-stellar core located close to the well-studied low-mass binary protostar IRAS 16293–2422. The physical structure of I16293E, constrained from Herschel continuum maps, shows very high central density (~3×10⁶ cm⁻³) and temperatures ranging from 11 to 15 K. The best model to reproduce our spectra of N₂H⁺ and N₂D⁺ hints at the presence of a static central core (~3000 au, 21′′) surrounded by an infalling gas layer. Furthermore, the high-J N₂H⁺ and N₂D⁺ lines that we observe are very sensitive to depletion because of their high critical densities (up to ~5×10⁶ cm⁻³), and our best model does not include any depletion towards the centre of the core. Because I16293E is embedded in a very active environment, it might undergo a local enhancement of cosmic rays that could partially counteract high densities and low temperatures that cause molecular depletion towards the centre of pre-stellar cores (Ceccarelli et al. 2014; Padovani et al. 2016). Additionally, because of the possible passage of a shock (Lis et al. 2002, 2016), the high-density nucleus in I16293E may have formed quickly, leaving not enough time for substantial freeze-out. Our tests cannot completely exclude depletion, especially at radii <1000 au (7"). High angular resolution observations to derive the physical structure of the inner part of the pre-stellar
core will be necessary in this regard. We could not reproduce the $N_2H^+$ and $N_2D^+$ lines using abundance profiles derived from chemical models. Nevertheless, the insights gained from these tests will be used to further constrain our models, especially concerning the radial dependence of the cosmic-ray ionisation rate.

I16293E exhibits some of the highest levels of deuterium fractionation observed in pre-stellar cores. This includes the detection of several multiply deuterated molecules. From the abundance profiles of our best models for $N_2H^+$ and $N_2D^+$, we derive a deuterium fraction of 0.44, among the largest values measured in similar objects.

The wide frequency coverage of the APEX data allowed us to detect the emission of 16 additional molecules, half of which are deuterated isotopologues. Given the complex morphology of the area and the limited angular resolution of our observations (13"–26"), we used the line profiles to disentangle different gas components in the line of sight. Some lines show multiple velocity components, and we suggest that some shocked gas from the outflows originating from the IRAS 16293-2422 A and B protostars contribute to the emission of molecules like CS and CH$_3$OH. This is very clearly seen in the brightest lines shown in Figure 2. There are, however, some lines that show a single-Gaussian profile and no hints of multiple velocity components. This is the case for HDO, HDS, NH$_2$D, and ND$_3$, as well as $N_2H^+$ and $N_2D^+$.

High spectral resolution observations coupled with radiative transfer modelling and supported by state-of-the-art chemical models are mandatory in order to dissect the contribution of different gas components towards I16293E, and finally to determine the chemical ingredients that are available for the future planetary systems that will form in the area. Future studies coupling observations and modelling should in particular explore how the enhanced deuteration affects different classes of molecules, with a particular focus on complex organic molecules.

**Acknowledgements.** The authors wish to thank the anonymous referee for the careful review of the manuscript. S.S., E.R., P.C. O.S., J.E.P. gratefully acknowledge the support of the Max Planck Society. The data was collected under the Atacama Pathfinder EXperiment (APEX) Project, led by the Max Planck Institute for Radio Astronomy at the ESO La Silla Paranal Observatory. This research has made use of data from the Herschel Gould Belt survey (HGBS) project (http://gouldbelt-herschel.cea.fr). The HGBS is a Herschel Key Programme jointly carried out by SPIRE Specialist Astronomy Group 3 (SAG 3), scientists of several institutes in the PACS Consortium (CEA Saclay, INAF-IFSI Rome and INAF-Arcetri, KU Leuven, MPIA Heidelberg), and scientists of the Herschel Science Center (HSC). This research has made use of spectroscopic and collisional data from the EMAA database (https://emaa.osug.fr and https://dx.doi.org/10.17178/EMAA). EMAA is supported by the Observatoire des Sciences de l'Univers de Grenoble (OSUG).

## References

Bacmann, A., Daniel, F., Caselli, P., et al. 2016, A&A, 587, A26

Bechtel, H. A., Steeves, A. H., & Field, R. W. 2006, ApJ, 649, L53

Brücken, S., Müller, H. S. P., Thorwirth, S., et al. 2006, J. Mol. Struct., 780, 3

Bocquet, R., Demaison, J., Cosléou, J., et al., 1999, J. Mol. Spectrosc., 195, 345

Bogey, M., Demuynck, C., & Destombes, J. L. 1986, Chem. Phys. Lett., 125, 383

Caselli, P., & Ceccarelli, C. 2012, A&A Rev., 20, 56

Caselli, P., & Dore, L. 2005, A&A, 433, 1145

Caselli, P., Pineda, J. E., Sipilä, O., et al. 2022, ApJ, 929, 13

Ceccarelli, C., Dominik, C., López-Sepulcre, A., et al. 2014, ApJ, 790, L1

Cazzoli, G., Dore, L., Puzzarini, C., et al. 2002, Phys. Chem. Chem. Phys. (Incorp. Faraday Trans.), 4, 35755

Clark, W. W., & De Lucia, F. C. 1976, J. Mol. Spectrosc., 60, 332

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

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

Daniel, F., Faure, A., Pagani, L., et al. 2016, A&A, 592, A45

Dixon, T. A., & Woods, R. C. 1977, J. Chem. Phys., 67, 3956

Dzib, S. A., Ortiz-León, G. N., Hernández-Gómez, A., et al. 2018, A&A, 614, A20

Helminger, P., & Gordy, W. 1969, Phys. Rev., 188, 100

Helminger, P., Cook, R. L., & De Lucia, F. C. 1971, J. Mol. Spectrosc., 40, 125

Herbst, E., & van Dishoeck, E. F. 2009, ARA&A, 47, 427

Hildebrand, R. H. 1983, QJRAS, 24, 267

Jørgensen, J. K., Favre, C., Bisschop, S. E., et al. 2012, ApJ, 757, L4

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

Kahle, K. A., Hernández-Gómez, A., Wyrowski, F., et al. 2023, A&A, 673, A143

Klapper, G., Surin, L., Lewen, F., et al. 2003, ApJ, 582, 262

Ladjeolate, B., André, P., Könyves, V., et al. 2020, A&A, 638, A74

Lin, S.-J., Pagani, L., Lai, S.-P., et al. 2020, A&A, 635, A188

Lique, F., Daniel, F., Pagani, L., et al. 2015, MNRAS, 446, 1245

Lis, D. C., Gerin, M., Phillips, T. G., et al. 2002, ApJ, 569, 322

Lis, D. C., Wootten, H. A., Gerin, M., et al. 2016, ApJ, 827, 133

Melosso, M., Bizzocchi, L., Dore, L., et al. 2021, J. Mol. Spectrosc., 377, 111431

Messer, J. K., De Lucia, F. C., & Helminger, P. 1984, J. Mol. Spectrosc., 105, 139

Müller, H. S. P., & Lewen, F. 2017, J. Mol. Spectrosc., 331, 28

Müller, H. S. P., Schlöder, F., Stutzki, J., et al. 2005, J. Mol. Struct., 742, 215

Ossenkopf, V., & Henning, Th. 1994, A&A, 291, 943

Padovani, M., Marcowith, A., Hennebelle, P., et al. 2016, A&A, 590, A8

Padovani, M., Ivlev, A. V., Galli, D., et al. 2018, A&A, 614, A111

Pagani, L., Belloche, A., & Parise, B. 2024, A&A, 691, A88

Pattle, K., Ward-Thompson, D. and Kirk, J. M., et al. 2015, MNRAS, 450, 1049

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

Pineda, J. E., Harju, J., Caselli, P., et al. 2022, AJ, 163, 294

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

Roberts, H., Herbst, E., & Millar, T. J. 2003, ApJ, 591, L41

Roy, A., André, P., Palmeirim, P., et al. 2014, A&A, 562, A138

Roueff, E., Loison, J. C., & Hickson, K. M. 2015, A&A, 576, A99

Saykally, R. J., Szanto, P. G., Anderson, T. G., et al. 1976, ApJ, 204, L143

Shirley, Y. L. 2015, PASP, 127, 299

Sipilä, O., Caselli, P., & Harju, J. 2019, A&A, 631, A63

Sokolov, V., Wang, K., Pineda, J. E., et al. 2017, A&A, 606, A133

Stark, R., Sandell, G., Beck, S. C., et al. 2004, ApJ, 608, 341

Xu, L.-H., Fisher, J., Lees, R. M., et al. 2008, J. Mol. Spectrosc., 251, 305

Zakharenko, O., Motiyenko, R. A., Margulès, L., et al. 2015, J. Mol. Spectrosc., 317, 41
<figure><img src="image_6.png" /><figcaption>**Fig. A.1.** Results of tests for the velocity radial profile. The green curve shows the results of the non-LTE modelling done assuming a static core (i.e. v=0 in the whole core). The blue dashed curve shows the results of the non-LTE modelling done assuming a static envelope (the velocity was set to -0.5 km s⁻¹ in the central part of the core and equal to zero at radii larger than 3000 au).</figcaption></figure>

## Appendix A: Velocity profile and abundance profile tests

Different velocity and abundance profiles have been tested to reproduce the N₂H⁺ and N₂D⁺ lines using the physical structure derived from Herschel observations (see Section 2.2).

The velocity radial profile that reproduces best our observations is represented by a static central core (v=0 at r< 3000 au) and an infalling outer gas layer (v=-0.5 km s⁻¹ at r> 3000 au). Our tests using different velocity radial profiles are shown in Figure A.1. The asymmetry in the N₂H⁺ 3-2 line cannot be reproduced with a static core (green curve in Figure A.1). Most of the lines, but particularly the high-J transitions of both N₂H⁺ and N₂D⁺, cannot be reproduced if we assume a static envelope and an infalling central part of the core (v=-0.5 km s⁻¹ at r< 3000 au, see the blue dashed curve in Figure A.1).

The abundance profile that reproduces best our observations is a step function with values of 4.5×10⁻¹⁰ for N₂H⁺ and 2×10⁻¹⁰ for N₂D⁺ in the inner 3000 au and a drop of 30% of the central values at radii larger than 3000 au. It is however expected for molecules to freeze-out in the central regions of pre-stellar cores because of the high densities and low temperatures (Caselli et al. 2022). Our tests assuming a complete depletion in the central 2000 and 1000 au are shown in Figure A.2. The high-J transitions of N₂H⁺ and N₂D⁺ are poorly reproduced if assuming a complete depletion in the central 2000 au (dashed blue curve in Figure A.2). It is important to note that these transitions have a high critical density (~10⁷ cm⁻³) and they are presumably emitting only in the dense gas at the centre of the pre-stellar core. A significant improvement can be seen if we assume a smaller radius for the complete depletion, 1000 au (green curve in Figure A.2), although the high-J transitions are still better reproduced without depletion (red curve in Figure 1). However, given the angular resolution of our Herschel observations (~2700 au), we would need a better resolved physical structure to make a more conclusive statement on the complete freeze-out radius for N₂H⁺ and N₂D⁺ in I16293E.

Figure A.3 shows the results of radiative transfer models using the radial abundance profiles from the chemical models described in Sipilä et al. (2019). We furthermore tested the effects of a radial dependence of ζ following the approach described in Padovani et al. (2018) and used in Redaelli et al. (2021). The abundance profiles that have been used are shown in Figure A.4.

## Appendix B: Dust temperature and H₂ column density maps of I16293E

The dust colour temperature (T_c) and H₂ column density (N(H₂)) maps of the I16293E core and the adjacent protostellar core IRAS 16293-2422 are shown in Fig. B.1. These are derived from Herschel/PACS and SPIRE/HiRes images at the wavelengths λ = 70, 100, 160, 250, 350, and 500 µm by fitting a grey body function to aligned surface brightness maps. Colour corrections recommended
<figure><img src="image_7.png" /><figcaption>**Fig. A.2.** Results of tests on the effects of depletion. The blue dashed curve shows the results of the non-LTE modelling done assuming a complete depletion in the inner 2000 au of the core, while the green curve assumes a complete depletion in the inner 1000 au of the core.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. A.3.** Results of tests using abundance profiles from chemical models. The blue curve shows the results of the non-LTE modelling using the abundance profiles of N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup> predicted from the chemical model presented in Sipilä et al. (2019). The green dashed and red dotted curves instead use a chemical model with a radial dependence of the cosmic-ray ionisation rate, with the high model (H) and with the low model (L) from Padovani et al. (2018), respectively.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. A.4.** Abundance profiles from chemical models. Left and central panels: Abundance profiles from different models for N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup>. Right panel: Abundance profiles for ND and NH from Sipilä et al. (2019) at 10<sup>5</sup>yr.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. B.1.** Dust colour temperature ($T_c$) and H<sub>2</sub> column density ($N(H_2)$) maps of the IRAS 16293 region derived using Herschel PACS images at 70, 100, and 160 µm, and deconvolved SPIRE HiRes images at 250, 350, and 500 µm. The position observed in the present work is indicated with a plus sign, and the D<sub>2</sub>H<sup>+</sup> peak discovered by Pagani et al. (2024) is indicated with a diamond. The protostellar core IRAS 16293-2422 dominates the right side of both images. The positions of the protostars IRAS16293-2422 A and B are shown with asterisks. The column density at the (0,0) is 2.8 × 10<sup>23</sup> cm<sup>-2</sup> and the T<sub>c</sub> there is 12.3 K. The effective beam of the SPIRE HiRes map at 500 µm is indicated in the bottom left corner of the image.</figcaption></figure>

in Herschel documentation<fnref n="4" /> were applied to the observed surface brightness values in the iterative solution, needed to take into account the dependence of the effective bandwidth of the instrument on the temperature of the source. We assumed that the dust opacity at 250 µm is $κ_{250}$ µm = 0.1 cm² g⁻¹, and that the dust emissivity index is $β = 2$. This assumption approximately corresponds to the dust opacity model of Ossenkopf & Henning (1994) for unprocessed dust particles with thin ice coatings. The positions of the dust emission peak in I16293E, the D<sub>2</sub>H<sup>+</sup> peak discovered by Pagani et al. (2024), and the protostars IRAS 16293-2422 A and B are indicated in the maps.

<footnote n="4">⁴ https://www.cosmos.esa.int/web/herschel/legacy-documentation</footnote>