© The Authors 2023

1 Introduction

Studying the physical and chemical processes involved in the early stages of low-mass star formation is crucial to understand the conditions that gave birth to the formation of stars such as our Sun. In particular, the analysis of the molecular inventory in low-mass protostars can shed light on the kinematics and the chemistry of these young stages.

One of the most studied sources that has been at the center of a plethora of molecular line studies is the archetypical Class 0 so-called protostar¹ IRAS 16293-2422, first discussed by Walker et al. (1986). This far-infrared source is a dust condensation located in the Lynds 1689 N dark cloud (L1689N) which is part of the Ophiuchus dark cloud and star-forming complex. It harbors a hierarchical system of low-mass young stellar objects (YSOs). A trigonometric parallax measured with very long baselines interferometry (VLBI) of the H₂O masers associated with one of the compact sources yielded a distance of ${141}_{-21}^{+30}$ pc(Dzib et al. 2018).

The composite source consists of two main condensations (Wootten 1989; Mundy et al. 1992) called IRAS 16293-2422 A and IRAS 16293−2422 B (sources A and B from now on) separated by 5" (~705 au). Source A has been found to be a double system itself (Wootten 1989; Maureira et al. 2020) which consists of the protostars A1 and A2 separated by 0.″3 (~42 au). On the other hand, the nature of source B still is not yet fully constrained, although several studies suggest that it is in a very early stage of star formation (Pineda et al. 2012; Hernández-Gémez et al. 2019a). The IRAS 16293-2422 A/B system is surrounded by an extended envelope with about a 6000−8000 au radius (Schöier et al. 2002; Crimier et al. 2010; Jacobsen et al. 2018) and has a mass of about 4–6 M_⊙ (Jacobsen et al. 2018; Ladjelate et al. 2020). 16293E, a very well-known prestellar core also embedded in L1689N, is located 1.′5 (~12700 au) east (see Fig. 1) of this extended envelope (Stark et al. 2004).

A multilobe molecular outflow system with scales up to 0.2 pc can be observed toward IRAS 16293−2422 (Walker et al. 1988; Mizuno et al. 1990). Two bipolar outflows have been reported to arise from source A (see Fig. 2), one of them extending in the east-west (E–W) direction (PA = 110°) and the other in the northeast-southwest (NE-SW) direction (PA = 60°; Castets et al. 2001; Hirano et al. 2001; Yeh et al. 2008). There is a third outflow with an extent of 0.01 pc in the northwest-southeast (NW-SE) direction (PA= 145°) that seems to be arising from source A and interacting with material around source B (Girart et al. 2014). While some authors suggest that source B could be indeed producing one of these outflows (Loinard et al. 2013; Oya et al. 2018), there is still debate on their true origin. Interestingly, one of these outflows arising from IRAS 16293−2422 A/B seems to be interacting also with 16293E, modifying the chemical abundances in this source (e.g., Lis et al. 2016).

The dust emission and the molecular content of both sources, A and B, have been studied in great detail with high angular resolution (up to ≈0.″1) with the Very Large Array (VLA; e.g., Chandler et al. 2005; Hernández-Gómez et al. 2019a), the SubMillimeter Array (SMA; Jørgensen et al. 2011) and the Atacama Large Millimeter/submillimeter Array (ALMA; Pineda et al. 2012; Jørgensen et al. 2016). In particular, the imaging spectral line surveys toward A and B (Jørgensen et al. 2011, 2016) with the SMA and ALMA, most recently the ALMA protostellar interferometric line survey (PILS; at 0.″5 resolution), reveal an immense number of lines from numerous species, including "complex organic molecules" toward both A and B (Manigand et al. 2020).

On larger scales, the molecules in the envelope have been studied in the course of the IRAS 16293−2422 millimeter and submillimeter spectral survey (TIMASSS; Caux et al. 2011) using pointed single-dish observations with the IRAM 30 meter telescope and the 15 meter James Clerk Maxwell Telescope (JCMT) at coarser angular resolution (10"–30"). Nevertheless, our knowledge of the spatial distribution of many molecules in the envelope and its surroundings is still quite limited as published studies deal only with single pointing observations (usually some centroid position between A1, A2 or B within a few arc second variation). In contrast, the few larger scale studies on the molecular environment further from the sources mainly focus on a small, selected number of lines (e.g., Castets et al. 2001; Lis et al. 2002). Indeed, a comprehensive systematic broad-band study of the more extended molecular environment around IRAS 16293−2422 A/B and E is still missing in the literature.

In this work, we present a study of the molecular environment of IRAS 16293−2422 A/B and E based on single-dish APEX observations. We have imaged a 3.′5 × 3'.5 sized region (0.12 × 0.12 pc²) centered between IRAS 16293−2422 and 16293E in a large number of molecular lines covering a significant part of the 870 μm atmospheric window.

The paper is organized as follows: In Sect. 2 we describe the details of the observations and the data reduction. In Sect. 3 we present the particulars of the identified molecules and show the maps with the molecular distribution in IRAS 16293−2422 A/B and 16293E. In Sect. 4, we present a detailed analysis of the physical properties of both sources and nearby related emission peaks based on temperature and column and volume density determinations from radiative transfer models. In Sect. 5, we discuss the consequences of our findings and compare the molecular line morphologies with dust continuum maps. Finally in Sect. 6, we summarize our results. Extensive appendices give detailed information on the observed line detections and derived quantities and present maps of the line emission.

2 Observations and data reduction

The observations were carried out with the Atacama Patfinder EXperiment 12 meter diameter submillimeter telescope (APEX) located on the Llano de Chajnantor in the Chilean High Andes at an altitude of 5107 m (Güsten et al. 2006).

2.1 Observations with LAsMA and FLASH+

The data were taken under project M-0102.F-9519A-2018 (P.I. Karl M. Menten) during several runs between 2019 July 19–28 with the Large APEX subMillimeter Array (LAsMA) multi-beam heterodyne receiver (Güsten et al. 2008) in the frequency range between 277 and 375 GHz. The LAsMA receiver is a 7 pixel multibeam sideband separating (2SB) receiver with the lower and upper sideband (LSB and USB, respectively), each 4 GHz wide and the centers of the bands separated by 12 GHz. Inspecting the data, we found calibration problems with the sixth LAsMA pixel. As consequence, we did not use data taken with this pixel.

The observations cover a region of 3.′5 × 3.′5 centered on a position offset (δ = 50'', α = −10'') from IRAS 16293 A/B at RA(J2000) = 16^h32^m22.78^s, Dec(J2000) = −24°28'38.7'', in order to include both cores. They were performed in the on-the-fly (OTF) mode and a total frequency coverage of 45.6 GHz has been recorded.

In addition, single pointing observations with the First-Light APEX Submillimeter Heterodyne instrument (FLASH+; Klein et al. 2014) were made toward IRAS 16293 A/B and were used only for line-identification purposes. These observations in total cover 17 GHz bandwidth, of which 9 GHz allow us to identify lines in frequency intervals between 476 and 493 GHz, while additional 2.3 GHz extend the frequency intervals covered by LAsMA between 276.5 GHz and 292.9 GHz. An overview of the observed frequency-bands and the observing dates is given in Appendix A. Both receivers were connected to fast Fourier transform spectrometers (FFTS; e.g., Klein et al. 2012). The FFTS modules provide a total of 2¹⁶ frequency channels per 4 GHz wide sideband, corresponding to a channel spacing (in frequency) of 61.04 kHz. This corresponds to 0.07 and 0.05 km s⁻¹ at our lowest and highest observed frequencies of 277 and 37 GHz, respectively. These numbers all are smaller than typical line widths observed (a few km s⁻¹).

During the observations, the LAsMA system temperatures varied between 150 and 250 K at frequencies below 310 GHz, between 200 and 400 K from 330 to 370 GHz, and between 400 and 600 K above 370 GHz. The conversion between antenna temperature $T_A^{*}$ and the main-beam brightness temperature T_MB is given by $T_{\text{MB}}=T_A^{*}\times {\eta }_{\text{FW}}\times {\eta }_{\text{MB}}$, where η_MB is the main-beam efficiency and η_FW the forward coupling efficiency. Based on Jupiter continuum pointings during the observations of the project, we estimate an average conversion factor of η_FW/η_MB = 1 / 0. 7 for the full array. The APEX half power beam width (HPBW) θ_B at frequency ν (GHz), in arcseconds is given by θ_B[″] = 7.″8 × (800/v[GHz]; Güsten et al. 2006) and thus varies at the observed frequencies between 17" and 23" (corresponding to 0.012 pc and 0.015 pc) respectively.

2.2 Observations with SEPIA660

Additional data were taken under project M-0102.F-9524C-2018 with the Swedish-ESO PI Instrument for APEX (SEPIA) at a tuning frequency of 691.473 GHz. The SEPIA660 receiver is one of three ALMA-resembling receiver cartridges installed at APEX and operates as single pixel heterodyne 2SB receiver with two polarizations (Belitsky et al. 2018). Upper and lower sideband each cover 8 GHz and the centers of the two bands are separated by 16 GHz.

The OTF observations with SEPIA660 cover a region of 200'' × 200'' centered on (50'', −10'') from IRAS 16293 A/B with a total integration time of 30 min. In addition, a larger region of 460" × 280" has been observed in the OTF mode with a total integration time of 20 min. System temperatures of the receiver hereby varied between 800 and 1200 K.

The SEPIA660 receiver was connected to 8 FFTS backends that each provide 2¹⁶ channels with a channel spacing of 61 kHz over 4 GHz of bandwidths. In this work, we present OTF maps of the CO (J = 6−5) transition at 691.473 GHz for which the velocity spacing is 0.03 km s⁻¹, much smaller that the observed line width. The APEX HPBW at this frequency is 9" which results in an increased angular resolution as compared with the LAsMA observations.

2.3 Data reduction

