© ESO 2021

1 Introduction

In the canonical picture of star formation, a dense core accretes material from a surrounding envelope, and, because of the conservation of angular momentum due to the accretion, a bipolar jet-outflow emerges. Prior to Class 0, a prestellar core provides the initial conditions for accretion and outflow-driving, and hence it is an important yet elusive initial component to study (Dunham et al. 2014; Padoan et al. 2014). The accretion and jet-launching are the active processes during the “protostellar” (Class 0-I) phase (Arce & Sargent 2006), as a protostellar disk also begins to form. Once the envelope dissipates and the outflow subsides (Class II), the disk is exposed around the young stellar object (YSO).

Several indicators capable of providing information about the evolutionary state of pre- and protostellar sources exist in the literature. The most used are the bolometric temperature L_bol, and the ratios between the masses of the disk and the envelope (M_disk/M_env) and between the bolometric and submillimeter luminosity (L_bol /L_submm). Furthermore, for sources with associated outflows, it is possible to assess their evolution through parameters such as momentum, energy, mass-loss rate, and momentum flux.

The Lupus complex, located near the Scorpious OB2 association, is known for its nearby (~150 pc) clusters of low-mass, premain sequence stars (Comerón 2008; Knude 2009). Within the 335° < l < 348° and 0° < b < 25° area, the Lupus complex was first found to have four subgroups named Lupus I, II, III, and IV. Around 70 T Tauri stars are found around the Lupus complex with estimated ages of about 1–3 Myr (e.g., Comerón 2008; Ansdell et al. 2016, 2018; Alcalá et al. 2015). Lupus were studied in ¹²CO, ¹³CO, and C¹⁸O in order to map the molecular gas (Hara et al. 1999; Vilas-Boas et al. 2000; Tachihara et al. 2001; Tothill et al. 2009). These maps have shown that Lupus is actually divided into nine subgroups, where only Lupus I, III, and IV present evidence of ongoing star formation. The C¹⁸O data were used to identify about ten cores in the most massive one, Lupus I (with a mass of ~1200 M_⊙), with column densities N(C¹⁸O) = (5–10) × 10¹⁴ cm⁻² indicating potential sites of star formation. Lupus III has a mass of about 300 M_⊙ and it hosts a rich cluster of T-Tauri stars, but not much star formation activity. Lupus IV is the third cloud of the complex to show evidence of star formation activity, including nine C¹⁸O dense cores with column densities of N(C¹⁸O) = (4–10) × 10¹⁴ cm⁻² and three H¹³CO⁺ cores (Tachihara et al. 2007).

The SOLA (“The Soul of Lupus with ALMA”, Saito et al. 2015) collaboration is composed by scientists with technical expertise in Atacama Large Millimeter Array (ALMA) and in infrared and optical techniques. The SOLA team is dedicated to a multiwavelength study of Lupus, including a catalog of all prestellar objects and YSOs (hereafter, “SOLA catalog”) which will be published elsewhere. Previously, using the SOLA catalog, 37 protostellar sources were targeted to search for binaries in Lupus in Cycle 2 (PI: Saito, 2013.1.00474.S). Among these sources, we have quantitative information about their submillimeter properties. The spectral setup of those Cycle 2 observations covered continuum, and (simultaneously) the ¹²CO(3–2)and HCO⁺(3–2) lines for evidence of outflow and infall at a single pointing for each source. The data presented in this work, observed during the ALMA Cycle 4, complements and extends the Cycle 2 study: we now aim to study the cores, envelopes, disks, and the full extent of outflows from sources at different evolutionary phases from prestellar through protostellar in Lupus. In this way, we want to characterize all these phenomena associated with the star formation process in a set of isolated sources at different stages of evolution in a similar nascent environment, and infer their physical and kinematic properties in relation to their stage of evolution. Several variables that can influence the evolutionary tracers challenging this task. To the intrinsic difficulty in measuring some indicators (such as choosing the distance for measuring the outflows opening angles), we must add the possible effects of precession and-or outflow interaction with the environment, and the uncertainty in the interpretation of the results depending on the outflow model adopted. However, we believe that the possibility of studying several sources in the same environment gives us the chance to evaluate how consistent the different indicators are. In addition, we have high-quality molecular outflow data that could shed light on the evolutionary state of the particular sources. In this way, we intend to evaluate the sources and test whether there is consistency between the evolutionary tracers commonly used.

Here we present the results of the study of seven sources. Six of them were selected from the SOLA catalog and one more was included since it was cataloged as a young Lupus protostar in previous studies. Although our sample is not statistically significant, we attempt to trace an evolutionary sequence in the sources of our sample by evaluating various evolutionary indicators from the literature. Spectral line and continuum observations from ALMA enable investigation of the molecular dust and gas environment surrounding the selected sources, on scales of 10³ AU. In Sect. 2 we describe the sample selection, and in Sect. 3 we describe the observations and explain the calibration and imaging procedure. We organized the Results section in three subsections: continuum emission (Sect. 4.1), single dish line emission (Sect. 4.2), and outflow emission (Sect. 4.3), where we describe the outflow morphology, physical parameters and the kinematics. In Sect. 5 we split the discussion in two subsections: in Sect. 5.1 we refer to the individual outflows and in Sect. 5.2 we discuss the characteristics of the source sample (spectral energy distribution, chemical, and outflow evolution). Finally we include an appendix, where we present a brief overview of each source, mention the particularities of some sources, and show the complete velocity cubes.

2 The sample

The SOLA catalog enacts a clump-finding algorithm and spectral energy distribution (SED) fitting to classify all sources in Lupus as prestellar or Class 0-III. The basis of this catalog is a set of millimeter and submillimeter images. Lupus I, II, III and IV clouds were observed at 1.1 mm with the AzTEC bolometer mounted on the ASTE telescope (Tsukagoshi et al. 2011), and the Lupus III cloud was also observed at 870 μm with the LABOCA bolometer mounted on the APEX telescope. In total these observations cover 7.5 deg² on the sky. A total of 251 sources were found in AzTEC and LABOCA images, utilizing CLUMPFIND as the structure search algorithm (Williams et al. 1994)¹. To classify these sources, we searched for optical to centimeter counterparts, estimating their spectral energy distributions (SEDs). We used cataloges and-or images from WFI, DENIS, 2MASS, AKARI, ATCA, SMA, SEST, WISE, IRAS, Spitzer, and Herschel telescopes, and correlated the positions at different wavelengths considering pointing errors from all telescopes. The AzTECand LABOCA objects are classified as 189 starless, four Class 0, seven Class I, 49 Class II, and two Class III objects.

We chose seven individual sources throughout the Lupus I, III, IV, and VI molecular clusters that span protostellar evolutionary phases in similar nascent environments: two each from the prestellar, Class 0, and Class I evolutionary phases, and an additional late-Class I source (see Table 1). These sources were chosen based on the following source criteria: (1) we selected sources based on SED shape. (2) For protostar candidate sources, we selected sources that are considered as single (non-multiple) sources. This was confirmed with the ALMA Cycle 2 Band 7 continuum observations with the resolution of 0.′′ 2, or 31 AU (PI: Saito, 2013.1.00474.S). This criterion results in five targets (two in each of the protostar categories Class 0 and Class I, and one classified as late-Class I), a complete sample of nonbinary protostars in Lupus. (3) Finally, we chose the two prestellar candidates that have a single, strong peak based on the concentration factor, C = F_peak/F_3σ, using AzTEC/ASTE 1.1 mm observations. These have a density greater than the critical density, to imply they are not transient clumps. A brief description of each source is provided in Appendix A.

Observing sources within Lupus allows us to assume similar environmental conditions (i.e., interstellar radiation fields, gas temperatures, chemical composition) and observing consistencies (i.e., beamsize, distance, spatial resolution, sensitivity). Recent calculations of Lupus I and III distances based on Gaia DR2 data have revealed that these clouds are located at 153.35 ± 4.64 and 154.75 ± 9.59 pc, respectively (Santamaría-Miranda, accepted), values which are in agreement with the distances derived by Sanchis et al. (2020). This, along with the distances reported by Zucker et al. (2020), leads us to adopt a mean distance of 150 pc for all the clouds in the complex. In any case, our conclusions will not be dramatically affected by this choice of distance, since the consistency will always be better than choosing sources from entirely different regions.

3 Data

This work is mainly based on dedicated observations using the total power and 7m array of ALMA, including single-field and mosaic modes. We have complemented these observations with data from the ALMA Science Archive; in order to study compact sources, we have supported these dedicated observations with archival data of continuum emission from the ALMA 12m array, and high-resolution molecular line data for those sources that drive molecular outflows. We have also included infrared data associated with the compact sources over a wide range of wavelengths.

In Sects. 3.1–3.4 we describe these four groups of data, respectively (our dedicated ALMA observations, archival continuum data, archival molecular data, and infrared data). In Sect. 3.5 we describe how the imaging and combination processes were carried out for each source.

3.1 SOLA observations

Using the ALMA 7m array and TP antennas, we observed the seven selected sources in Band 6 (1.3 mm observations. Project code: 2016.1.01141.S, P.I: Takahashi, Satoko). Three of the sources were observed with a single field, and four using small mosaics (see Table 1). The data were taken on 12-20 Dec. 2016; 15-30 Apr. 2017; and 19-21 May 2017. The observationswere made with baselines between nine and 45 m, from which we achieved angular resolutions of ~8′′ (1240 AU).

The correlator configuration included four spectral windows with 4096 channels each. The two spectral windows dedicated to observe ¹²CO(J = 2–1) (ν_rest = 230.538 GHz) and SiO(5–4) (ν_rest = 217.10498 GHz) covered 500 MHz each at a frequency resolution of 122.1 kHz (velocity resolution 0.16 km s⁻¹); the two spectral windows dedicated to observe N₂ D⁺(3–2) (ν_rest = 231.321828 GHz) and C¹⁸O(J = 2–1) (ν_rest = 219.560358 GHz) covered 125MHz each at frequency resolution of 30.5 kHz (velocity resolution 0.04 km s⁻¹). With respect to the weather conditions, the precipitable water vapour (PWV) values ranged between 0.25 and 2.16 mm during the total power observations, and approximately 1.8 mm for the interferometric observations. Calibration of the raw visibility data was performed using the standard reduction script for the Cycle 4 data provided by the ALMA Observatory. This pipeline ran within the Common Astronomical Software Application (CASA 5.1.0, McMullin et al. 2007) environment. The integration time and the calibrators used to correct for instrumental and atmospheric disturbances (flux, phase, and bandpass) are listed in Table 1.

In the total power observations, the frequency setup was exactly the same as in the case of the interferometric observations. In this work we directly use the total power images delivered by ALMA.

The calibrated interferometric data were cleaned in CASA to produce continuum images. The spectral line cubes were produced by subtracting the continuum and applying a standard cleaning with primary beam correction. Imaging and data combination methods for the data are explained in greater detail in Sect. 3.5.

3.2 ALMA archive continuum data

Aiming for a better understanding of the finer structure of the continuum sources, we complemented the observations with data products from the ALMA archive. Specifically, we searched the archive for all continuum data available for our sample, and retrieved high-angular resolution (taken with the so-called ALMA 12m array) continuum images of five of the sources: IRAS 15398-3359, IRAS 16059-3857, Merin 28, Sz 102, and J160115-41523. The observational parameters, such us frequency band, central frequency, number of antennas, shortest and longest baseline, time on-source, and flux, phase and bandpass calibrators are compiled in Table 2.

For some of these sources, we self-calibrated and cleaned the data, and we found that the quality of the continuum images (signal-to-noise ratio; S/N) did not improve considerably with respect to the archival product images (since the continuum sources are quite compact and the on-source exposure times are fairly short). Moreover, the small flux variations from nonself-calibrated and self-calibrated continuum emission was always within the 10% uncertainty of the flux measurements. Hence in this case we decided to use the pipeline delivered images.

3.3 ALMA archive molecular data

We also retrieved high-angular resolution data cubes from the ALMA archive for the sources presenting outflow activity (IRAS 15398-3359, IRAS 16059-3857, and J160115-41523).

The 12m array observations toward IRAS 15398-3359 (2013.1.00879.S) were carried out during ALMA Cycle 2 on 2014 April 30 with 34 antennas, and on 2014 May 19 and June 6 with 36 antennas. Since the extent of the outflow is small, these observationswere only a single pointing, enough to observe the entire outflow. The spectral lines observed were ¹²CO(2–1), C¹⁸O(2–1), and SO(J_N = 6₅, 5₄).

Additional data from IRAS 16059-3857 (also named Lupus 3 MMS) were taken using the 7m array, the 12m array and the total power (TP) from ALMA over the whole outflow extent (2017.1.00019.S). The total power data were observed during the period of 2018 April 7 – 2018 June 9. The 7m data were observed during the period 2018 April 2 – 2018 May 11, and 31 pointings were used to make a mosaic map. The 12m data were observed on 2018 June 4 and 5 and 105 pointings were used to make the mosaic map.

The high-angular resolution observations associated with J160115-41523 (2015.1.01510.S) consist of ¹²CO(2–1) observationswith a single pointing toward the central source, and were taken on 2016 April 2. The Band 6 continuum data used in this work does not correspond to this project, since we found in the archive this source observed in Band 6 with an angular resolution similar to that found in Band 3 and 7 (see Sect. 3.2).

The observational parameters, such as configuration, achieved angular resolution, number of antennas, shortest and longest baseline, time on-source, and flux, phase and bandpass calibrators are compiled in Table 3.

3.4 IR data

In order to fit the spectral energy distribution (SED) of the sample’s sources we use data from IR to optical wavelengths (Lopez, SOLA Collaboration et al., 2020 in prep.), along with the fluxes we measure in ALMA Band 6. A table detailing the data (band, wavelengths, FWHM, and fluxes) has been included in Appendix B.

The association between ALMA and optical-IR counterparts was first based on a visual inspection. In the case that the association was not clear for any source, we compare the astrometric accuracy/pointing error of each telescope with the distance to the source position. Finally, we checked the angular resolution (FWHM) of each instrument if there was any doubt about the association. We followed the same procedure by Santamaría-Miranda et al. (accepted), where a more extensive list of telescopes, instruments, and errors can be found.

3.5 Imaging and data combination

