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

# Orbital and Mass Constraints of the Young Binary System IRAS 16293-2422 A

María José Maureira¹, Jaime E. Pineda¹, Dominique M. Segura-Cox¹, Paola Caselli¹, Leonardo Testi², Giuseppe Lodato³, Laurent Loinard⁴,⁵, and Antonio Hernández-Gómez⁴,⁶,⁷

¹Max-Planck-Institut für extraterrestrische Physik (MPE), Gießenbachstr. 1, D-85741 Garching, Germany; maureira@mpe.mpg.de

²European Southern Observatory, Karl-Schwarzschild-Strasse 2, D-85748 Garching, Germany

³Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria 16, I-20133 Milano, Italy

⁴Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México Morelia, 58089, México

⁵Instituto de Astronomía, Universidad Nacional Autónoma de México, Apartado Postal 70-264, Ciudad de México 04510, México

⁶Max-Planck-Institute für Radioastronomie Auf dem Hügel 69, D-53121 Bonn, Germany

⁷IRAP, Université de Toulouse CNRS, UPS, CNES, Toulouse, France

Received 2020 March 30; revised 2020 May 22; accepted 2020 May 22; published 2020 July 2

## Abstract

We present 3 mm ALMA continuum and line observations at resolutions of 6.5 au and 13 au, respectively, toward the Class 0 system IRAS 16293-2422 A. The continuum observations reveal two compact sources toward IRAS 16293-2422 A, coinciding with compact ionized gas emission previously observed at radio wavelengths (A1 and A2), confirming the long-known radio sources as protostellar. The emission toward A2 is resolved and traces a dust disk with an FWHM size of ~12 au, while the emission toward A1 sets a limit to the FWHM size of the dust disk of ~4 au. We also detect spatially resolved molecular kinematic tracers near the protostellar disks. Several lines of the *J* = 5–4 rotational transition of HNCO, NH₂CHO, and t-HCOOH are detected, with which we derived individual line-of-sight velocities. Using these together with the CS (*J* = 2–1), we fit Keplerian profiles toward the individual compact sources and derive masses of the central protostars. The kinematic analysis indicates that A1 and A2 are a bound binary system. Using this new context for the previous 30 yr of Very Large Array observations, we fit orbital parameters to the relative motion between A1 and A2 and find that the combined protostellar mass derived from the orbit is consistent with the masses derived from the gas kinematics. Both estimations indicate masses consistently higher (0.5 ≲ *M*₁ ≲ *M*₂ ≲ 2 *M*₄) than previous estimations using lower-resolution observations of the gas kinematics. The ALMA high-resolution data provides a unique insight into the gas kinematics and masses of a young deeply embedded bound binary system.

*Unified Astronomy Thesaurus concepts: Star formation (1569); Protostars (1302); Close binary stars (254); Circumstellar dust (236); Circumstellar gas (238)*

## 1. Introduction

IRAS 16293-2422 (hereafter IRAS 16293) is a well-known bright ($L_{bol} \approx 21 L_{⊙}$; Jacobsen et al. 2018) low-mass protostellar system located in the Ophiuchus molecular cloud at a distance of 141 pc (Dzib et al. 2018). It has been widely studied because it was one of the first protostellar systems consistent with being in the Class 0 stage (Andre et al. 1993) and also because of its active “hot corino” chemistry with numerous complex molecules tracing compact regions with a high excitation (≥100 K) temperature (Ceccarelli et al. 2000; Caux et al. 2011; Jørgensen et al. 2016). Further, IRAS 16293 was one of the first Class 0 sources recognized as a multiple system. It was first resolved into two objects in centimeter observations showing source A to the south and source B to the north, separated by ~700 au or 5" (Wootten 1989). Later, Mundy et al. (1992) confirmed the protostellar nature of A and B by detecting compact dust thermal emission at 3 mm, coincident with the location of the centimeter sources. Source B shows a single peak from centimeter to submillimeter wavelengths, and it is believed to be at a very early stage of evolution (Pineda et al. 2012; Jørgensen et al. 2016; Hernández-Gómez et al. 2019). On the other hand, Wootten (1989) 2 cm Very Large Array (VLA) observations revealed early on two compact sources, A1 and A2 within A, separated by ~50 au and aligned approximately perpendicularly to the line connecting A and B. Later, Chandler et al. (2005) also resolved source A into two submillimeter peaks named Ab and

Aa separated by ~90 au and aligned similarly to A1–A2. The peak of Aa was located in between A1 and A2. Similar results were obtained with ALMA 1.3 mm observations with a resolution of 0.′′25 or 35 au (Sadavoy et al. 2018). By imaging the longest baselines in their observations, they recovered peaks oriented similar to the submillimeter peaks Ab and Aa identified in Chandler et al. (2005), with an additional weaker peak next to Aa, named Aa*. The lack of clear correspondence between the submillimeter and centimeter sources within source A prevented the confirmation of the nature of A1 and A2 as protostellar. For instance, Chandler et al. (2005) argued that A1 was a shock feature due to a precessing jet. This claim was based on the large proper motions of A1 with respect to A2 (Loinard 2002), and the shift in P.A. of A1 with respect to A2. On the other hand, Hernández-Gómez et al. (2019) recently argued that A1 is the location of a protostar, due to the nearly constant flux of A1 over time, inconsistent with the expectation of a shock feature. Despite this debate, the observation of possibly three outflows powered within A (Mizuno et al. 1990; Stark et al. 2004; Girart et al. 2014; van der Wiel et al. 2019) further supported the multiple protostellar nature of source A.

Sources A and B are embedded within a dense core with a mass of 4–6 $M_⊙$ enclosed within several 1000 au (Jacobsen et al. 2018; Ladjelate et al. 2020). Source B has a face-on configuration, and its mass has been constrained from 0.1 $M_⊙$ up to about 1 $M_⊙$ using interferometric observations at a resolution of 70–140 au (Pineda et al. 2012; Jacobsen et al. 2018; Oya et al. 2018). Source A shows a fattened disk-like
structure with a radius of about 100 au in observations at a resolution of 30–140 au (Jørgensen et al. 2016; Sadavoy et al. 2018). The same observations reveal a velocity gradient along the major axis of the disk-like structure, and the velocity profile was used to constrain the mass for source A, resulting in values around 1 $M_⊙$ (Oya et al. 2016; Jacobsen et al. 2018). On the other hand, Loinard et al. (2007) estimated a mass of 2–3 $M_⊙$ from assuming that A1 and A2 were two protostars whose relative motion was a circular orbit in the plane of the sky.

Here, we present ALMA Band 3 continuum observations with a resolution of 0.′046 (6.5 au) that reveal for the first time two compact sources at wavelengths tracing dust thermal emission, coincident with the location of the centimeter compact sources A1 and A2, thus confirming IRAS 16293 A1–A2 as a binary. Further, we present 3 mm line emission at a resolution of 13 au used to study the gas kinematics. This paper is organized as follows. In Section 2, we describe the observations and data reduction. In Section 3, we present the results of the continuum and line analysis. In Section 4, we discuss the implications of our results and analyze the positions of the new ALMA epoch and previous VLA observations and derived orbital parameters. Section 5 presents the summary and conclusions.

## 2. Observations

IRAS 16293 (A and B) was observed on 2017 October 8 and 12 using the most extended Cycle 5 configuration of ALMA (41.4 m–16.2 km baseline range) in Band 3 (single pointing) with a total time on-source of 1.25 hr. The maximum recoverable scale is ~0.′5 (or 70 au),<fnref n="8" /> and the pointing center of the observations was ICRS 16:32:22.63-24:28:31.8. The bandpass/flux calibrator and phase calibrator were J1517-2422 and J1625-2527, respectively. The observations were part of the cycle 5 project ID:2017.1.01247.S. (PI: G. Dipierro).

The spectral setup consisted of one continuum window centered at 99.988 GHz with 128 channels and a total bandwidth of 2 GHz and four windows of 960, 1920, 960, and 1920 channels with widths of 22.070 kHz centered at the frequencies of ¹³CO (1–0), C¹⁷O (1–0), C¹⁸O (1–0), and CS (2–1), respectively. CASA 5.4.1 (McMullin et al. 2007) was used for both calibration and imaging. Calibration of the raw visibility data was done using the standard pipeline. When imaging the continuum, we iteratively performed phase-only self-calibration with a minimum solution interval of 9 s. Afterwards we performed two amplitude self-calibration iterations, with a minimum solution interval of 60 s. The final continuum data set after phase+amplitude self-calibration was imaged using the tclean task with the multi-scale deconvolver and a robust parameter of 0.5. We tried different scales for the multi-scale imaging, and for the final image, we adopted those that resulted in the minimal residuals and no significant artifacts. The adopted values are four scales of 0, 8, 24, and 72 pixels (pixel size of 6 mas). These scales correspond to approximately 0 (point source), 1, 3, and 9 times the beam size. The beam size, beam P.A., and noise of the final continuum image are 0.′048 × 0.′046 (6.5 au), 79.3, and 15 µJy beam⁻¹, respectively.

The continuum self-calibration solutions were applied to all four line windows, to improve the signal-to-noise ratio (Brogan et al. 2018). Continuum subtraction was performed with CASA task uvcontsub by selecting line-free channels. We identified CS (J = 2–1) and several lines of the J = 5–4 rotational transition of HNCO, NH₂CHO, and t-HCOOH. These, along with ¹³CO (J = 1–0) were the brightest emission lines in the data. One extra bright line was also observed, but due to uncertainties with their identification, we do not use it for our analysis, and it is not shown here. C¹⁷O (J = 1–0) was undetected, and C¹⁸O (J = 1–0) showed only extended, very weak emission and thus was not used for the analysis. We cleaned CS (2–1), the J = 5–4 rotational transition of HNCO, NH₂CHO, and t-HCOOH, and ¹³CO (1–0) cubes using natural weighting, the multi-scale deconvolver with scales of 0, 5, and 15 pixels, a pixel size of 0.′018, a channel width of 0.38 km s⁻¹, and a uvtaper parameter of 0.′06 to better recover moderately extended features. Those scales correspond to 0 (point source), 1, and 3 times the beam. As with the continuum image, we picked those scales based on examination of the residuals. The cube of the CS (2–1) and ¹³CO (1–0) shows strong imaging artifacts (e.g., stripes) due to missing flux in the velocity range 3.24–5.9 km s⁻¹ (VLSR of the large-scale cloud). For the t-HCOOH, we produced an average cube using four sublevels of the J = 5–4 transitions that were isolated. Table A1 summarizes the lines used with their corresponding levels, sublevels, frequencies, and upper-level energies. The average final beam size, beam P.A., and noise for the lines are 0.′106 × 0.′084, 55.7, and 1.1 mJy/beam per channel.