For the data reduction, we used the CLASS program that is part of the GILDAS² software package developed by the Institut de Radioastronomie Millimétrique (IRAM). For the mapping of each spectral line, the local standard of rest (LSR) velocity-scale was first modified to correspond to the frequency of the considered transition. In order to improve the signal-to-noise ratio in the spectra, each two adjacent channels were averaged, resulting in a velocity-resolution of about 0.1 km s⁻¹ per channel, which is adequate for resolving the observed spectral lines.

A linear spectral baseline was then calculated for the line surroundings, usually considering channels covering velocities below −10 km s⁻¹ and above 15 km s⁻¹, such that the line emission around the systematic cloud velocity of 4 km s⁻¹ is properly excluded. After the baseline-subtraction, the resulting spectrum was extracted in a velocity-range between −25 km s⁻¹ and 30 km s⁻¹. The windows used for baseline-calculation and extraction were adjusted for the processing of broader spectral lines such as the CO (J = 3−2) transition.

2.4 Mapping

The maps of the IRAS 16293−2422 A/B+E area were created by extracting all the observed lines of a given transition from the APEX data and combining and gridding them using the CLASS-module xy_map. Since the edges of the resulting maps are not covered by all pixels, they show a higher noise and are therefore not considered during the analysis. In order to determine the noise of a particular map, an integrated-intensity map was created from nearby line-free channels.

Broad line wings indicate outflows, which are visualized by creating integrated intensity maps from channels that correspond to the line wing emission. Contours of these maps were drawn in a way that regions with signal to noise ratio higher than 3 (3σ) are highlighted, for which the mean noise over the whole map was considered. The mapping of cloud cores and outflows was performed for all transitions present in the LAsMA data.

3 Results

To identify the detected molecular lines, we have used the Jet Propulsion Laboratory (JPL) line catalog³ (Pickett et al. 1998) and the Cologne Database for Molecular Spectroscopy (CDMS)⁴ (Müller et al. 2001, 2005; Endres et al. 2016). The line identification was carried out at the positions of IRAS 16293 A/B and E, since occurring transitions are expected to be brightest in the immediate surroundings of the cloud cores. While we included frequency windows observed with FLASH+ for the line identification at IRAS 16293 A/B, we note that these data were obtained as single pointings and therefore the corresponding transitions could not be mapped.

In our observations, we identify a total of 144 spectral lines from 36 different species, which subdivide into 20 distinct molecules and 16 of their isotopologues. 33 of these species show a total of 126 transitions in the frequency range covered by LAsMA, which was used to map the larger-scale emission around the YSOs. A full list with all the observed transitions can be found in Appendix C. Maps of the emission distribution of each detected transition were created and are shown in Appendix J. We show only the brightest lines of species with many transitions, since the fainter transitions do not provide us with any additional information about the spatial distributions of molecules.

3.1 Cloud cores in L1689N

In general, the molecular line maps reveal a complex morphology with larger-scale extended emission and multiple peaks, the two most prominent of them toward 16293E and the location of the protostars IRAS 16293 A and B. Since the spatial resolution in our observations is about 20", the two protostars, which are separated by 5", cannot be resolved individually. Instead, the observations aim to examine the large-scale molecular structure of L1689N and interactions of outflows with the cloud. Since the prestellar core 16293E is a cold source (e.g., Stark et al. 2004), it is best traced by the N₂H⁺ (J = 3−2), N₂D⁺ (J = 4−3), and DCO⁺ (J = 5−4) spectral lines, as shown in Fig. 1. In contrast, the extended envelope surrounding the A/B sources is traced by almost all molecules for which emission is detected in the data. Both cloud cores are slightly resolved in the observations and a number of molecular lines trace larger scale emission in which the cores are embedded.

3.2 The outflow structure

In order to trace the general gas content of the cloud, CO is used as a tracer for molecular hydrogen. Its high abundance results in strong line emission with extended line wings that help us to identify molecular outflows based on the observed gas velocities relative to the cloud along the line of sight. Since the systemic cloud velocity is about 4 km s⁻¹, velocities below this value trace blueshifted emission, while higher velocities were measured in case of redshifted emission from gas moving away from us.

Figure 2 shows integrated intensity (moment 0) maps of the CO (6⁻5) transition, which reveal the complex large-scale structure in L1689N. In previous studies, Walker et al. (1988) identified several emission peaks based on CO (2−1) observations of this region. Later, Hirano et al. (2001) and Castets et al. (2001) renamed these emission peaks as E1, E2, W1, W2, and HE2. We have included the position of these emission peaks (see Table B.1) in Fig. 2. In comparison to these previous works, we choose slightly different coordinates to describe the emission peaks based in our formaldehyde maps (see Sect. 3.3). In addition, the emission peak SW can be seen in the southwest of the CO (6−5) maps, which was not covered by our LAsMA observations. The overall CO morphology is dominated by two bipolar outflows which both originate in the vicinity of IRAS 16293 A/B.

One of these large-scale outflows is the east-west outflow, which contributes to emission peak W1 with the western red lobe. The blue lobe of this outflow extends broadly across the region east of IRAS 16293 A/B and appears to split around 16293E, which was already noted by Stark et al. (2004). As indicated by the blue arrows and velocity contours of CO (3−2) in the bottom panel of Fig. 2, this split results in a broader outflow region with higher velocities north of the prestellar core and a more compact region with lower relative velocities south of it. The latter corresponds to the emission peak HE2.

The second outflow extends in northeast-southwest direction and is the main contributor to emission at the peaks E1 and E2. These positions show redshifted emission with high relative velocities at the position E2. The southwestern counterpart of this outflow is not covered by our broad bandwidth LAsMA observations but instead can be seen in the further extended SEPIA maps of the CO (6−5) transition. Emission peaks originating from these outflows are seen in the maps of numerous species, which are analyzed in Sects. 3.3, 3.4, and 4.2.

3.3 Formaldehyde emission peaks

The formaldehyde (H₂CO) and methanol (CH₃OH) molecules are known to be good tracers of molecular outflows (Bachiller & Pérez Gutiérrez 1997). Since formaldehyde is important for the analysis in Sect. 4, we discuss the observational results from the H₂CO lines in more detail.

Figure 3 presents the spatial distribution of the H₂CO $J_{K_a,K_c}=4_{1,4}-3_{1,3}$ transition using different velocity ranges for integrating the emission. While the moment 0 map shown in the upper panel includes emission from the full H₂CO line core, the middle and bottom panels of Fig. 3 respectively capture the blue and redshifted emission seen across the region. We note that the distribution of the emission is similar to that of CO (6−5), concentrated around the previously defined peaks. While these maps do not show significant emission at 16293E, a narrow line can be seen in the spectra of Fig. 4 at this position.

In addition, these maps clearly show a separation of the emission peak W2 north of IRAS 16293 A/B from the remaining outflow structure. Since we do not detect a counterpart at a southern location, it is unlikely that the emission originates from a bipolar outflow. In contrast, SCUBA 450 μm maps shown in Sect. 5.1 reveal emission at this position, suggesting a cold dust source to be present.

Analogous to CO, the lines of formaldehyde also show profiles with extended line wings. In order to clarify differences of the line profiles across the L1689N region, Fig. 4 displays the formaldehyde spectra which occur at each of the emission peaks. While the line profile of the H₂CO $J_{K_a,K_c}=4_{1,4}-3_{1,3}$ transition is very narrow at the position of 16293E, all other peak positions show a broader line profile. Additionally, line profiles at A/B and W2 show signs of self-absorption due to high optical depths.

Especially the redshifted emission at the emission peaks W1 and E2 extends to high velocities of about 15 km s⁻¹. This is in agreement with our observations of the CO (3−2) transition, for which we detect outflow velocities of up to 20 km s⁻¹. Similar velocities were also detected by Hirano et al. (2001) in SiO emission, who also found velocities up to 20 km s⁻¹. These positions therefore are likely to be dominated by shocked gas as a consequence of outflow activities from IRAS 16293 A.

3.4 Distribution of different molecules

In order to obtain an overview of the molecular content in the region, we examined all emission peaks for occurring transitions. The results of this process are summarized in Table 1 and Appendix C.

In general, a majority of the considered molecules is detected around the protostars IRAS 16293 A and B, as this is one of the positions at which the line identification was carried out. Exceptions are the transitions of doubly deuterated ammonia and the 372.4 GHz ortho-H₂D+ line, which, at our sensitivity, can only be detected near 16293E.

Most of the sulfur bearing species are present in the surroundings of the protostars A/B and in the E1 outflow peak. The absence of strong emission in the other outflows is unusual, as those species are characteristic for warm gas and shocks, as described by Wakelam et al. (2004). A reason for this might be that the observed SO and SO₂ transitions arise from high energies above the ground states and have small Einstein A coefficients, such that they are only seen toward positions with high excitation and column densities.

Previous works (e.g., Hirano et al. 2001; Castets et al. 2001) have shown that SiO emission traces very well the outflows in the L1689N region, which is true for fast shocks in general (Schilke et al. 1997). The observed frequency windows only cover the SiO (8−7) transition, which is detected mostly in the envelope of IRAS 16293−2422 A/B. Weak SiO emission is also detected toward the E1 emission peak. Again, the excitation conditions at the other positions might not be sufficient for the relatively high-J line of SiO to be detected.

The prestellar core 16293E shows emission of many deuterated species such as DCO⁺, N₂D⁺, and NHD₂. This confirms it to be a quiescent cold source, as the enhanced deuteration is a result of strong exothermic reactions that proceed in cold material (see e.g., Roberts & Millar 2000). In addition to these deuterated species, other molecules are detected and mapped for the first time in this source, such as CN, NO, and weak emission of CS and CH₃OH.