In this section we describe the imaging and data combination procedure we carried out for the different data sets, which includes observationsfrom the SOLA project (hereafter 7m and total power observations, see Sect. 3.1), and, for some sources, the data obtained from the ALMA archive (hereafter 12m observations) described in Sect. 3.3.

Due to the physical size of the dishes and the minimum separation (shortest baseline) between any two, interferometers never sample the central region of the(u, v)-plane. The incompleteness of the (u, v)-coverage at low spatial frequencies (know as the short-spacing problem) makes interferometers insensitive to the emission on large angular scales. The effect of missing short-spacings is negligible for objects that are small in comparison to the extent of the primary beam, but for objects with large extended structures, the lack of sensitivity toward low spatial frequencies is a severe shortcoming. To overcome the short-spacing problem the data need to be combined with those of a single-dish telescope, which are capable of measuring the total power. For this reason we have combined interferometric and single-dish observations from those data sets where we encounter large-scale extended emission (outflows) in single-dish observations.

For AzTEC-lup1-2, AzTEC-lup3-5, Sz 102, and Merin 28, with only SOLA project data, the 7m calibrated visibilities were Fourier transformed and cleaned with the CASA task tclean. We set the Briggs weighting parameter robust = 0.5 for the continuum and ¹²CO(J = 2–1) images as a compromise between angular resolution and S/N (typical beam of 7.′′ 0 × 4.′′ 5, PA = 86°), and we set robust = 2 – which ensures higher S/N – for the C¹⁸O(J = 2–1), N₂ D⁺(3–2) and SiO(4–3) images, where the emission is fainter (typical beam of 8.′′ 0 × 4.′′8, PA = 82°). The rms noise level in the continuum images is around 2–4 mJy beam⁻¹, and ~100–150 mJy beam⁻¹ in the 0.17 km s⁻¹ channels of the line velocity cubes. The images produced from the total power single dish have an angular resolution of 29′′. For these sources the total power data only shows cloud emission, therefore we did not perform the combination with the 7m array data.

In the case of IRAS 15398-3359 the 7m array visibilities (¹²CO and C¹⁸O) were concatenated with the 12m visibilities from the archive data and cleaned together using the CASA task tclean, setting a Briggs weighting parameter robust = 0.5 for the ¹²CO(J = 2—1) and robust = 2 for the C¹⁸O(2-1). The FWHM of the clean beam after imaging the concatenated 7m+12m data is 0.′′ 6 × 0.′′5, with a position angle 54. °0 and a rms noise level of 7 mJy beam⁻¹ in a spectral channel (velocity resolution 0.16 km s⁻¹). The interferometric (7m and 12m) and total power images were combined using the CASA task feather. This algorithm converts each image to the gridded visibility plane, combines them, and reconverts them into a combined image. The combined image preserves the detailed structures revealed by the interferometric data while also recovering the extended emission filtered by them. We also note that the high-resolution archival data of this source includes the SO line, which was not observed under the SOLA project. We re-imaged these observations and obtained data cubes with an angular resolution of 0.′′ 6 × 0.′′5 with a position angle 36. °5 and a velocity resolution of 0.16 km s⁻¹, reaching an rms of 6.7 mJy beam⁻¹.

Combining the interferometric data (7m and 12m) available for IRAS 16059-3857 we have obtained, after imaging with a Briggs weighting parameter robust = 0.5, a ¹²CO(2–1) data cube with an angular resolution of 1.′′ 7 × 1.′′1, a beam position angle 86. °5, and a S/N of 13 mJy beam⁻¹ at a spectral channel width of 0.32 km s⁻¹. Since we have total power data of the full outflow associated with this source, we combined it with the interferometric data via the CASA task feather, obtaining images that recover the extended emission flux while preserving the high angular resolution.

In the case of J160115-41523 we do not have the complete outflow mapped by the 12m array, but only one single pointing ¹²CO(2–1) observation toward the central source. The angular resolution of this high resolution data is 0.′′ 9 × 0.′′8 with a beam position angle 69. °1 and a rms noise level of15 mJy beam⁻¹ at a velocity resolution of0.16 km s⁻¹, obtained using a value of Briggs weighting parameter robust = 0.5. The 7m observations of the full outflow were cleaned, setting the Briggs weighting parameter robust = 0.5 for the continuum and ¹²CO(J = 2–1) images, and combined with the total power data by using the CASA task feather, recovering the full outflow emission. The angular resolution of this combined images is 7.′′ 5 × 4.′′7 (PA = −80°) and the velocity resolution 0.17 km s⁻¹, achieving an rms 75 mJy beam⁻¹.

4 Results

4.1 Continuum emission

Millimeter continuum thermal dust emission was detected with ALMA toward five sources in our sample. As shown in Table 4, the five sources detected as point-like in the dedicated 7m array data are resolved in high resolution ALMA archive data, except IRAS 15398-3359 which is detected as a point-like source within low-level extended emission, even in 12m ALMA data.

In the upper panels of Fig. 1 we show the continuum emission associated with each source. Compact emission from AzTEC-lup1-2 and AzTEC-lup3-5 is not detected by ALMA, so we display the AzTEC/ASTE 1.1 mm continuum emission which traces the clumps where the sources are embedded. In IRAS 15398-3359, IRAS 16059-3857, J160115-41523, Sz 102 and Merin 28, the ALMA Band 7 continuum emission traces the disks, while extended emission is also detected associated with IRAS 15398-3359. These high resolution images, show that the sizes of the sources range between 0.′′ 2 and 0.′′ 8 (~40 to 120 AU).

In order to estimate the continuum emission parameters (major and minor axis sizes, position angles, integrated fluxes, and peak fluxes)and their uncertainties, we fit 2D-Gaussians by using the CASA task imfit. Although this model may not be the best one for all our sources, we found that the maximum values of the residual maps are below 4σ, so they are a good approximation. Assuming that the sources inherently have a circular shape, we estimate that the inclination angles of the disks are between 50° and 70°². Table 5 lists for each source the Gaussian fitting parameters: position, deconvolved FWHM of the major and minor axes, position angle, size of the source, and inclination. In the case of IRAS 15398-3359, since the source is marginally resolved we can not be sure if the inclination calculated from the Gaussian fit is the real one, however it is in good agreement with that calculated previously by other authors based on H₂CO(5₁₅–4₁₄) ALMA observations, assuming that the outflow cavity has a parabolic shape and its velocity is proportional to the distance to the protostar (Oya et al. 2014). We propose that the continuum emission detected in high resolution images corresponds mostly to the disks (typical disk sizes <150 pc), and the emission detected in the ALMA 7m observations is emitted by the dust envelope+disk, as the angular resolution of these observations (~7′′) at the Lupusdistance corresponds to ~1000 AU.

In the central panel of Fig. 1 we show the Spectral Energy Distributions (SEDs) of the seven sources. Filled red circles correspond to the IR fluxes described in Sect. 3.4, and to mm measurements taken from the ALMA 7m continuum images. Open blue circles represent the fluxes measured in the high-resolution ALMA images. The fluxes were fitted by a modified black-body model, shown with a black dashed line.

We estimate the dust temperature T_d from fitting the IR SED of the sources (see Fig. 1). To obtain the temperature of the cold dust, we fit a black-body model to the emission at millimeter and FIR wavelengths. The fit was optimized using the Python scipy (Virtanen et al. 2020) task curve_fit, which uses a nonlinear least squares algorithm to fit a function to data. The errors come from the uncertainty of the black body fit, and are obtained from the covariance matrix provided by the task. Usually this emission is fitted by a gray-body model with a variable submillimeter slope, however due to the lack of information in the submillimeter spectrum and the low resolution of some infrared data we opted to fit a simple black body model.

In fitting, we do not include high angular resolution millimeter emission (open blue circles) as the fluxes could be affected by emission filtering, and in most cases are likely lower limits. A black body fits well with the data in the FIR for AzTEC-lup1-2, AzTEC-lup3-5, IRAS 15398-3359, IRAS 16059-3857, and J160115-41523, and the derived temperatures show the presence of cold dust (between 8.3 and 32.0 K), while Sz 102 and Merin 28 fit with a black body at a higher temperature (144.3 and 275.0 K, respectively). However, for the latter source, the lack of FIR data makes this value unreliable. The estimated T_d values for each source are listed in Table 6.

The SED corresponding to AzTEC-lup1-2, which only includes FIR data, is typical of prestellar cores, and emission appears to pertain to isothermal dust with ~10 K, peaking at λ > 100 μm. The protostar seems to be completely covered by gas and dust and is obscured with a large optical depth by the dust envelope. No conclusion can be reached from the stellar-black body radiation. The source has no compact emission and only extended millimetric emission is detected, so we are probably observing the moment when the dust cloud starts to collapse and the disk has not yet formed.

In the case of AzTEC-lup3-5, we do not have a mm detection either; however, this source emits over a wide wavelength range and the SED has the appearance of a less embedded object. On the other hand, the black body fit yields a dust temperaturevalue close to 10 K, so it would be one of the coldest sources in our sample. As we cannot see the disk because there is no compact mm emission we consider two possible scenarios: (1) that we are actually seeing two objects in the line of sight (a prestellar core and an extinct reddened star), or (2) indeed it is a protostar that has already collapsed and has not yet formed the disk around it, so we do not see the compact emission (this would be rare).

The SEDs we obtained for IRAS 15398-3359, IRAS 16059-3857, and J160115-41523 correspond to Class 0-I objects, which are dominated by infalling material from the parent molecular cloud to the central object with the eventual presence of material flowing from the object toward the surroundings (outflows). When part of the envelope is dissipated by the outflows, the object becomes detectable in the near infrared. This scenario is in accordance with IRAS 15398-3359, IRAS 16059-3857, and J160115-41523, which are sources that have outflows associated. The FIR emission in these sources can be fit by black bodies with temperatures ~30 K.

In contrast to any other source of our sample, the SEDs of Sz 102 and Merin 28 show peak flux at a wavelength near to 22 μm, rather than in the FIR. In both cases, also in contrast to the other sources, the mm points are above the black body curve, so the fitting is probably not appropriate, especially in the case of Merin 28, where no IR observations are available at wavelengths larger than 22 μm.

Furthermore, we estimate the mass associated with the submillimeter continuum emission. Assuming isothermal dust emission, well-mixed gas and dust, and optically thin emission, the dust masses are given by $$M_{\text{d}}=\frac{S_{\nu }{D}^2}{B_{\nu }(T_{\text{d}}){\kappa }_{\nu }},$$(1)

where S_ν is the continuum flux density, D is the distance to the source, B_ν is the Planck function, T_d is the dust temperature, and κ_ν is the mass dust opacity coefficient (Hildebrand 1983). The dust opacity was estimated from κ_ν = 0.1(ν/10¹²Hz)^β cm⁻² gr⁻¹ (e.g., Beckwith & Sargent 1991). The value of β was estimated by fitting the flux densities at different millimetric wavelengths with a power law, F_ν ∝ ν^2+β, based on measuresof the present observations and the literature. Since we are assuming optically thin dust emission and not taking into account possible scattering effects (Zhu et al. 2019), the derived masses should be considered as lower limits. We derive the dust temperature by fitting the SED of every source. For this purpose we also use ALMA archival data from the sources in different spectral bands (see Table 2). The measured fluxes are listed in Table 6. The error in β was derived from the fit assuming a 10% uncertainty in the ALMA observations due to calibration errors. For sources with insufficient measurements to make the fit, we adopted a value of β = 1.8 ± 0.2 a typical opacity spectral index value for the ISM at submillimetric wavelengths (Draine 2006). In all cases, the β values obtained from the fitting (see Table 6) range between 0 and 1, according to the values of β reported by other authors (Draine 2006; Beckwith & Sargent 1991; Ribas et al. 2017; Ansdell et al. 2018).

The values of M_d, listed in the last three columns of the Table 6 (expressed in jovian masses), were calculated from fluxes measured in the dedicated observations and in high resolution continuum archive images, showing similar results, which are within the errors. The errors in the dust masses were calculated by propagating the errors in fluxes, dust temperatures, and β. Since for both AzTEC-lup1-2 and AzTEC-lup3-5 we cannot detect compact emission, the estimated dust mass values were calculated by using a flux equal to three times the rms and therefore represent an upper limit. For three of the sources where compact emission was detected (IRAS 15398-3359, IRAS 16059-3857, and J160115-41523) we note that the dust mass calculated from the high resolution observations is lower than that calculated from the 7m observations. This could have two explanations: the flux loss produced by the extended emission filtering, or to the fact that at high resolution we are measuring the mass of the dust associated to the disk and at low resolution the mass associated with the envelope and the disk. This second hypothesis seems to work well for relatively large samples (e.g., Tobin et al. 2020). We make the caveat that unresolved edge-on disks with large envelopes can be confused with disks with high inclination and the 12m observations may include emission from these envelopes. In the case the case of Sz 102 and Merin 28 (where the masses calculated from the measured fluxes with the 7m and 12m arrays are practically the same) we consider that the envelope mass would be negligible, and we are only measuring the disk mass. Particularly in the case of Sz102, other studies conducted with ALMA observations have shown that no traces of a massive envelope are observed (Louvet et al. 2016).

Fig. 1 Upper panels: continuum emission toward AzTEC-lup1-2, AzTEC-lup3-5 (ASTE, 1.1 mm), IRAS 15398-3359 and IRAS 16059-3857 (ALMA, 0.8 mm). Contours: 4, 10, 20 and 30 σ. Middle row panels: spectral energy distribution. A black-body could be fitted considering only the FIR and mm emission (black dashed line). The open blue dots corresponding to emission detected in the high resolution mm images are not included in the fit. The dust temperature derived from the fitting is shown in the top right corner. Bottom panels: line emission detected toward the sources at the position of the protostar.

Fig. 1 (Cont.) Same as Figure 1 for J160115-41523, Sz 102 and Merin 28.

4.2 Line emission

The lower panels of Fig. 1 show the line emission associated with the continuum sources described inSect. 4.1. The spectra shown were obtained from the total power data. The interferometric data alone, with the achieved sensitivity, show the detection of only a few lines. Although it is not possibleto spatially resolve the line sources in the total power data, we extract the spectra within a 2″ box at the position where the continuum sources are located. This procedure is sufficient to achieve our goal, which is simply to identify the molecules present in the source envelopes. The Gaussian fit performed to the lines is very simple, and although it is not perfect for lines with more complex structure (double peaked or asymmetric profiles), it allows us to determine a systemic velocity andan overall line width. The fit was optimized using the python task curve_fit, and the parameters are shown in Table 7 along with their uncertainties.