## 3. Results

### 3.1. ALMA 3 mm Continuum Emission

Figure 1 shows the high-resolution 3 mm ALMA observations toward IRAS 16293. The 3 mm counterparts of the radio sources A1 and A2 are clearly detected. The compact emission toward A2 is resolved into an elongated disk-like structure while the compact and brighter emission toward A1 appears unresolved. The projected separation of A1 and A2 is 0.′38 or 54 au. The further improvement in sensitivity achieved through self-calibration resulted in further detection of weaker extended structures. A circumbinary disk-like structure with a semimajor axis of about ~0.′7 (100 au) and P.A. ~50° is observed, in agreement with previous lower-resolution observations (Jørgensen et al. 2016; Oya et al. 2016; Sadavoy et al. 2018). Overlaid on the circumbinary disk-like structure, there are also newly seen complex narrower features.

We performed 2D Gaussian fits on the image plane of the bright and compact continuum emission toward A1 and A2. We report here the results of fits with background subtraction since the residuals for this fit compared with one without background subtraction are substantially better (Figure B1). See Appendix B for details of the procedure and comparisons. Table 1 lists the results of the fit. Using the sizes in Table 1, the inferred inclinations assuming circular geometry are 59° ± 4 and 74° ± 1 from the plane-of-sky for A1 and A2, respectively. To obtain an estimation of the inclination of the circumbinary disk-like structure, we performed another 2D Gaussian fit. For this, we used three components, one for each of the compact emissions toward A1 and A2, and one for the extended structure. Assuming circular geometry, we inferred an inclination of 64° ± 1 for the circumbinary disk-like structure

<footnote n="8">⁸ From the ALMA Cycle 5 proposer’s guide. We note that this value is close to that obtained using 0.983× (wavelength/L<sub>s</sub>), with L<sub>s</sub> as the 5th percentile baseline in the observations. This is an empirically determined relation published in the ALMA Cycle 8 Technical Handbook.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** ALMA high angular resolution images of the multiple protostellar system IRAS 16293. Left: 3 mm continuum observations showing the triple nature of IRAS 16293. Source B to the north, hosts a single protostar, embedded within a 40 au across close to close-on disk-like structure. The flux from the two compact 3 mm sources to the south remains significant after removing contamination from free-free, unambiguously confirming the radio sources A1 and A2 as a binary protostellar system. Right: zoom-in view toward source A. The bright and compact (major axis <12 au) sources, separated by 54 au, likely correspond to two individual circumstellar disks. Extended dust structures surrounding the circumstellar disks are also revealed.</figcaption></figure>

with a P.A. of 50° ± 1. Although the inclinations inferred for both of the compact sources and the extended material agree within 10°, their P.A. are misaligned. The P.A. of the compact resolved emission toward A2 is about 138°, resulting in an almost ~90° misalignment with that of the circumbinary disk-like structure (P.A. ~ 50°). Although a similar misalignment is derived for A1 according to the results from the fit, the compact A1 emission is unresolved, making this measurement uncertain.

In the context of previous radio and submillimeter observations, it is likely that the two compact emissions come from small circumstellar disks. Previous VLA observations provide a spectral index at the low-frequency end, which, for sources A1 and A2, are consistent with free-free emission from an ionized jet. Hernández-Gómez et al. (2019) reported a spectral index from VLA observations of 0.5 ± 0.2 and 0.7 ± 0.2 for A1 and A2, respectively. From the most recent 7 mm observations with the VLA (epoch 2013) where A1 and A2 are clearly resolved (Hernández-Gómez et al. 2019), we estimated fluxes of 1.83 ± 0.08 and 2.24 ± 0.13 mJy for A1 and A2 at 3 mm. This results in free-free contamination of 37% ± 2% and 53% ± 3% for A1 and A2, respectively. These are conservative upper estimates of the free-free fluxes, since attempts to fit the compact sources with the sum of a point source and a Gaussian result in even less flux coming from the point source (taken as the unresolved free-free) than the above extrapolation (Appendix B). Thus, at least half of the flux in both sources comes from thermal dust emission tracing the location of two protostars.

The 3 mm compact emission toward A2 is perpendicular to the bipolar ejecta observed at centimeter wavelengths within 100 au (Loinard et al. 2007; Pech et al. 2010). Recent ~100 au resolution water maser observations also show blueshifted emission from A2, moving along the bipolar ejecta direction

(Dzib et al. 2018). The bipolar ejecta is also aligned with a 0.1 pc scale CO molecular outflow (Mizuno et al. 1990; Stark et al. 2004), but this molecular outflow has no clear counterpart below ~700 au (Yeh et al. 2008; Girart et al. 2014; van der Wiel et al. 2019), and the blueshifted lobe is at the opposite side of the blueshifted water maser emission (Dzib et al. 2018). Despite the confusion with the molecular outflow, the evidence from the ejecta and the flux coming from thermal dust emission support that the compact structure revealed in the 3 mm observations toward A2 is tracing a dust circumstellar disk with an approximate radius of 12 au (from the FWHM in Table 1).

The emission toward A1 is not well resolved. There are no previous observations of ejecta from A1, and this source also has not been unambiguously associated with an individual molecular outflow yet. However, besides the CO outflow mentioned above and possibly related to A2, there are two other CO outflows powered within A (Mizuno et al. 1990; Kristensen et al. 2013; Girart et al. 2014; van der Wiel et al. 2019): one oriented east-west and one pointing toward B. We note that although it is expected than one of them is being driven by A1, none of these CO outflows are aligned with the inferred (although unresolved) minor axis of the A1 dust disk. Our observations constrain the size of the dust circumstellar disk around A1 to ≲3.5 au (from the FWHM in Table 1).

Finally, we note that the previous super-resolution images at frequencies >200 GHz that revealed the peaks Ab and Aa, not matching the location of A1 and A2 (particularly A2; Chandler et al. 2005; Chen et al. 2013; Sadavoy et al. 2018), were likely affected by the high optical depth of the surrounding material, which prevented the clear detection of the embedded compact sources. Our observations show that the location of Ab and Ab from Sadavoy et al. (2018) is near the substructures around A1, while Aa*, the weakest additional peak next to Aa identified in Sadavoy et al. (2018), is near A1 (see Figure B2).
<table><caption>**Table 1**

Fit to the 3 mm Compact Sources A1 and A2</caption><thead><tr><th>Source</th><th>R.A.<br>(J2000)</th><th>Decl.<br>(J2000)</th><th>Deconvolved Size<br>(arcsec)</th><th>P.A.<br>(°)</th><th>Peak Flux Density<br>(mJy beam<sup>-1</sup>)</th><th>Integrated Flux<br>(mJy)</th></tr></thead><tbody><tr><td>A1</td><td>16:32:22.878</td><td>-24:28:36.684</td><td>(24.6 ± 0.8, 12.8 ± 1.2)</td><td>119.8 ± 3.6</td><td>4.18 ± 0.02</td><td>4.91 ± 0.04</td></tr><tr><td>A2</td><td>16:32:22.851</td><td>-24:28:36.647</td><td>(87.3 ± 2.5, 23.7 ± 2.0)</td><td>137.9 ± 1.0</td><td>1.73 ± 0.03</td><td>4.23 ± 0.11</td></tr></tbody></table>

### 3.1.1. *Masses from 3 mm Continuum Emission*

Gas masses from the 3 mm continuum are commonly estimated using:

$$M = 100 \frac{d^2 S_{\nu}}{B_{\nu}(T_d) \kappa_{\nu}} \quad (1)$$