As can be seen in Table 1 and Appendix J, the distribution of methanol is similar to that of formaldehyde, as both molecules can be observed in all of the positions. The reason for their comparable distribution is the formation of methanol from formaldehyde through hydrogenation on CO-rich ices (Fuchs et al. 2009, 2020).

3.5 Rotation of the envelope

The protostellar objects A/B are surrounded by a large gas and dust envelope. This structure has a mass of about 4–6 M_⊙ (Jacobsen et al. 2018; Ladjelate et al. 2020) and radius of about 6000–8000 au (Schöier et al. 2002; Crimier et al. 2010; Jacobsen et al. 2018). Several studies have focused in determining the large-scale properties of this structure based on the kinematics of a number of molecules. For example, the extended envelope was suggested to have infall motions based on the blue-skewed profile seen in CS (7−6) and (5−4) molecular maps at spatial scales of about 39 (Narayanan et al. 1998). In addition, based on CS (3−2) and (2−1) observations, Menten et al. (1987) suggested the presence of a velocity gradient that could be associated with the rotation of the envelope.

To investigate such rotational motions, we have produced velocity field maps for some species of the close environment around IRAS 16292−2422 A/B. In Fig. 5 we show the moment 0 and moment 1 maps of the C¹⁷O (3−2), H¹³CO⁺ (4−3), N₂D⁺ (4−3), and N₂H⁺ (3−2) transitions. These maps include emission of the line core, covering a velocity interval of 2 − 6 km s⁻¹ for C¹⁷O (3−2), H¹³CO⁺ (4−3). In case of N₂D⁺ (4−3), emission between 2.5−5 km s⁻¹ is integrated in order to exclude nearby HDCO emission. For N₂H⁺ (3−2), we only consider emission from the main hyperfine structure components, by choosing an integration range of 2.3−5 km s⁻¹. We note that the N₂D⁺ (4−3) transition and the main component of the N₂H⁺ (3−2) transition both consist of multiple unresolved components as a consequence of the ¹⁴N hyperfine structure. We therefore assume that possibly occurring changes in relative intensities of these components due to optical depth effects do not affect the observed velocity field.

As a reference, a yellow circle of 50 radius is included to Fig. 5 which indicates the approximate extent of the envelope. A large-scale velocity gradient is observed in the NE–SW direction in the C¹⁷O moment 1 map. For H¹³CO⁺, the velocity gradient is mostly aligned in the N–S direction, while for N₂D⁺ it is distributed in the E–W direction, although the molecular gas distribution follows the same trend as in the C¹⁷O map. Less emission is shown due to the clipping limits used to create the map (10σ for N₂D⁺ and 4σ for C¹⁷O). For N₂H⁺, the velocity field is more complex, but the velocity spatial distribution resembles overall that of the C¹⁷O as well. Closer to IRAS 16293−2422, the velocity field in the N-S direction prevails for most of the maps, which could be interpreted as rotation of the envelope itself. This would agree with the interpretation of Menten et al. (1987). However, based on the full velocity field morphology seen in our maps, we cannot disentangle a pure rotation motion associated with the envelope and the larger cloud scale velocity field that dominates the whole studied region.

4 Analysis

4.1 Estimation of the gas temperature in L1689N

The formaldehyde molecule has a long-standing history as thermometer and density probe of molecular clouds. This is because the H₂CO molecule is a slightly asymmetric rotor and the transitions between levels with different K_a values are dominated by collisions. Therefore, line ratios involving different K_a ladders can be used to determine the kinetic temperature of the gas (e.g., see Mangum & Wootten 1993). To derive the excitation temperature T_ex, which is equal to the kinetic temperature under the assumption of a local thermodynamic equilibrium (LTE) conditions, we made use of the H₂CO rotational transitions detected in our observations.

In the LTE approximation, the excitation temperature between two energy levels u and l can be written as $$T_{\text{ex}}=-\frac{E_{\text{u}}-E_{\text{l}}}{k}{\left[\ln \left(\frac{N_{\text{u}}}{N_{\text{l}}}\frac{g_{\text{l}}}{g_{\text{u}}}\right)\right]}^{-1},$$(1)

where E_u and E_l are the energies of the upper and lower levels, respectively, g_u and g_l are the degeneracies and N_u and N_l are the corresponding column densities (see Appendix D for more details).

The observed formaldehyde transitions allow the derivation of excitation temperatures for ortho- and para-H₂CO separately. Since only three transitions per species are available, these temperatures were obtained from the respective line ratios. For para-H₂CO was therefore possible to use the line ratios (4_2,2–3_2,1)/(4_0,4–3_0,3) and (4_2,3–3_2,2)/(4_0,4–3_0,3). In case of ortho-H₂CO, the (4_3,2–3_3,1) and (4_3,1–3_3,0) lines are blended. For this reason, the ratio [(4_3,2–3_3,1) + (4_3,1–3_3,0)]/[2 × (4_1,4−3_1,3)] was considered. In the Table D.1, we show the main parameters of the considered H₂CO transitions based on the CDMS and JPL catalogs.

The temperatures were calculated at each position, at which both of the considered transitions are detected with a 3σ significance. In order to include most of the emission, main beam temperatures were integrated in a velocity-range between 2 km s⁻¹ and 6 km s⁻¹ for computing the intensities of para-H₂CO. To properly include the emission of both blended ortho-H₂CO lines, the integrations for this species were carried out on the H₂CO (4_3,2−3_3,1) spectra in a velocity interval from −5 km s⁻¹ to 10 km s⁻¹.

Based on the created line ratio and temperature maps shown in Fig. 6, it was possible to derive the temperatures at the previously discussed positions, by including all mapped temperatures in a 10" radius for computing a weighted average. Statistical uncertainties were estimated based on Gaussian error propagation of the map RMS (see Table C.1). The resulting temperatures and uncertainties are shown in Table 2. These uncertainties do not include the calibration uncertainty, which is conservatively estimated to be about 20%. It can be expected that this uncertainty is additionally transferred to the temperature values.

The temperatures of para-H₂CO were calculated by averaging the values from both para-H₂CO maps at the associated positions. As the transitions of formaldehyde with high upper level energies are too weak for temperature estimates in the direction of 16293E, a temperature of 12 K based on the results of Stark et al. (2004) is assumed for further calculations.

It is worth noticing that the derived temperatures in the ortho-H₂CO map are about 20-30 K higher than derived from the para-transitions. The reason for this deviation might be the upper level of the ortho-H₂CO transitions. With an upper level energy of 140.9 K, these transitions trace a population which requires higher temperatures to be excited, compared to the para-H₂CO transitions.

Previously, van Dishoeck et al. (1995) studied the emission of H₂CO toward IRAS 16293−2422 A/B based on observations obtained with the Caltech Submilimeter Observatory (CSO) and the JCMT. They used the rotational diagram method and found a rotational temperature of 80 ± 10 K. Also, from H₂CO line ratios within the same ladder, they derived values for the kinetic temperature between 60 and 140 K. Later, Ceccarelli et al. (2000) used the H₂CO data from van Dishoeck et al. (1995, among others) and confirmed that the emission in IRAS 16293−2422 A/B originates from two regions: A hot core region with temperatures above 100 K and an outer and colder layer with lower H₂CO abundance. Using p-H₂CO, we obtained an average temperature value of 62.4 ± 0.7 K (see Table 2), which is somewhat low compared with the values from previous studies. The o-H₂CO temperature of 90.1 ± 4.4 K is in better agreement with both van Dishoeck et al. (1995) and Ceccarelli et al. (2000). For a more detailed comparison of the temperatures for different layers in IRAS 16293−2422, see Appendix G.

Furthermore, Castets et al. (2001) studied the H₂CO and SiO emission based on observations with the IRAM 30 meter telescope, the Swedish-ESO Submillimetre Telescope (SEST), and the Infrared Space Observatory (ISO) and derived values for the other emission peaks between 80 and 150 K. In this work, we derive lower p-H₂CO temperature values (about 35 K on average) for the rest of the emission peaks. Since H₂CO is more reliable thermometer than SiO, we argue that these lower values reported in Table 2 are more accurate.

In order to test the validity of the assumption of LTE, we compared the calculated temperatures and line ratios from ortho-and para-H₂CO with non-LTE models. These models are computed with the RADEX radiative transfer code (van der Tak et al. 2007) considering different H₂ volume densities and using collision rates from Wiesenfeld & Faure (2013). The detailed RADEX computation can be found in Appendix E. We find that RADEX and LTE models are in better agreement for higher volume densities, while the temperature values are underestimated by about 20% for lower values of the density. These effects of non-LTE are partially compensated by opacity effects at the position of the protostars A/B.

4.2 Column densities

To derive the column densities for all the species, we have used the LTE radiative transfer CLASS extension weeds. This software allows us to compute a synthetic spectrum based on a set of physical parameters, namely column density, temperature, source size, velocity shift and line width (Maret et al. 2011). For the computation, weeds solves the radiative transfer equation by assuming that the emission originates from a single layer. The synthetic spectrum emulates the individual line profiles for each species and can be directly compared with the observations.

Applying the formalism introduced in Appendix D, the total column density N_tot can be related to the column density of an upper energy level E_u by $$\frac{N_{\text{tot}}}{N_{\text{u}}}=\frac{Q_{\text{rot}}}{g_{\text{u}}}\exp \left(\frac{E_{\text{u}}}{k{T}_{\text{ex}}}\right),$$(2)

where Q_rot(T) is the rotational partition function which describes the sum over all rotational energy levels in a molecule for a given excitation temperature (Mangum & Shirley 2015).

Inserting the derived temperature estimates in Eq. (2) would, in combination with Eq. (D.2), also allow the computation of column densities based on the observed line profiles. The advantage of determining the column densities with weeds is that the software additionally considers line blending with other molecular species and optical depths effects in the calculation of the line profiles, thus giving us a more reliable estimate.