It is possible to see that the chemical complexity varies from one source to another. All sources exhibit ¹²CO(2–1) and C¹⁸O(2–1) emission at velocities ranging 4–5 km s⁻¹, but AzTEC-lup1-2, IRAS 15398-3359, and IRAS 16059-3857 are chemically more complex: all three show N₂ D⁺(3–2) and DCN(3–2) emission, and IRAS 15398-3359 presents weak SiO(5–4) and CH₃OH(6(1,5)–7(2,6)) lines, the latter one at a different velocity than the rest (at –0.4km s⁻¹). In all sources the CO line width is many times wider than the other lines. The systemic velocities were derived from the peak velocities of the molecules which trace highest densities (C¹⁸O, N₂ D⁺, and DCN), considering an error equal to half the channel width. These values coincide with the velocity of the absorption dips in the optically thick lines.

Both AzTEC-lup1-2 and AzTEC-lup3-5 spectra show multiple components. Since these sources have not been spatially resolved using infrared observations nor the millimeter emission from ASTE, and they do not show any compact emission in the ALMA data, we could not rule out the possibility that there could be more than one source projected in the plane of the sky toward these positions. However, the AzTEC-lup1-2 spectrum is chemically more complex than that of AzTEC-lup3-5, showing N₂ D⁺ and DCN emission. These molecules are high density tracers: N₂ D⁺ shows the presenceof cold gas and sometimes traces quiescent clumps, and DCN survives at low temperatures where other molecules freeze out (see discussion in Sect. 5.2.2).

The IRAS 15398-3359 spectrum shows prominent ¹²CO wings, which are a typical feature of outflow sources. This source is the only one that shows SiO emission (with v_peak = 4.9 km s⁻¹ and Δ v = 1.4 km s⁻¹) that is usually associated with young energetic outflows. Furthermore, it presents a methanol line at a different velocity (v = –0.4 km s⁻¹), which does not trace the velocity of the cloud, and it is probably generated in a shock.

The spectrum in IRAS 16059-3857 shows an asymmetric double peak of ¹²CO with large broadening wings, which are probably associated with its outflow.

The three remaining sources (J160115-41523, Sz 102, and Merin 28) show only ¹²CO and C¹⁸O emission, with multiple peaks and asymmetric lines. We note that the two C¹⁸O emission peaks in J160115-41523 have similar amplitudes, while those detected in Sz 102 are asymmetric. Merin 28 has a single C¹⁸O peak and the maximum of ¹²CO is blue-shifted, with v_peak at 3.3 km s⁻¹ compared to the v_lsr at 4.2 km s⁻¹.

4.3 Molecular outflows

In this section we introduce the outflows associated with three of the sources of our sample (IRAS 15398-3359, IRAS 16059-3857, and J160115-41523), by showing the CO(2–1) integrated intensity emission from their red and blue lobes (Fig. 2) which delineate and describe well their overall morphology. We estimate the physical parameters of each outflow (Table 8), and study their kinematics via the moment 0, 1 and 2 maps (Figs. 4, 6, and 8) and the position-velocity (PV) diagrams (Figs. 5, 7, and 9). We present the velocity cubes of all three outflows in Figs. 10, 12, and 13. We also report the emission associated with Sz 102 and Merin 28, in which we can observe gas at different velocities than the systemic velocity.

4.3.1 First glimpse at the CO outflow emission

Figure 2 shows the integrated ¹²CO(2–1) intensity images associated with the three observed sources in which outflows have been detected. In Table 8we list the systemic velocity, the integration ranges used to build the redshitfed and blueshifted images (Δ v_rad, determined by channels with emission above 3σ), the position angle (PA), the opening angle (θ_op), the size, and the dynamical timescales (t_dyn) of the outflows, among other measurements obtained from a homogeneous analysis of the moment images for all three outflows. We derive these quantities for the red- and blue-shifted lobes of each outflow separately. Δ v_rad = ∣v_sys − v_max∣ is the outflow spread in radial velocity, where v_sys is the systemicvelocity and v_max is the velocity up to which it is possible to detect outflow emission over the 3σ threshold. The outflow opening angles are derived from the width of the outflow at a projected distance of 1000 AU from the originating source. The outflow position angles are taken as the bisectors of the angles formed by the lines used to define the opening angles. The outflow inclinations are obtained by assuming that they are perpendicular to the disks (Table 5). The sizes are estimated in projection in the plane of the sky (considering the emission over 3 σ threshold), and then corrected by the inclination. The errors in the measured sizes are given by the angular resolution of each image (Δsize = beam/2). We estimate the dynamical timescales using the sizes and the radial velocity spreads, Δ v_rad, deprojected by each outflow inclination angle (which is the complementary angle of the disk inclination i): t_dyn = size/(v_rad cos(90-i)).

A common characteristic of all three outflows is that the outflow opening angle at a given distance decreases as the radial velocity increases. Also, in general, the opening angles are smaller if measured at larger projected distances from the originating source. This speaks of the intrinsic difficulty in measuring the outflow opening angle and compare measurements from different authors in the literature. Here we have opted to take the angles from inspecting the integrated intensity images at a certain projected distance of 1000 AU, since at that distance all three outflows show a structure in which it is possible to perform this measurement, while at greater distances they become irregular.

The outflow associated with the compact continuum source linked with IRAS 15398-3359 (top panel in Fig. 2) has been very well studied (see Sect. 5.1.1 and Appendix A). This outflow is very collimated (θ_op ~ 30°) and shows adistinct bipolar morphology. It shows brighter CO emission at the cavity walls of the blue-shifted southwest lobe with some arc bridges crossing it, and a more knotty or irregular morphology in the red-shifted lobe. Faint red-shifted emission is detected at the tip of the blue-shifted lobe. This outflow spreads over 27 km s⁻¹ in radial velocity, with the blue and red lobes extending 15.3 and 11.6 km s⁻¹ from the systemic velocity, respectively. Its size (deprojected using an inclination of 66. °14 with respect to the line of sight) is 4350 ± 94 AU. The dynamical time derived from the size and velocity Δv_rad of each lobe,is of about 300 ± 20 yr, revealing that it would be the shortest of the outflow sources in our sample.

A compact continuum millimeter source is coincident with IRAS 16059-3857 and we consider it the driving source of the outflow shown in the central panel of Fig. 2. This outflow is less collimated than that associated with IRAS 15398-3359 (θ_op ~ 75°) and although it looks bipolar, its two lobes are not exactly opposed, but make a 148° angle. The CO outflow extends for 23 km s⁻¹, with the blue and red lobes extending for 9.6 and 13.2 km s⁻¹ respectively. The measured outflow size corrected by its inclination (53. °2) is 39 250 ± 320 AU, with the redlobe being almost two times larger than the blue one. The red lobe size estimate represents a lower limit since it is possible that this outflow may extend beyond the observed field of view. We also estimate its average dynamical time to be 5500 ± 460 yr.

Close to the central source, a V-shaped structure in both lobes traces the walls of the outflow cavity. As pointed out before by Yen et al. (2017), the blue-shifted lobe is oriented due west and the red-shifted toward the northeast. Away from the central source, the emission from the red-shifted lobe is mostly seen in its southern wall, which extends at least up to the limits of the field of view. The blue-shifted lobe presents two singular features: (i) about 3400 AU west from the source there are several gas shreds that extend due south, (ii) about 7300 AU west from the source lies the center of a shell-like structure. This structure spatially coincides with gas moving at red-shifted velocities and a previously reported Herbig-Haro object (HH 78). The possible causes of this peculiar structure will be discussed in Sect. 5.1.2.

The bottom panels of Fig. 2 show the CO emission distribution of the molecular outflow associated with the millimeter continuum source linked to J160115-41523. The outflow velocity spread of the blue and red lobes are 5.6 and 2.6 km s⁻¹, respectively, with a total spread of 8.2km s⁻¹, and the total size is 40 550 ± 1214 AU. Furthermore, its average dynamical time is >18 000 ± 700 yr, making it the largest outflow in our sample.

Close to the central source (bottom right panel of Fig. 2), the outflow shows a bipolar morphology with the red and blue-shifted lobes pointing in a southeast and west direction. It comprises two V-shaped structures associated with the outflow cavity walls of both lobes, and we measure an average opening angle of ~60°.

The presence of blue-shifted emission eastward from the central source (i.e., toward the red-shifted lobe) could be due to the closeness of the outflow main axis to the plane of the sky (it is inclined about 24° with respect to the plane of the sky), and that the wide-angle flow is expanding toward both blue- and red-shifted directions. Further from the central source, the CO emission from the outflow comes mainly from the walls of the cavity dragged by the wind and two bow-like structures at the tips of both lobes. Noteworthy, at a distance of about 5500 AU west from the central source, the blue-shifted lobe deviates about 20° north (from a PA of 266° close to the young star to a PA of 282° at the tip of the blue-shifted lobe of the outflow), turning into a northwest direction.

In two other sources of our sample we have detected signs of the presence of outflows, although we have not been able to delineate them clearly. In Fig. 3 we show moment 0 maps of the blue- and red-shifted emission associated with Sz 102 and Merin 28.

Sz 102 is associated with the Herbig-Haro objects HH 228 W and HH 228 E1, E2, E3, and E4, that extend far beyond the observed field. The dedicated 7m array observations are able to detect blue- and red-shifted emission in the east-west direction, which is in agreement with the position angle of the jet (first reported by Heyer & Graham 1989) associated with Sz 102 (see Appendix A). The velocity ranges in which the blue- and red-shifted emission is detected are 1.0 and 0.7 km s⁻¹, respectively. Probably we are detecting the remnants of the molecular envelope in which the source formed.

The red-shifted emission associated with Merin 28 is an arc-shape structure located at 82′′ of the central source subtending an angle between PA 55° and 75°. This structure is only present in a few channels covering a velocity range of only 0.5 km s⁻¹. The blue-shifted emission is detected barely over the noise level and only in two velocity channels. We probably are not able to detect it since it is contaminated with cloud emission. Both maps show strong emission at the central source position.

Fig. 2 12CO(2–1) integrated intensity images showing the blue-shifted and red-shifted emission for the outflows associated with IRAS15398-3359 (top panel, 12m+7m data), IRAS 16059-3857 (middle panel, 12m+7m data), and J160115-41523 (bottom panels, 7m data at left and 12m data at right). Each panel is labeled with the name of the source. The minimum and maximum velocities over which the emission is integrated are chosen based on visual inspection of the velocity channels (Table 8). Left panels: contours correspond to 10, 20, and 40 times the rms noise level (the rms noise for the full outflow moment 0 blue and red maps are 0.03–0.04, 0.04–0.05, and 0.2–0.08 Jy beam−1 km s−1 for IRAS 15398-3359, IRAS 16059-3857, and J160115-41523, respectively). Bottom right panel: 12CO high-angular resolutionintegrated intensity emission of the blue- and red-shifted lobes toward the central region of J160115-41523, indicated with a box in the bottom left panel. The contours correspond to 3, 5, and 7 times the rms noise levels (rms of the blue and red maps are 0.04 and 0.1 Jy beam−1 km s−1, respectively). The locations of the driving sources match the corresponding millimeter continuum emission, and are marked with black stars in the images. The estimated PA of each outflow lobe is indicated by black arrows. The beam sizes are shown as black ellipses in the lower left of each panel.

Fig. 3 Blue- and red-shifted emission associated with Sz 102 (blue contours: 2, 4, 6, 8 and 10 rms; red contours: 2, 3, and 4 rms) and Merin 28 (contours: 4, 8, and 12 rms). Plus signs in the top panel display the position of the Herbig-Haro objects HH 228 (over the DSS2 IR image), that lie out of the region mapped by the dedicated 7m array observations (black circle). The positions of the continuum emission are indicated with a black star.

4.3.2 Mass, momentum, and energy from the outflows

In order to recover the whole emission from the outflows and therefore to better estimate their masses, we combined the interferometric and total power data of IRAS 15398-3359, IRAS 16059-3857, and J160115-41523, as mentioned in Sect. 3.5 above.

We usethe expressions from Mangum & Shirley (2015) adapted to the ¹²CO(2–1) transition to estimate the column densities. From these, we derive the masses, the momenta and kinetic energies of each lobe separately, and the corresponding values for the entire outflows. In particular, we use the following equation to derive the column density of the gas (this expression is valid in the optically thin case only): $$N_{\text{tot}}(\text{CO})=\frac{3h}{8{\pi }^3{\mu }^{2}J}\frac{Q_{\text{rot}}{\text{e}}^{E_u/k{T}_{\text{ex}}}}{{\text{e}}^{h\nu /k{T}_{\text{ex}}}-1}\frac{{\displaystyle \int T_b}\text{dv}}{[J_{\nu }(T_{ex})-J_{\nu }(T_{\text{bg}})]},$$(2)

where J_ν(T) = (hν∕k)∕(e^hν∕kT − 1) is the Rayleigh-Jeans equivalent temperature. For the ¹²CO(2–1) transition, it can be expressed as: $$\begin{array}{ll}{N}_{\text{tot}}(\text{CO})=\hfill & 1.196\times {10}^{14}\frac{(T_{\text{ex}}+0.921){\text{e}}^{16.596/T_{\text{ex}}}}{{\text{e}}^{11.065/T_{\text{ex}}-1}}\hfill \\ \hfill & \times \frac{T_{\text{B}}{\Delta }_{\text{v}}}{[J_{\nu }(T_{\text{ex}})-J_{\nu }(T_{\text{bg}})]}.\hfill \end{array}$$(3)

To get this expression we use the line strength S = J/(2J + 1) = 2/5, the dipole moment μ = 1.1 × 10¹⁹C cm, the partition function Q_rot = k T_ex /(hB₀)+1/3 with a rigid rotor rotation constant B₀ = 57.635968 GHz, and the degeneracy g = 2J+1 = 5. We put the T_B in K and the Δv velocity interval in km s⁻¹.

We measure the average brightness temperature for the outflow (T_B) for each velocity channel assuming a background temperature of 2.7 K. We avoid the channels contaminated with cloud emission since it could not be disentangled from that of the outflow (the velocity ranges of cloud emission are 4.2–6.8, 3.0–6.8, and 3.14–5.0 km s⁻¹ for IRAS 15398-3359, IRAS 16059-3857, and J160115-41523, respectively). In addition, we estimate the density for two different values of T_ex: 25 and 100 K. This range of temperatures have been used to estimate masses of outflows associated with other youngstars (Arce & Sargent 2006, van Kempen et al. 2009a, Arce et al. 2013). Using these T_ex values we find that some channels are optically thick at the emission peak, which sets a lower limit for the densities.