where $S_ν$ is the integrated flux density, $B_ν$ is the Planck function, $κ_ν$ is the dust opacity, and $d$ is the distance. Equation (1) assumes optically thin emission (Hildebrand 1983). $S_ν$ is taken from Table 1 and corrected by the free-free contamination. We assume a dust temperature of 100 K (Oya et al. 2016; van't Hoff et al. 2020) and a gas-to-dust ratio of 100. For the dust opacity, we adopt a range of values between 0.23 cm$^2$ g$^{-1}$, appropriate for dense material at a very young evolutionary state (Ossenkopf & Henning 1994; Dunham et al. 2016), and 0.81 cm$^2$ g$^{-1}$, corresponding to dense material at a more evolved state (Agurto-Gangas et al. 2019). For obtaining an estimate of the mass of the extended emission toward source A, we use the integrated flux of all emission above $3\sigma$ after subtracting the contribution from the compact sources A1 and A2. We obtained a total integrated flux from the extended emission of ~73 mJy. The range of masses calculated with the two opacities for A1, A2 and the extended emission correspond to $(1-3) \times 10^{-3} M_⊙$, $(1-3) \times 10^{-3} M_⊙$, and 0.03–0.1 $M_⊙$, respectively. There are important additional uncertainties on these values arising from uncertainties in the dust temperature, optical depth, and as previously mentioned, the true contribution from free-free emission. All of these factors contribute independently and can decrease or increase the reported values within a factor of a few (Ballering & Eisner 2019). Further, given that at the scales of the observed circumstellar disks, dust scattering can decrease the intensity at millimeter wavelengths, which are likely also optically thick, our estimates for the mass of the compact sources should be taken as conservative lower limits (Liu 2019; Ueda et al. 2020). The order of magnitude of the lower limits reported here is comparable to other compact ($\lesssim$40 au) circumstellar disks in Class I multiple systems observed at high ($\lesssim$25 au) resolution, derived in a similar fashion (Takakuwa et al. 2017; Alves et al. 2019; Cruz-Sáenz de Miera et al. 2019).

### 3.2. Molecular Lines

Figures 2 and 3 show the moment 0 and moment 1 maps for CS (2–1), HNCO (5–4), NH$_2$CHO (5–4), and t-HCOOH (5–4) molecular lines. The moment 0 and 1 maps are integrated over the same velocity ranges. In addition, only pixels with emission >3$\sigma$ per channel were considered. The velocity range for CS was split into two to avoid channels with artifacts (due to missing flux). The two ranges correspond to [−7.02, 2.86] km s$^{-1}$ and [6.28, 12.74] km s$^{-1}$. For HNCO and t-HCOOH, the

velocity range is [−7.02, 12.72] km s$^{-1}$, while for NH$_2$CHO, we restricted the range to [−0.18, 12.74] km s$^{-1}$ since at blueshifted velocities, the channels also showed emission from an unrelated line. Similar maps for $^{13}$CO(1-0) are shown in Appendix Figure F1.

The moment 0 maps for HNCO (5–4), NH$_2$CHO (5–4), and t-HCOOH (5–4) show a similar distribution (Figure 2). Their integrated intensities are both brighter around the location of A1, and more compact than the CS (2–1) emission. Except for the CS emission toward A1, none of the tracers seem to peak at the location of A1 or A2. The moment 1 maps in Figure 3 show that the ALMA observations presented here resolve the known velocity gradient toward source A (Pineda et al. 2012; Oya et al. 2016; van't Hoff et al. 2020), approximately aligned with the major axis of the circumbinary disk-like structure. High-velocity blueshifted and redshifted components are clearly revealed in the CS (2–1) map toward A1 and A2, respectively. Higher-velocity components can also be identified in the other tracers, except for NH$_2$CHO, where the restricted velocity range used to avoid line contamination resulted in a removal of the high-velocity blueshifted component (which is also present in this line).

Interestingly, the location where the integrated intensity of our molecular tracers is enhanced (particularly HNCO, NH$_2$CHO, and t-HCOOH) coincides with the location of the dust substructures detected in the continuum emission (black contours in Figure 2). Similarly, outside the region in between the sources, the location of the higher-velocity material in the moment 1 maps coincides with the location of the dust substructures. The region covered by the dust substructures and the most compact molecular lines extends up to 20–40 au from the protostar positions. This is around the lower limit of the centrifugal barrier radius inferred for source A in Oya et al. (2016), from observations with a resolution of ~70 au. In their model, the centrifugal barrier corresponds to a radius at which the infalling of material suddenly stops due to the conservation of angular momentum, within which a rotationally supported disk is expected to form. The infall velocity peaks at the centrifugal radius (twice the size of the centrifugal barrier). From there, rotational motions dominate, with increasing velocity, down to the centrifugal barrier. Thus, higher velocities due to a combination of infall and rotation are expected in the region marking the transition from the inner envelope to the circumstellar disk. Although in a binary system the gas kinematics are expected to be more complex than in that simple model, accretion shocks and structures related to the transition from the circumbinary material to the circumstellar disks are also expected. In simulations, spiral structures connecting the circumbinary material to the circumstellar disks are observed (Matsumoto et al. 2019). Further, Mösta et al. (2019) find that these spiral structures can take the form of complex tightly wound features depending on the eccentricity of the system, the mass ratio, and the specific orbital phase. That scenario has been given as a possible explanation for the curved filamentary features observed within the circumbinary
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Moment 0 maps toward IRAS 16293 A (color). The velocity range for the CS (2–1) was split into two to avoid channels with artifacts due to missing flux. The two ranges were [−7.02, 2.86] km s⁻¹ and [6.28, 12.74] km s⁻¹. For HNCO and t-HCOOH, the velocity range was [−7.02, 12.72] km s⁻¹, while for NH₂CHO, we restricted the range to [−0.18, 12.74] km s⁻¹ since at blueshifted velocities, the channels also showed emission from a different line. The black contours show the 3 mm continuum emission at levels 124, 320, and 448 µJy to identify the circumbinary disk-like structure, the compact sources A1 and A2, and the smaller-scale substructures around them. The crosses mark the peak location for A1 and A2 (Table 1). The beam is shown in the bottom left corner of each panel.</figcaption></figure>

material of the binary Class I system BHB2007-11 (Alves et al. 2019), where the filamentary substructures are revealed at the same scales as the substructures in IRAS 16293 A. Thus the spatial correlation between the dust substructures and the enhancement of the integrated intensity and higher velocities of our observed lines is consistent with these features tracing shocks, spirals, or more complex features associated with the transfer of material from the circumbinary ambient into the circumstellar disks.

Figure 3 also shows a clear almost 90° misalignment between the major axis of the circumstellar disk around A2 and the direction of the velocity gradient of the circumbinary material. Although such misalignment between the rotation axis of the circumbinary material and that of the region close to the protostars was expected given the known ~90° difference in the P.A. of the bipolar outflows/ejection and the rotation axis toward source A down to 70 au scales (Mizuno et al. 1990; Loinard et al. 2007; Pineda et al. 2012; Girart et al. 2014; Oya et al. 2016; van der Wiel et al. 2019), the higher-resolution observations show that this misalignment persists down to the smallest scales resolved in our line observations (~13 au). This type of misalignment has been seen before for close protostellar binaries, for example, the Class I system IRS 43 with a separation of ~74 au (Brinch et al. 2016) in which the individual circumstellar disks were found to be significantly misaligned (>60°), in inclination and P.A. Further, the orbital plane of the binary was constrained to be oriented close to face-on, while the circumbinary material was oriented close to edge-on. Misaligned configurations for the rotation axis of individual circumstellar disks, the circumbinary material, and the orbital motion naturally arise in simulations where turbulence is included in the star-forming cloud (Offner et al. 2010; Bate 2018; Lee et al. 2019). For instance, members of a multiple system might form a few thousands of astronomical unit apart, from gas with different angular momentum, and later move closer to form a bound tight binary (or higher multiplicity) system (Offner et al. 2016; Bate 2018; Kuffmeier et al. 2019; Lee et al. 2019). Misalignment can also be the product of binary formation in an elongated structure whose minor axis is misaligned with the initial rotation axis (Bonnell et al. 1992). Finally, subsequent accretion of material with a misaligned angular momentum can also explain the misalignment between the compact dust disks and the surrounding rotating material, as observed toward IRAS 16293 A (Bate 2018).

### 3.3. Line-of-sight Velocities

We use the HNCO (5–4), NH₂CHO (5–4), and t-HCOOH (5–4) lines, which do not show significant artifacts from missing large-scale structures to extract the velocity of A1 and A2 along the line of sight. To get a robust estimation of the source velocities, we use position–velocity (p–v) diagrams at a direction along the velocity gradient (Figure 3) around each source. The adopted P.A. corresponds to 65° and 30° for A1 and A2, respectively. The width of the p–v cuts was set to match the beam (5 pixels). See Figure E2 for a diagram of the cuts overlaid on the moment 1 maps. The procedure to obtain the line-of-sight velocity consisted of fitting the p–v diagrams at intermediate velocities with a linear gradient and extracting the velocity at the position of the source. For the linear gradient fit, in each channel along the p–v cut, we fit a Gaussian to the emission and then fit the peaks of the Gaussians. The velocity range that was fit with a linear gradient corresponding to 0.9–4 km s⁻¹ and 4.7–8 km s⁻¹ for A1 and A2, respectively. The selected range covers the region where the compact source is located in the p–v diagram as well as the region in which a single linear velocity gradient is observed. Figure 4 shows the p–v diagrams for A1 and A2. The colored lines show the gradient fit for each molecular line. The final line-of-sight velocity is given as the average and its associated error among the molecular lines for each source. Table 2 summarizes the results. We obtained line-of-sight systemic velocities of 2.1 ± 0.1 km s⁻¹ and 5.8 ± 0.1 km s⁻¹ for A1 and A2, respectively. We repeated the procedure using p–v diagrams with P.A. values differing from the previous one by ±10°. The resultant velocities are in agreement within uncertainties.

Our estimation of the line-of-sight velocity for A2 has the caveat that it uses the emission from tracers that were brighter around A1. If these tracers are not tracing closely the material near A2, the line-of-sight velocity might be different from the estimated value. Since the high-velocity redshifted component approaching A2 in Figure 3 starts at about 6 km s⁻¹, similar to the redshifted CO emission from outflow lobes toward source A (Mizuno et al. 1990), the calculated line-of-sight velocity of 5.8 km s⁻¹ is consistent with an upper limit. On the other hand, Dzib et al. (2018) observed a water maser at a line-of-sight velocity of 2.1 km s⁻¹ at a location consistent with A2 and moving along the direction of one of the CO molecular outflows. Thus we can consider the velocity of this water maser
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Moment 1 maps toward IRAS 16293 A (color). The velocity range for the CS (2–1) was split into two to avoid channels with artifacts due to missing flux. The two ranges were [−7.02, 2.86] km s⁻¹ and [6.28, 12.74] km s⁻¹. For HNCO and t-HCOOH, the velocity range was [−7.02, 12.72] km s⁻¹, while for NH₂CHO, we restricted the range to [−0.18, 12.74] since at blueshifted velocities, the channels also showed emission from a different line. The black contours show the 3 mm continuum emission at levels 124, 320, and 448 µJy to identify the circumbinary disk-like structure, the compact sources A1 and A2, and the smaller-scale substructures around them. The crosses mark the peak location for A1 and A2 (Table 1). The beam is shown in the bottom left corner of each panel.</figcaption></figure>

emission as a strict lower limit to the line-of-sight velocity for A2. This lower limit would result in both sources having the same line-of-sight velocity.

### 3.4. Mass Constraints from Gas Kinematics

To investigate the velocity profile of the material near A1 and A2, we use *p*-v diagrams along the gradient observed in the moment 1 maps. Although these gradients are misaligned with the major axis of the circumbinary disk-like material (P.A. ~ 50°) and that of the circumstellar disk around A2 (P.A. ~ 140°), such misalignment does not rule out that the high-velocity components closest to the protostars might be tracing the Keplerian rotation of the gas circumstellar disks with misaligned axes (e.g., Jensen & Akeson 2014), a scenario that we test in the following sections.

Figure 5 shows CS (2–1) and HNCO (5–4) *p*-v diagrams toward A1 and A2. The P.A. of the cuts are the same as those used for obtaining the line-of-sight velocities (see Figure 4) and correspond to direction along the velocity gradient at the location of each source (Figure F2). These directions are 15° and 20° different from the major axis of the circumbinary disk-like material, for A1 and A2, respectively. The HNCO (5–4) structures in the *p*-v maps are similar to those of NH₂CHO (5–4) and t-HCOOH (5–4; Figure 4). In this work, we focused on the gas kinematics close to the individual sources and located outside of the region in between the sources. This is because the region between the sources is likely more affected by the interaction between them. The CS (2–1) emission in the outer regions shows a resolved profile of increasing velocity toward A1 and A2 (in the bottom left quadrant for A1 and upper right quadrant for A2 in Figure 5). In the following, we analyze individually the velocity profiles observed in Figure 5 for A1 and A2.

#### 3.4.1. A1

The CS (2–1) high-velocity blueshifted component toward A1 can be identified down to a velocity of ~−8 km s⁻¹ with the velocity peak located within 1 beam (≈0.09 or 13 au) of the protostar position. On the other hand, the emission from the other lines represented in Figure 5 by the HNCO (5–4) show a linear gradient around the protostellar location. The morphology of HNCO (5–4) can be explained by emission arising from rotating gas concentrated at a narrow range of distances from the protostar (e.g., Lindberg et al. 2014; Yen et al. 2014).

This is because if $v_{\text{rot}}$ is the rotational velocity of the rotating region at a radius $R$, then we expect the velocity across the major axis to change as $V_{\text{rot}} \times \Delta_{\text{offset}}/R$, with $\Delta_{\text{offset}}$ the distance to the source location. As discussed in Section 3.2, the emission from HNCO might be tracing shocks or spiral features located at the transition from the circumbinary ambient into the circumstellar disks. The HNCO blueshifted velocity peak is located farther away from the protostar than that of the CS. The latter could then be tracing the Keplerian motion of the circumstellar disk around A1. This is also in agreement with the interpretation of Keplerian motion of the H₂CS (7–6) in Oya et al. (2016). This tracer showed similar (although unresolved) high-velocity emission to the CS (2–1) in this work. This high-velocity emission was inconsistent with the velocities of the farther away infalling and rotating material traced by other species.

Given the resolved velocity profile of the blueshifted emission, we use the “upper edge” method (Seifried et al. 2016) to extract the line-of-sight velocity of the gas as a function of the distance from A1. Seifried et al. (2016) used simulations and showed that having a resolution of about 15 au is important to be able to identify a Keplerian profile using the “upper” edge of the emission in a *p*-v diagram. Thus, our observations are adequate for applying this method to determine if the profiles are consistent with Keplerian rotation. We find that the velocity profile is indeed consistent with Keplerian rotation (i.e., $v \propto r^{-0.5}$) up to a distance of 20 au from A1. The extracted points beyond this distance depart from a Keplerian power law, in agreement with a steeper profile closer to $v \propto r^{-1}$. The data points that are consistent with a Keplerian power law and the corresponding fit are shown in Figure F3. The Keplerian profile fit assumes the systemic line-of-sight velocities in Table 2. The resultant Keplerian curve is overlaid on the *p*-v diagram in Figure 5. The solid part of the curve shows the regions used for the Keplerian fit, while the dashed lines are an extrapolation of the fit. From the Keplerian fit to the gas, we obtain a mass of $0.8 \pm 0.04 M_{\odot}$ for A1. This value is already corrected by the inferred inclination of the extended structure that the CS (2–1) line (and extended dust) is tracing (i.e., using $i = 64^{\circ}$), as we do not resolve and/or detect line emission directly associated with the small dust circumstellar disk. We note that the inclination is only slightly different from that inferred from the compact dust emission ($i = 59^{\circ}$) and would change the resultant mass only by 10%.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** P-v diagrams of observed molecular lines centered around A1 (top) and A2 (bottom). The orientation of the cut for each source is indicated in the bottom right corner (see Figure F2 for the orientation). The vertical black lines mark the position of A1 (top) and A2 (bottom), respectively. The inclined solid lines show the resultant linear gradient fit to each molecular line with which we obtain a line-of-sight velocity for each molecular line and source. The horizontal black lines show the average line-of-sight velocity of 2.1 and 5.8 km s⁻¹ calculated using the three molecular lines for A1 (top) and A2 (bottom). All of the panels show the same velocity range and identical spatial scales.</figcaption></figure>

<table><caption>**Table 2**
Line-of-sight Velocities</caption><thead><tr><th></th><th>HNCO (5-4)<br>(km s<sup>-1</sup>)</th><th>NH<sub>2</sub>CHO (5-4)<br>(km s<sup>-1</sup>)</th><th>t-HCOOH (5-4)<br>(km s<sup>-1</sup>)</th><th>Mean<br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>A1</td><td>2.12 ± 0.07</td><td>2.30 ± 0.03</td><td>1.86 ± 0.04</td><td>2.09 ± 0.10</td></tr><tr><td>A2</td><td>5.90 ± 0.26</td><td>5.84 ± 0.21</td><td>5.70 ± 0.16</td><td>5.81 ± 0.05</td></tr></tbody></table>

As mentioned above, the method that we used to extract the velocity profile results in larger reported uncertainties of the masses (of a few tens of percent), when tested with simulations (Seifried et al. 2016). Thus, we adopt here a larger, 30%, uncertainty (based on Table 6 of Seifried et al. 2016), which results in an individual mass of 0.8 ± 0.3 $M_\odot$. By using a cut along the velocity gradient, we minimize the contamination by infall motions, which nevertheless are likely present because the material is not oriented edge-on. We expect that the contamination by infall motions would change the mass within the large adopted uncertainty (Seifried et al. 2016). As with the

line-of-sight velocity analysis with *p*-v diagrams, we repeated the Keplerian fit procedure using *p*-v diagrams with P.A. values differing from the previous one by ±10°. The resultant masses were in agreement within the uncertainties of 30%. We note that although the ¹³CO (1–0) also shows a profile of increasing velocities toward A1 (Figure E1), the observations do not recover emission close to the source, and thus, ¹³CO (1–0) does not trace the regions where the profiles were consistent with Keplerian rotation for the CS (2–1).

Although we find that the velocity profile of the CS emission within 20 au is consistent with a Keplerian power law, we cannot rule out that this high-velocity material is tracing gas that is infalling from the edge of the circumbinary disk-like structure (traced by the HNCO) due to, for example, a loss of angular momentum. In this case, the mass inferred by the Keplerian fit would be overestimating the mass. Using a simple infall and rotation model in Section 3.4.2 for A2, we show that the overestimation factor is similar to the adopted 30% uncertainty in the Keplerian mass.
<figure><img src="image_6.png" /></figure>

**Figure 5.** P–v diagrams of the CS (2–1) centered around A1 (left) and A2 (right). The orientation of the cut for each source is in the top left corner and coincides with the orientation of the straight lines in Figure F2. The vertical black lines mark the position of the sources, while the horizontal black lines mark the systemic velocity for each source (Table 2). The blue and red curves show the Keplerian rotation fit to the gas using only the bottom left quadrant and upper right quadrant for A1 and A2, respectively. The solid portion of the curve indicates the range of distances from the central source that were consistent with Keplerian rotation and hence used for the Keplerian profile fit. The dashed lines are an extrapolation of the solid curve Keplerian fit. The resultant masses are displayed in each case, corrected by inclination with *i* = 64° (inclination of the extended emission surrounding the circumstellar disks, traced by the molecular lines and dust). Contours show the HNCO (5–4) counterpart emission (same as in Figure 4). In both panels, the contours start at 3*σ* and end at 7*σ*, increasing in steps of 1*σ* = 1.4 mJy beam⁻¹. Both panels have the same velocity range and spatial scales.

3.4.2. A2

Figure 5 shows that the high-velocity redshifted component can be identified up to a velocity of ~13 km s⁻¹ for the CS, with the velocity peak located ~0.′′14 or 20 au from the protostar. The emission from HNCO (5–4) is weaker near A2 compared with A1, but a linear gradient can also be identified (Figure 4). High-velocity redshifted gas, reaching velocities similar to the CS (2–1) are also observed in the HNCO (5–4). We follow the procedure in Section 3.4.1 to investigate if the profile could be consistent with a Keplerian power law. In addition, since the velocity peaks of CS, HNCO, and the location of the dust substructure to the southwest of A2 coincide, the high-velocity peak in A2 might be associated with gas that is infalling and rotating from the circumbinary disk-like structure to the circumstellar disk. Thus, to estimate the mass in that scenario, we also use a simple model of infall and rotation (Appendix C), similar to Oya et al. (2016).

Following the procedure done for A1 (Section 3.4.1), we
find that the velocity profile of the CS toward A2 in Figure 5 is
consistent with Keplerian rotation (i.e., *v* ∝ *r*⁻⁰.⁵) up to a
distance of 50 au. The extracted points beyond these distances
depart from a Keplerian power law, in agreement with a steeper
profile closer to *v* ∝ *r*⁻¹. The data points that are consistent
with a Keplerian power law and the corresponding fit are
shown in Figure F3. The resultant Keplerian curve is overlaid
on the *p*–*v* diagram in Figure 5. We obtain a mass of 1.4 ±
0.4 *M*☉. Similar to A1, we corrected by an inclination of
*i* = 64° and assumed an uncertainty of 30% (Seifried et al.
2016). We also repeated the Keplerian fit procedure using *p*–*v*
diagrams with P.A. values differing from the previous one by
±10°. The resultant masses were in agreement within the

uncertainties of 30%. For this source, ¹³CO (1–0; Figure F1)
shows a velocity profile similar to that of CS (2–1).

For the combination of rotation and infall scenario, we
compared the velocity profile with a curve of a simple model of
rotating material undergoing infall with conservation of angular
momentum described in detail in Appendix C. This simple
model is similar to the one used to provide the size of a
centrifugal barrier for IRAS 16293 A by Oya et al. (2016), as
well as to describe in general the kinematics outside of the
Keplerian rotation region in other protostellar sources (Sakai
et al. 2014; Yen et al. 2014). By assuming that the rotation axis
of the material is perpendicular to the major axis of the
circumbinary disk-like structure, we find that a centrifugal
barrier of 20 au and mass of 0.9 M<sub>☉</sub> can reasonably reproduce
the velocity profile along the same cut used for the Keplerian fit
(see Figure 6). This profile assumes the same inclination as that
used for the Keplerian fit (i = 64°). A smaller centrifugal
barrier can also reproduce the profile if the central mass is
increased, while a larger centrifugal barrier cannot reproduce
emission that is as close to the protostar as the one observed in
Figure 6. We note that the values of the masses inferred in the
two cases agree under the adopted uncertainty of a few tens of
percent.

**4. Discussion**

An important conclusion from our kinematics analysis is that
the point-source masses are consistently larger than previous
estimations with five to six times lower resolution, in which the
location of the point sources was not resolved (Pineda et al.
2012; Oya et al. 2016; Jacobsen et al. 2018). They derived a
mass of ~1 M☉, assuming a single source. This is in agreement
with our results within the uncertainties. However, recent work
<figure><img src="image_7.png" /><figcaption>**Figure 6.** P-v diagrams of the CS (2-1) and HNCO (5-4) centered around A2, the same as the right panel in Figure 5. The purple curve shows the velocity obtained considering infall with rotation with conservation of angular momentum around a central mass of 0.9 M<sub>☉</sub>, and a centrifugal barrier of 20 au (where the innermost part of the purple profile ends). The curved is corrected by inclination and projection effects considering that the rotation axis is perpendicular to the major axis of the circumbinary material (P.A. = 50°). The kinematics of an infalling and rotating flow of gas can also (in addition to pure Keplerian motion) explain the velocity profile toward A2. In both cases, the central mass agrees within the uncertainties.</figcaption></figure>

comparing observations and synthetic observations of the kinematics of the circumbinary and circumstellar disks for the binary system L1551 NE further support our finding of higher protostellar mass for IRAS 16293 A than previously estimated. L1551 NE is a Class I system at a similar distance, showing a similar separation and inclination as IRAS 16293 A. Takakuwa et al. (2017) found maximum molecular line velocity differences (with respect to the protostellar system velocity) of about 4 km s<sup>-1</sup>, while the protostars' line-of-sight velocity difference is 1.3 km s<sup>-1</sup>, compared to about 8 km s<sup>-1</sup> and 3.7 km s<sup>-1</sup>, respectively, for IRAS 16293 A. Takakuwa et al. (2017) find that a simulation with a total protostellar mass of 0.8 M<sub>☉</sub> is in agreement with the observed kinematics. By comparison, the higher velocity differences observed toward IRAS 16293 A would imply a higher mass, consistent with our results. Future comparisons with synthetic observations as in Takakuwa et al. (2017) will help to further constrain the individual masses.

### 4.1. Are A1 and A2 Bound?

For observed embedded multiple systems with separations <100 au, the usual assumption is that they are bound (Tobin et al. 2016). For this close binary pair, we have measured their motions, and thus we can estimate the required total mass for the pair to be bound. We use the most recently published proper motions for A1 and A2 (Hernández-Gómez et al. 2019) that used data spanning almost 30 yr, along with our new line-of-sight velocity to obtain the total velocity for each source. The measured individual line-of-sight velocities in Section 3.3 lead to a relative line-of-sight difference between A1 and A2 of 3.7 km s<sup>-1</sup>. The total relative velocities (A1–A2) on the line of sight and plane-of-sky are -3.7 ± 0.2 km s<sup>-1</sup> and 5.2 ± 0.6 km s<sup>-1</sup>, respectively, yielding a total relative velocity magnitude of 6.4 ± 0.5 km s<sup>-1</sup>. The

minimum mass for A1 and A2 to be a bound pair follows from the following condition:

$$E_{\text{kin}} + E_{\text{pot}} \le 0 \qquad (2)$$

where the kinetic and potential energy are given by:

$$E_{\text{kin}} = \frac{1}{2} \sum_{1,2} M_i (v_i - v_{\text{com}})^2 \qquad (3)$$

$$E_{\text{pot}} = \frac{-GM_1 M_2}{|\mathbf{r}_1 - \mathbf{r}_2|} \qquad (4)$$

$M_i$ corresponds to the mass of source i, $v_i$ is the velocity of source i, $v_{\text{com}}$ is the velocity of the center of mass, and $|\mathbf{r}_1 - \mathbf{r}_2|$ is the separation between the sources. The minimum total mass can be thus calculated as:

$$\frac{1}{2} \frac{(v_1 - v_2)^2}{G/|\mathbf{r}_1 - \mathbf{r}_2|} \le M_1 + M_2. \qquad (5)$$

Assuming that the point sources lie in the plane of the circumbinary disk-like material for which we calculated a P.A. of 50° and inclination with respect to the plane-of-sky of 64° (Section 3.1), the deprojected distance between A1 and A2 is 0.678 or 95.6 au. This results in a minimum mass of 2.2 ± 0.3 M<sub>☉</sub>. As discussed in Section 3.3, the line-of-sight velocity difference between A1 and A2 might be smaller. If we consider the line-of-sight velocity difference to be zero, we obtain a strict minimum mass of 1.7 ± 0.3 M<sub>☉</sub> for the system to be bound. In Section 3.4, we used the gas kinematics and estimated a combined A1+A2 mass of 2.2 ± 0.5, suggesting that the pair is indeed bound.

### 4.2. Motion of the Protostars

Further insight into the dynamics of this now bound binary system can be obtained through the study of their proper motion. We examined all VLA positions preceding the ALMA measurement and reported in Hernández-Gómez et al. (2019). We selected from this work (and references therein) only those observations for which by eye the two sources were resolved and free from ejecta contamination. These criteria resulted in 12 VLA observations with the epochs listed in Table D1. Ejecta emission from A2 has routinely been seen in VLA < 22 GHz observations since 2006 (Loinard et al. 2007). As a result of this, all selected observations after 2005 correspond to frequencies >33 GHz. The left panel in Figure 7 shows our selected absolute positions for A1 and A2 through time, overlaid on the most recent ALMA 3 mm continuum observation. A1 has moved faster in the plane of the sky than A2, resulting not only in a change in orientation but also in the separation between the sources. The separation has changed from 47 au in the first observation to 54 au in the most recent one. These changes are clearly seen in the right panel in Figure 7 showing the relative position of A1 with respect to A2 through time. We note that the relative positions are not affected by differences in the absolute astrometry accuracy of the observations. The change in time of the relative positions is important since if A1 and A2 are orbiting around their center of mass, the relative trajectory of A1 with respect to A2 is also described by a Keplerian orbit, the parameters of which also provide the total mass of the system, including contributions from both protostars+compact disks, and the gas mass enclosed within the orbit in the case of this young embedded system. The relative motions in Figure 7 further suggest that
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Motions of A1 and A2 over a 30 yr period. All observations prior to 2017 correspond to VLA observations recently compiled in Hernández-Gómez et al. (2019), and in which by visual inspection the two sources were resolved and free from ejecta contamination. Left: absolute motions of A1 and A2 overlaid on the ALMA 3 mm observation. Right: relative motions of A1 with respect to A2 overlaid on the ALMA 3 mm observation. When errors are not visible, it is because they are smaller than the symbol.</figcaption></figure>

they are tracing a bound orbital trajectory. In agreement with this expectation, Hernández-Gómez et al. (2019) tried a quadratic fit to the relative A1–A2 positions. Although that parabolic trajectory provided a better fit to their data (particularly their most recent data), our ALMA observation is in disagreement with the prediction of their quadratic fit. The quadratic fit predicts a reversal in the decl. motion of A1 with respect to A2 around 2015, indicating that A1 started to move north (see their Figure 7), while our new ALMA data point shows that A1 kept moving south (Figure 7). Such discrepancy might be due to the use of jet contaminated and/or low-resolution data, which strongly overwhelms the subtle astrometric variations due to the orbital motions and that were discarded in our analysis. Despite our careful selection, we note that the errors in Figure 7 are only those from the Gaussian fit to the compact sources and thus do not take into account errors produced by unresolved ejecta emission that can temporarily displace the observed center from the true center. As data earlier than 2010 had typically lower resolution (Table D1), it is likely that the differences between positions at closely separated epochs or among different epochs in the early VLA data (e.g., epoch 1989 in Figure 7) are due to this temporary “wobble” from unresolved ejecta.

To further investigate the relative A1–A2 trajectory in light of the new ALMA observations, we first calculated the best ballistic trajectory using our selected VLA epochs. Then, we compared the predicted position along this trajectory with the observed ALMA position. Figure 8 shows the relative positions and the calculated ballistic trajectory (straight line). The new ALMA epoch is in disagreement with the prediction from the ballistic trajectory. Given this and our previous conclusion that the system is likely bound (Section 4.1), we also fit all relative positions (VLA+ALMA) with Keplerian orbital trajectories using a Monte Carlo approach (see Appendix D for details of the fitting procedure and posterior distributions). Figure 8 shows the most likely Keplerian orbital trajectories along with the ballistic line. The predictions for all epochs along both trajectories (ballistic and Keplerian orbit) are also marked. Both trajectories have similar predicted positions before the ALMA

observation. If we recalculate the best ballistic trajectory using all epochs, we obtain slightly better residuals for the Keplerian orbital trajectory. Yet, the Keplerian orbital trajectory has seven free parameters, while the ballistic trajectory only has four. Thus, although the current data do not allow us to securely rule out the ballistic trajectory, the results from the gas kinematics provide independent evidence in support of the Keplerian trajectory. Future monitoring observations are key to further strengthen this conclusion while also allowing for improved accuracy on the Keplerian orbit determination.

The full posterior distributions, medians, and confidence intervals of the orbital parameters are given in Appendix D and displayed in Figure D1. The resultant orbits have a period of $P = 362_{-73}^{+133}$ yr and semimajor axis of $a = 80.26_{-8.81}^{+14.60}$ au. The inclination also seems to be well constrained ($i = 58.69_{-3.82}^{+3.39}$) and is similar to the inclination of the circumbinary gas and dust emission. The other angles needed to fully define the orientation of the orbit will require further observations in order to obtain accurate constraints. This is because tests in which a couple of epochs were not considered (2005.2 and 2003.5) provided different values for these parameters, and also for the eccentricity ($e = -0.19_{-0.06}^{+0.09}$), while the rest remained consistent. The total mass derived from the Keplerian orbital trajectory is $M_{\text{tot}} = 4_{-1}^{+1} M_{\odot}$. This derived total mass from the orbit is in agreement with previous estimations assuming a simplified plane-of-the-sky circular orbit (Chandler et al. 2005; Loinard et al. 2007; Pech et al. 2010); although, in these studies the protostellar nature of A1 had not been confirmed. Since the presence of gas is also contributing to the total derived mass, we can take this value along with its large uncertainty as an upper limit to the combined mass of the point sources. In addition, we can also provide an independent upper limit using the luminosity of source A. This source has been modeled by Jacobsen et al. (2018), although only as a single source, and they found that a luminosity of $18 L_{\odot}$ resulted in good agreement with their observations. Since a $2 M_{\odot}$ pre-main-sequence star has an approximate luminosity of $\sim 10 L_{\odot}$ at the birthline (Stahler & Palla 2005), we can set this value as the upper limit to the most massive protostar. Thus the combined
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Relative motions of A1 with respect to A2 (black markers with errors). The straight line shows the best ballistic trajectory using only the VLA data. The square symbols overlaid on the straight line show the predicted positions for all observations along the ballistic trajectory. The purple curves show 100 bound Keplerian orbital solutions drawn from posterior distributions (Figure B1). The dot symbols along the orbital trajectory show the predicted positions for all observed positions along a randomly selected Keplerian orbital trajectory. Lines connecting the observed positions with the predicted ones along the randomly selected Keplerian orbit are also shown for clarity. When errors are not visible, it is because they are smaller than the symbol.</figcaption></figure>

evidence, from the gas kinematics, stellar kinematics, and luminosity, results in individual protostellar masses reasonably constrained in the range $0.5 \lesssim M_1 \lesssim M_2 \lesssim 2 M_\odot$. The order of the previous relation takes into account that relative to source B, A1 has moved significantly more than A2 (Loinard 2002; Chandler et al. 2005; Pech et al. 2010). These constraints result in a mass ratio between ~0.3 and 1.

Further constraints on the mass ratio could be obtained by also deriving the Keplerian orbital parameters of the center of mass of source A with respect to source B, which are also consistent with being bound (see Appendix D), forming a triple hierarchical system. However, the separation between A and B is an order of magnitude larger than the separation between A1 and A2, resulting in possible orbits between the center of mass of A relative to B of several $10^3-10^4$ yr. Thus we cannot constrain this orbit with the current 30 yr of observations.

**5. Summary and Conclusions**

We use high-resolution continuum (6.5 au resolution) and line (13 au resolution) 3 mm ALMA observations toward the Class 0 multiple system IRAS 16293-2422. In this work, we analyzed the southern source in this system (IRAS 16293 A). In addition, we use observations from the VLA covering a period of 30 yr (Wootten 1989; Hernández-Gómez et al. 2019) to review the motion of the compact sources within IRAS 16293 A. Our results can be summarized as follows:

1. The two radio sources A1 and A2 are unambiguously detected in the 3 mm continuum observations, and a considerable fraction of the flux in both sources is

consistent with thermal dust emission. Thus, the 3 mm observations confirm the protostellar nature of both sources (in particular of A1), which remained debated due to conflicting results between previous >200 GHz and <43 GHz VLA observations. The peaks (Aa, Aa*, and Ab) previously seen in >200 GHz observations were likely affected by optical depth, which prevented the clear detection of the embedded compact sources (particularly A2). Some of these peaks are tracing substructures in the extended emission instead.

2. The compact emission toward A2 is resolved and is consistent with a dust circumstellar disk with an FWHM size of 12 au, oriented perpendicular to the previously observed bipolar ejecta at centimeter wavelengths as well as perpendicular to the disk-like circumbinary dust emission. The compact emission toward A1 is unresolved, setting a limit to the FWHM size of the dust circumstellar disk of 3.6 au.

3. Complex substructures extending from 20 to 40 au from the protostars are also observed. They are associated with regions where the emission of several lines of the $J = 5-4$ rotational transition of HNCO, NH$_2$CHO, and t-HCOOH is enhanced. Similarly, they are associated with regions where these tracers as well as the CS (2-1) show higher velocities. Thus, these substructures might be tracing shocks or spiral features at the transition from the circumbinary structure into the circumstellar disks.

4. We use the compact emission from HNCO (5-4), NH$_2$CHO (5-4), and t-HCOOH (5-4) to estimate individual line-of-sight velocities for A1 and A2, yielding
a line-of-sight velocity difference of 3.7 km s⁻¹. The CS (2–1) traces clear high-velocity emission associated with the positions of the protostars. The velocity profiles from the locations of the sources and up to a few tens of astronomical unit toward the outskirts are consistent with a Keplerian power law. The Keplerian power law implies individual masses of 0.8 ± 0.3 M☉ and 1.4 ± 0.4 M☉ for A1 and A2, respectively. The velocity profiles can also be explained by material that is rotating and infalling from the circumbinary disk-like structure to the circumstellar disks, resulting in smaller although comparable masses.

5. The most recent reported proper motions from VLA observations, our newly measured line-of-sight velocities, and protostellar masses from the gas kinematics indicate that the binary system A1–A2 is bound.

6. The new positions from the ALMA observations depart from the predicted position along a ballistic trajectory inferred from the VLA observations, suggesting the observation of an orbital trajectory. We fit orbital parameters to the relative positions of the VLA+ALMA observations resulting in orbital solutions with a period of 362<sup>+133</sup><sub>-73</sub> yr, semimajor axis of 80.26<sup>+14.60</sup><sub>-8.81</sub> au, and inclination consistent with that of the extended circumbinary material. The results also indicate a low eccentricity (e = 0.19<sup>+0.09</sup><sub>-0.06</sub>), but future observations are needed to better constrain the geometry of the orbit. The total mass derived from the orbital fit is M<sub>tot</sub> = 4<sup>+1</sup><sub>-1</sub> M☉. The independent mass constraints from the gas kinematics and the stellar kinematics are in agreement within the uncertainties, which, when added to luminosity restrictions, result in individual masses reasonably constrained in the range 0.5 ≲ M<sub>1</sub> ≲ M<sub>2</sub> ≲ 2 M☉.

The range of protostellar masses inferred from the orbital analysis and the gas kinematics is consistently higher than previous estimations using lower-resolution observations of the gas kinematics or models with a single source. Given the current mass of the IRAS 16293 A and B envelope of 5 ± 1 M☉ at scales of a few 1000 au, the binary system A and single source B are also likely bound, forming a triple hierarchical system. Future monitoring observations, as well as detailed modeling with simulations, will help to further constrain the dynamics and individual masses of this deeply embedded triple system.

The authors thank Bo Zhao and Kedron Silsbee for insightful theoretical discussions on binary formation and orbits. L.L. acknowledges the financial support of DGAPA, UNAM (project IN112417), and CONACYT, México. This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 823823 (DUSTBUSTERS). The authors also thank Emanuel Caux for discussions about the use of CASSIS.

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.01247.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.

<table><caption>**Table A1**
Properties of the Observed Transitions</caption><thead><tr><th>Molecule</th><th>Transition</th><th>Frequency (GHz)</th><th><i>E</i><sub>up</sub> (K)</th></tr></thead><tbody><tr><td>CS</td><td>2–1</td><td>97.980953</td><td>7</td></tr><tr><td>HNCO</td><td>5<sub>3,2</sub>–4<sub>3,1</sub></td><td>109.833487</td><td>391</td></tr><tr><td>HNCO</td><td>5<sub>3,3</sub>–4<sub>3,3</sub></td><td>109.833487</td><td>391</td></tr><tr><td>NH<sub>2</sub>CHO</td><td>5<sub>1,4</sub>–4<sub>1,3</sub></td><td>109.753578</td><td>19</td></tr><tr><td>t-HCOOH</td><td>5<sub>2,4</sub>–4<sub>2,3</sub></td><td>112.287144</td><td>29</td></tr><tr><td>t-HCOOH</td><td>5<sub>4,1</sub>–4<sub>4,0</sub></td><td>112.432319<sup>a</sup></td><td>67</td></tr><tr><td>t-HCOOH</td><td>5<sub>3,3</sub>–4<sub>3,2</sub></td><td>112.459621</td><td>45</td></tr><tr><td>t-HCOOH</td><td>5<sub>3,2</sub>–4<sub>3,1</sub></td><td>112.467007</td><td>45</td></tr><tr><td><sup>13</sup>CO</td><td>1–0</td><td>110.201354</td><td>5</td></tr><tr><td>C<sup>17</sup>O<sup>b</sup></td><td>1–0</td><td>112.358982</td><td>5</td></tr><tr><td>C<sup>18</sup>O<sup>c</sup></td><td>1–0</td><td>109.782173</td><td>5</td></tr></tbody></table>

**Notes.** All of the cubes were imaged with a channel width of 0.38 km s⁻¹. Frequencies and transitions were based on the full CDMS and JPL catalogs available within the CASSIS software (developed by IRAP-UPS/CNRS http://cassis.irap.omp.eu).

<sup>a</sup> There is an additional transition (t-HCOOH 5<sub>4,2</sub>–4<sub>4,1</sub>) separated by 0.07 km s⁻¹, and thus unresolved in our observations.

<sup>b</sup> Undetected.

<sup>c</sup> Detected, emission is very weak and extended and thus was not used in this study.

## Appendix A
Molecular Line Properties

Table A1 summarizes the frequencies and upper energy levels of the transitions analyzed in this work. This list does not include all of the lines in the data, as there are several weaker lines, which were not suitable for a kinematical analysis and whose identification was beyond the scope of this study.

## Appendix B
Fits to the Compact Emission toward A1 and A2

We carried out 2D Gaussian fits in the image plane of the bright and compact continuum emission toward A1 and A2. We performed fits with and without background subtraction. For the fit without any background subtraction, we use the CASA task imfit, restricting the fitting region to a square of size 0"/18. For the fit with background subtraction, we extracted a square sub-image around each source with a size of 0"/20. We then made a spline fit of the background with the SciPy task bisplrep, masking out the compact emission. The mask for source A1 consisted of a circular region of radius 0"/07, while A2 consisted of an ellipse following the orientation of the compact emission with semimajor and semiminor axes of 0"/11 and 0"/07, respectively. We checked that the background fit shows no hints of compact emission at the center. We then subtracted the background fit from the data and performed a 2D Gaussian fit with imfit without any restriction. Figure D1 shows the data, the fit with background subtraction, and the residuals of the fit with and without background subtraction for both sources. The positions of the peak and the peak flux do not change with the type of fit while the size and integrated flux are 1.5–2× smaller when the fit is done after background subtraction. Since the residuals for the fit with background subtraction are substantially better, we list the results of the background subtraction fit in Table 1, and we use these fluxes for the mass estimation. The positional errors are ~1 × 10⁻⁴ arcseconds and ~7 × 10⁻⁴ arcseconds for A1 and A2, respectively.
<figure><img src="image_10.png" /><figcaption>**Figure B1.** Observed 3 mm compact emission and comparison of residuals from the Gaussian fit with and without background subtraction for A1 (top) and A2 (bottom). From left to right, the panels show the observed 3 mm emission, the Gaussian fit to the compact emission plus the background fit, the residuals between the observed emission and the Gaussian plus background fit, and the residuals when doing the Gaussian fit without background subtraction. In all panels, the x-axis and y-axis are in units of pixels. Contours in the fourth and fifth columns correspond to σ levels of −7, −5, −3, 3, 5, and 7. Negative σ levels are shown as dashed contours.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure B2.** Zoom-in view toward source A in the ALMA 100 GHz observations. The bright and compact emission correspond to the protostellar sources A1 and A2 (epoch 2017.9). The crosses mark the positions of the previously detected peaks using lower-resolution (>$18$ au) observations at frequencies >200 GHz (Chandler et al. 2005; Chen et al. 2013; Sadavoy et al. 2018). The emission at frequencies >200 GHz is optically thick and does not let one see the peaks tracing the compact sources embedded (particularly A2), which are clearly revealed in the 100 GHz observations, even when degrading the resolution to match previous observations. The positions for Aa, Ab, and Aa* correspond to those reported in Sadavoy et al. (2018; epoch 2015.5) recalculated with respect to the peak of source B to account for proper motion between the epochs (although small ~0''/01). While Ab and Aa are tracing the substructures on the extended disk-like structure, Aa* (the weakest among the three) peaks near A1. However, the Aa* peak location (to the south of A1 in the figure) is inconsistent with the trajectories discussed in Section 4.2 and the motions in Figure 7, and thus not a reliable tracer of A1.</figcaption></figure>

We also attempted a two-component fit (Gaussian+point source) to the compact 3 mm emission toward A1 and A2, after background subtraction. For A1, all of the properties of the Gaussian except for the peak flux and integrated flux remain consistent with the previous single Gaussian fit. The peak and

integrated flux are reduced by 25% and 30%, respectively. This is similar to the free-free contamination (~37%) estimated using the A1 radio spectral index obtained from VLA observations (Hernández-Gómez et al. 2019). For A2, the resultant peak and integrated fluxes for the point source were consistent with zero, while the free-free contamination using the A2 radio spectral index (Hernández-Gómez et al. 2019) was ~53%. Figure B2 shows the location of compact submillimeter sources Aa, Aa*, and Ab taken from Sadavoy et al. (2018) overlaid on the ALMA 3 mm observations presented in this work.

## Appendix C
### Simple Velocity Profile of Infalling and Rotating Gas

We use the simple model in which the gas motion is approximated by particles moving in a plane under the influence of a central object of mass M. The particles are considered to have a constant angular momentum (perpendicular to the plane of motion) along the trajectories, resulting in a centrifugal barrier of radius:

$$r_{CB} = \frac{j^2}{2GM}, \qquad (C1)$$

where j is the specific angular momentum. This is equivalent to the trajectories calculated by Ulrich (1976), if we consider only the motion in the plane perpendicular to the rotation axis (equatorial plane). The infall $v_{inf}$ and rotational $v_{rot}$ velocities are given by:

$$v_{inf} = \sqrt{\frac{2GM}{r} - \frac{j^2}{r^2}}, \qquad (C2)$$

$$v_{rot} = \frac{j}{r}, \qquad (C3)$$

where r is the distance to the protostar calculated as

$$r = \sqrt{x^2 + y^2} \qquad (C4)$$

$$x = d \times \Delta_{\text{off}} \cos(\text{P.A.} - \text{P.A.'}) \qquad (C5)$$

$$y = d \times \Delta_{\text{off}} \sin(\text{P.A.} - \text{P.A.'}) / \cos(i), \qquad (C6)$$
where d is the distance to the source, $\Delta_{off}$ is the offset in arcseconds along the p-v cut with a position angle P.A., and i is the inclination with respect to the plane-of-sky. P.A.' corresponds to the position angle of the major axis of the equatorial plane. Then, the final velocity profile $v_{los}$ on the p-v map will be given by:

$$v_{los} = v_{inf} \frac{y}{r} + v_{rot} \frac{x}{r} \quad (C7)$$

## Appendix D: Orbit Fitting

We use the open-source software package orbitize! (Blunt et al. 2020), which uses a parallel-tempered Markov Chain Monte Carlo algorithm (Foreman-Mackey et al. 2013; Vousden et al. 2016) to fit orbits using both positional and line-of-sight velocity observations (Table 2). We selected from Hernandez-Gomez et al. (2019 and references therein) only those observations for which by eye the two sources were resolved

<figure><img src="image_12.png" /><figcaption>Figure D1. Corner plot showing the posterior distribution of the orbital parameters from top to bottom: semimajor axis in au, eccentricity $e$, inclination $i$ in degrees, argument of periastron $\omega$ in degrees, longitude of ascending node $\Omega$ in degrees, epoch of periastron passage $\tau$ (measured as fraction of orbit compared to a reference date; see Appendix D), and total mass $M_{tot}$ in $M_{\odot}$. The red lines indicate the median for each parameter, and the dashed lines correspond to symmetric confidence intervals around the median and enclosing 68% of the probability.</figcaption></figure>
<table><caption>**Table D1**
Selected VLA Observations for Analysis of Protostellar Motion</caption><thead><tr><th>Epoch (year)</th><th>Frequency (GHz)</th><th>Synthesized Beam</th><th>References</th></tr></thead><tbody><tr><td>1987.7</td><td>15</td><td>0"19 × 0"09</td><td>(1), (2), (3)</td></tr><tr><td>1989.1</td><td>8</td><td>0"34 × 0"19</td><td>(4), (2), (3)</td></tr><tr><td>1989.1</td><td>22</td><td>0"18 × 0"09</td><td>(4), (2), (3)</td></tr><tr><td>1994.3</td><td>8</td><td>0"34 × 0"16</td><td>(5), (2), (3)</td></tr><tr><td>2003.5</td><td>43</td><td>0"09 × 0"05</td><td>(2), (3), (6)</td></tr><tr><td>2003.7</td><td>8</td><td>0"39 × 0"19</td><td>(2), (3), (7)</td></tr><tr><td>2005.2</td><td>43</td><td>0"30 × 0"17</td><td>(7), (3)</td></tr><tr><td>2011.2</td><td>41</td><td>0"30 × 0"14</td><td>(3)</td></tr><tr><td>2011.4</td><td>41</td><td>0"13 × 0"10</td><td>(3)</td></tr><tr><td>2011.4</td><td>41</td><td>0"08 × 0"05</td><td>(3)</td></tr><tr><td>2012.9</td><td>33</td><td>0"10 × 0"04</td><td>(3)</td></tr><tr><td>2013.0</td><td>41</td><td>0"09 × 0"04</td><td>(3)</td></tr></tbody></table>

**References.** (1) Wootten (1989); (2) Chandler et al. (2005); (3) Hernández-Gómez et al. (2019); (4) Mundy et al. (1992); (5) Loinard (2002); (6) Rodríguez et al. (2005); (7) Loinard et al. (2007).

and free from ejecta contamination by visual inspection. Table D1 summarizes the epoch, frequency, resolution, and references for the selected VLA observations. We input the 12 VLA positions and the single ALMA positions and radial velocity measurements. To assess possible priors, we first explore the parameter space with a grid fitting method (Köhler et al. 2008). The final implemented priors are Gaussian with mean 100 au and standard deviation of 50 au for the semimajor axis *a*, uniform between 0 and 0.85 for eccentricity *e*, LogUniform with a minimum of 2 $M_\odot$ and a maximum of 6 $M_\odot$ for the total mass $M_{\text{tot}}$, uniform between 40° and 80° for the inclination *i*, uniform between 0° and 360° for both the argument of periastron $\omega$ and the longitude of ascending node $\Omega$, and Gaussian with a mean of 0.5 and standard deviation of 0.3 for the epoch of periastron passage $\tau$. The latter is expressed in orbitize! as a fraction of the orbital period past a specified reference date $t_{\text{ref}}$ (default 2020 January 1) and thus with possible values between 0 and 1. These priors are only weakly informed. We note that changing all priors to uniform within comparable ranges does not significantly affect our results. The distance was fixed to 141 pc. We ran 20,000 steps with 1000 walkers per temperature with 20 temperatures.

We removed 10,000 steps as burn-in. The resultant median values and confidence intervals from the posterior distribution are $a = 80.26_{-8.81}^{+14.60}$ au, $e = 0.19_{-0.06}^{+0.09}$, $i = 58.69_{-3.82}^{+3.39}$, $\omega = 214.90_{-57.15}^{+41.27}$, $\Omega = 315.20_{-3.61}^{+2.66}$, $\tau = 0.29_{-0.13}^{+0.73}$, and $M_{\text{tot}} = 3.93_{-0.80}^{+1.09} M_\odot$. The full posterior distributions, medians, and confidence intervals are displayed in Figure D1.

## Appendix E
### Bound Analysis A–B

Following the method in Section 4.1, we first use the most recently published proper motions for B (Hernández-Gómez et al. 2019) and previous determinations of its line-of-sight velocity (Pineda et al. 2012), to estimate a minimum total mass $A + B$ for the triple system (A1, A2, and B) to be bound. The proper motions alone show that A1 has moved significantly with respect to A2 and B ($5 \pm 1$ km s⁻¹ in both cases). On the other hand, A2 has not moved significantly with respect to B ($1.5 \pm 1$ km s⁻¹), suggesting that the center of mass of A1–A2 is located closer to A2 (i.e., mass ratio A1/A2 < 1). Using a mass ratio between 0 and 1 and the analysis in Section 4.1 results in a total minimum mass between ~2 and 7 $M_\odot$, for IRAS 16293 to be a bound triple. Source B has a mass close to 1 $M_\odot$ between its circumstellar disk and protostellar masses<fnref n="9" /> (Pineda et al. 2012; Oya et al. 2018), and source A has a combined mass of at least 1 $M_\odot$, while the mass in the large-scale envelope around the three sources is about 4–6 $M_\odot$ (Jacobsen et al. 2018; Ladjelate et al. 2020). Thus it is reasonable to conclude that B is also bound to A1 and A2, forming a hierarchical triple system.

## Appendix F
### Additional Figures

This section contains additional figures discussed in the main text. Figure F1 shows the moment 0 and 1 of ¹³CO(1–0). Figure F2 shows the direction and approximate width of the p–v cuts used to obtain the individual line-of-sight velocities for A1 and A2 (Figure 4), as well as an estimation of their masses (Figures 5 and 6). Figure F3 shows the data points extracted from the p–v maps in Figure 5, using the “upper-edge” method (Seifried et al. 2016). Overlaid are the resultant Keplerian power-law fits (solid lines).

<footnote n="9">⁹ We note that source B is face-on and extremely optically thick within the recoverable scales of our data; thus, we could not perform the same kinematic analysis as in source A.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure F1.** $^{13}$CO (1–0) moment 0 (left) and 1 (right) maps toward IRAS 16293 A (color). The velocity range for the $^{13}$CO (1–0) was split into two to avoid channels with artifacts due to missing flux. The two ranges were [−7.02, 2.86] km s$^{-1}$ and [6.28, 12.74] km s$^{-1}$, the same as CS (2–1) in Figures 2 and 3. The black contours show the 3 mm continuum emission at levels 124, 320, and 448 µJy to identify the circumbinary disk-like structure, the compact sources A1 and A2, and the smaller-scale substructures around them. The crosses mark the peak location for A1 and A2 (Table 1). The beam is shown in the bottom left corner of each panel.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure F2.** Moment 1 maps toward IRAS 16293 A (color), the same as in Figure 3. The black straight lines show the direction and width of the *p*–*v* cuts used for the gas kinematic analysis (Sections 3.3 and 3.4). The P.A. of the cuts correspond to 65° and 30° for A1 and A2, respectively.</figcaption></figure>
<figure><img src="image_15.png" /></figure>

**Figure F3.** CS (2–1) velocity profiles extracted using the upper edge method (Seifried et al. 2016) for A1 (left) and A2 (right). The symbols correspond to the data points, in agreement with a Keplerian power-law profile. Errors are equal to the channel width (0.38 km s⁻¹). The lines show the Keplerian power-law fit to the data, resulting in the mass constraints discussed in Section 3.4.

## ORCID iDs

María José Maureira https://orcid.org/0000-0002-7026-8163

Jaime E. Pineda https://orcid.org/0000-0002-3972-1978

Dominique M. Segura-Cox https://orcid.org/0000-0003-3172-6763

Paola Caselli https://orcid.org/0000-0003-1481-7911

Leonardo Testi https://orcid.org/0000-0003-1859-3070

Giuseppe Lodato https://orcid.org/0000-0002-2357-7692

Laurent Loinard https://orcid.org/0000-0002-5635-3345

Antonio Hernández-Gómez https://orcid.org/0000-0001-7520-4305

## References

Agurto-Gangas, C., Pineda, J. E., Szűcs, L., et al. 2019, A&A, 623, A147

Alves, F. O., Caselli, P., Girart, J. M., et al. 2019, Sci, 366, 90

Andre, P., Ward-Thompson, D., & Barsony, M. 1993, ApJ, 406, 122

Ballering, N. P., & Eisner, J. A. 2019, AJ, 157, 144

Bate, M. R. 2018, MNRAS, 475, 5618

Blunt, S., Wang, J. J., Angelo, I., et al. 2020, AJ, 159, 89

Bonnell, I., Arcoragi, J.-P., Martel, H., & Bastien, P. 1992, ApJ, 400, 579

Brinch, C., Jørgensen, J. K., Hogerheijde, M. R., Nelson, R. P., & Gressel, O. 2016, ApJL, 830, L16

Brogan, C. L., Hunter, T. R., & Fomalont, E. B. 2018, arXiv:1805.05266

Caux, E., Kahane, C., Castets, A., et al. 2011, A&A, 532, A23

Ceccarelli, C., Castets, A., Caux, E., et al. 2000, A&A, 355, 1129

Chandler, C. J., Brogan, C. L., Shirley, Y. L., & Loinard, L. 2005, ApJ, 632, 371

Chen, X., Arce, H. G., Zhang, Q., et al. 2013, ApJ, 768, 110

Cruz-Sáenz de Miera, F., Kóspál, Á, Ábrahám, P., Liu, H. B., & Takami, M. 2019, ApJL, 882, L4

Dunham, M. M., Offner, S. S. R., Pineda, J. E., et al. 2016, ApJ, 823, 160

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

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Girart, J. M., Estalella, R., Palau, A., Torrelles, J. M., & Rao, R. 2014, ApJL, 780, L11

Hernández-Gómez, A., Loinard, L., Chandler, C. J., et al. 2019, ApJ, 875, 94

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

Jacobsen, S. K., Jørgensen, J. K., van der Wiel, M. H. D., et al. 2018, A&A, 612, A72

Jensen, E. L. N., & Akeson, R. 2014, Natur, 511, 567

Jørgensen, J. K., van der Wiel, M. H. D., Coutens, A., et al. 2016, A&A, 595, A117

Köhler, R., Ratzka, T., Herbst, T. M., & Kasper, M. 2008, A&A, 482, 929

Kristensen, L. E., Klaassen, P. D., Mottram, J. C., Schmalzl, M., & Hogerheijde, M. R. 2013, A&A, 549, L6

Kuffmeier, M., Calcutt, H., & Kristensen, L. E. 2019, A&A, 628, A112

Ladjelate, B., André, P., Könyves, V., et al. 2020, arXiv:2001.11036

Lee, A. T., Offner, S. S. R., Kratter, K. M., Smullen, R. A., & Li, P. S. 2019, ApJ, 887, 232

Lindberg, J. E., Jørgensen, J. K., Brinch, C., et al. 2014, A&A, 566, A74

Liu, H. B. 2019, ApJL, 877, L22

Loinard, L. 2002, RMxAA, 38, 61

Loinard, L., Chandler, C. J., Rodríguez, L. F., et al. 2007, ApJ, 670, 1353

Matsumoto, T., Saigo, K., & Takakuwa, S. 2019, ApJ, 871, 36

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, CASA Architecture and Applications, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Mizuno, A., Fukui, Y., Iwata, T., Nozawa, S., & Takano, T. 1990, ApJ, 356, 184

Mösta, P., Taam, R. E., & Duffell, P. C. 2019, ApJL, 875, L21

Mundy, L. G., Wootten, A., Wilking, B. A., Blake, G. A., & Sargent, A. I. 1992, ApJ, 385, 306

Offner, S. S. R., Dunham, M. M., Lee, K. I., Arce, H. G., & Fielding, D. B. 2016, ApJL, 827, L11

Offner, S. S. R., Kratter, K. M., Matzner, C. D., Krumholz, M. R., & Klein, R. I. 2010, ApJ, 725, 1485

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

Oya, Y., Moriwaki, K., Onishi, S., et al. 2018, ApJ, 854, 96

Oya, Y., Sakai, N., López-Sepulcre, A., et al. 2016, ApJ, 824, 88

Pech, G., Loinard, L., Chandler, C. J., et al. 2010, ApJ, 712, 1403

Pineda, J. E., Maury, A. J., Fuller, G. A., et al. 2012, A&A, 544, L7

Rodríguez, L. F., Loinard, L., D'Alessio, P., Wilner, D. J., & Ho, P. T. P. 2005, ApJL, 621, L133

Sadavoy, S. I., Myers, P. C., Stephens, I. W., et al. 2018, ApJ, 869, 115

Sakai, N., Sakai, T., Hirota, T., et al. 2014, Natur, 507, 78
Seifried, D., Sánchez-Monge, Á, Walch, S., & Banerjee, R. 2016, MNRAS, 459, 1892

Stahler, S. W., & Palla, F. 2005, The Formation of Stars (New York: Wiley-VCH)

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

Takakuwa, S., Saigo, K., Matsumoto, T., et al. 2017, ApJ, 837, 86

Tobin, J. J., Looney, L. W., Li, Z.-Y., et al. 2016, ApJ, 818, 73

Ueda, T., Kataoka, A., & Tsukagoshi, T. 2020, ApJ, 893, 125

Ulrich, R. K. 1976, ApJ, 210, 377

van der Wiel, M. H. D., Jacobsen, S. K., Jørgensen, J. K., et al. 2019, A&A, 626, A93

van't Hoff, M. L. R., van Dishoeck, E. F., Jørgensen, J. K., & Calcutt, H. 2020, A&A, 633, A7

Vousden, W. D., Farr, W. M., & Mandel, I. 2016, MNRAS, 455, 1919

Wootten, A. 1989, ApJ, 337, 858

Yeh, S. C. C., Hirano, N., Bourke, T. L., et al. 2008, ApJ, 675, 454

Yen, H.-W., Takakuwa, S., Ohashi, N., et al. 2014, ApJ, 793, 1