For deriving column densities at each of the seven emission peaks listed in Table B.1, we extracted and averaged the spectra in a 10" radius around the respective positions. This radius is motivated by the beam size of the telescope and the fact that the cores are slightly resolved in the observations. In order to obtain comparable results for all transitions, we adopted a constant source size of 20" in the computation of the synthetic spectra.

Since the highest upper energy level for the majority of the observed transitions is about 40 K, hence much lower than the upper level energy of 140.9 K of the higher energy ortho-H₂CO lines, the temperatures derived from para-H₂CO were used in the computations of synthetic spectra.

Some of the presented lines, for example the formaldehyde transitions, show self-absorption due to nonnegligible optical thickness. In these cases, the fit is guided by the line wings instead of the self-absorbed line core which may lead to an underestimation of the derived column densities. While fitting the line wings gives reasonable results for most of the self-absorbed lines, it is not possible to compute satisfactory models for the optical thick transitions of CO (3−2), HCO⁺ (4−3), and HCN (4−3). Fitting these lines would require at least a two layer model, which goes beyond the scope of this work.

Instead, the column densities of these molecules were estimated based on the column densities of their less abundant isotopologues C¹⁷O, H¹³CO⁺, and H¹³CN, by converting the derived column densities according to the isotope ratios in the interstellar medium. For this we adopted a ratio of ¹²C/¹³C = 68 from Milam et al. (2005) and a ratio of ¹⁶O/¹⁷O = 1790 from the values suggested by Wilson & Rood (1994). The observed H¹³CN (4−3) transition is blended with the SO₂ (13_2,12−12_1,11) transition, due to which the derived SO₂ model has been subtracted from the spectrum before estimating the H¹³CN column density.

In cases where we did not detect emission stronger than three times the RMS, upper limits on the column densities were constrained instead. For these, the median line width of the observed transitions at a respective position was applied to determine the column density at which the resulting computed line profile would fulfill the detection criterion.

Appendix F presents an overview of the derived molecular abundances and simulated line profiles. Based on a conservative estimate on the calibration uncertainty of about 20%, it can be expected that this uncertainty is transferred to the column density values. As an example, the radiative transfer model for the H₂CS (8_1,7−7_1,6) transition at 278.9 GHz at the position of IRAS 16293 A/B is shown in Fig. F.1.

Figure 7 visualizes the column densities and molecular abundances relative to CO at the considered emission peaks in L1689N. Figure 7a displays the column densities in a mosaic plot with logarithmically scaled tiles, in order to give a quick overview of the presence and absence of species at the respective positions. While this scaling emphasizes the presence of less abundant species, the order of magnitudes differences in the molecular abundances are seen more clearly in the bar plot shown in Fig. 7b. These figures show that the molecular species are widely distributed across the whole cloud. The few species detected exclusively at the protostellar core IRAS 16293 A/B are mainly less abundant isotopes of otherwise abundant species and molecules for which we only cover fainter transitions. The complexity of the distribution of the species across the cloud is discussed further in Sect. 5.4.

Since the column densities of H₂CO are important for the further derivation of the H₂ volume densities (see Sect. 4.3), we also tested the plausibility of the derived values by creating multicomponent non-LTE models for the molecule using the CASSIS⁵-RADEX software. The details of this line modeling can be found in Appendix G.

We find that the self-absorbed line profiles at IRAS 16293 A/B and W2 can be reproduced with a three and two layer model respectively. The three physical components toward A/B hereby consist of a hot corino in addition to a warm and an extended cold envelope. The derived LTE column density at this position is in between the values for the hot corino and warm envelope. A similar result is obtained by the line modeling at the W2 position, in which the column density derived with LTE assumptions is also in between the values of both considered physical components. Single layer non-LTE models at the positions E and HE2 are in agreement with the LTE models.

4.3 H₂ volume densities

To determine the H₂ volume densities n present in L1689N, a similar approach as above can be used for the line ratios and excitation temperatures derived from the H₂CO transitions (5_1,5− 4_1,4)/(4_1,4−3_1,3). In order to compare these transitions without the need to consider different filling factors, the maps were smoothed to the same resolution prior to the further analysis.

Since the ratios of lines with different J are sensitive to density deviations (Mangum & Wootten 1993), the excitation temperature derived with these transitions is much lower than the actual kinetic temperature in the cloud. Due to this, it is possible to get a crude estimate of the H₂ volume density by comparing RADEX models for the (5_1,5⁻4_1,4)/(4_1,4⁻3_1,3) ratio with the temperatures derived in Sect. 4.1. In order to explain this process, the H₂ volume density is derived in detail for the position of IRAS 16293 A/B.

Assuming LTE conditions, the kinetic temperature at this position is 62.4 K, as derived from para-H₂CO transitions with J_up = 4 (see Sect. 4.1). Mapping of the H₂CO (5_1,5−4_1,4)/(4_1,4− 3_1,3) line ratios results in a value of 1.02 at A/B, hence an excitation temperature of about 38 K, much lower than the kinetic temperature. By comparing RADEX models with a range of H₂ volume densities as displayed in Fig. 8a, it can be seen that a volume density estimate of log₁₀(n(H₂)) = 6.7 predicts the kinetic temperature of 62.4 K for the observed line ratio of 1.02. Based on the temperature uncertainties, this method eventually leads to uncertainties in the H₂ volume density of a factor on the order of 1.5.

RADEX models with an ortho-H₂CO column density of 2.9 × 10¹⁴ cm⁻² were computed for the position of IRAS 16293 A/B, whereas a common column density of 3.9 × 10¹³ cm⁻² was applied for the other positions, on which we observe the transitions to be mostly optically thin. The RADEX models are shown next to the datapoints for the H₂CO J_up = 5/J_up = 4 ratios in Fig. 8. The resulting H₂ volume densities can be found in Table 3.

While we observe high volume densities of 5.0 × 10⁶ cm⁻³ at the positions of the protostars IRAS 16293 A/B, the emission peaks E2 and HE2 which relate to outer outflow positions show lower values of 5.0 × 10⁵ cm⁻³ and 6.3 × 10⁵ cm⁻³. The H₂ volume density estimations at the remaining emission peaks E1, W1, and W2 result in a common intermediate value of 1.6 × 10⁶ cm⁻³. These volume densities are in agreement with the values considered to test the assumption of LTE conditions in Fig. E.1.

5 Discussion

5.1 Dust continuum maps of L1689N

In order to compare the results of Sect. 3 with the dust continuum, archival data of the region surrounding L1689N were obtained at wavelengths between 4.5 μm and 450 μm, which were taken with the following instruments/telescopes: IRAC2/Spitzer⁶, PACS/Herschel Space Observatory⁷ and the SCUBA-2 JCMT⁸.

A Spitzer 4.5 μm continuum map is shown in Fig. 9a. Besides continuum emission, lines from H₂ and vibrational excited CO are included in the 4.5 μm band, both of which can be excited by shocks (Reach et al. 2006). This band shows emission near all outflow-positions except for E2, whereas the emission does not peak at IRAS 16293 A B, E, and W2. In addition, the emission near 16293E indicates a splitting of the eastern blue lobe, of which the northern part is followed up by two bow-shocks in the east, slightly outside of our map coverage with LAsMA. Nevertheless, the CO (3−2) submillimeter line (see bottom panel of Fig. 2) shows strong emission in the direction of these bow-shocks, indicating that they likely originate from the E-W-outflow.

Another strong continuum source exists south of IRAS 16293 A/B, which Pagani et al. (2015) identify as a young stellar object driving an outflow pointing toward the position HE2. This source can also be seen in the Herschel 70 μm map (Fig. 9b) which shows the dust continuum map together with contours of blue-shifted CO (3−2) emission. Since CO clearly shows the influence of the E–W-outflow on HE2 (also discussed in Sect. 3.2), it is likely that the outflows of both sources contribute to the conditions in this region.

As longer wavelengths can be attributed to colder dust, the maps of 160 μm and 450 μm (Figs. 9c, d) emission show the bulk of the cloud material. As it can be seen in these maps, the cold core 16293E is more clearly detected at these wavelengths. In addition, more extended emission is observed in all directions, particularly to the southeast of 16293E at 160 μm.

Finally, the coldest dust components are probed by continuum maps with wavelengths of 450 μm and above. As shown in Fig. 9d, this long wavelength emission is mainly seen in the close vicinity of the cloud cores IRAS 16293 A/B and E. Additional emission can be seen at the position W2 north of IRAS 16293 A/B, which will be discussed further in Sect. 5.5.

5.2 Spatial distribution of selected species at the emission peaks

The maps shown in Appendix J show the morphology of the molecules in the IRAS 16293 complex. In order to visualize more easily the differences in the spatial distribution among the emission peaks, we have overlaid the emission for some species that have strong large-scale emission for a more direct comparison.

In Fig. 10a, we show the moment 0 maps of the H₂CO (4₁₄−3₁₃) transition (grayscale), overlaid with the N₂D⁺ (4−3), DCO⁺ (5−4), and N₂H⁺ (3−2) emission in contours. From this figure, we see that all four species have the same peak around IRAS 16293 A/B. Indeed, the emission of all considered molecules in all panels peaks at the exact position of the embedded protostars A and B. In contrast, Fig. 10a shows that the emission peaks in the prestellar core 16293E differ for the presented species. The peak of DCO⁺ emission is located 5.3" toward the south of the reference position for 16293E (marked with a white point as in previous figures), while the N₂H⁺ peak is located 1.1" to the north and the N₂D⁺ is 7.2" offset toward the east of 16293E respectively. The offsets between some species such as N₂D⁺ and ND3 were noted by Lis et al. (2016), suggesting that they might be due to either intrinsic abundance gradients or optical depth effects.