Using all this, we estimate the mass as: $$M_{\text{out}}=\mu mh\Omega N_{\text{tot}}/X_{\text{CO}},$$(4)

where μ is the mean molecular weight, which is assumed to be equal to 2.76 after allowing for a relative helium abundance of 25% by mass (Yamaguchi et al. 1999), m is the hydrogen atom mass (~1.67 × 10⁻²⁴ g), Ω is the area, and a CO abundance of X_CO = 10⁻⁴ (e.g., Lacy et al. 1994).

We further calculate the momentum and the kinetic energy of the outflows correcting by each of the outflow inclinations and using: $$P_{\text{out}}=M_{\text{out}}\left(\frac{V_r-V_{\text{LSR}}}{\text{cos}(i)}\right),$$(5) $$E_{\text{out}}=\frac{1}{2}{M}_{\text{out}}{\left(\frac{V_r-V_{\text{LSR}}}{\text{cos}(i)}\right)}^2.$$(6)

The outflow masses, along with their momenta and energies, are listed in the seventh, eighth and ninth columns in Table 8. We note the caveat that these values should be treated as lower limits, since we do not correct for the missing flux at cloud velocities, nor the opacity of the emission. The masses corresponding to a T_ex of 25 are typically a factor of 2.5 smaller than those estimated using 100 K, and this factor is a fair representation of the uncertainty in the mass measurements. In our outflow sample, the masses are in the range of 10⁻⁴ and 10⁻⁵M_⊙, the momenta are roughly a few times 10⁻³M_⊙ km s⁻¹, and the energies are between 10⁴⁰ and 10⁴¹ erg.

In the case of IRAS 15398-3359 the derived mass, momentum, and energy agree, in order of magnitude, with those derived by Bjerkeli et al. (2016b) and Dunham et al. (2014), based on ¹²CO(2–1) data observed with the Submillimeter Array (SMA) and the JCMT, respectively.

We estimate lower limits for the average mass loss rate (Ṁ_out = M_out∕t_dyn) and the average rate of linear momentum injected by each outflow, also known as the flux force (Ṗ_out = P_out∕t_dyn). These values are also listed in Table 8 (Cols. 10 and 11), resulting in mass loss rates ranging 10⁻⁹−10⁻⁷ M_⊙ yr⁻¹, and rates of linear momentum ranging 10⁻⁸ − 10⁻⁶ M_⊙ yr⁻¹ for the outflow lobes. Finally, a rough estimate of the accretion rate toward each young star can be made by assuming it is 10% of the outflow mass loss rate (Ṁ_acc = 0.1 Ṁ_out, Pudritz & Banerjee 2005, Ellerbroek et al. 2013). The average numbers for Ṁ_acc are on the order of 10⁻⁸−10⁻⁷ M_⊙ yr⁻¹.

Fig. 4 IRAS 15398-3359. Left: red lobe moment maps (from 5.9 to 16.7 km s−1). Right: blue lobe moment maps (from –10.2 to 4.8 km s−1). The systemic velocity is ~5.1 km s−1. The positionof HH 185 matches with the 2MASS Ks band peak (15:43:01.311, –34:09:14.81) (see also Heyer & Graham 1989) which is shown in black contours in the bottom right panel. Moment 1 and 2 maps only show the emission over 4 rms.

4.3.3 Kinematics

This section contains figures showing the moment maps of first and second order (integrated weighted radial velocity and integrated weighted velocity dispersion, respectively) separated by the cloud velocity, along with several position-velocity (PV) diagrams along and across the outflows of IRAS 15398-3359, IRAS 16059-3857, and J160115-41523 (Figs. 4–9). We analyze the observations and comment on our findings for each individual outflow. We discuss these findings andinterpret them with respect to the nature of the outflows in the next section.

The moment 1 of the ¹²CO emission of IRAS 15398-3359 (Fig. 4, middle panels) shows a nonuniform velocity pattern on each lobe of the outflow. For instance, in the redshifted lobe there are several high-velocity spots, while in the blueshifted lobe, the velocity changes with distance from the protostar from about +1 km s⁻¹ (cyan) to +3 km s⁻¹ (yellow) and then again to +1 km s⁻¹ (v_cloud = 5.1 km s⁻¹). A characteristic high-velocity spot is distinguished close to the tip of the blueshifted lobe. All these high-velocity spots are identified in the moment 2 images of IRAS 15398-3359 (lower panels), as well as three more spots close to the protostar position in the blueshifted side of the outflow. It is worth noting that some of these spots positionally coincide with the emission from the HH 185 object, which is comprised by two sources in the 2MASS images: one unresolved, at the tip of the blueshifted lobe and the other, extending along the outflow body. The faint emission detected southeast of the outflow is discussed in Appendix C.

Figure 5 shows a cut along the IRAS 15398-3359 outflow axis (upper panel) revealing the presence of four knots (which we define as sudden increases of radial velocity at an almost fixed position) in the red lobe and six additional knots in the blue lobe within 14′′ from the driving source. These are labeled as R1, R2, R3, R4, B1, B2, B2, B4, B5, and B6, and most coincide with the spots traced in the moment maps. Red and blue-shifted high-velocity emission is also detected close to the disk position. In addition, the cuts made at 4′′ and 8′′ from the protostar (600 and 1200 AU at the considered distance in projection) across the outflow axis (see lower panels of Fig. 5) reveal two or even three (see e.g., C2 cut) structures with large radial velocity dispersion extending from the cloud velocity to high velocities almost parallel. Most of these structures are probably related to the walls or the tips of the outflow.

Figure 6 shows the moment maps for IRAS 16059-3857 (moment 0 map have already been shown in contours in Fig. 2). A large difference between the two lobes is quickly noticed. Regarding the kinematics of the outflow, the redshifted lobe is smooth and shows a change of velocity from the outer parts of the outflow cavity walls (close to the cloud velocity of 4.7 km s⁻¹) to the inner parts (> 10 km s⁻¹). Unlike in the case of IRAS 15398-3359, the redshifted lobe of IRAS 16059-3857 does not show any high-velocity or high-dispersion spot or knots. The blueshifted side of the outflow in the vicinity of the source (<3500 AU) is v-shaped, but it loses regularity abruptly and continues further west in a blueshifted intricate emission at the position of the HH 87 (16:09:12.8, –39:05:02), surrounded by gas nearly at the cloud velocity. Evidence of the presence of shocked gas can be found in theincrease of velocity dispersion (moment 2 images in Fig. 6) and, of course, its perfect positional match with HH 87. In addition, redshifted emission is detected matching spatially with the blue lobe. We discuss in the Appendix C the possible origin of this structure.

Since the redshifted side of the IRAS 16059-3857 outflow shows a more regular behavior, we only perform PV analysis along this lobe (Fig. 7), rather than including the blueshifted lobe. We identify three different parabolic structures marked with the red, blue, and green dashed lines in the top panel of the figure. There may be more parabolas, but we think our identification may serve well as a proof of concept of the existence of several overlapping kinematic structures. In the PV diagrams transversal to the outflow axis (middle and bottom panels of Fig. 7) we identify various ring-like structures that could be associated with the parabolas found in the longitudinal PV cut. The top of these structures, approximately flat in radial velocity, are indicated by red, blue, and green arrows. The rings move from lower to higher radial velocities as the cut is taken further from the protostellar position. In addition, these transversal PV cuts clearly show the smooth velocity gradient between the edges of the cavity (lower velocities) and the inner regions of the outflow.

The moment maps of the J160115-41523 outflow apparently trace its cavity walls well in the blueshifted side (Fig. 8). The cavity structure ends in a strong blueshifted arc of emission, possibly formed by shocked gas at about 116′′ and a PA of −78° from the protostar position (there are other high-velocity dispersion spots at 117′′ and PA = − 83°, 121′′ and − 87°, and at 104′′ and − 93°). Another arc-like feature with higher velocity dispersion is detected in the path of this side of the outflow, at about 80′′ from the protostar with a PA of −85°. It is also worth noting that close to the protostar, the most prominent side of the cavity is the southernmost. A faint cavity wall is also detected in the left side lobe, tracing an X shape centered on the protostar (see Fig. 13). In the red-shifted side of the outflow (left panels of Fig. 8) the observations show gas at different position angles from the protostar (from 80° to 110°). The most prominent features (higher velocity and dispersion velocity) are a point-like spot at about 77′′ and PA = 108°, and an arc-like complex structure subtending an angle between PA 90° and 100° at about 92′′. There is also another arc further from the protostar and closer to the cloud emission at 114′′ with a PA ~ 96°. There are still another two CO cloudlets at PA (78′′, 70°) and (63′′, 80°). We also detected red-shifted emission southwest of the source (60′′, − 117°). All these features are indicated with black arrows in the middle and bottom panels in Fig. 8.

Figure 9 shows the position-velocity diagrams obtained from the low- and high-angular resolution data (full outflow and close to the position of the source J160115-41523, respectively), following the outflow axis. In the upper panel the PV diagram obtained from the low-resolution data (cut width of 30′′, PA = 98°) shows features related with high velocity gas suggesting the presence of shocks, even though (except in the terminal region of the outflow) theemission is not very clear. In the middle panel, the PV diagram of the central region (cut width of 1′′, PA = 89°) shows high-velocity emission from the disk of the young star detected at the central position. Also, a parabolic structure crossing the cloud velocity is detected in the red lobe. This structure is related with the cavity walls of the redshifted side of the outflow, which should be oriented close to the plane of the sky in order to show such a crossing of the cloud velocity. It has no correspondence in the blueshifted side.

Fig. 5 IRAS 15398-3359. Upper panel: position-velocity diagram of the 12CO(2–1) emission along the outflow axis with a cut width of 1′′. The arrowsshow the presence high velocity gas in both lobes. Bottom panels: position-velocity diagrams of the 12CO(2–1) emission across to the outflow axis at an angular distance of 4′′ and 8′′ from the central source in each lobe.

Fig. 6 IRAS 16059-3857. Left: red lobe moment maps (from 4.9 to 17.9 km s−1). The red-shifted emission spatially coinciding with the blue lobe. Right: blue lobe moment maps (from –4.9 to 4.28 km s−1). We note the emission that extends southward perpendicular to the lobe. The systemic velocity is ~4.7 km s−1. At the end, a bubble-shaped structure is detected. Over the moment 2 map the black dashed contours show theDSS2 emission (6300–6900 Å) that reveals the presence of the Herbig-Haro object HH 87, matching a region of high velocity dispersion. Moment 1 and 2 maps only show the emission over 4σ.

Fig. 7 IRAS 16059-3857. Up: position-velocity diagram of the 12CO(2–1) emission along the red lobe outflow axis with a cut width of 1′′. Bottom: position-velocity diagrams of the 12CO(2–1) emission along 1′′-wide cuts perpendicular to the outflow axis at distances of 1000 AU (C1), 1500 AU (C2), 2000 AU (C3), 2500 AU (C4), 3000 AU (C5), 3500 AU (C6), 4000 AU (C7), and 4500 AU (C8).

5 Discussion

We split this section in two parts. First we analyze the three individual outflows observed in our sample (IRAS 15398-3359, IRAS 16059-3857, and J160115-41523) to discern their nature. In second place we discuss the physical characteristics of the whole sample of young protostars, trying to classify them in an evolutionary series based on their properties. We also discuss the links between evolutionary stage, the molecular observations, and the outflow characteristics.

5.1 Individual outflows

The CO outflows are generally seen as low-velocity shell structures around jets that produce high-velocity shocks along their paths. The higher velocity emission is usually found farther out from the source (Lee et al. 2000)and the position-velocity diagrams usually show one of two distinctive kinematic features: a parabolic structure originating at the driving source, or a convex spur structure with the high-velocity tip near known H₂ bow shocks or HH objects. The former parabolic PV structures could be produced by wide-angle winds (Raga & Cabrit 1993), while the latter are generally modeled under a jet-driven bow shock scenario (Shu et al. 1991).

Regarding the ejection, it is supposed to be intimately related with accretion (Hartigan et al. 1995). A variety of instabilities in the accretion disk can cause the accretion of material from a circumstellar disk of a forming star to be episodic (e.g., Dunham & Vorobyov 2012), and this would in turn correspond to an episodic outflow rate (Vorobyov et al. 2018). In Class 0 and I sources, which are still embedded in their parental cores, the most significant evidence of episodic accretion comes from jets that show a series of knots along their axis (Santiago-García et al. 2009; Hirano et al. 2010; Plunkett et al. 2015).

We adopt the wide-angle wind and the jet-driven wind models to discuss the nature of the outflows of our sample through interpreting the details of the ALMA observations. However, in this study we do not pretend a direct identification between the outflows and these models. We are aware that they are simplified models and that more sophisticated ones would be needed to explain in detail the ejection mechanisms.

Fig. 8 J160115-41523: moment maps of the red (left, from 4.1 to 6.7 km s−1) and blue (right, from –1.5 to 3.6 km s−1) lobes. The systemic velocity is ~4.1 km s−1. Moment 1 and 2 maps only show the emission over 4rms. The black dashed lines in the right upper panel show a 27° deflection in the blue lobe. The arrows point to the features described in the main text, where high velocity or dispersion emission is detected.

5.1.1 IRAS 15398-3359

Figure 10 displays some selected ¹²CO(2–1) velocity channels of the outflow associated with IRAS 15398-3359 at high and low velocities (see Appendix D for the complete velocity cubes of the three outflows). A close inspection of the redshifted low-velocity channels (v ~ 6.5–6.8 km s⁻¹) led us to notice the presence of some CO loops or bow-shocks at different position angles measured from the protostar position. These outflow lobes seem to trace the expanding shell of separate ejections. Moreover, this possibility is strengthened by the fact that the position of the bullet-like shocks shown in the position-velocity diagram in Fig. 5, are not associated with gas in the central lane of the outflow, but at different position angles from the protostar (see high-velocity channels in Fig. 10). The loops and hot-spots seen at different velocities in the CO velocity cube agree well with an scenario with ejections in different directions from the central source. As a proof of concept, we try to identify these outflows using four highly eccentric ellipses (shown with black lines and labeled as R1, R2, R3, and R4 in the Fig. 10). These four structures have different sizes and position angles, and we have drawn them so all four share the same position in one of their apices. The semi-major axes of the ellipses are larger the greater the position angle is. By looking carefully at the rest of the redshifted velocity channels we notice evidence of the presence of these structures along the whole velocity range reached by the outflow. The CO emission appears sometimes delineating the border of the ellipses and in other cases approximately at the intersections (hot-spots) of two or more of them, as happens in the higher velocity channels (see lower panels of the Fig. 10). We further extend this analysis to the more spread-out blueshifted lobe identifying four ellipses with semi-major axes longer than those of the redshifted lobe by a factor of 1.85 (labeled as B1, B2, B3, and B4 in the Fig. 10). These ellipses are not perfectly counter-aligned with their redshifted counterparts by an angle of about 1°–3°. The geometric parameters of all these ellipses are listed in Table 9. The two sets of ellipses describe qualitatively well the emission observed at both low and high velocities. The set of ellipses of one of the outflow’s sides have various sizes. They seem to change their sizes and, perhaps more interestingly, their orientations in intervals of ~ 2000 AU and ~10°, respectively, from smaller to larger.

Although these structures were identified manually, and the total number and position of the ellipses may not be exactly as presented here, this analysis allows us to have a good first approximation, and gives us the opportunity to discuss possible ejection mechanisms of the source. Our analysis leads us to identify each pair of counter-aligning ellipses with an episodic bipolar ejection of material from IRAS 15398-3359. In principle, if the features were ejected with a similar velocity and inclination with respect to the line of sight, the largest ellipses would be associated with the oldest ejections. We estimated the dynamical time for each ejection taking the inclination effect into account (Table 9) and the younger ejection may appear to be more spaced in time than the older for both sides of the outflow (about 50 and 80 yr from the oldest and the youngest in the western side). Due to their smaller size, the dynamical times of the redshifted (eastern) ejections are ~1.4 times shorter. Since the ejections appear to be mostly bipolar in the plane of the sky, it would also be expected that each pair of counter-aligned ejecta share the same inclination with respect to the line of sight. Therefore, the differences in the sizes of two opposite ejecta should be caused by differences in the density of the environment or in the ejection velocity. The former scenario seems more plausible, given the high degree of alignment of the paired ejections and the symmetry of the whole set of ellipses. Nevertheless, it is true that the four pairs of ejections are not exactly bipolar by a few degrees. The deviations from strict bipolarity could be explained by an intrinsic difference in the ejection angle on both sides or, more possibly the displacement of the protostar in the plane of the sky (proper motions) or even the orbital motion of a multiple system. In fact, outflows with a similar appearance have been found associated with high-mass protostars; see for instance the cases of Cepheus A HW2 and S140 (Cunningham et al. 2009; Zapata et al. 2013; Weigelt et al. 2002), where the multiple ejections with different orientations have been explained as producedby a tilt of the ejecting system during the periastron passage of a companion in a very eccentric orbit.

The difference in the position angles of the ejections in IRAS 15398-3359 suggests that the outflow is ejecting material episodically and that the originating source is perhaps precessing, as already proposed by Bjerkeli et al. (2016b). These authors proposed that the different PA between the two outflow sides could possibly be due to precession of the ejection axis, originated by the tidal interaction of a binary companion. However, an inspection of the 3.6, 4.5, and 8.0 μm Spitzer images show that the infrared emission is not directly associated with IRAS 15398-3359 but may probably stem from shocks of the blueshifted outflow close to the star (see Fig. 11). Regardless of the cause for the tilt on the system, Jørgensen et al. (2013) found what seems to be the chemical footprint of an accretion burst in the past 100 to 1000 yr (a ring of H₁₃CO⁺ which has been destroyed at the center by the water vaporized during the luminosity bump created by the accretion burst). This time estimate coincides quite well with the dynamical time of the youngest of the ejections (400–600 yr for the redshifted and the blueshifted side respectively).

Another clue that makes us think that the outflow could have an episodic behavior is the pattern detected in the position-velocity diagram (Fig. 5). Inspecting this figure we detect four high-velocity spur-like features (jumps in radial velocity) in the red lobe and six in the blue lobe. These spur-like features in a PV diagram are reminiscent of the jet-driven wind model described in (e.g., Lee et al. 2000). In addition they match several large-density knots positions on the velocity cubes. The spurs are spaced quite evenly, suggesting the existence of an episodic outflow (e.g., Plunkett et al. 2015). They probably trace bow-shocks formed by variations in the mass loss rate or jet velocity, which likewise can be caused by variations in the accretion rate, or side-shocks that are produced when a new outflow ejection is launched in a different orientation and sweeps partially the trail of dragged material left by a previous ejection. As seen in the upper and lower panels of the Fig. 10, some high-velocity spurs could be originated in the intersections of two or more ellipses (i.e., the collision of different ejections). The shape of each of the high-velocity signatures in the PV matches with a jet-driven model, showing a spur structure in the PV diagram along the jet axis, as described by Lee et al. (2000).

Finally, if we consider the fact that IRAS 15398-3359 may be somehow undergoing bursts of accretion and driving a new outflow ejection in different orientations, it is likely that the estimate of our previous outflow’s dynamical time is wrongly reported. The different inclinations with respect to the plane of the sky of each ejection would have implicationsin the calculated dynamical times. It is likewise worth to note that knowledge about the nature of an outflow is crucial for correlating this dynamical time estimate against any other relevant measurement of the protostar.

Fig. 9 J160115-41523. Upper panel: position-velocity diagram of the 12CO(2–1) emission along the outflow axis with a cut width of 30′′ (low resolution, full outflow). Arrows point out high velocity emission. Middle panel: position-velocity diagram of the 12CO(2–1) emission along the outflow axis with a cut width of 1′′ (high resolution, central region). We indicate the emission from the disk. A parabolic structure, indicated with a red dashed line,crossing the cloud velocity is detected in the red lobe. Bottom panels: position-velocity diagrams of the 12CO(2–1) emission across to the outflow axis at an angular distance of 2′′ and 5′′ from the central source in the red lobe, and 2′′ in the blue lobe.

Fig. 10 IRAS 15398-3359: 12CO(2–1) emission velocity channel maps. The star indicates the position of the compact continuum source. Black ellipses show the four bipolar structures.

Fig. 11 Spitzer IRAC emission at 3.5 μm (top right), 4.5 μm (bottom left), and 8.0 μm (bottom right) associated with IRAS 15398-3359. Top left pannel: blue-shifted and red-shifted 12CO(2–1) emission. The black star indicates the position of the compact source. The infrared emission extending southwestward may stem from shocks of the blueshifted outflow close to the star. The spot in RA, Dec(J2000) = (15:43:01.3; −34:09:15.4) seen at 3.4 and 4.5 μm coincides with the terminal part of the blue lobe.

5.1.2 IRAS 16059-3857

The outflow associated with IRAS 16059-3857 presents two lobes with very different characteristics: the red lobe, previously analyzed in Sect. 4.3.3, shows the typical kinematic features of an episodic wide-angle outflow, while the blue lobe shows a peculiar emission distribution. The presence of more than one parabolic structure apparently sharing a similar origin in the PV diagram of Fig. 7 could indicate an episodic outflow that has had multiple ejections (e.g., Zhang et al. 2019). Interestingly, unlike in the case of IRAS 15398-3359 these parabolic structures in the PV diagrams correspond better with the wide-angle wind model (Lee et al. 2000). The episodic ejections of IRAS 16059-3857 do not show abrupt spurs in velocity (i.e., no strong velocity shocks), but a continuous velocity increase as they move farther from the driving system. In this case, the different ejections seem to be launched with a similar orientation, but small perturbations may be hindered due to the fact that they seem to be launched almost isotropically, but with different thrusts depending on the polar angle (see wide-angle wind model details in the literature, e.g., Lee et al. 2000).

Figure 12 shows the channel maps of the redshifted lobe from 5.88 to 8.42 km s⁻¹. After analyzing the images, we could identify several elliptical structures. The parameters of these ellipses (central position, major and minor axis, and position angle) are listed in Table 10. We realized that there are at least two elliptical features in each velocity channel, making two distinct sets that apparently evolve with velocity. By moving away from the systemic velocity (4.7 km s⁻¹), the central position of the ellipses moves along the axis of the outflow, away from the continuum source. This Hubble-law-like behavior is again expected when a wide-angle wind is launched into an environment with a density distribution sinusoidally stratified along the polar axis (more dense in the equatorial plane than in the polar caps). Hence we hypothesize that in the case of IRAS 16059-3857 a wide-angle wind is episodically ejected. This different nature compared with that of the jet-driven wind in IRAS 15398-3359 may be due to environmental conditions (different density distributions) or to intrinsically different ejection mechanisms (isotropic vs polar wind), but present observations do not allow us to discern between both scenarios.

Fig. 12 IRAS 16059-3857: 12CO(2–1) emission velocity channel maps. The star indicates the position of the compact continuum source. We identified elliptical structures labeled as e1 and e2.

5.1.3 J160115-41523

The outflow from J160115-41523 is apparently not as complicated as the other two in our sample. It shows at first sight more or less classicalcavity walls, probably excavated by a high-velocity jet. Close to the protostar, the high-angular resolution ALMA data apparently shows the typical biconical shape of outflows (see Fig. 13). Although most of the evidence indicates that the outflow coincides with a jet-driven model, there is other evidence that agrees with a wide-angle wind model, as the parabolic structure in the CO emission of the PV diagram of its redshifted lobe (Fig. 9). The eastern and western tips of the outflow end up both in bow-shock structures that show larger velocities and dispersions.

There is however the possibility that this outflow is not so, let us say, well-behaved. The upper panels of Fig. 8 show the moment 0 maps of the outflow associated with J160115-41523. We already noted that the morphology of each lobe is different, and a change of direction of the blueshifted lobe is noticeable at about 6000 AU from the central source (see dashed lines in Fig. 8). In this same blueshifted lobe, between theoutflow deviation and the outflow end there are various quasi-parallel stripes of gas oriented close to north-south and seen at the cloud’s velocity. Also, as noted in previous sections, the redshifted lobe shows CO emission at position angles north of the main outflow lobe. We wonder if these features at both sides of the protostar position are related, and if it is possible that the outflow from J160115-41523 drove several ejections with different orientations as IRAS 15398-3359. If this is the case, the north-south stripes of gas in the blueshifted side of the outflow, could be related with the interaction of different outflow ejections and the cloud surface (which emission is more prominent in the western side of the FoV). These hypotheses are right now mere speculations, but the data show hints that the quiet and typical scenario may not be the final answer for this source. New high-angular resolution and better sensitivity observations along the outflow may help to interpret the nature of the outflow.

Fig. 13 J160115-41523: blue (up) and red (down) high resolution velocity channel maps.

5.2 Characteristics of the source sample

In the following section we discuss the results obtained from the SEDs, chemistry, and outflow parameters with the purpose of determining the evolutionary state of the sources in our sample.

5.2.1 SED evolution

As premain sequence stars evolve, their envelope mass is expected to decrease due to accretion onto the disk-protostar system and-or dispersion by outflowing winds. Thereby, the ratio of the disk/envelope mass can be regarded as an age indicator of a protostar. The envelope material also affects the characteristic spectrum of the escaping light. It absorbs light from lower wavelengths, and reprocesses it into IR and (sub)mm wavelengths, providing the characteristic two-bump Spectral Energy Distribution of many protostars: the submm bump is explained by a T_eff< 100 K blackbody which corresponds to the cold dust of the envelope emission (T_d in Table 6, hereafter T_cold), and the second near-infrared bump, with larger effective temperature, corresponds to the warmer dust emitting in the inner disk, as well as the scattered and reprocessed light from the protostar (we characterize it as T_warm). The SEDs of premain sequence stars and the ratio between the bolometric luminosity and the luminosity from radio up to submillimeter wavelengths only (L_bol - L_submm), have been thus utilized as age indicators, based on the relative importance of these two bumps. Hence the spectral index between 2 and ~25 μm (Adams et al. 1987; Greene et al. 1994; Young & Evans 2005), and an average color measurement translated into a blackbody temperature: the bolometric temperature, T_bol (Myers & Ladd 1993), are commonly invoked to classify the premain sequence stars into evolutionary stages or classes (Andre et al. 1993). However, there are some known problems with such indicators due to the non-isotropic nature of the accretion-outflow processes and the different inclinations of the protostellar objects (e.g., Crapsi et al. 2008; Evans et al. 2009, and references therein). Sometimes various indicators are needed to obtain a more realistic evolutionary situation, although the bolometric temperature has been the most trusted age proxy, and it is the base upon which some age estimates are constructed (e.g., Young & Evans 2005).

In this work we collect some commonly used protostellar age indicators for our sampled objects. We summarize some of them in Table 11 along with certain physical characteristics. A quick inspection of the SEDs (Fig. 1) shows that only one of the sources does not have emission bluewards of ~25 μm (AzTEC-lup1-2) suggesting this source is probably in the prestellar phase, only surrounded by cold dust. AzTEC-lup1-2 has a one-bump SED while the rest have emission at near infrared wavelengths and most of them show a second bump peaking bluewards of 25 μm as well. We estimate the black-body temperatures of these two bumps for every SED (T_cold and T_warm, see Table 11). At the warmer end of the SEDs we have Sz 102 and Merin 28. For these two sources, T_warm was estimated by fitting the fluxes at the shortest wavelengths, since the separation between the two bumps is not so evident, suggesting that they are probably in a more evolved stage in which the surroundings of the young stars is becoming more transparent.We should keep in mind that stars with face-on or low-inclined disks can also show a SED with two ill-separated bumps (e.g., Crapsi et al. 2008; Enoch et al. 2009), but in this case, the disks of Sz 102 and Merin 28 are closer to edge-on. The rest of the sources display SEDs with two clear bumps. We see in the following that fine-tuning the age classification of our young stellar sample is not straightforward, even after using several of the age indicators at hand.

We use the SEDs to estimate T_bol as: $$T_{\text{bol}}=1.25\times {10}^{-11}\langle \nu \rangle \text{K},$$

where the mean frequency is $$\langle \nu \rangle =\text{ frac{}{\displaystyle \int \nu }S\text{{\_}\nu \text{}d}\nu \text{}{}{\displaystyle \int S}\text{{\_}\nu \text{}d}\nu \text{}}.$$