In Fig. 10b, we compare the spatial distribution of H₂CO (4₁₄−3₁₃) transition (grayscale), together with the SO (7₆−6₅), CS (6−5), and the C¹⁷O (3−2) emission peaks. We see a more pronounced difference around E1 for CS and SO, which have an offset from E1 of 12.1" and 6.8" respectively. The C¹⁷O emission around E1 is weak and does not extend toward E2. It seems there is also an offset between these four species at HE2 but it is more complicated to see as the emission is not very strong at that position.

In Fig. 10c, we show a similar comparison between HNC (4−3), CS (6−5), and HCN (4−3). At HE2, there is also an offset between CS and HCN. In contrast, these species are distributed quite similar around E1, while HNC which seem to peak about 30" to the south of E1. In Fig. 10d, the HCO⁺ (4−3), H¹³CO⁺ (4−3), and SiO (8−7) in contours are shown. These 3 species peak at about the same position south of E1. Of these species, only H₂CO and HCO⁺ are present at HE2 and seem to peak at about the same position.

Hirano et al. (2001) found an anticorrelation between SiO (2−1) and H¹³CO⁺ (1−0), (their Fig. 4), suggesting that the SiO emission arises from the region in which the outflow from IRAS 16293 A/B is interacting with the prestellar core 16293E.

We observe higher energy transitions of these two species in our data (SiO (8−7) and H¹³CO⁺ (4−3)) but we do not see such anticorrelation, probably because our observed transitions have larger critical densities.

It is worth noticing that the emission of H₂D⁺ (1−0) behaves differently to other species in the prestellar core 16293E. In Fig. 11, we show the continuum emission at 450μm obtained with the JCMT SCUBA-2 (grayscale), overlaid with emission of H₂D⁺, NH₂D, NHD₂, and N₂H⁺ in contours. The deuterated species show an emission peak in the northeast with respect to the reference position of 16293E as traced by N₂H⁺ (see Fig. 10). Interestingly, the NH₂D emission also does not peak exactly at the position of the A/B protostars, in contrast to the majority of the molecular tracers. The H₂D⁺ emission seems to be located exclusively in the eastern part of the prestellar core, without expanding further toward the N₂H⁺ peak. Pagani et al. (2015) also note this behavior of H₂D⁺ when comparing their data with the distribution of the continuum emission. Based on the finding that the position of the dust emission peak changes with wavelength, they suggest a temperature gradient in this region with H₂D⁺ tracing the coldest and highest column density gas.

The offset between the continuum emission and the H₂D⁺ peaks indeed is likely related to the difference in the physical conditions (temperature and density) between these two spots. In particular, determining the temperature from H₂CO line ratios at these peaks is not possible since we do not detect emission of the excited H₂CO transitions (4₂₃−3₂₂,4₂₂−3₂₁,4₃₁−3₃₀,4₃₂−3₃₁) at these positions. In contrast, N₂H⁺ is well detected in the vicinity of 16293E.

Figure 12 shows a map of the line-ratio between the N₂H⁺ (4−3) and N₂H⁺ (3−2) transitions, where we also display the H₂D⁺ emission in blue contours as a reference. Several pixels close to the H₂D⁺ peak show lower line ratios as compared to the continuum peak position. Since the N₂H⁺ (4−3)/N₂H⁺ (3−2) ratio probes to first order the pressure of the gas, Fig. 12 indicates that H₂D⁺ peaks at a position of significant change in pressure. This might also be indicative of lower temperatures in this region, which would agree with the interpretation of having the coldest part of the prestellar core close to the H₂D⁺ peak. Unfortunately, it is not straightforward to assess if there are also important density differences between the two peaks. A density gradient in this region would imply variations in the observed line ratios, which prevents us from discriminating between the effects of the two physical quantities on this region.

Note that the W2 position shows large values for the N₂H⁺ (4−3)/N₂H⁺ (3−2) line ratio. The spectra at the corresponding pixels show an artifact close to the N₂H⁺ (4−3) transition, which is likely the cause of these seemingly high line ratios. The actual N₂H⁺ (4−3) line emission at W2 is very faint and not detected with a 3σ significance.

Based on the contour plots shown in Fig. 10, it is clear that we observe the spacial distribution of different molecular species to differ across the whole eastern cloud. This seems to be in contrast of what we observe at the position of the protostars A/B, for which we observe almost all species to peak at the exact position of the protostellar core. Also, we do not see significant differences between species in the peak positions at W1. In contrast, species do not show a clear emission peak at W2, but rather broadly distributed weak emission. The observed morphological differences between the studied species could be related to the interaction between the outflows from IRAS 16293–2422 and the cold core 16293E. Therefore, we study the kinematics in 16293E with more detail in the next section.

5.3 The interaction between 16293E and the outflows from IRAS 16293–2422 A/B

As previously discussed, one of the outflows arising from IRAS 16293–2422 AB is proposed to be interacting with the prestellar core 16293E (e.g., see Stark et al. 2004; Lis et al. 2016). As earlier works were limited to observations of a small number of species, we aim to test this scenario with our unbiased survey of the region.

5.3.1 Velocity offsets between the deuterated and nondeuterated species at 16293E

Since both deuterated and nondeuterated species have been studied in the prestellar core 16293E, a difference in the velocities between them with respect to the cloud velocity (4.0 km s⁻¹) has been observed. In particular, deuterium-bearing species such as ND, N₂D⁺, DCO⁺, D₂H⁺, and ND₂H show peak velocities closer to 3.5−3.7 km s⁻¹ (Lis et al. 2002; Vastel et al. 2004, 2012; Gerin et al. 2006; Bacmann et al. 2016), while other tracers such as H¹³CO⁺ and C¹⁷O show peak velocities closer to 3.8−4.0 km s⁻¹ (Lis et al. 2002; Stark et al. 2004). Vastel et al. (2012) suggest that such behavior could be due to the interaction between the outflows from A/B and the prestellar core 16293E, introducing a difference between the deuterated and nondeuterated molecules. Since this effect could also be introduced by the different beam sizes corresponding to the different observations of these species, we tested if this trend can be confirmed from our observations, which have the same beam.

Based on Gaussian fits to estimate the integrated intensity of all the observed species reported in Table F.2, it can be concluded that we do see a similar trend in our data. Indeed, the deuterated species H₂D⁺, DNC, NH₂D, NHD₂, and N₂D⁺ have a peak velocity of 3.6 km s⁻¹, while the DCO⁺ and N₂H⁺ have a peak velocity of 3.7 km s⁻¹. On the other hand, the rest of the species have larger peak velocities between 3.8 and 4.0 km s⁻¹, closer to the cloud velocity (except for CN and NO). Interestingly, this velocity offset trend for deuterated species is not followed by the HDCO, D₂CO, and DCN molecules which appear to be at the cloud velocity.

5.3.2 Kinematics in the vicinity of 16293E

Lis et al. (2016) study the emission of NH₂D based on HIFI Herschel Space Observatory observations at 494.454 GHz. They extracted the spectrum of this transition at three different positions in 16293E and estimated the peak velocities and line widths of the individual spectra (see their Fig. 5). They conclude that the change in the peak velocity and line width of NH₂D in the vicinity of 16293E is an evidence of the outflowing gas interacting with this source.

To investigate if such behavior holds for other molecular species in this source, we made use of our unbiased survey to see if this is valid for other molecular species. As seen in Figs. 10 and 11, different species are distributed differently along 16293E. Therefore, we have selected three positions where we have clear emission for all deuterated species that peak near 16293E: H₂D⁺, DCO⁺, DNC, N₂D⁺, NHD₂, and NH₂D. The positions are shown in Fig. 11 and are located at the coordinates: NE: (16^h32^m30.6^s, −24°28'46.7"), NW: (16^h32^m28.9^s, −24°28'46.7"), SW: (16^h32^m28.9^s, −24°29'13.7"). Although these positions are different to the ones chosen by Lis et al. (2016), they are still well separated. We have extracted the spectra at each of these positions and fit a Gaussian to the line profiles using the CLASS software to determine their velocity peak and line widths. The results from the fits are shown in Table I.1. It is worth mentioning that both NHD₂ and NH₂D present hyperfine structure in their spectra, resulting in three lines close in frequency for each species. We fit the visible HFS components under the assumption of same line widths for all components using the component-separations from the CDMS. The fit spectra are shown in Fig. 13 and the Appendix I.

Most species show a similar trend as Lis et al. (2016) find for NH₂D. More precisely, we see that the peak velocity increases from the NE to the SW position of the core, which is indicative of the interaction with the outflows powered by IRAS 16293–2422. Although, this trend was also observed to be present on larger scales, as shown in Fig. 5 of Sect. 3.5. Due to this larger scale velocity field, the peak velocity alone is not a clear indicator of motions from the proposed core outflow interactions. Nevertheless, a similar trend of line widths increasing from NE to SW positions is observed near 16293E, with the exception of DCO⁺ and the weak lines of NHD₂. This may be a good indicator for a core-outflow interaction, as the observed velocity gradient should not influence the line widths. Since the emission of DCO⁺ peaks further southwest compared to the other considered species, it may trace a differently effected layer of the gas. As compared to Lis et al. (2016), we extracted the analyzed spectra further to the inner part of the cold cloud core, which may be located at the NH₂D emission peak. Such trends therefore might indicate that the interaction of the outflows is propagating toward the inner part of the core.

5.4 Chemical differences at the selected positions

In the previous sections, we analyze the distribution of molecular species across the L1689N cloud and derived column densities of these species at the visible emission peaks. Table 1 gives an overview about the detected species at all peaks, while the derived molecular abundances given in Table F.8 are visualized in Fig. 7. To emphasize further the chemical differences between pairs of sources, we created Venn diagrams that cover the occurrences of different species (see Appendix H).

These diagrams show that E1 not only contains a larger number number of molecules than E2, but all the species observed in E2 are also found in E1. In a previous study, Castets et al. (2001) suggested that the simultaneous presence of H₂CO and SiO emission in E1 might indicate that E1 is actually younger than E2, given that this source was observed to only emit SiO and this molecule is expected to last longer than H₂CO. From our maps, we see that E1 is indeed chemically richer than E2, while we do observe the presence of H₂CO emission at E2. Interestingly, both E1 and E2 show SO emission, while only E1 shows SO₂ emission.