We also estimate the bolometric luminosity, L_bol: $$L_{\text{bol}}=4\pi d^2{\displaystyle \int S_{\nu }}\text{d}\nu ,$$

where d is the distance and S_ν is the flux at the frequency ν. We define the ratio L_bol/L_submm, in which L_submm is obtained integrating the SED at λ > 350 μm. Their valuesand uncertainties are in Table 11 and were obtained integrating by two different methods (Simpson and trapezoidal rules), before further averaging (final estimate) and subtracting (uncertainty) them both. The infrared and (sub)millimeter spectral indices (α_IR and β_mm, respectively) of the source sample are another typical evolutionary indicators shown in Table 11. The table contains the ratio between the disk and envelope mass, M_disk/M_env as well, calculated from the masses estimated by the continuum emission observed using the 12m-only (M_disk) and the substracting between the 7m-only (M_env+disk) and the 12m-only array. This ratio cannot be calculated for Sz 102 and Merin 28, since, as mentioned in Sect. 4.1, the envelope masses would be zero.

AzTEC-lup1-2 and AzTEC-lup3-5 are the two sources with no disk detection (M_disk/M_env = 0) as well as a L_bol/L_submm ratio lower than 2; therefore, we consider that they are probably in the prestellar phase. While AzTEC-lup1-2’s T_bol is, as expected, very low (10 K), AzTEC-lup3-5’s is over 100 K due to the second bump in the SED with T_warm = 1000 K. Higher spatial resolution infrared observations are needed to decide whether this second warm dust component is truly associated with the cold dust component or it is just spatially overlapped. In the meantime, it seems that for AzTEC-lup3-5 T_bol is not a robust tracer of its evolutionary stage.

On the upper end of the M_disk/M_env ratio (which diverges) Sz 102 and Merin 28 have a similar luminosity (L_bol ~ 0.1L_⊙), and show the two largest T_bol values (>450 K) and L_bol/L_submm ratios (>130). Following Table 1 of Enoch et al. (2009), these two sources may be classified as late Class I protostars, since 300 K < T_bol < 650 K. Their infrared spectral indices however are not flat or negative as expected for Class I stars, indicating the problems of α_IR as a reliable age proxy in these early evolutionary stages. We notice too that Sz 102 has a warm temperature, larger than 4000 K, which may indicate that the protostar has already taken off the envelope veil, maybe aided by a more face-on disk which is in addition much smaller and less massive than Merin 28’s disk. Sz 102 has a lower T_cold than Merin 28 though. These two sources are likely transiting a similar evolutionary stage, but the T_bol and the ratio L_bol/L_submm do not agree about which is at an earlier stage. We note at this point that these SEDs are not corrected by extinction, which may produce differences in the results.

The remaining three sources (IRAS 15398-3359, IRAS 16059-3857, and J160115-41523) show prominent CO emission associated with their outflows, suggesting that they entered into a prominent accretion phase recently (outflow dynamical times less than a few 10⁴ yr). Following again the classification scheme of Enoch et al. (2009) based on T_bol, IRAS 16059-3857 may be the youngest of the three (early Class 0), while IRAS 15398-3359 and J160115-41523 may be more evolved (both in the late Class 0 or early Class I stage) (see also Chen et al. 1995). We obtain the exact same order when taking the ratio L_bol/L_submm into account. Regarding the disk to envelope mass ratio, the three sources are arranged in the opposite order as set by the bolometric luminosity. If we take as valid the trend based on the T_bol, then the ratio between disk and envelope mass should be taken with care when dating objects in a similar evolutionary stage. We should make the caveat that IRAS 15398-3359 is one order of magnitude brighter than the other two, though. It should be noted that since the disks are near edge on, they may actually be older than their SED appears (Crapsi et al. 2008).

Overall, the most used evolutionary indicators (T_bol, and the ratios L_bol/L_submm and M_disk/M_env) seem to agree in diagnosing which are generally more or less evolved sources, but they do not agree completely, and sometimes the different indicators lead to contradictions when arranging sources within a similar stage. The infrared or submillimeter spectral indices α and β are not robust tools for this sample.

5.2.2 Chemical evolution

Regarding the detected molecular lines in the present ALMA observations, ¹²CO and C¹⁸O are present in all the observed sources. They are observed as part of the cloud emission, the envelope, the disk, and the outflow. Only AzTEC-lup1-2, IRAS 15398-3359, and IRAS 16059-3857 were associated with emission from other lines (see Tables 4 and 7, and Fig. 1). In particular, these three sources show emission from the deuterated species N₂ D⁺ and DCN. Their hydrogenated counterparts are known to trace dense gas and low temperature regions. N₂ H⁺ is mainly associated with quiescent gas from envelopes and interstellar medium filaments (e.g., Tobin et al. 2013), while HCN has also been observed at the cavity walls and head-shocks of outflows (e.g., Busquet et al. 2017). At very low temperatures, CO and other heavy and abundant molecules freeze-out onto dust grains. As temperatures increase due to the growing of a protostar, the CO molecules are released into the gas-phase again. It is worth to note that N₂ H⁺ is thought to be destroyed as CO is released during this process (e.g., Roberts et al. 2003; Busquet et al. 2011; Stephens et al. 2015). This is most important in HII regions or PDRs (Photo Dissociation Regions), but it is also a process seen when protostars evolve and increase their temperature (e.g., L1157, Kwon et al. 2015). Therefore, N₂ H⁺ is expected mainly in cold (and therefore less evolved) dusty envelopes where CO is still depleted in the dust grain icy mantles, since nitrogen-bearing species are relatively unaffected by freeze-out (Charnley 1997). Besides, larger deuterium fractionation (or D-enrichment) has been principally detected in prestellar and Class 0 sources (Ceccarelli et al. 2014). Again, relatively larger abundances of deuterated species are expected due to the very cold temperatures on the order of 10 K and the freeze-out of abundant neutral species like CO. DCN had been reported in low-velocity emission associated with outflows (e.g., Codella et al. 2012), while N₂ D⁺ is found in dense envelopes of young protostars, mainly in their colder outskirts (e.g., Tobin et al. 2013). We can use these molecules as indicators of stellar evolution in the first stages of star formation. The presence of DCN and N₂ D⁺ in three of the sampled sources (AzTEC-lup1-2, IRAS 15398-3359, and IRAS 16059-3857) place them at the earliest stages of protostellar evolution, while the non-detection of the nitrogen-bearing deuterated species in the rest of the sources suggests they are warmer and thereby more evolved. Interpretation of the chemistry of disks and envelopes is complicated, even when large samples of objects in the same evolutionary phase are observed with high-angular resolution by sensitive molecular surveys (e.g., Miotello et al. 2019), therefore our results should be taken with care.

For IRAS 15398-3359 the observations show the detection of other lines such as CH₃OH, SO, and SiO. These molecules typically trace grain-mantle released into the gas phase, and their abundances in the shocked gas are enhanced with respect to the quiescent gas (Jiménez-Serra et al. 2005; Burkhardt et al. 2016). Some of these are known to trace shocks in outflows, and their presence is usually interpreted as a sign of youth in protostellar outflows (the presence of these molecules associated with IRAS 15398-3359 has already been discussed by Jørgensen et al. 2013). However, in this case, the presence of fresh gas enhancing the shock chemistry could be not only a matter of the protostellar youth, but also due to the ejections at different position angles of the outflow. The youngest ejections may encounter fresh gas (from a supply of non-previously shocked icy-coated grains) each time the system changes its orientation, thus producing a rich chemistry regardless of the protostellar age. This situation in IRAS 15398-3359 should be taken as a cautionary call when using shocks or chemistry as an outflow signpost of youth. The nature of the individual systems can modify the interpretation of the observations in at least one of every ten sources.

5.2.3 Outflow evolution

Acknowledging the complications imposed by the inaccuracies of the measurements of the physical characteristics of the outflows, and despite realizing that they have different intrinsic natures (wind-driven vs jet-driven) or peculiarities (precession, pulsation), in this section we compare our findings with similar efforts in the literature. We also test some commonly referenced correlations used for evolutionary trends of outflows.

In order to compare the outflow physical characteristics we picked up two similar works on a sample of outflows from low-mass stars: Arce & Sargent (2006) and Dunham et al. (2014). These two papers show a list of outflows whose originating sources are categorized in different evolutionary stages. We have split the sources into three main groups (Class 0, Class I, and Class II sources), and derived median values for the kinetic energy, the mass loss rate and the flux force (see Table 12). We left out three sources classified either as Class 0 or Class I in Dunham et al. (2014).

Molecular outflows have been presented in the literature as a possible evolutionary tracer for their corresponding protostellar system (e.g., Curtis et al. 2010; Arce & Sargent 2006; Plunkett et al. 2013; Y"i"ld"i"z et al. 2015; Arce et al. 2007). Given their ubiquity, morphological and kinematic (velocity) extent, sufficient emission in order to be detected, and the fact that they do not need high-angular resolution observations to be analyzed, the outflows are indeed a promising observational probe. However, their measured dynamical times do not trace well the evolutionary stage, and molecular outflows from Class I and Class II sources can be found with similar dynamical times as outflows from Class 0 protostars. Hence dynamical times of molecular outflows are not good age tracers and should be taken as lower limits for the age of their corresponding originating sources. Measurements of the energy and momentum of outflows resulted in the general consideration that powerful outflows are younger (e.g., Bontemps et al. 1996; Saraceno et al. 1996). This has been derived from the correlation of the outflow momentum flux (analogous of the rate of momentum, Ṗ_out, reported in this work) with the bolometric luminosity and-or the envelope mass. As previously mentioned, one of the best age proxies for young stellar objects is the bolometric temperature (Ladd et al. 1998), and therefore it would be good to test the correlation of the outflow parameters against T_bol when seeking for an age estimator. Although this should be tested in a statistically significant sample, we compare the three studied outflows under the view of previous statistical studies to see whether they fit in with general trends.

In this work, we have detected classical molecular outflows from three sources: IRAS 15398-3359, IRAS 16059-3857, and J160115-41523 (see Table 4). Based on the SOLA catalog SED fitting, these sources are classified as Class 0, Class 0, andlate-Class I, respectively (see Table 1) but, in the analysis of Sect. 5.2.1 based on T_bol and the L_bol/L_submm ratio, we classified them as late Class 0, early Class 0, and late Class 0, respectively (all three may be Class 0 protostars). We use the latter arrangement in the following discussion.

Regarding the bolometric temperature, although the measurements are somewhat uncertain due to the typical scarcity of experimental data, the results for the three outflows reflect that J160115-41523 might be older than IRAS 15398-3359, and both older than IRAS 16059-3857 (90, 68, and 42 K, respectively). We note however that if we use the ratio of disk mass to envelope mass (another typical protostellar age proxy, that becomes larger with time) the reverse order is obtained, with M_disk/M_env ratios of 0.3, 0.4, and 3.6 for J160115-41523, IRAS 15398-3359, and IRAS 16059-3857.

One method of characterizing outflows is based on their emission structure in spatial-velocity space, specifically length, opening angle, and velocity extent. Curtis et al. (2010) found that the longest outflows in their sample were Class 0, compared with later classes; while reporting scatter, they indicate that the class-dependent length is significant. They found average lobe lengths of (140 ± 20) and (88 ± 11) arcsec (which at their assumed distance of Perseus, d = 250 pc, corresponds to 3.5 × 10⁴ and 2 × 10⁴ AU) for Class 0 and I sources respectively. According to Table 8 in the present work, we report the outflow lobe sizes (length, corrected for inclination) to be between 2 × 10³ and 2.5 × 10⁴ AU. IRAS 16059-3857 shows the longest outflow (also likely extending beyond the FOV), consistent with the finding from Curtis et al. (2010) that the youngest Class 0 objects drive longer outflows. The outflow of J160115-41523 (also a late Class 0 or early Class I) is consistent with the Class I average outflow length from Curtis et al. (2010). While for these two outflows the trend seems to qualitatively apply, on the other hand the shortest outflow in our sample (IRAS 15398-3359’s) is also a Class 0. Hence, the lobe size apparently is not a consistent evolutionary indicator for all sources in our sample, provided T_bol is.

Outflow opening angle (θ_op) is also presented in Table 8. Arce & Sargent (2006) showed that, even while outflow morphology varies from source to source, a significant trend holds such that outflow opening angle appears to broaden with time. They fit the outflow angle as a function of age for 17 sources with the relation: $$\log (\theta /\text{deg})=(0.7\pm 0.2)+(0.26\pm 0.4)\log (t/\text{yr}).$$

The ages were determined from the T_bol-age relation from Ladd et al. (1998). Among the outflows in our sample, IRAS 15398-3359 has a significantly smaller opening angle than IRAS 16059-3857 and J160115-415235. The latter have opening angles more consistent with Class I sources in Arce & Sargent (2006). The trend of outflow opening angle is consistent with that of lobe size, previously discussed. In both cases, it seems that IRAS 15398-3359 is the least evolved, while IRAS 16059-3857 and J160115-415235 have similar characteristics.

Curtis et al. (2010) also presents maximum outflow velocity as an evolutionary indicator; velocity is apparently correlated with length, so that the Class 0 sources drive longer and faster outflows than Class I sources. The velocity extent is larger in IRAS 16059-3857 than in J160115-41523 (apparently decreasing with age). The velocity extent of IRAS 15398-3359 is the largest, although thought to be a late Class 0. Although the velocity has not been corrected for inclination angle, either in Table 8, or in Curtis et al. (2010), the lack of a trend in our sample holds for speed, and not merely observed radial velocity.

Momentum, energy, mass-loss rate, and momentum flux are kinematic properties of outflows that have also been linked to age in the literature. Momentum and energy are thought to decrease from Class 0 to I. Momentum has been argued to be a more representative quantity because outflows are thought to be momentum-driven (Shu et al. 1991, 2000; Masson & Chernin 1993; Lee et al. 2000). Curtis et al. (2010) found average momentum for Class 0 and Class I sources to be ⟨P₀⟩ = 0.7 M_⊙ (km s⁻¹)⁻¹, ⟨P_I⟩ = 0.1 M_⊙ (km s⁻¹)⁻¹; and average energy to be ⟨E₀⟩ = 1.4 × 10⁴⁴ erg, ⟨E_I⟩ = 1.0 × 10⁴³, respectively (omitting from the means SVS 13, which may be an anomalously strong Class I outflow).