It can be noticed that the emission peak W1 and W2 share the same species, although W2 shows SO, HNC, and CS less abundant, while being richer in N₂H⁺. The emission peaks that show the largest number of detected species besides IRAS 16293 A/B correspond to the prestellar core 16293E and HE2. Similar to previous studies, we mostly see probes of cold and dense gas at the E position. We find comparable molecules toward HE2 and E, with the exception of H₂D⁺ and NHD2 being absent in the HE2 emission peak. While 16293E is a quiescent cloud core, the line profiles at the position HE2 (see Fig. 4) show rather blueshifted shocked gas. As can be seen in Fig. 7, all species detected at IRAS 16293 A/B and E are also observed at HE2, which is in favor of an interaction between the E–W outflow and the prestellar core.

Based on the abundant molecules, E1 can be interpreted as prime region of core-outflow interactions, since it shows moderate line wings and both, dense gas tracers as well as likely shock tracers such as methanol and formaldehyde. Shock dominated positions are W1 and E2, both with high velocity red-shifted wings and a lack of dense gas tracers, while being bright in methanol and formaldehyde. These two peaks exhibit basically the same number and type of observed species. The only difference is that W1 has CS emission, while E2 does not, although the CS maps suggest that a contribution from the envelope around A/B sources can affect W1.

5.5 On the nature of W2

The emission at the position of the W2 peak was already reported in early studies of the L1689N region (e.g., Mizuno et al. 1990). However, the origin of this emission has not yet been completely clarified. On the one hand, Hirano et al. (2001) interpreted the peak as northern boundary of the western red lobe of the EW outflow where the outflowing gas encounters the ambient medium. On the other hand, Castets et al. (2001) suggested that W2 could be either an additional outflow of IRAS16293 A/B with a deeply embedded and therefore invisible counterpart, or an embedded source itself.

As can be seen in Fig. 9d, the SCUBA-2 map at 450 μm reveals long wavelength continuum emission at the W2 emission peak, which suggests the presence of a cold dust source at this position. This conjecture is supported by the line profiles at this position (see Fig. 4), which show moderate line wings alongside self-absorption. Compared to the shock dominated regions W1 and E2, this region therefore seems to be relatively quiescent while the high optical depth suggests the presence of dense gas. Based on the H₂ densities given in Table 3, the embedded object is likely to be less massive then the IRAS 16293 A/B cores.

Comparing the occurring species at this position, we note a similarity between the chemistry of W1 and W2 (see Sect. 5.4). While we detect emission from a large number of species toward W2, the distribution of this emission (see Fig. 10 and Appendix J) only shows a clear peak for H₂CO and CH₃OH, while other species show rather broadly distributed weak emission in the region. Emission from other species observed at W2 therefore may be attributed to the extended envelope of IRAS 16293 A/B instead of the cold source embedded at W2.

It still remains unclear if the embedded object seen at W2 interacts with the outflows of IRAS 16293 A. The CO (3–2) maps suggest that W2 indeed is located at the northern boundaries of the E–W outflow, while the narrow line profiles at W2 do not suggest a strong interaction, as it is the case at E1. There may be a scenario in which the small scale NW-SE outflow, which Girart et al. (2014) observe to be redirected by IRAS 16293 B, extends on larger scales up to W2. Nevertheless, since we do not observe a southern counterpart of this outflow with a spacing comparable to W2, it is unlikely that this outflow extends to larger scales.

5.6 Molecular abundances in L1689N as compared to the literature

Adapting a canonical CO/H₂ ratio of 10⁻⁴, the derived column densities of HCO⁺ and DCO⁺ species at the position of IRAS 16293 A/B (see Sect. 4.2 and Table F.9) are comparable to the estimates from Stark et al. (2004), who derive HCO⁺/H₂ and DCO⁺/H₂ ratios on the order of 1 × 10⁻⁹ and 2 × 10⁻¹¹. In contrast, the ratio of N₂H⁺/H₂ = 2.2 × 10⁻¹⁰ at IRAS 16293 A/B is about a magnitude larger than the value of 3 × 10⁻¹¹ estimated by Stark et al. (2004), although it is in agreement with the value of 1.4 × 10⁻¹⁰ found by Jørgensen et al. (2004).

Also, the SiO column density at the position E1 is about an order of magnitude smaller than expected based on the results of Castets et al. (2001), who derived a value of 5.8 × 10¹³ cm⁻² for the total SiO column density toward this position. This deviation originates in the different analyzed transitions, as Castets et al. (2001) observe the SiO (2−1), SiO (3−2), and SiO (5−4) transitions, which have much smaller critical densities than the SiO (8−7) transition that we use in our estimations.

In general, most of the derived abundances at IRAS 16293 A/B and E are comparable with the mean values for Class 0 protostars and prestellar objects found in the sample of Jørgensen et al. (2004). Major differences to the sample of Jørgensen et al. (2004) exist for N₂H⁺, which is an order of magnitude less abundant in IRAS 16293 A/B then in other Class 0 protostars, where they observe a N₂H⁺ /H₂ ratio of 2.5 × 10⁻⁹. Also the observed DCO⁺ abundance of 9.8 × 10⁻¹¹ with respect to H₂ is almost an order of magnitude larger in 16293E than in the compared prestellar objects, that show abundances of 2.3 × 10⁻¹¹. The latter might be caused by the computation of our abundances relying on CO, while other authors used millimeter and submillimeter dust continuum emission to derive H₂ column densities. In cold regions such as 16293E the CO might be frozen out onto grains which would enhance molecular abundances relative to CO.

For Class 0 protostars in the Orion region specifically, Johnstone et al. (2003) derived H₂CO and CH₃OH column densities of 4 × 10¹³ cm⁻² and 2 × 10¹³ cm⁻² respectively. While these values are in agreement with values we derive for the outer layers of the protostar envelope in Appendix G, our derived single-layer LTE column densities in Table F.8 are about an order of magnitude larger. Therefore, the observed transitions are likely to include emission originating from the hot corino of the protostar.

Besides that, the ³²SO/³⁴SO ratio at IRAS 16293 A/B suggest an ³²S/³⁴S isotope ratio of 13.9, only about half of the value 24.4 ± 5.0 derived by Chin et al. (1996) which is commonly cited for comparable sources. Interestingly, this value is closer to the ³²S/³⁴S isotope ratio of ${16.3}_{-1.7}^{+2.1}$ observed toward the Galactic center (Humire et al. 2020). While the SO line profile appears mostly Gaussian, a reason for the deviation might be a high opacity of the considered SO transitions even without showing features of major self-absorption.

Finally, the deuterium fractionation as derived from the DCO⁺ /HCO⁺ ratio is enhanced at 16293E, where a fractionation of 0.05 is detected. The prestellar objects considered by Jørgensen et al. (2004) show on average deuterium fraction of 0.03. However, the high deuterium fraction at 16293E is in agreement with the findings of Lis et al. (2002), who observed a D/H ratio of about 10% in this source. In contrast, a smaller fractionation of 0.006 is observed at IRAS 16293 A/B, which is in agreement with the average deuterium fractionation of 0.007 in Class 0 protostars, as found in the sample of Jørgensen et al. (2004). The high deuterium fractionation in addition to the presence of several deuterated species make 16293E an interesting target for deeper integrating observations with the goal of investigating the variety of deuterated species in prestellar objects.

6 Summary

In this work, we have studied the large-scale spatial distribution of the molecular content in the environment of IRAS 16293–2422. This work fills an important information gap, since many previous studies toward IRAS 16293–2422 and its environment were based on pointed observations, therefore missing the complex spatial molecular morphology. The analysis is based on LAsMA and FLASH+ APEX observations covering frequency ranges between 277-375 GHz and 476–493 GHz respectively. The results of this study are summarized as follows:

• We identify a total of 144 transitions from 36 different molecular species in the observations. This is the first time that maps of the emission of such a large number of molecules have been presented for this region, which brings new information on the morphology of the molecular environment of IRAS 16293–2422.

• The produced large-scale molecular maps of L1689N show the two cores of IRAS 16293–2422 A/B and 16293E well separated and embedded in the molecular cloud, with extended envelopes surrounding them. Also, the emissions trace the two known outflows driven by the Class 0 protostar of IRAS 16293–2422 A. Outflow velocities up to 20 km s⁻¹ are detected in the northeastern outflow of IRAS 16293–2422 A.

• To easily compare and visualize the chemical wealth of each of the emission peaks, we have produced Venn diagrams, which allow us to directly analyze the molecular content. We find that the emission peaks related with the outflows show emission from other molecules besides the typical shock tracers (SiO, CO) such as HCN, CS, and NO. In addition, emission from deuterated species also is observed at these positions.

• A large-scale velocity gradient is observed for some species. Part of this gradient might be associated with the envelope around the A/B sources as it was suggested in previous studies, although we cannot disentangle the pure envelope rotation motion from the large-scale cloud kinematics.

• We have estimated new kinetic temperatures from para-H₂CO line ratios in L1689N at all considered positions, except for the prestellar core 16293E which shows only very weak H₂CO emission. The derived temperatures range from 31.3 K at the W2 position to 62.4 K as measured at IRAS 16293 A/B. These values are in agreement with previous H₂CO temperature estimations in the literature.

• We have computed new column densities for all the detected species in all emission peaks by producing synthetic spectra using the LTE radiative transfer CLASS module weeds. Upper limits were computed where appropriate and the corresponding abundances relative to CO were calculated.

• We have produced a complex non-LTE radiative transfer model with multiple physical components to reproduce the self-absorption H₂CO line profiles observed toward all sources, except for those were the wings produced by the outflows have an important contribution. We find that the results from this RADEX model are in agreement with the average values derived from our LTE modeling.