The outflows in our sample appear to follow a trend with age, as their energies are arranged in the same order as the T_bol. The momentum almost follows the same trend, with a reversal between IRAS 15398-3359 and J160115-41523.

Momentum flux (or rate of injected momentum, Ṗ_out) may be another age indicator, since it appears to correlate with the envelope mass in both Class 0 and I sources (Bontemps et al. 1996; Saraceno et al. 1996). For IRAS 15398-3359, J160115-41523, and IRAS 16059-3857, the envelope masses are 2.3 M_J, 20.2 M_J, and 36.0 M_J, respectively, and do not correlate with the measured momentum fluxes. However, the measured momentum fluxes inversely correlate with the bolometric temperatures (i.e., the lower the temperature, the larger the momentum flux), hence this could be another useful indicator for the age of outflows.

Our assessment of outflow characteristics, including morphology and energetics as evolutionary indicators, reveals that the nuances of the outflows themselves render these trends unreliable on a source-by-source comparison. In summary, if we assume that IRAS 16059-3857 is the least evolved (“youngest”), and J160115-41523 is the most evolved (“oldest”), with IRAS 15398-3359 intermediate, then most trends for evolution found in the literature do not consistently hold for all three sources. Based on limited statistics in our sample, we do not negate the statistical significance of overall trends such as outflow opening angle increasing with time, or length decreasing with evolution. On the other hand, the age of our sample does apparently correlate with energy and momentum flux. Nonetheless, we point out that specific characteristics of individual sources may render one (observational) trend more robust than another. Observational limitations may include, but are not limited to, uncertainties in distance determination or inclination of the outflow. Additionally, a thorough investigation into the outflows known to be launched from sources undergoing episodic accretion (Y"i"ld"i"z et al. 2015), may reveal additional revisions to the understanding of observational outflow trends.

6 Summary

In this paper we present the analysis and results of the millimeter continuum and molecular gas emission associated with seven prestellar and protostellar objects located in the Lupus molecular clouds. Based on intermediate angular resolution dedicated observationsand high angular resolution data from the ALMA archive, we have analyzed the dusty disks and envelopes associated with each source and the molecular outflows observed for some of them. We present for the first time fully interferometrically mapped data of the outflows associated with IRAS 16059-3857 and J160115-41523. Based on our results we have attempted to delineate an evolutionary sequence for our source sample by testing different evolutionary indicators (see Table 13).

From the SEDs and the continuum emission we obtained the bolometric temperature T_bol, the ratios L_bol/L_submm and M_disk/M_env, and the IR and millimetric spectral indices (α and β). We found that the first three indicators seem to be good at determining whether a source is inside one of the typical protostellar Classes, but they are not as effective in arranging (by age) sources in similar evolutionary stages. The spectral indices are not good diagnostic tools for our sample. Based on these indicators we found that AzTEC-lup1-2 and AzTEC-lup3-5 are the youngest sources (no disk detection) being probably in the prestellar phase, while Sz 102 and Merin 28 seem to be the most evolved (they are probably in the Class I) since the surroundings of the young stars appear have become more transparent.The three remaining sources (IRAS 15398-3359, IRAS 16059-3857, and J160115-41523) show prominent ¹²CO emission associated with their outflows and appear to be in an intermediate evolutionary phase in between the two previous groups (they are probably in the Class 0), having recently entered into an accretion phase. For these three sources we found distinct ordering results when using different age indicators.

We detected ¹²CO and C¹⁸O emission associated with all the sources in our sample. Three of them (AzTEC-lup1-2,IRAS 15398-3359, and IRAS 16059-3857), chemically more complex, also present emission of deuterated species (DCN and N₂ D⁺) suggesting that they would be younger. Only one source in our sample (IRAS 15398-3359) presents CH₃OH and SiO emission, usually interpreted as a sign of youth in protostellar outflows. However, as the observations show, IRAS 15398-3359 (the most luminous protostar in our sample too) is possibly associated with a precessing outflow, which may produce shocks in the interactions between the ejections and the still fresh unperturbed gas in different directions.

As we mentioned, the parameters obtained from the SEDs do not help us to discern between the evolutionary state of the three sources presenting clearly detected outflows. In an effort to find additional indicators, we have tested several evolutionary indicators obtained from the outflow emission such as the dynamical time, opening angle, maximum outflow velocity, momentum, energy, and momentum flux. Although most of the evolutionary indicators do not consistently hold for the three sources, it is most likely that IRAS 16059-3857 is the youngest, followed by IRAS 15398-3359 and J160115-41523.

In addition to this, we have characterized the molecular outflows from their morphology and kinematics and contrasted them with the wide-angle and jet-driven models. We identified four pairs of bipolar elliptical structures in the outflow associated with IRAS 15398-3359, which have different sizes and position angles. We propose that we are possibly detecting an episodic precessing outflow where every ellipse corresponds to a bipolar ejection. The PV diagram supports the hypothesis of an episodic outflow, showing high-velocity spur-like features quite evenly spaced in both lobes, reminiscent of a jet-driven wind model. On the other hand, IRAS 16059-3857 shows two lobes with very different morphologies. The redshifted lobe exhibits the typical kinematic features of an episodic wide-angle outflow, evidenced by the presence of more than one parabolic structure in the PV diagram. The elliptical structures in the redshifted lobe agree with the hypothesis of a wide-angle wind episodically ejected, rather than a jet-driven outflow. Otherwise, the terminal part of the blueshifted lobe shows a peculiar bubble-shaped structure which coincides with HH 78. We propose this is caused because the outflow is emerging from the primordial cloud where the protostar forms. Finally, J160115-41523 also shows some peculiarities. The blueshifted lobe suffers a change of direction and various quasi-parallel stripes of gas appear close to the cloud velocity, whereas the red-shifted lobe shows CO emission at position angles north of the main outflow lobe. We wonder if these features on both sides of the protostar position are related, or if it is possible that the outflow has driven several ejections with different orientations as in IRAS 15398-3359. Data with higher angular resolution and more sensitivity will be needed to find an answer.

Acknowledgements

We are very grateful to the anonymous referees for all the comments and suggestions that helped a lot to improve this work. This paper makes use of the following ALMA data: ADS/JAO.ALMA 2016.1.01141.S, ADS/JAO.ALMA 2011.0.00777.S, ADS/JAO.ALMA 2013.1.00879.S, ADS/JAO.ALMA 2011.0.00628.S, ADS/JAO.ALMA 2013.1.00244.S, ADS/JAO.ALMA 2015.1.00306.S, ADS/JAO.ALMA 2013.1.00474.S, ADS/JAO.ALMA 2016.1.00571.S, ADS/JAO.ALMA 2016.1.00459.S, ADS/JAO.ALMA 2016.1.01239.S, ADS/JAO.ALMA 2016.1.01239.S, ADS/JAO.ALMA 2015.1.01510.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The authors are grateful to the visiting program to JAO and ESO that facilitated this collaborative effort. M.M.V. is funded by post-graduate scholarships from the Argentine Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET). Id.G.-M. is partially supported by MCIU-AEI (Spain) grant AYA2017-84390-C2-R (co-funded by FEDER).

Appendix A Summary of individual sources

In this appendix we briefly describe some respects of the bibliography of each source.

A.1 AzTEC-lup1-2 and AzTEC-lup3-5

The sources AzTEC-lup1-2 and AzTEC-lup3-5 are located in Lupus I and Lupus III clouds, respectively, and were classified in the SOLA catalog as prestellar cores. The peak intensity from ASTE/AzTEC 1.1 mm continuum are I_peak = 455 mJy beam⁻¹ and I_peak = 169 mJy beam⁻¹, respectively. AzTEC-lup3-5 is associated with the IR source J160859.7-390737 detected in all the Spitzer bands. The cores where the sources are immersed were cataloged by Benedettini et al. (2018) from Herschel images. The authors calculate a column density N(H₂) = 29.61 × 10²¹ cm⁻³ and N(H₂) = 7.59 × 10²¹cm⁻³, for AzTEC-lup1-2 and AzTEC-lup3-5 respectively, and classified both as prestellar.

A.2 IRAS 15398-3359

IRAS 15398-3359 is located in Lupus I and is embedded in the dense core B228 (Barnard et al. 1927). B228 has a mass of 0.3–0.8 M_⊙ (Reipurth et al. 1993; Shirley et al. 2000) and is located at a rest velocity of 5.1 km s⁻¹ (Mardones et al. 1997; van Kempen et al. 2009b). The IRAS 15398-3359 protostar was first identified by Heyer & Graham (1989), who also discovered HH 185, a Herbig-Haro object driven by this source, later on cataloged as Class 0 protostar by Evans et al. (2009). The outflow driven by IRAS 15398-3359 was discovered by Tachihara et al. (1996) and mapped in HCO⁺(4–3) and various CO transitions by van Kempen et al. (2009b), who determined a temperature of ~100–200 K for the outflowing gas, based on the line ratios in different transitions.

From H₂CO and CCH ALMA observations, Oya et al. (2014) detected the outflow extending northeast-southwest, with an inclination angle of 20° (almost along the line of sight). They estimated an upper limit to the outflow mass of 0.09 M_⊙. Bjerkeli et al. (2016b) detected theoutflow in ¹²CO, ¹³CO, HCO⁺, and N₂ H⁺ using SMA (Submillimeter Array) data, and calculated the mass (2 × 10⁻⁴ M_⊙), momentum (4 × 10⁻⁴ M_⊙ km s⁻¹), momentum rate (10⁻⁶ M_⊙ km s⁻¹ yr⁻¹), mechanical luminosity (10⁻³ L_⊙), kinetic energy (5 × 10⁴⁰ erg), and mass-loss rate (4 × 10⁻⁸ M_⊙ yr⁻¹). Using C¹⁸O observations, Yen et al. (2017) identified some signatures of infall. Their kinematical models allowed them to estimate the protostellar mass of 0.01 M_⊙, and the diskradius of 20 AU.

A.3 IRAS 16059-3857

IRAS 16059-3857 (also named Lupus 3 MMS) is an embedded protostar in the Lupus III Molecular Cloud (Comerón 2008). It was first discovered by Tachihara et al. (2007) using H¹³CO⁺(1–0) and 1.2 mm continuum observations. They point out the diffuse and faint cometary shape IR emission pointing toward the source, suggesting that this emission is a jet coming from a embedded young protostar and physically associated with IRAS 16059-3857 and the Herbig-Haro object HH 78. Lupus 3 MMS has been detected in the mid-infrared at 3.6–70 μm with Spitzer (Tachihara et al. 2007; Chapman et al. 2007; Merín et al. 2008) from which it has been classified as a Class 0 protostar (Tachihara et al. 2007; Dunham et al. 2008; Evans et al. 2009).

A.4 J160115-41523

The source SSTc2d J160115.6 415235.3 is a deeply embedded protostar in the Lupus IV Molecular Cloud at a distance of 150 pc (Comerón 2008), and is associated with 2MASS J16011549-4152351 located near the center of the globular filament GF 17 (Chapman et al. 2007). Through the infrared detections, it was classified as a flat-spectrum source (Merín et al. 2008). The source was also identified as a transition disk candidate by Ansdell et al. (2018), and the gas and dust masses of the surroundings of this young star have been determined in an ALMA disk survey carried out in Lupus previously (0.007 and 6 × 10⁻⁵ M_⊙ for its envelope and disk respectively, Ansdell 2017).

A.5 Sz 102

The source Sz 102 (Schwartz 1977), also known as Th 28 (The 1962), is associated with the Herbig-Haro object HH 228, which consists of an east–west jet emanating from the driving source, and the emission of the Herbig–Haro objects HH 228 E, HH 228 W, HH 228 E2, HH 228 E3, and HH 228 E4 (Krautter 1986; Graham & Heyer 1988; Comerón & Fernández 2011). Proper motions of >260 km s⁻¹ have been measured for HH 228 W, E, E2, E3, and E4 (Wang & Henning 2009; Comerón & Fernández 2011). The system of HH objects associated with Sz 102 appears to nearly lie on the plane of the sky (Comerón & Fernández 2011), which is also consistent with the faintness of the central source (probably due to the obscuration by a circumstellar disk seen nearly edge-on, Hughes et al. 1994). Finally, the jet presents hints of rotation (Coffey et al. 2004, 2007) and a remarkably high mass–loss rate in the outflow (~1.2 × 10⁻⁸ M_⊙ yr⁻¹, Coffey et al. 2008).

A.6 Merin 28

Merin 28 is associated with the 2MASS source J 16245177-3956326 and was classified as a flat spectral source (Spezzi et al. 2011). This object is not detected at optical wavelengths as expected for young objects still partially embedded, and appears point-like in all the 2MASS and Spitzer infrared images. The lack of detection in the HCO⁺(3–2) line observations casted doubts on the existence of a residual dense envelope seen in other sources with flat infrared spectral indices (van Kempen et al. 2009c).

Appendix B Tables

The table inthis appendix include information referred to the infrared data mentioned in Sect. 3.4

Appendix C Peculiarities of the individual outflows

In this appendix we describe some peculiar features of the outflows associated with IRAS 15398-3359 and IRAS 16059-3857, and intend to explain them.

Fig. C.1 Left: 2MASS Ks band image (gray-scale) and 12CO(2–1) emission associated with the blue and red lobes of IRAS 15398-3359 (blue and red contours, respectively). The source 2MASS 15430576-3410004, located in the left bottom corner in the image, could be driving the 12CO(2–1) emission detected southeast of IRAS 15398-3359. The red ellipses indicate possible ejection interacting with the gas in the vicinity of IRAS 15398-3359. The black circle indicates the region zoomed in the right panels. Right: velocity channel maps of the 12CO(2–1) emission (color scale) near the cloud velocity. SO(65, 54) and C18O(2–1) moment 0 maps are shown in black and pink contours, respectively.

C.1 IRAS 15398-3359

The source IRAS 15398-3359 has been very well studied (e.g., Jørgensen et al. 2013, Oya et al. 2014, Bjerkeli et al. 2016a, Yen et al. 2017). However we have detected molecular emission which is not related with the source nor with its previously-identified outflow.

In the right panels of Fig. C.1 we show the ¹²CO(J = 2–1) emission in three channels close to the cloud velocity. We also show the C¹⁸O(J = 2–1) and SO(J_N = 6₅, 5₄) moment 0 emission (integrated from 3.24 to 4.52 km s⁻¹). At these velocities a structure can be seen almost perpendicular to the main outflow associated with IRAS 15398-3359. This structure is also revealed by the SO(J_N = 6₅, 5₄) and SiO(4–3) emission (the latter is not shown here), indicating the possible presence of shocked gas. The latter two molecules are known as good tracers of physical shocks in the interstellar and circumstellar material, and are frequently detected associated with bipolar outflows of Class 0 and I protostars (e.g., Tafalla et al. 2010, Podio et al. 2015, Sakai et al. 2016). The C¹⁸O(J = 2–1) emission, which shows the disk associated with the central source and traces parts of the internal IRAS 15398-3359 outflow walls, also presents an elongation in the direction of this filament.

One possibility to explain this feature is the presence of another source located southeast of IRAS 15398-3359, throwing an outflow thatimpacts and excites the dense environment of IRAS 15398-3359. In the panels of Fig. C.1 the red ellipses indicate possible ejections coming from the southeast. The variations in the position of these shocks may be due to the source being precessing and therefore outflowing in different directions. These bullets could be even crossing the IRAS 15398-3359 outflow, being responsible for the lower CO emission intensity at the source position in the 4.28 km s⁻¹ velocity channel (bottom left panel in Fig. C.1), if the gas in these southern ejections had a lowertemperature. Moreover, the crossing of the southern outflow may explain the C¹⁸O elongation perpendicular to the IRAS 15398-3359 outflow. Although there may have been some contamination of the outflow originating from IRAS 15398-3359 our main analysis of that outflow is not affected.

Fig. C.2 IRAS 16059-3857. Up: 12 μm WISE emission. The star symbol indicates the position of the compact source associated with IRAS 16059-3857. The black contours shows the outflow emission. Bottom: moment maps of the shell-shape structures detected in the blue lobe. Black ellipses show the structures detected from the moment 0 map. Blue contours show the dss2 emission associated with HH 78.

We have also identified a possible source responsible for launching this second outflow: 2MASS 15430576-3410004, located about 69′′ south east of IRAS 15398-3359 (exactly at RA, Dec(J2000) 15:43:05.76, –34:10:00.4). The ellipses describing the tips of the observed shocks can all have the same origin at the position of this source.

C.2 IRAS 16059-3857

The blue lobe of the outflow associated with IRAS 16059-3857 in the vicinity of the protostar (<3500 AU) is V-shaped,but at greater distances it has two very particular features: on the one hand the shreds of CO emission south of the outflow cone, and on the other hand the bubble-shaped structures at the terminal part of the outflow. The bubble-shaped structures spatially coincide with some red-shifted emission (see upper left panel of Fig. 6) and with the Herbig–Haro object HH 78 (RA, Dec(J2000) = 16:09:12.8, –39:05:02). This HH object, can be observed in the Spitzer bands as a small (~4′′) faint nebulosity, mostly at 3.6 and 4.5 μm.

Comerón & Fernández (2011) proposed that HH 78 originates from a jet emerging from the original molecular cloud, since deep infrared observations of the region (Nakajima et al. 2003) showed that the visible portion of HH 78 can be traced as a well collimated jet up to the position of IRAS 16059-3857 (Tachihara et al. 2007), which is almost surely at the origin of this jet. The near-infrared part of the jet is bright in the K_S band, probably because of shock-excited H₂ emission coincident with the molecular outflow detected by Tachihara et al. (2007). This suggests that the jet would originate from a location deeply embedded in the cloud that hosts IRAS 16059-3857 and moves toward the side of the cloud that faces our direction, showing itself as the visible counterpart of HH 78 as it emerges. This scenario is in line with what is revealed by the present observations. The data suggests that the CO bubble-shaped structure seen at the western tip of the outflow is showing the outflow escaping from the molecular cloud in which the source is immersed.

In the upper panel of Fig. C.2 (WISE 12 μm emission), it can be seen that the outflow associated with IRAS 16059-3857 is immersed in a dark region or molecular pocket, at the center of which is forming the protostar. We see that the blue lobe tip coincides with the border of that dark region so we think that the structures detected there could be produced when the outflow reaches the limit of the dense pocket of gas and dust. The lower panels of the figure, show two elliptical structures clearly revealed by the moment 0 map. These structures spatially coincide with the DSS2 optical emission (shown in blue contours), produced at the position of HH 78. The two ellipses match the regions of highest velocity and dispersion, as observed in the moment 1 and 2 maps. The complex arc structures detected in these images could be due to the interaction of the molecular outflow with the edge of the dense pocket in which the system is embedded.

In Fig. C.3 we display the same ellipses over the velocity channel maps from 0.8 to 5.6 km s⁻¹, covering the blueshifted channels, the cloud velocity (4.6 km s⁻¹), and also part of the redshifted emission. The edges of the ellipses are detected in the extent of this velocity range, along with a very complex structure comprising arcs and filaments.

Fig. C.3 IRAS 16059-3857: velocity channel maps showing the structures in the blue lobe. Shell-like structures are detected coinciding with the position of the Herbig Haro object HH 78 (dss2 emission shown in magenta contours). The ellipses are the same as those in Fig. C.2.

We suggest that the high velocity gas emission (both at blue- and red-shifted velocities) at thes position of HH 78 is possibly due to physical interactions between different ejections from the protostellar system or between the jet and the walls of the molecular cloud.

Appendix D Complete velocity cubes

In this appendix we include the complete velocity cubes with 0.32 km s⁻¹ width channels, for IRAS 15398-3359 and IRAS 16059-3857, and 0.17 km s⁻¹ width channels for J160115-41523.

Fig. D.1 IRAS 15398-3359: blue lobs channel maps.

Fig. D.2 IRAS 15398-3359: red lobe channel maps.

Fig. D.3 IRAS 16059-3857: blue lobe channel maps.

Fig. D.4 IRAS 16059-3857: red lobe channel maps.

Fig. D.5 J160115-41523: blue and red lobes channel maps.

References

All Tables

All Figures

Fig. 1 Upper panels: continuum emission toward AzTEC-lup1-2, AzTEC-lup3-5 (ASTE, 1.1 mm), IRAS 15398-3359 and IRAS 16059-3857 (ALMA, 0.8 mm). Contours: 4, 10, 20 and 30 σ. Middle row panels: spectral energy distribution. A black-body could be fitted considering only the FIR and mm emission (black dashed line). The open blue dots corresponding to emission detected in the high resolution mm images are not included in the fit. The dust temperature derived from the fitting is shown in the top right corner. Bottom panels: line emission detected toward the sources at the position of the protostar.

In the text

	Fig. 1 (Cont.) Same as Figure 1 for J160115-41523, Sz 102 and Merin 28.
In the text

Fig. 2 12CO(2–1) integrated intensity images showing the blue-shifted and red-shifted emission for the outflows associated with IRAS15398-3359 (top panel, 12m+7m data), IRAS 16059-3857 (middle panel, 12m+7m data), and J160115-41523 (bottom panels, 7m data at left and 12m data at right). Each panel is labeled with the name of the source. The minimum and maximum velocities over which the emission is integrated are chosen based on visual inspection of the velocity channels (Table 8). Left panels: contours correspond to 10, 20, and 40 times the rms noise level (the rms noise for the full outflow moment 0 blue and red maps are 0.03–0.04, 0.04–0.05, and 0.2–0.08 Jy beam−1 km s−1 for IRAS 15398-3359, IRAS 16059-3857, and J160115-41523, respectively). Bottom right panel: 12CO high-angular resolutionintegrated intensity emission of the blue- and red-shifted lobes toward the central region of J160115-41523, indicated with a box in the bottom left panel. The contours correspond to 3, 5, and 7 times the rms noise levels (rms of the blue and red maps are 0.04 and 0.1 Jy beam−1 km s−1, respectively). The locations of the driving sources match the corresponding millimeter continuum emission, and are marked with black stars in the images. The estimated PA of each outflow lobe is indicated by black arrows. The beam sizes are shown as black ellipses in the lower left of each panel.

In the text

Fig. 3 Blue- and red-shifted emission associated with Sz 102 (blue contours: 2, 4, 6, 8 and 10 rms; red contours: 2, 3, and 4 rms) and Merin 28 (contours: 4, 8, and 12 rms). Plus signs in the top panel display the position of the Herbig-Haro objects HH 228 (over the DSS2 IR image), that lie out of the region mapped by the dedicated 7m array observations (black circle). The positions of the continuum emission are indicated with a black star.

In the text

Fig. 4 IRAS 15398-3359. Left: red lobe moment maps (from 5.9 to 16.7 km s−1). Right: blue lobe moment maps (from –10.2 to 4.8 km s−1). The systemic velocity is ~5.1 km s−1. The positionof HH 185 matches with the 2MASS Ks band peak (15:43:01.311, –34:09:14.81) (see also Heyer & Graham 1989) which is shown in black contours in the bottom right panel. Moment 1 and 2 maps only show the emission over 4 rms.

In the text

	Fig. 5 IRAS 15398-3359. Upper panel: position-velocity diagram of the 12CO(2–1) emission along the outflow axis with a cut width of 1′′. The arrowsshow the presence high velocity gas in both lobes. Bottom panels: position-velocity diagrams of the 12CO(2–1) emission across to the outflow axis at an angular distance of 4′′ and 8′′ from the central source in each lobe.
In the text

Fig. 6 IRAS 16059-3857. Left: red lobe moment maps (from 4.9 to 17.9 km s−1). The red-shifted emission spatially coinciding with the blue lobe. Right: blue lobe moment maps (from –4.9 to 4.28 km s−1). We note the emission that extends southward perpendicular to the lobe. The systemic velocity is ~4.7 km s−1. At the end, a bubble-shaped structure is detected. Over the moment 2 map the black dashed contours show theDSS2 emission (6300–6900 Å) that reveals the presence of the Herbig-Haro object HH 87, matching a region of high velocity dispersion. Moment 1 and 2 maps only show the emission over 4σ.

In the text

	Fig. 7 IRAS 16059-3857. Up: position-velocity diagram of the 12CO(2–1) emission along the red lobe outflow axis with a cut width of 1′′. Bottom: position-velocity diagrams of the 12CO(2–1) emission along 1′′-wide cuts perpendicular to the outflow axis at distances of 1000 AU (C1), 1500 AU (C2), 2000 AU (C3), 2500 AU (C4), 3000 AU (C5), 3500 AU (C6), 4000 AU (C7), and 4500 AU (C8).
In the text

Fig. 8 J160115-41523: moment maps of the red (left, from 4.1 to 6.7 km s−1) and blue (right, from –1.5 to 3.6 km s−1) lobes. The systemic velocity is ~4.1 km s−1. Moment 1 and 2 maps only show the emission over 4rms. The black dashed lines in the right upper panel show a 27° deflection in the blue lobe. The arrows point to the features described in the main text, where high velocity or dispersion emission is detected.

In the text

Fig. 9 J160115-41523. Upper panel: position-velocity diagram of the 12CO(2–1) emission along the outflow axis with a cut width of 30′′ (low resolution, full outflow). Arrows point out high velocity emission. Middle panel: position-velocity diagram of the 12CO(2–1) emission along the outflow axis with a cut width of 1′′ (high resolution, central region). We indicate the emission from the disk. A parabolic structure, indicated with a red dashed line,crossing the cloud velocity is detected in the red lobe. Bottom panels: position-velocity diagrams of the 12CO(2–1) emission across to the outflow axis at an angular distance of 2′′ and 5′′ from the central source in the red lobe, and 2′′ in the blue lobe.

In the text

	Fig. 10 IRAS 15398-3359: 12CO(2–1) emission velocity channel maps. The star indicates the position of the compact continuum source. Black ellipses show the four bipolar structures.
In the text

Fig. 11 Spitzer IRAC emission at 3.5 μm (top right), 4.5 μm (bottom left), and 8.0 μm (bottom right) associated with IRAS 15398-3359. Top left pannel: blue-shifted and red-shifted 12CO(2–1) emission. The black star indicates the position of the compact source. The infrared emission extending southwestward may stem from shocks of the blueshifted outflow close to the star. The spot in RA, Dec(J2000) = (15:43:01.3; −34:09:15.4) seen at 3.4 and 4.5 μm coincides with the terminal part of the blue lobe.

In the text

	Fig. 12 IRAS 16059-3857: 12CO(2–1) emission velocity channel maps. The star indicates the position of the compact continuum source. We identified elliptical structures labeled as e1 and e2.
In the text

	Fig. 13 J160115-41523: blue (up) and red (down) high resolution velocity channel maps.
In the text

Fig. C.1 Left: 2MASS Ks band image (gray-scale) and 12CO(2–1) emission associated with the blue and red lobes of IRAS 15398-3359 (blue and red contours, respectively). The source 2MASS 15430576-3410004, located in the left bottom corner in the image, could be driving the 12CO(2–1) emission detected southeast of IRAS 15398-3359. The red ellipses indicate possible ejection interacting with the gas in the vicinity of IRAS 15398-3359. The black circle indicates the region zoomed in the right panels. Right: velocity channel maps of the 12CO(2–1) emission (color scale) near the cloud velocity. SO(65, 54) and C18O(2–1) moment 0 maps are shown in black and pink contours, respectively.

In the text

	Fig. C.2 IRAS 16059-3857. Up: 12 μm WISE emission. The star symbol indicates the position of the compact source associated with IRAS 16059-3857. The black contours shows the outflow emission. Bottom: moment maps of the shell-shape structures detected in the blue lobe. Black ellipses show the structures detected from the moment 0 map. Blue contours show the dss2 emission associated with HH 78.
In the text

	Fig. C.3 IRAS 16059-3857: velocity channel maps showing the structures in the blue lobe. Shell-like structures are detected coinciding with the position of the Herbig Haro object HH 78 (dss2 emission shown in magenta contours). The ellipses are the same as those in Fig. C.2.
In the text

	Fig. D.1 IRAS 15398-3359: blue lobs channel maps.
In the text

	Fig. D.2 IRAS 15398-3359: red lobe channel maps.
In the text

	Fig. D.3 IRAS 16059-3857: blue lobe channel maps.
In the text

	Fig. D.4 IRAS 16059-3857: red lobe channel maps.
In the text

	Fig. D.5 J160115-41523: blue and red lobes channel maps.
In the text