• The derived H₂CO temperatures were used in conjunction with additional line ratios of H₂CO lines with different J_up to estimate the H₂ volume densities in L1689N, resulting in values up to 5 × 10⁶ cm⁻³ at IRAS 16293–2422 A/B and values on the order of 1 × 10⁶ cm⁻³ for the other considered positions.

• We tested the proposed scenario of an interaction of one of the outflows arising from IRAS 16293–2422 A/B with the prestellar core 16293E. In this process, we were able to confirm the velocity offset between deuterated and nondeuterated species reported in literature and the trend of increasing velocities and line widths along the NE–SW axis across 16293E. These results, in combination with the high number of molecular species observed in the eastern cloud region, confirm the presence of an interaction between prestellar core and outflow.

• We show dust continuum maps obtained with Spitzer, Herschel, and JCMT telescopes. These maps reveal emission at the position W2 north of IRAS 16293 A/B, which appears to be unrelated to the outflow-structure of the protostars. A comparison between the continuum and molecular maps suggests that the origin of the emission at W2 might be due to a colder dust source embedded in the L1698N cloud, as major dust-continuum emission is observed for long wavelengths at this position. It is unclear if the embedded source is influenced by the outflows of IRAS 16293–2422 A.

In conclusion, combining all the results obtained in this work allows us to give a more complete and detailed view on the complex molecular distribution in the environment of IRAS 16293–2422. We have have been able to determine new values of physical parameters such as kinetic temperatures, column densities, and volume densities, and present a detailed description of the chemical content found at each of the cloud cores and emission peaks, including the interaction region between IRAS 16293–2422 and 16293E.

Acknowledgements

The authors are grateful to the anonymous referee for the useful comments and suggestions that helped improve this paper. We thank Arnaud Belloche for an early reading of the manuscript and for valuable suggestions. This publication is based on data acquired with the Atacama Pathfinder Experiment (APEX) under programme ID [M-0102.F-9519A-2018]. APEX is a collaboration between the Max-Planck-Institut für Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory. This work was partially funded by the Collaborative Research Council 956 "Conditions and impact of star formation" funded by the Deutsche Forschungsgemeinschaft (DFG).

Appendix A Observed frequency-bands

Appendix B Positions of the studied emission peaks

Appendix C Full list of the detected molecular line transitions

Table C.1 provides a full list of the spectral lines that were identified in our data. The list includes the name of the species, the quantum numbers of the transition, the frequency in MHz, the upper level energy, E_up, in K and the Einstein A coefficient, A_ij, in s⁻¹.

The second column describes the transition according to the values stated by the CDMS (Usually J or J_F for linear molecules, J_K for symmetric tops and $J_{K_a,K_c}$ for asymmetric tops). An exception is methanol, for which we give the symmetry state and, as subscript, K for the E-type and K for the A-type species. NO, CN, and C₂H are described as N_J,F, where we additionally give the parity of NO in the subscript.

The sixth column shows which of the receivers used in our observations cover the frequency of the respective line, where it should be noted that OTF-maps were only observed with LAsMA. The RMS noise is given in the seventh column of Table C.1, the values correspond to the RMS of the respective spectrum at A/B (in Kelvin) or the average RMS over each map (in K km s⁻¹) for mapped transitions indicated with an L in the sixth column. The analyzed LAsMA spectra and cubes have velocity resolutions between 0.14 and 0.10 km s⁻¹, which correspond to the lowest and highest observed frequencies of 277 and 375 GHz. The spectra taken with FLASH+ were smoothed to a velocity resolution of 0.08 km s⁻¹.

The last column lists which of the considered positions show line emission with a signal to noise ratio above three. In order to determine the signal to noise ratio of a spectral line, we averaged all spectra in a 10" radius around the respected position, motivated by the beam size of the telescope.

Appendix D Determination of the H₂CO excitation temperatures under LTE conditions

In the LTE approximation, the population between two energy levels n_u and n_l is described by the Boltzmann distribution $$\frac{n_{\text{u}}}{n_{\text{l}}}=\frac{g_{\text{u}}}{g_{\text{l}}}\exp \left(-\frac{E_{\text{u}}-E_{\text{l}}}{k{T}_{\text{ex}}}\right),$$(D.1)

where g_u and g_l are the degeneracies of the upper and lower levels, respectively, E_u and E_l are their energies and T_ex is the excitation temperature; k is the Boltzmann constant. The degeneracy g_u of an upper level can be computed as g_u = g_J × g_K × g_I, where g_J = 2 J + 1 is the rotational degeneracy and g_K = g_I = 1 applies as we consider H₂CO transitions in the same symmetry state (Mangum & Shirley 2015).

Assuming optically thin emission, the column density N_u of a transitions upper level is proportional to the integrated line intensity I of this transition (e.g., see Goldsmith & Langer 1999) and is given by $$N_{\text{u}}=\frac{8\pi k{v}^2}{h{c}^3{A}_{\text{ul}}}I,$$(D.2)

where ν is the frequency, c is the speed of light, A_ul is the spontaneous emission Einstein (A) coefficient of a respective transition. I is computed by integrating the main beam temperature T_mb over the considered velocities according to $$I={\displaystyle \int T_{\text{mb}}\text{\hspace{0.17em}dv}}.$$(D.3)

For this analysis, the ratio of column densities N_u/N_l as derived from the integrated main beam temperatures is used in Eq. D.1 instead of the relative level populations. The excitation temperature can therefore be written as $$T_{\text{ex}}=-\frac{E_{\text{u}}-E_{\text{l}}}{k}{\left[\ln \left(\frac{N_{\text{u}}}{N_{\text{l}}}\frac{g_{\text{l}}}{g_{\text{u}}}\right)\right]}^{-1}.$$(D.4)

The observed formaldehyde transitions allow the derivation of excitation temperatures for ortho- and para-H₂CO separately. Since only three transitions per species are available, these temperatures were obtained from the respective line ratios. For para-H₂CO was therefore possible to use the line ratios (4_2,2−3_2,1)/(4_0,4−3_0,3) and (4_2,3−3_2,2)/(4_0,4−3_0,3). In case of ortho-H₂CO, the (4_3,2−3_3,1) and (4_3,1−3_3,0) lines are blended. For this reason, the ratio [(4_3,2−3_3,1) + (4_3,1−3_3,0)]/[2 × (4_1,4−3_1,3)] was considered. Table D.1 lists the applied parameters for these calculations.

Appendix E Temperatures from RADEX non-LTE models

In order to test the validity of the assumption of LTE conditions, the correlation between calculated temperatures and line ratios from ortho- and para-H₂CO are compared with non-LTE models that were computed with the RADEX radiative transfer code (van der Tak et al. 2007) considering different H₂ volume densities n and using collision rates from Wiesenfeld & Faure (2013).

Figure E.1 shows several RADEX models for the kinetic temperature as function of line ratio alongside datapoints from our maps, which correspond to the individual pixels displayed in Fig. 6. The uncertainties on these points were calculated by Gaussian error propagation of the mean RMS from the associated maps.

The RADEX models were computed assuming a formaldehyde line width of 3.0 km s⁻¹ and considering only the cosmic microwave background for the background temperature of 2.73 K. Based on the computations in Sect. 4.2, a H₂CO column density of 3.8 × 10¹⁴cm⁻² was derived at at the position of IRAS 16293 A/B and 5.2 × 10¹³cm⁻² on average at the other positions. Since the weeds computations did not distinguish between ortho- and para-formaldehyde, the individual ortho- and para-column densities for the RADEX modeling are determined using a ratio of ortho/para = 3 (Kahane et al. 1984).

A comparison between the results from the LTE temperature derivation and the low column density non-LTE RADEX models in Fig. E.1b shows that the excitation temperatures computed from the data agree best with RADEX models for high volume densities, as the assumption of LTE holds in these cases. In contrast, the temperatures are underestimated for lower volume-densities by up to 20%. High column densities of H₂CO lead to an increase in optical depth for the observed transitions, which affects the (4_0,4−3_0,3) and (4_1,4−3_1,3) transitions more than the fainter upper energy state transitions. This effectively leads to an increase of the observed line ratios and therefore to an overes-timation of the derived excitation temperatures when assuming optically thin emission. In contrast, RADEX accounts for these optical depth effects, resulting in the temperature offset for high volume densities seen in Fig. E.1a. Due to this, temperatures derived at A/B are likely to be in better agreement with non-LTE temperatures than values derived at the other emission peaks.

Appendix F Derived parameters from the radiative transfer modeling with the CLASS module weeds

We performed a number of LTE radiative transfer models with the CLASS-weeds module in order the reproduce the observed molecular line profiles on each of the individual sources. The derived parameters from such models are shown in this Appendix. The derived values are the column density in cm⁻², the full width at half maximum (FWHM) in km s⁻¹ and the offset (δ) in km s⁻¹ with respect to the local standard of rest (LSR) velocity of the source (4 km s⁻¹). We have used the derived temperatures reported in Table 2 for each source. In addition to these values, a constant source size of 20" was assumed during the modeling and the velocities derived in Sect. 4.1 were applied.

In Fig. F.1 we show an example of a synthetic spectrum obtained with weeds for the H₂CS (8_1,7⁻7_1,6) transition detected in IRAS 16293–2422 A/B. Parameters used for the computation of simulated line profiles can be found in Table F.1-F.7. Based on a conservative estimate on the calibration uncertainty of about 20%, it can be expected that this uncertainty is transferred to the column density values.

An overview of the derived column densities and upper limits is presented in Table F.8, which states values derived from self-absorbed or distorted line profiles in parenthesis. Column densities that are based on optically thin isotopes are given in the last three rows of the table.

For comparing these values more easily with common literature, abundances of all detected transitions relative to CO are given in Table F.9. For this, the CO column density was derived from the C¹⁷O column density by assuming a ratio of ¹⁶O/¹⁷O = 1790 based the values suggested by Wilson & Rood (1994).

Appendix G Multilayer non-LTE radiative transfer model for H₂CO

In this part of the Appendix we describe in detail the multicomponent line modeling of H₂CO toward IRAS 16293–2422. The H₂CO temperatures and column densities derived with CLASS-weeds and RADEX in the main body of this work already describe very well most of the line profiles. Nevertheless, some of the H₂CO lines present self-absorption features. To model these lines properly, a more sophisticated model is needed. For that reason, we have used a non-LTE multilayer radiative transfer model to reproduce such line profiles. In this case, we have used the CASSIS⁹ software (Vastel et al. 2015) to perform such models. We have also used the collision rates for the ortho- and para-H₂CO - H₂ system derived by Wiesenfeld & Faure (2013).

In CASSIS, it is possible to create a model that consists of any number of physical components. Each of these physical components is defined by six parameters: Column density, temperature, linewidths (FWHM), LSR velocity, source size and H₂ volume density. CASSIS makes use of the Monte Carlo Markov Chain (MCMC) method (Hastings 1970) to explore the space of parameters and to find the best combination between them to reproduce the line profiles. This is done by means of a χ² minimization.

Not all the emission peaks could be modeled, since some of them (E1, E2, W1) present strong outflow wings in the line profiles (see Fig. 4). In contrast, the sources IRAS 16293–2422 A/B, 16293E, W2, and HE2 exhibit (approximately) Gaussian line profiles, although in some cases it is clear that there is contribution from the outflows.

Appendix G.1 IRAS 16293-2422 A/B

Some of the H₂CO lines (such as the 4_1,4−3_1,3 transition) in this source present self-absorption features (see Fig. G.1). Indeed, a number of molecules detected toward this source have been reported to have self-absorption features. In some cases, this feature is associated with the presence of an absorbing foreground cloud at a velocity of 4.2km s⁻¹ (e.g., Coutens et al. 2012; Bottinelli et al. 2014). To reproduce such absorption in our model, we need an extended and cold component. For our full radiative transfer model, we have considered a total of three physical components. The first component is compact and warm, associated with the hot corino in this source. The second component is associated with the warm envelope, while the third component correspond to the outermost extended and cold layer. The presence of these 3 physical components in IRAS 16293–2422 have been also identified in other molecules such as HNCO (Hernández-Gómez et al. 2019b). Some of the parameters were fixed during the modeling. First, the size of the hot corino has been fixed to 1" as the individual sources A and B have angular sizes below this limit. Since our observations cannot resolve between A and B, we assume that the bulk of the emission associated with the hot corino comes from a single warm source with this size. The temperature of the hot corino has been fixed to 90 K motivated by the results from the H₂CO line ratios presented in the main body of this work. For the warm envelope, we assume that the emission fills the beam, and therefore its size was fixed to 20". The temperature was also fixed to 50 K as preliminar models show that this values does not change significantly during the modeling. Since the extended component is assumed to be very extended, we fixed its size to 100". We have noted that is not easy to reproduce the absorption observed in the 4_1,4−3_1,3 transition, unless a very low temperature is considered. When having this parameter as free during the modeling, the absorption is easily destroyed during the modeling when increasing the temperature. For this reason, we have fixed the temperature of this cold component to 10 K. To determine the best choice for the velocity of each component, we have made some Gaussian fits to the line profiles and choose the best values based on them. In particular, we need to fix the velocities of all sources to reproduce the self-absorption feature correctly.

Finally, we have used the H₂ density radial profile determined by Crimier et al. (2010) to choose a suitable H₂ average value for each of the components depending on their sizes. The best parameters derived with CASSIS-RADEX are shown in Table G.1, while the produced synthetic spectrum is shown in Fig. G.1. As it can be seen, the non-LTE radiative transfer model is in good agreement with the observations. Interestingly, the velocity of the absorption is seen at 4.05 km s⁻¹, closer to the systemic velocity of the source.

A very useful utility of CASSIS-RADEX is that both forms of ortho-para H₂CO can be minimized at the same time during the modeling by using an extra parameter called iso, which corresponds to the ortho/para ratio. As a first approach, we fixed the ortho/para ratio to 3.0. However, we have noticed that changing this value between 1–3 does not affect the final result.

It is worth mentioning that Ceccarelli et al. (2000) studied the H₂CO emission in IRAS 16293–2422 and proposed a two components model that consisted of an evaporation region that corresponds to the hot corino in this source and a cold outer region. Their derived H₂CO abundances in these components are lower by a factor of about 2–4 compared with our model. These differences might be due to the different radial H₂ density profiles used for the modeling. In addition, we have use different number of components and different component sizes, which can also affect the computed abundance values.

Appendix G.2 16293E

The prestellar core 16293E exhibits narrow Gaussian H₂CO line profiles. A total of four formaldehyde transitions are detected toward this source (see Fig. G.2). Since these line profiles show a simpler morphology compared with IRAS 16293–2422 A/B, only one physical component was used during the modeling. Previous works have estimated the temperature in this source (12 K for the gas (Lis et al. 2002); 16–20 K for the dust (Stark et al. 2004); 11 K for the dust (Bacmann et al. 2016)). Also, the H₂ density has been recently estimated in different ranges: 3.3 × 10⁷ cm⁻³ (Pattle et al. 2015); (1.1–1.9)×10⁷ cm⁻³ (Lis et al. 2016); 1.4 × 10⁷ cm⁻³ (Bacmann et al. 2016). For our model, we have left the temperature to be a free parameter varying between 12–20 K. The H₂ density was varying also between 1 × 10⁶–1 × 10⁸ cm⁻³. The ortho/para ratio was fixed to 3.0. We verified that varying this ratio did not change the final results. The best parameters obtained with CASSI-RADEX are shown in Table G.2 and Fig. G.2, respectively. A temperature of about 12 K and an approximate H₂ density of 3 × 10⁷ cm⁻³ resulted from our model. We note that although these values give a better reduced χ² value during the minimization, there is some emission lacking for the 4_1,4−3_1,3 and 4_0,4−3_0,3 lines in our model.

Appendix G.3 W2

For the emission peak W2, four H₂CO transitions are detected in the data. In particular, the line at 281 GHz (4_1,4−3_1,3) presents self-absorption features. In this case two physical components are need to reproduce the data. It is worth mentioning that some of these lines are clearly contaminated by the outflows. Nevertheless, the radiative transfer model can satisfactorily reproduce most of the line profiles. The results from such modeling are presented in Table G.3 and Fig. G.3 respectively. Note that in the case of the self-absorption, some emission is lacking, even when the absorption is well reproduced. Compared with the temperatures derived from the H₂CO line ratios, we need much smaller values. Such discrepancy might be due to the degeneracy of our model. This is because we have seen that different combinations of the parameters can give similar results, although in most of these cases the self-absorption is destroyed. Therefore they are not representative of the physical conditions in this emission peak that seems to be associated with a cold dust source. While we require an extended and cold component as in the case of IRAS 16293–2422, lower values for the H₂ density in this component are needed.

Appendix G.4 HE2

The morphology of the H₂CO line profiles in the emission peak HE2 is not completely Gaussian. Nevertheless, they can still be modeled with CASSIS-RADEX. Four transitions are clearly detected in this source and since there is no clear self-absorption detection, we used only one physical component to model the lines. The results obtained with CASSIS-RADEX are summarized in Table G.4 and Fig. G.4 respectively.

Appendix H Venn diagrams for the occurrence of molecules

In this section we show Venn diagrams to visually compare the occurrence of different molecules on each of the emission peaks E, E1. E2, W1, W2, and HE2.

Appendix I Line profiles of deuterated species in the vicinity of 16293E

We extracted spectra of three positions where we see clear emission for all deuterated species that peak near 16293E: H₂D⁺, DCO⁺, DNC, N₂D⁺, NHD₂, and NH₂D. The exact coordinates are the following: NE: (16^h32^m30.6^s, -24°28'46.7"), NW: (16^h32^m28.9^s, -24°28'46.7"), SW: (16^h32^m28.9^s, -24°29'13.7"). In order to test the existence of a trend in line width and velocity shift of these lines, we conducted a Gaussian fit to these line profiles. The resulting fits to the spectra are shown in Figs. I.5, the derived peak velocity and FWHM values are shown in Table I.1.

Appendix J Integrated intensity maps for the observed molecules

In this part of the Appendix we show all the velocity-integrated intensity (moment 0) maps for all the detected molecules in the observations. The figures are sorted based on the frequencies of the corresponding transitions (see Table C.1).

The name of the corresponding molecule and its transition are indicated in red color in the upper left side of each map, while the velocity range considered for the integration is given the lower right corner. In case of blended lines the integration range was chosen such that all lines are covered, the contributing transitions are given in the caption of the corresponding figures. The beam is shown in the lower left corner of the maps. The color bar in the right side has units of K km s⁻¹. The contour levels are drawn in black color at 3σ, in steps of 2σ between 4σ and 10σ and in steps of 10σ afterwards. The markers show the position of the emission peaks introduced in Sect. 3. Given that some methanol transitions exhibit similar spatial morphologies, we have not mapped them all (see Table C.1). The intensity is given in terms of antenna temperatures $T_A^{*}$. To convert to brightness temperature T_MB, a multiplicative factor of 1/0.7 needs to be considered.

Some transitions only show weak emission, such that the signal to noise ratio of a single pixel is not sufficient for a significant detection. In these cases, we averaged the spectra in a 10" radius around the respective emission peaks to confirm a detection. The averaged spectra are shown next to the associated maps. Some of these spectra are smoothed to a lower velocity resolution in order to improve the signal to noise ratio for the visualization.

References

---

All Tables

All Figures