© The Authors 2024

1 Introduction

Stellar mass is perhaps the most fundamental stellar property. Nevertheless, measuring stellar masses for the youngest protostars remains difficult because these objects are deeply embedded in an infalling core and envelope and cannot be observed directly in optical or near-infrared wavelengths (e.g. Di Francesco et al. 2007; Tychoniec et al. 2021). Detections of Keplerian rotation in the protostellar disks around protostars are paramount for constraining the stellar masses of these systems. Historically, Keplerian rotation has been difficult to detect (1) because bright molecular lines such as CO tend to become optically thick within the envelope, obscuring the disk, and (2) due to a lack of high-spatial resolutions needed to disentangle the disk from the envelope (e.g. Ohashi et al. 1997; Schreyer et al. 2006; Jørgensen et al. 2007; Murillo et al. 2013; van’t Hoff et al. 2018).

The Atacama Large Millimeter/submillimeter Array (ALMA) can provide the high-resolution molecular line observations that are necessary to make detections of Keplerian rotation in young disks possible. Moreover, the significant improvements in sensitivity with ALMA have opened the possibility of observing fainter and optically thin isotopologues, which allows us to detect disks through the envelope (e.g. Murillo et al. 2013; Yen et al. 2014; Ohashi et al. 2014; Aso et al. 2015; Ginsburg et al. 2018; Reynolds et al. 2021). ALMA is even capable of exploring deviations from pure Keplerian rotation that could be due to planets (e.g. Pinte et al. 2018, 2020; Teague et al. 2019, 2021), though primarily in older protoplanetary disks where the envelope is less of a confounding factor. Nevertheless, despite the improvements afforded by ALMA, well-constrained stellar masses have only been measured for a limited number of low-mass protostars to date (e.g. Murillo et al. 2013; Ohashi et al. 2014, 2023; Aso et al. 2017, 2023; Flores et al. 2023; Thieme et al. 2023).

VLA 1623 is in the ρ Oph A star-forming region at 140 pc (Ortiz-León et al. 2018) and is one of the rare protostellar systems with detected Keplerian rotation. Specifically, VLA 1623 is the canonical Class 0 protostellar object (André et al. 1993) and a hierarchical multiple system, with sources VLA 1623A and VLA 1623B (separated by ~140 au) located towards the centre of the core and VLA 1623W roughly 1500 au west of the pair¹ (Bontemps & Andre 1997; Looney et al. 2000; Murillo & Lai 2013). VLA 1623A is itself a compact binary, with its Aa and Ab components separated by ~14 au (Harris et al. 2018) and surrounded by a large (R ≈ 180 au) circumbinary disk (Murillo et al. 2013). Molecular line observations of the circumbinary disk show evidence of Keplerian rotation and yield a combined stellar mass of between 0.2 and 0.5 M_⊙ (Murillo et al. 2013; Hsieh et al. 2020). Measuring the masses of the other stellar components, however, has been challenging due to confusion with the circumbinary disk and insufficient sensitivity to robustly detect and resolve the disks in gas tracers.

Constraining the stellar masses is an important step towards examining the origins of the VLA 1623 system. Observations suggest that the B and W components may be more evolved (non-coeval) due to a lack of outflow and the shape of their spectral energy distributions (SEDs), which imply a later evolutionary stage (Murillo et al. 2018a). Harris et al. (2018) examined the proper motion of the stars and propose that VLA 1623W may have been ejected. An ejection scenario could explain the lack of envelope detected towards VLA 1623W (Michel et al. 2022) as well as the large velocity streamers in the system (Hsieh et al. 2020; Mercimek et al. 2023). Conversely, previous molecular line observations have been used to estimate the other stellar masses (Ohashi et al. 2022; Mercimek et al. 2023), and those estimations, even if not concrete measurements, indicate that VLA 1623B and VLA 1623W may be relatively massive compared with VLA 1623A. If confirmed, these masses would pose difficulties for the ejection scenario (Reipurth et al. 2010).

Here we present new molecular line observations of C¹⁷O (3–2) of the VLA 1623 system from ALMA. We used these data and radiative transfer models of rotating disks to constrain the stellar masses for the stellar components in VLA 1623. C¹⁷O is a rarer species that is less optically thick within the envelope than the line tracers that have previously been observed, thereby making it easier for us to disentangle the Keplerian rotating disks of the protostars from each other. The paper is organized as follows: Sect. 2 details the observations, Section 3 shows that the C¹⁷O traces the disk and the other lines trace outflows or infall, Sect. 4 describes the radiative transfer models, Sect. 5 discusses the stellar masses and the possible origins of the VLA 1623 system, and finally Sect. 6 gives our conclusions.

2 Data

ALMA observed VLA 1623 in Band 7 on 20 and 22 April 2019 for project 2018.1.01089.S in the C-4 configuration with baselines of 15 m – 740 m. The phase calibrator was J1650-2943, and the flux/bandpass calibrator was J1924-2914 (20 April) and J1517-2422 (22 April). The data phase centre was at 16:26:26.35 -24:24:30.55, between VLA 1623A and VLA 1623B, meaning that VLA 1623W is outside of the primary beam full width at half maximum (FWHM)². The time on source was 52.74 min in total. The observations were set up with spectral windows on C¹⁷O (3−2) at 337.061 GHz, H¹³CO⁺ (4−3) at 346.998 GHz, SO₂ (8_2,6−7_4,4) at 334.673 GHz, and H¹³CN (4−3) and SO₂ (13_2,12−12_1,11) at 345.34 GHz. There was also one continuum window with a central frequency of 348.2227 GHz and bandwidth of 1.85 GHz.

We applied self-calibration to the continuum window and the three line windows using line-free channels. We used three rounds of phase-only self-calibration, starting with a shallow clean and a long solution interval and progressing to deeper cleans with shorter solution intervals. Due to the high (> 1000) source signal-to-noise ratio (S/N), we also used two rounds of phase and amplitude self-calibration using long integration times. We used solution intervals of 3 min, 60.6 s, and 30.3 s for the phase-only self-calibration and 3 min and 90.9 s for the amplitude and phase self-calibration. The clean depths for the five rounds were 0.92, 0.45, 0.20, 0.17, and 0.13 mJy beam⁻¹. The final continuum map including all line-free channels have a sensitivity of 95 µJy beam⁻¹ at 341.8541 GHz for a beam of 0.3″ × 0.24″ with Briggs weighting and robust = −0.5. The maximum recoverable scale for the data is 3.8″ using the 5th percentile baseline.

The self-calibration gain solutions were applied to each of the line spectral windows. We applied continuum subtraction using a fit order of 1. The final line cubes are made using tclean interactively with robust = 0.5. Table 1 gives the rest frequency, channel width, sensitivity per channel, and beam size for each detected line. For C¹⁷O (3−2), we report the rest frequency of the main hyperfine component only. The H¹³CN (4−3) and SO₂ (13_2,12−12_1,11) lines may be blended. Hereafter, we refer to these data as (blended) H¹³CN to distinguish these data from the other detected lines. Table 1 also lists the velocity range used to make all moment maps (unless a different range is specified).

3 Results

3.1 Overview

Figure 1 shows the dust continuum data for VLA 1623 with the main components labelled. We recover the circumbinary disk around the protobinary, Aa and Ab, but we do not resolve the protobinary itself (labelled as A). Hereafter, we refer to the protobinary as VLA 1623A. The circumstellar disks for B and W are also labelled.

Figure 2 shows the Moment 0 data for each line with contours of dust continuum. For the sake of clarity, the Moment 0 maps have not been primary-beam-corrected. C¹⁷O and H¹³CO⁺ have extended emission, where both molecules trace the circumbinary disk and inner envelope around A and B. At the locations of the circumstellar disks for A and B, we see a relative deficit (absorption) in C¹⁷O and H¹³CO⁺ that is consistent with the circumstellar disks shadowing the brighter circumbinary disk material or self-absorption (e.g. Murillo et al. 2013; Hara et al. 2021). The circumstellar material is bright in SO₂ and (blended) H¹³CN, and both lines also show compact emission within the circumbinary disk towards the south-west.

Finally, although faint, we also find hints (~2.5Σ) of the vertical streamer passing through VLA 1623W in C¹⁷O (3−2) that matches what Mercimek et al. (2023) found in C¹⁸O (2−1) observations (labelled S in Fig. 2). VLA 1623W is only detected in C¹⁷O (see Sect. 3.3), with a system velocity of about 1.8 km s⁻¹. The system velocities of VLA 1623A and VLA 1623B are roughly 4 km s⁻¹ and 2.3 km s⁻¹, respectively.

3.2 Outflows and infall

The C¹⁷O data show a large X-shaped pattern centred on the VLA 1623A/B system (see Fig. 2). Weaker emission along this pattern is also seen in SO2 and (blended) H¹³CN, indicating that this emission may be associated with the outflow and not the inner envelope. VLA 1623A has overlapping outflows with position angles of ~125º (e.g. Murillo et al. 2018b; Hara et al. 2021; Ohashi et al. 2022) in good agreement with the orientation of the extended C¹⁷O data. We find no C¹⁷O gas within the outflow itself, however, suggesting that C¹⁷O is tracing the outflow cavity walls only, likely due to higher densities in the walls over the outflow itself and higher temperatures leading to more sublimation of CO isotopologues. We note that Hsieh et al. (2020) instead identify similar structures in C¹⁸O (2−1) as gas streams in the envelope. To distinguish between these cases, Fig. 3 compares the C¹⁷O integrated intensity data with contours of c-C₃H₂ at 217.8 GHz from Murillo et al. (2018b). Hydrocarbons like c-C₃H₂ are excellent tracers of outflow cavities (e.g. Sakai et al. 2014; Tychoniec et al. 2021; Ohashi et al. 2022). Broadly, we find good agreement between c-C₃H₂ and the X-shaped extended emission in C¹⁷O, which is why we attribute this material to the outflow wall rather than a gas streamer.

Emission from SO₂ has also been associated with jets and outflows of protostellar sources (e.g. Wakelam et al. 2005; Feng et al. 2020). Figure 4 shows the high velocity red- and blueshifted gas from SO₂ (8_2,6−7_4,4) for VLA 1623B (a similar trend is seen for the (blended) H¹³CN data). Both line data show a gradient aligned with the disk minor axis, which matches expectations for high velocity material associated with the jet or outflow and not a rotating disk. Due to confusion with the circumbinary disk, we cannot identify a clear gradient for VLA 1623A.

Alternatively, the SO₂ emission could be tracing infalling gas rather than outflowing gas. Sakai et al. (2014) find that sulfur-bearing species can trace the centrifugal barrier, a transition layer in disks with weak shocks due to infalling gas hitting the disk.

In these cases the sulfur-bearing gas is confined near the disk (e.g. Sakai et al. 2017). Codella et al. (2024) also find a large streamer extending south from VLA 1623B in C¹⁸O (2−1) and several transitions of SO, which they interpret as an accretion flow. While we do not see this streamer in our SO₂ gas, it is detected faintly in the (blended) H¹³CN data towards the southern part of the primary beam. Therefore, the SO₂ and (blended) H¹³CN data may be tracing a mix of both infall and outflow.

Fig. 1 Dust continuum map at 341.854 GHz (0.877 mm) zoomed in on the VLA 1623 sources. Protostars are labelled A (for the unresolved Aa and Ab binary), B, and W.

3.3 Circumstellar disks

Circumstellar disks around protostars are often well detected using rare isotopologues of CO like C¹⁷O (e.g. Tychoniec et al. 2021). While C¹⁷O appears to trace the outflow cavity walls on larger scales (see Fig. 3), we find that the high velocity gas towards the circumstellar disks have velocity gradients consistent with Keplerian rotation indicating that the compact emission indeed arises from the circumstellar disk.

Figure 5 shows C¹⁷O (3−2) Moment 1 maps of all three disks. To ensure a good signal-to-noise, we spectrally smoothed the data to 0.4 km s⁻¹ channels and we calculated the Moment 1 velocities using immoments in CASA with strict channel selections to avoid lower velocity gas that could be confused with emission outside of the disks. The velocity ranges used were 0.8−3.2 km s⁻¹ and 5.6−7.2 km s⁻¹ for VLA 1623A, (−5.2)−(−2.4) km s⁻¹ and 6.8−9.6 km s⁻¹ for VLA 1623B, and (−2.0)−5.2 km s⁻¹ for VLA 1623W. Moreover, we masked each map using thresholds in Moment 0 and dust continuum to isolate the Moment 1 data for each disk.

All three disks show blue-red gradients indicative of rotation. For VLA 1623A, the gradient is for the circumbinary disk as we lack the spatial resolution to resolve any gradients in the circumstellar disks of Aa and Ab. For VLA 1623B, the C¹⁷O emission is complicated by gas emanating from the circumbinary disk (see Fig. 2). Nevertheless, Fig. 5 shows that we can use high velocity (e.g. >4 km s⁻¹ from the systemic velocity) C¹⁷O emission to help isolate its circumstellar material. The velocity gradient seen in VLA 1623B aligns well with the major axis of the disk. Ohashi et al. (2022) find a similar gradient in VLA 1623B with high-velocity CS (5−4) data. While they found that the CS data traced both the disk and outflows, the high velocity CS gas formed a gradient that was perpendicular to the outflow direction. Thus, we conclude that the high velocity C¹⁷O (3−2) emission is similarly tracing rotation in the disk.

Figure 5 also shows a clear velocity gradient for VLA 1623W going from the northern to the southern part of the disk along the major axis as expected with disk rotation. This gradient agrees well with the C¹⁸O (2−1) gradient presented by Mercimek et al. (2023). We lack sensitivity in the C¹⁷O data to detect a velocity gradient in the streamers like Mercimek et al., but the broad centroid velocities of the two streamers in C¹⁷O (3−2) match those seen in C¹⁸O (2−1).

In the next section, we model the C¹⁷O (3−2) emission using disk models with Keplerian rotation. Broadly, we found that the emission from the highest velocity channels (greatest deviation from the systemic velocity) follow a r^−0.5 profile, as expected for Keplerian rotation in disks. Nevertheless, other dynamical processes such as infall from the envelope or streamers can produce gradients. While such features can be tested using rotation curves (e.g. following Ohashi et al. 2014; Aso et al. 2015; Maret et al. 2020), the disks are sufficiently resolved or confused with circumbinary material such that the centroid position becomes unreliable. For simplicity with our modelling, we assumed that the C¹⁷O (3−2) emission solely traces the rotationally supported disks and that any potential contributions from the surrounding envelope or gas streamers are negligible. To ensure this assumption is appropriate, we employed velocity cuts similar to Fig. 5 to exclude the low-velocity emission that could be dominated by these components (see Sect. 4 for details).

Fig. 2 Moment 0 maps (in Jy beam−1 km−1 s−1) of C17O (3−2), H13CO+ (4−3), SO2 (8−7), and (blended) H13CN (4−3). The velocity ranges are given in Table 1. Black contours show the Stokes I continuum at 7 and 49 mJy beam−1. The beam for each map is in the lower-left corner. The upper panels use log scaling to highlight the extended emission, whereas the bottom panels have linear scaling. The line maps are not primary-beam-corrected, for the sake of visualization. The first panel shows labels for VLA 1623W (W) and the velocity streamers (S) from Mercimek et al. (2023) for clarity.

Fig. 3 C17O (3−2) data from Fig. 2 with contours of c-C3H2 (217.8 GHz) integrated intensity from Murillo et al. (2018b) to show the outflow. The c-C3H2 integrated intensity data are evaluated over a velocity range of 2.65−5.36 km s−1, and the contours go from 0.02 to 0.05 Jy beam−1 km s−1 in steps of 0.005 Jy beam−1 km s−1. The map resolutions are given in the bottom-left corner.

Fig. 4 Red- and blueshifted SO2 (8−7) emission overlaid on the dust continuum data. The integrated intensities were generated using the velocity ranges indicated in the legend. Contours correspond to 10, 15, 20, 25, 30, 35, and 40 σ, where σ = 2.5 mJy beam−1 km s−1 for both the blue- and redshifted data. The SO2 map resolution is in the lower-left corner.

Fig. 5 C17O Moment 1 maps for all three disks using restricted velocity channels (see the main text for details). Black contours show the Stokes I continuum at 7 and 49 mJy beam−1. The beam for each map is in the lower-left corner. The Moment 1 maps are masked using thresholds in Moment 0 and dust continuum to avoid noisy pixels and cleanly show the gradients in the disks.

4 Modelling

We modelled the C¹⁷O (3−2) data for VLA 1623A, VLA 1623B and VLA 1623W using the pdspy code, which fits Keplerian rotating disk models to spectral line datasets directly in the uυ plane (Sheehan et al. 2019). pdspy uses radiative transfer codes and Bayesian sampling to identify a best-fit disk model and the full posterior distributions for the parameters given a set of observations. Full details of how the code works can be found in Sheehan & Eisner (2018), but we describe the technique briefly here, including any differences from what has previously been documented. In short, pdspy generates a 2D axisymmetric model of a protostellar disk. The surface density is described by a power-law function, $$\text{Σ}={\text{Σ}}_0{\left(\frac{R}{R_{\text{disk }}}\right)}^{-\gamma },$$(1)

where Σ₀ is the surface density normalization, R is the radial distance from the star in cylindrical coordinates, R_disk is the radius where the disk is truncated, and γ is the surface density power-law exponent. For VLA 1623B and VLA 1623W, we fixed R_in, the inner radius where the disk is truncated, to 0.1 au. As VLA 1623A is a close separation binary with evidence of a cavity (Harris et al. 2018), we left R_in as a free parameter. Rather than fitting for Σ₀, we instead integrated the surface density to calculate the total (gas and dust) mass, M_disk, and used that as a free parameter in our fits. We assumed the temperature follows a power law profile in radius, T(R) = T₀ (R/1 au)^−q, and that the disk is vertically isothermal. The scale height of the disk as a function of radius could then be derived from the temperature at a given radius and the equations of hydrostatic equilibrium, and similarly the volume density follows from the surface density and scale height. We modelled the C¹⁷O (3−2) emission (no dust) assuming a constant C¹⁷O abundance with a factor of 1/2240 relative to CO and a fixed CO abundance of 10⁻⁴ relative to molecular hydrogen (e.g. following Czekala et al. 2015). We also adopted a constant turbulent velocity, a_turb, throughout the disk but leave its value as a free parameter. pdspy then uses the radiative transfer modelling code RADMC-3D and ray tracing to generate synthetic channel maps of the model disk with a given inclination (i) and position angle (PA), and then Fourier transforms the synthetic channel maps into the uυ plane using the galario code (Tazzari et al. 2018). This last step ensures that the synthetic model data match the spatial scales covered by the observations. We note that the position angle is defined based on the direction of the angular momentum vector of the disk, east of north, following Czekala et al. (2015).

The model parameter values as well as the source system velocity, υ_sys, and location in the field (x₀, y₀) are optimized to the observations. Unlike in Sheehan et al. (2019), we use the dynesty code (Speagle 2020) instead of a Monte Carlo Markov chain (MCMC) to find the best-fit values. dynesty uses nested sampling to integrate the likelihood function and calculate the Bayesian evidence and the posteriors for each parameter. We opted to use dynesty here because we find nested sampling to be more robust with multi-modal posteriors and it also has well-defined stopping criteria.

Since the C¹⁷O (3−2) observations include emission from VLA 1623A, B, and W, we applied additional constraints on the model and data to ensure that we fit each source separately. First, we provided pdspy with relatively strict priors on the location of the source and allowed that position to vary minimally (by <±0.3″ in either direction) so that the code is forced to fit the emission from the target of interest. Second, we restricted the channels and uυ ranges to select emission where the target source is dominant or at least clearly separable from all other sources as well as from the envelope and outflows. For the channel restrictions, we excluded velocities between 3.25 and 5.5 km s⁻¹ from the fit to A and between −2.5 and 6.75 km s⁻¹ from the fit to B. For the uυ restriction, we excluded baselines with uυ < 100 kλ (≳2″) for B and uυ < 200 kλ (≳ 1″) for W. The uυ cuts are less severe for B because there is less extended emission in the vicinity of B once the channel cuts previously described are taken into account. For VLA 1623A, we are fitting the extended circumbinary disk and as such, the velocity channel cuts are sufficient to exclude the envelope-scale without additional restrictions in the uυ range. To ensure that these choices do not affect the fit quality, we carefully check the results of the models.

Table 2 lists the best-fit model parameters for VLA 1623A, B and W. The best-fit values are determined from the maximum likelihood models from the posterior and the uncertainties are calculated as the range around the best-fit values that include 68% of the posterior samples. We also added a 10% uncertainty in quadrature to most parameters to represent additional uncertainties on those values imparted by systematic errors in, for example, flux calibration and source distance. We excluded this additional uncertainty from x₀, y₀, i, PA, and υ_sys as they are primarily geometric and are therefore less likely to be impacted by the aforementioned systematics. There may be additional systematic effects such as the choice of model that lead to larger errors on the measured values than those presented here (e.g. Premnath et al. 2020). We further note that several of the best-fit parameters are not well constrained or have unusually low or high values. For example, some of the best-fit models have very low power-law exponents for q and γ, implying flat profiles for temperature or surface density. We caution against over-interpreting the value of these disk parameters given that we had to impose strict velocity and uυ limits. These cuts will limit the spatial extent of the disk, which will affect our ability to constrain the radial profiles (e.g. power-law indices). We focus instead on the results for the stellar mass, which we can reasonably constrain with the high-velocity gas, but we present the best-fit values for the other disk parameters for completeness.

Figs. 6–8 compare the channel maps from the observations with the synthetic channel maps of the best-fit model to demonstrate the fit quality for each source. We find that the models fit the data well, with little or no contours above the 5σ level appearing in the residual channel maps associated with the target disk. Emission that exceeds 5σ in the residuals generally arises from structures not included in the model but are present in the field of view. For example, Fig. 7 shows significant redshifted residual emission above the 5σ level in the channel maps of VLA 1623B, but this residual emission comes from the circumbinary disk, which is not included in the model of VLA 1623B. This broad agreement between the observations and the model indicates that our simple model of Keplerian rotation was indeed reasonable and that other dynamical processes such as infall are either negligible or limited to larger scales than that of the disk.

Figure 9 compares the observed position-velocity data along the disk major axis (background colour) with the synthetic position-velocity data (contours) from the best-fit models of VLA 1623A, B, and W. Though position-velocity diagrams represent a reduction of the dimensionality of the data, they can be a useful tool for by-eye comparison with Keplerian rotation. Figure 9 shows simple Keplerian rotation curves for fixed stellar masses to demonstrate that the masses derived from our fitting are consistent with what might be obtained from such a simple comparison.

Fig. 6 Comparison of the channel maps for VLA 1623A with the best-fit Keplerian disk model. The top two rows show (blue- and redshifted) channel maps from the observations, the middle two rows show the best-fit model, and the bottom two rows show the residuals. The central velocities of the channel maps are indicated in the top two rows. For the model and residual channel maps, we used galario to match baselines and then cleaned the synthetic data with similar imaging parameters as the data. The contours start at 5σ and continue in increments of 20σ, with dashed contours showing negative levels. We excluded channels between 3−5.75 km s−1 as these channels were not fit for VLA 1623A.

5 Discussion

5.1 Stellar masses

Previously, Keplerian rotation modelling had only been done for VLA 1623A. Murillo et al. (2013) determine a (combined) stellar mass of M_★ ~ 0.2 M_⊙ using models of Keplerian rotation with infall and early (Cycle 0) ALMA C¹⁸O (2−1) data. With more sensitive data, Hsieh et al. (2020) measure a combined stellar mass of M_⊙ ~ 0.3−0.5 M_⊙ for VLA 1623A. In both cases, the mass values for VLA 1623A are a total mass for the tight binary, VLA 1623Aa and VLA 1623Ab (Harris et al. 2018). The remaining stellar components, VLA 1623W and VLA 1623B, have estimated masses from the velocity gradient assuming Keplerian rotation, rather than fitting the Keplerian motion itself. For VLA 1623B, Ohashi et al. (2022) measure a mass of M_★ ≲ 1.7 M_⊙ from CS (5−4) data with υ_rot = 7.8 km s⁻¹ at R = 27 au, and for VLA 1623W, Mercimek et al. (2023) estimate a mass of M_★ ≈ 0.45 M_⊙ from C¹⁸O (2−1) data with υ_rot ≈ 3 km s⁻¹ at R ≈ 50 au.

Broadly, our measured stellar masses from the radiative transfer models agree well with the previously inferred dynamical masses. We find best-fit stellar masses of M_⋆ = 0.64 ± 0.06 M_⊙ for VLA 1623W and $M_{\star }={1.9}_{-0.2}^{+0.3}{M}_{\odot }$ for VLA 1623B, and a combined stellar mass of M_⋆ = 0.27 ± 0.03 for VLA 1623A (see Table 2). Since these masses are measured from C¹⁷O (3–2), we have a rarer isotopologue of carbon monoxide than C¹⁸O and less confusion with the infalling envelope and streamers (Murillo et al. 2013; Mercimek et al. 2023). Moreover, our modelling process fits the full 3D (position, position, velocity) dataset with physically motivated models rather than estimating the mass from the peak line emission in a few channels at a radius offset from the disk centre. In addition, pdspy fits visibility data. The visibility plane encodes information at several times smaller spatial scales than is recovered by deconvolution (Jennings et al. 2020), and as such we are not limited by the beam resolution and can better constrain the models on sub-beam scales.

These measurements put constraints on the stellar masses for all the components (with VLA 1623A combined) and indicate a hierarchy. VLA 1623B has the highest mass, with more than twice the mass of the other two stellar components combined. VLA 1623W is second in mass with VLA 1623A in third. These more accurate stellar masses yield valuable insights into the physical conditions of VLA 1623, and its formation and evolution. We discuss individual scenarios in Sect. 5.2.

Fig. 7 Same as Fig. 6, but for VLA 1623B. Contours start at 5Σ and continue in increments of 3Σ. Channels between −2 and 6.5 km s−1 and baselines <100 kλ are excluded.

Fig. 8 Same as Fig. 6, but for VLA 1623W. Contours start at 5Σ and continue in increments of 3Σ. Baselines <200 kλ are excluded from the imaging, and the data are binned into 0.5 km s−1 channels to improve image signal-to-noise.

Fig. 9 Position-velocity (PV) diagrams for VLA 1623A (left), VLA 1623B (centre), and VLA 1623W (right). The background images show the observed data, and contours show the synthetic data from the best-fit model for each source. The PV diagrams were extracted with an aperture width of ~0.35″, i.e. about one beam-width, perpendicular to the extracted spatial direction for VLA 1623B and VLA 1623W as the emission for both sources is marginally resolved. As VLA 1623A is much better resolved, we instead used a width of ~0.7″ to extract its PV diagram. Velocities between 3.25 km s−1 and 5.5 km s−1 were masked for VLA 1623A to avoid channels with significant envelope emission, and velocities between −2.5 km s−1 and 6.75 km s−1 were masked for VLA 1623B to avoid confusion from emission associated with VLA 1623A. We note that these cuts lead to asymmetries in the number of channels available on either side of the line centre. We also show example Keplerian rotation curves for stars with a range of stellar masses as a simple visual comparison.

5.2 Origin of the VLA 1623 system

The detection of Keplerian rotation towards both VLA 1623W and VLA 1623B indicates that these sources are genuine pro-tostars (see also Ohashi et al. 2022; Mercimek et al. 2023). Nevertheless, the VLA 1623 system is a complex environment and its origin is not straightforward. Here, we consider different scenarios to explain the formation and evolution of VLA 1623.

5.2.1 Dynamical interactions

One theory to explain the VLA 1623 system is that VLA 1623W was ejected. An ejection scenario has been proposed because the proper motions of VLA 1623W indicate that it points back to the VLA 1623A/B protostars (Harris et al. 2018). A dynamical interaction could also explain the counter-rotating disks between VLA 1623A and VLA 1623B (Takaishi et al. 2021; Ohashi et al. 2022), the lack of envelope emission around VLA 1623W, and the large gas streamers (Murillo et al. 2013; Mercimek et al. 2023).

Nevertheless, there are challenges to explaining VLA 1623W as an ejected object. First, multi-body interactions generally result in the ejection of the lowest-mass object (Reipurth et al. 2010). From our stellar mass measurements, VLA 1623A has the lowest mass by more than a factor of two even if the circumbi-nary disk is included. Second, a dynamical interaction leading to ejection is expected to produce a spiral structure in the circumbi-nary disk (Takaishi et al. 2021), which is not observed, although a spiral structure caused by dynamical interactions could dissipate quickly (Cuello et al. 2023). Third, VLA 1623W does not appear to get close enough to VLA 1623A/B for a dynamical ejection to be likely. In Appendix A, we extrapolate the proper motion of VLA 1623W and VLA 1623B back in time to find a closest plane-of-sky encounter of about 650 au. Reipurth et al. (2010) simulated dynamical interactions of triple systems of protostars in dense cores at separations between 50 au and 400 au and found that ejections were weaker and less common at 400 au than at 50 au (see also the results from stellar flybys from Cuello et al. 2023). Our analysis shows that the closest separation only becomes <100 au at the 6th percentile, whereas the closest separation is ≳400 au at the 25th percentile (see Appendix A). For such wide projected separations, we would expect that the stars would be too distant to have a substantial gravitational effect on their relative orbits (Sadavoy & Stahler 2017). As a result, a dynamical ejection scenario appears unlikely to explain the VLA 1623 system.

5.2.2 Chance alignment

Since VLA 1623W has a different system velocity relative to VLA 1623A and B, it may be an unrelated young stellar object (YSO) that happens to have a chance alignment with the VLA 1623A/B system. Murillo & Lai (2013) identify VLA 1623 as having non-coeval YSO components based on their SEDs (see also Murillo et al. 2018a). VLA 1623W also lacks envelope emission and an outflow (Murillo et al. 2013; Santangelo et al. 2015; Hara et al. 2021), which could indicate that it is a more evolved YSO that formed separately.

Following Tobin et al. (2022), we estimated the probability that VLA 1623W is a companion source using Bayes theorem, $$P(c\mid d)=\frac{P(d\mid c)P(c)}{P(d)},$$(2)

where P(d|c) is the probability we would detect a true companion source, P(c) is the probability of having a companion at a given separation, and P(d) is the probability of a detection of any YSO. We adopted P(d|c) = 0.75, similar to Tobin et al. (2022). Tobin et al. (2022) measured their companion detection probability based on a dust mass sensitivity of ~ 1 M_⊕ for dusty disks. Since we have higher sensitivity observations of ~0.01 M_⊕ for source detection (e.g. Sadavoy et al. 2019), our adopted value of P(d|c) should be considered a lower limit. We also used P(c) = 0.2 from Tobin et al. (2022), which is based on companion statistics from Perseus and Orion for a separations ≲ 1000 au. Finally, we calculated the probability of having a detection, P(d), with $$P(d)=0.75P(c)+\left[1{\text{e}}^{-0.75\text{Σ}\pi r^2}\right][1-0.75P(c)],$$(3)

where Σ is the stellar surface density and r is the separation being considered. Equation (3) gives the probability of detecting any source at the observed sensitivity. The first term is the likelihood of detecting a true companion, assuming 75% are detectable at our sensitivity, and the second term is the likelihood of detecting an unrelated source in a specified area for a given YSO stellar density (for further details, see Appendix A of Tobin et al. 2022). We used r = 10″ for the area. To estimate the YSO stellar density, we calculated the 11th nearest neighbour at the position of VLA 1623 from the Gaia-corrected YSO catalogue of Grasser et al. (2021). The 11th nearest neighbour gives Σ₁₁ = 1470 YSO pc⁻². This stellar density is slightly below the peak YSO stellar density for L1688 of 2000 YSOs pc⁻² from Gutermuth et al. (2009), but VLA 1623 is located off the cluster centre and so the 11th nearest neighbour should be more representative.

With our assumed values, we calculate a probability that VLA 1623W is a companion protostar of VLA 1623 of P(c|d) = 0.55, which suggests that VLA 1623W has nearly equal chances to being a member or an unrelated YSO. Nevertheless, this probability of 55% is likely a lower limit. First, VLA 1623W appears protostellar in nature (see Sect. 5.2.3), and the stellar surface densities from both Gutermuth et al. (2009) and Grasser et al. (2021) are dominated by Class II sources. Since protostars (embedded YSOs) are less common than Class II YSOs (e.g. Gutermuth et al. 2009; Dunham et al. 2015), the probability that VLA 1623W is a companion object should be higher than our current estimate. Second, the streamers towards VLA 1623W indicate that it is interacting with the dense core itself and therefore is likely within the dense core and not completely unrelated to the VLA 1623 system. These factors increase the likelihood that VLA 1623W is a true companion source, so we consider VLA 1623W to be a member YSO.

5.2.3 In situ formation

Here, we define objects that formed out of the same natal core environment as forming in situ. In general, there are two main mechanisms that have been proposed to produce multiple stellar systems: disk fragmentation and turbulent fragmentation (Offner et al. 2023). Briefly, disk fragmentation occurs when a circum-stellar disk becomes unstable to gravity and fragments to form additional stars, primarily at small (≲100 au) separations (e.g. Bonnell & Bate 1994; Kratter et al. 2010). Turbulent fragmentation occurs when density perturbations within the natal dense core become massive enough for self-gravity, allowing them to collapse and form independent objects at a wide range of separations (e.g. Offner et al. 2010, 2016; Kuffmeier et al. 2019). Although the separation distributions can broadly indicate the formation mechanism, simulations show that wide binary systems can shrink to smaller orbits on timescales of ≲0.1 Myr (Offner et al. 2012), which is less than the protostellar lifetime. As a result, differences in angular momentum vectors may be more illuminating, as disk fragmentation predicts aligned vectors and turbulent fragmentation implies random vectors (see Offner et al. 2023, and references therein).

VLA 1623A and VLA 1623B have a projected separation of roughly 200 au, whereas VLA 1623W is over 1000 au away from the pair. The wide separation for VLA 1623W would indicate that it formed via turbulent fragmentation, but from separation alone, it is difficult to conclude the origins of A or B, especially given the size of the circumbinary disk. Nevertheless, VLA 1623A and VLA 1623B have misaligned velocity gradients indicative of counter-rotation (see Fig. 5 and Ohashi et al. 2022) and the disk inclinations and positions angles do not agree as well (see Table 2). Both of these factors are consistent with turbulent fragmentation. Moreover, since VLA 1623B is more massive than VLA 1623A and its circumbinary disk, it seems unlikely that it was formed via disk fragmentation and then was later perturbed through dynamical interactions to have a misaligned rotation axis.

Turbulent fragmentation does not require that all stellar components formed at the same time. Murillo & Lai (2013) and Murillo et al. (2018a) compared the SEDs and circumstellar material for the VLA 1623 components and concluded that VLA 1623W appeared more evolved than VLA 1623A and VLA 1623B. Indeed, there is a strong, collimated outflow coming from the VLA 1623A/B system, whereas VLA 1623W has no detected outflows or envelope (Murillo et al. 2013; Santangelo et al. 2015; Hara et al. 2021; Michel et al. 2022), in agreement with this protostar being older. Thus, VLA 1623W may have formed first through turbulent fragmentation and then VLA 1623A/B formed later.

Although the protostars in VLA 1623 may not be entirely coeval, this does not mean that VLA 1623W must be an evolved YSO (e.g. Flat or pre-main sequence). First, VLA 1623W is detected in C¹⁸O (Mercimek et al. 2023) and C¹⁷O. Protostel-lar disks are generally warmer, which favour the detection of rare isotopologues in the gas phase, and these disks also tend to be more massive and optically thick, which will shield the molecules from being selectively photo-dissociated as seen in older protoplanetary disks (e.g. Miotello et al. 2017; van’t Hoff et al. 2018; Zhang et al. 2020; Artur de la Villarmois et al. 2019). Some protoplanetary disks have been detected in C¹⁷O (e.g. Zhang et al. 2021), but these observations are generally towards the inner radii, whereas we detect C¹⁷O towards the entire disk of VLA 1623W (e.g. beyond the dust disk; see Fig. 5). Younger protostars tend to be brighter with warmer disks, which means that rare CO isotopologues will be detected in the gas phase out to larger radii (e.g. van’t Hoff et al. 2018). Second, VLA 1623W appears to be flared at 0.87 mm and 1.3 mm (Michel et al. 2022), indicating that larger dust grains have not yet had time to settle to the midplane. Protoplanetary disks tend to be geometrically thin at those wavelengths (e.g. Villenave et al. 2020), whereas flaring has been detected towards some protostellar disks (e.g. Sheehan et al. 2022; Ohashi et al. 2023). Finally, VLA 1623W is an optically thick, relatively massive disk (Harris et al. 2018; Sadavoy et al. 2019; Michel et al. 2022), and disk mass tends to decrease with evolutionary stage for a given star-forming region (e.g. Tobin et al. 2020). Table 2 gives an equivalent disk dust mass of 6 M_⊕, assuming a gas-to-dust mass ratio of 100. This dust mass places VLA 1623W in the upper half of disk masses from the Ophiuchus disk survey, ODISEA, and is consistent with Class I (or Flat) disks (Williams et al. 2019). Although there are more massive disks around evolved YSOs in other clouds, the ODISEA survey shows that disk masses in Ophiuchus tend to be lower than other nearby systems, making the mass of VLA 1623W significant relative to the other Ophiuchus disks.

Even though it has a youthful disk, VLA 1623W lacks an envelope component typical of protostellar sources. Since VLA 1623W is located along the collimated outflow axis (in projection), its envelope could have been stripped (e.g. de Gouveia Dal Pino 2005; Ladd et al. 2011). In this case, we would observe a less embedded SED and the protostar would have less accretion to drive an outflow. Currently, there is no evidence of shocked gas towards VLA 1623W (e.g. it is not detected in SO₂), although we could lack the sensitivity. While we cannot clearly classify VLA 1623W given the uncertainties from its SED and challenges with its high disk inclination, we believe that VLA 1623W is more likely to be a protostellar source and less likely to be a more evolved (e.g. Class II) object.

For the origin of the VLA 1623 protostellar system, we propose that the main stellar components, VLA 1623A, VLA 1623B, and VLA 1623W, formed via turbulent fragmentation, whereas the tight VLA 1623A binary system (Aa and Ab) at the centre of a large circumbinary disk may have formed via disk fragmentation given their small separation, although we cannot rule out that all the stellar components formed via core fragmentation and migrated. A full origin picture must account for the wide range of stellar masses for each of the stars, explain why the circumbinary disk formed around the lowest-mass stellar companion, and determine the stability of the stellar system and the circumbinary disk.

5.3 VLA 1623 stability

In this section, we consider the gravitationally stability of the protostellar disks and the VLA 1623 system. For disk stability, we used the disk-to-star mass ratio to assess the gravitational stability of the circumstellar and circumbinary disks. Systems with low disk-to-star mass ratios should be more stable to fragmentation, whereas higher ratios would be unstable (e.g. Vorobyov 2010). Mercer & Stamatellos (2020) find that typical ratios of ≳ 30% favour fragmentation in disks based on a sample of low-mass stars. This cutoff matches what is seen in the L1448 IRS3 system, where the fragmenting circumbinary disk of L1448 IRS3B has a large disk-to-star mass ratio of 25% and the gravitationally stable disk around L1448 IRS3A has a smaller disk-to-star mass ratio of 3% (Reynolds et al. 2021). Gravitational stability in L1448 IRS3 was independently assessed using a Toomre Q analysis and the presence or absence of spiral structures and fragments in the disks themselves (Tobin et al. 2016; Reynolds et al. 2021).

From our best-fit models in Table 2, we calculate ratios of M_disk/M_⋆ ~ 20% for the circumbinary disk around VLA 1623A and M_disk/M_⋆ ~ 0.3% for both VLA 1623B and VLA 1623W. These results imply that the stellar disks are likely stable under their own self-gravity, even with the uncertainties on the disk masses from the models. The circumbinary disk, however, has an intermediate ratio, which could indicate that it is close to being unstable and may undergo fragmentation in the future. Nevertheless, there is no evidence of ongoing spiral structure or fragmentation in the circumbinary disk (e.g. Sadavoy et al. 2018; Harris et al. 2018), both of which are considered signposts of gravitational instabilities (e.g. Kratter et al. 2010) and are detected in L1448 IRS3B (Tobin et al. 2016). We cannot rule out that the circumbinary disk fragmented in the past (e.g. to form the tight VLA 1623Aa and VLA 1623Ab binary) and has since reached a more stable point.

In the case of the circumbinary disk, spiral structure could also arise due to dynamical interactions with its close companion, VLA 1623B, rather than gravitational instabilities. It is interesting that we see no evidence of perturbations in the circumbinary disk given the close (in projection) separation between VLA 1623A and VLA 1623B. Specifically, VLA 1623B is the most massive YSO in the system (e.g. its stellar mass exceeds the estimated combined mass of VLA 1623A and the circumbinary disk by a factor of six) and the projected separation of VLA 1623A and VLA 1623B is approximately equal to the minimum value given their disk sizes. The circumbinary disk, however, is well fit by an isolated disk model with no significant residuals indicating deviations from Keplerian rotation. Moreover, both disks have low gas temperatures (Murillo et al. 2018b), which suggest no dynamical interactions or perturbations are taking place. Both factors imply that VLA 1623B may have a large line-of-sight separation and is not physically close enough to VLA 1623A to perturb or heat the gas in the circumbi-nary disk. A large line-of-sight separation could also indicate that VLA 1623B is an unrelated YSO, but given their close separation in projection, they are more likely to be true companion. Using the same probability definitions as Sect. 5.2.2, we find P(c|d) = 0.95 for VLA 1623A and VLA 1623B to be companion YSOs, assuming P(c) = 0.14 for separations < 500 au (Tobin et al. 2022) and a distance of r = 1.5″. With such a high probability, VLA 1623B should be considered a companion protostar.

Finally, we considered whether VLA 1623W is itself grav-itationally bound to the VLA 1623A/B YSOs, assuming it is a true companion source (see Sect. 5.2.2). For simplicity, we ignored any contributions from the surrounding dense core, and calculated the potential energy between VLA 1623W and VLA 1623A/B using the star and disk masses in Table 2 and a projected separation of r = 10″, and we calculated a kinetic energy based on the relative 3D motions of the A/B and W from their system velocities (assuming 3 km s⁻¹ for A/B and 1.75 km s^-1 for W) and proper motions (see also Appendix A). The resulting energies are Ω = −GM_ABM_W/r ~ −0.2 × 10³⁷ J for the gravitational energy and $K=\frac{1}{2}{M}_{\text{W}}{\sigma }_{3D}^2\approx 1\times {10}^{37} \text{J}$ for the kinetic energy. These back-of-the-envelope calculations imply that the turbulent kinetic energy is roughly a factor of five larger than the gravitational energy from the YSOs alone. Moreover, the calculated gravitational energy is an upper limit, since we only have a projected separation between VLA 1623A/B and VLA 1623W, and the true 3D distance could be much larger. Thus, the kinetic energy appears several times larger than the gravitational energy, indicating that VLA 1623W is unlikely to be bound to the VLA 1623A/B system and could disperse, although the source may still be bound to the dense core.

6 Conclusions

We present new ALMA molecular line observations of the VLA 1623 system. We primarily focused on C¹⁷O (3–2) observations, which trace the disks of the protostars and show velocity gradients consistent with Keplerian rotation. We used the radiative transfer modelling code pdspy to model the C¹⁷O (3–2) emission for VLA 1623A, VLA 1623B, and VLA 1623W, obtaining constraints for their stellar masses. Our main conclusions are:

• 1.

  We measure stellar masses of 0.27 M_⊙, 1.9 M_⊙, and 0.64 M_⊙ for VLA 1623A (Aa and Ab combined), VLA 1623B, and VLA 1623W, respectively. These masses are in good agreement with previous estimates that used different tracers and techniques;

• 2.

  Based on the new mass measurements and an analysis of the proper motion of the stars, we disfavour a scenario where VLA 1623W was ejected from the central core;

• 3.

  Following Tobin et al. (2022), we estimate a probability of VLA 1623W being a companion source of nearly 55%, which means we cannot rule out that it is an unrelated YSO along the line of sight. Nevertheless, based on the apparent youth of the disk and the presence of gas streamers connecting VLA 1623W to the VLA 1623A/B protostars, we favour it being a genuine companion source;

• 4.

  We propose a scenario where VLA 1623A, VLA 1623B, and VLA 1623W formed initially from turbulent fragmentation, although VLA 1623A may have undergone disk fragmentation to produce the tight binary, Aa and Ab. The protostars may not be entirely coeval with each other, but differences in age should not exceed the protostellar lifetime;

• 5.

  The disks around VLA 1623B and VLA 1623W appear gravitationally stable based on very low disk-to-star mass fractions. The circumbinary disk has an intermediate fraction of 20%, which could indicate instability and future fragmentation. Nevertheless, we see no evidence of spiral structure in the circumbinary disk, either from gravitational instability or dynamical interactions.

We also find that VLA 1623W appears to be unbound to the other protostars, suggesting that it may disperse. As such, these observations represent a rare snapshot in time of a multiple-protostellar system prior to its diffusion or strong dynamical interactions. While the presence of streamers towards VLA 1623W (Mercimek et al. 2023) indicates that some interactions between the stars and envelope may be occurring, we do not see disturbed gas motions in the envelope or circumbinary disk. Moreover, given the mass of VLA 1623B and its projected separation, we would expect it to influence the circumbinary disk, although we detect no evidence of perturbations with the current data. Future analyses that model the circumbinary disk may yet discover deviations from Keplerian rotation due to gravitational perturbations from the more massive VLA 1623B, which would show the onset of dynamical interactions in action.

Acknowledgements

We thank the anonymous referee for valuable comments that improved the publication. The authors thank the NAASC for support with the ALMA observations and data processing. The authors would also like to thank Ewine van Dishoeck for valuable discussions. SIS acknowledges the support for this work provided by the Natural Science and Engineering Research Council of Canada (NSERC), RGPIN-2020-03981, RGPIN-2020-03982, and RGPIN-2020-03983. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.01089.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. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.

Appendix A Proper motion and closest approach

This appendix outlines the calculations for the proper motions for VLA 1623B and VLA 1623W. We excluded VLA 1623A from these calculations since the combination of it being a tight binary and having a circumbinary disk can confuse its centroid position. We took positions of VLA 1623B and VLA 1623W from the literature that present high-resolution data (≲ 1″) and where the epoch of observations could be identified. Table A.1 lists the corresponding archival data, where Column 1 gives the reference, Column 2 gives the epoch of observation, Column 3 gives the observing frequency, Columns 4 and 5 give the source right ascension and its adopted error, Columns 6 and 7 give the source declination and its adopted error, Column 8 gives the geometric mean size of the synthesized beam, and Column 8 gives the phase calibrator for the observations. We follow Sadavoy et al. (2018) in using the map resolution scaled by the source peak S/N to estimate the position uncertainties. For the present study and those of Sadavoy et al. (2019) and Harris et al. (2018), both sources had peak S/N > 100, so we adopted a minimum error of 5 mas.

Figure A.1 shows the relative position of VLA 1623W compared to the epoch in Sadavoy et al. (2019). We show the VLA 1623W results only as an example as the corresponding plot for VLA 1623B is very similar. The dashed line shows our adopted proper motion, which we measured using the data from Sadavoy et al. (2019), Harris et al. (2018), Murillo et al. (2013), and Chen et al. (2013) only. We excluded the data point from Maury et al. (2012), because they did not provide sufficient precision for a proper motion analysis (Sadavoy et al. 2018). We further excluded our new measurement presented here due to a noticeable offset of roughly 15 mas in RA and 30 mas in Dec compared to the general trend (dashed line), which we found for both VLA 1623W and VLA 1623B. As such, this 15 – 30 mas offset appears to be systematic for the entire map.

A systematic offset in position could be explained by the present study using different calibrators. Table A.1 lists the phase calibrator for each observation. Most studies used J1625-2527, but this work used J1650-2943. From the ALMA Technical Handbook, we can expect a 1 σ pointing accuracy of ~15 mas, which is of comparable order to the systematic offsets. To address possible astrometry uncertainties from the calibrators themselves, we added in quadrature a 15 mas position error to the position errors from Table A.1 before fitting for the proper motion using linear least squares.

We find a proper motion for VLA 1623W of $${\mu }_{\alpha }\cos \delta =-12.4\pm 5.1\text{\hspace{0.17em}mas\hspace{0.17em}}{\text{yr}}^{-1}$$(A.1) $${\mu }_{\delta }=-25.6\pm 6.3\text{\hspace{0.17em}mas\hspace{0.17em}}{\text{yr}}^{-1}\text{, }$$(A.2)

which are slightly different from Harris et al. (2018) but within errors. Harris et al. (2018) measured a proper motion of VLA 1623W using their data and those of Murillo et al. (2013) only. There is a possibility that the position of VLA 1623W may also have differences due to wavelength since VLA 1623W is viewed nearly edge-on and shows evidence of flaring (Michel et al. 2022). There is a slight difference in RA position between Harris et al. (2018) at 870 µm and the 1.3 mm measurements from Murillo et al. (2013) and Sadavoy et al. (2019), whereas their Dec positions are in good agreement. VLA 1623W appears nearly vertical in Dec, which would make any temperature variations with disk scale height appear only in RA. Nevertheless, Murillo et al. (2013) have a different phase calibrator, which could also induce an offset. Therefore, we cannot conclude that the shift in RA from Harris et al. (2018) is due to temperature stratification.

Fig. A.1 Proper motion fits for VLA 1623W using the data from Table A.1. Offsets are shown relative to Sadavoy et al. (2019). Best-fit proper motions are given by the dashed lines. Only the 1.3 mm (230 GHz) data are used for the fit. Error bars include a 15 mas pointing uncertainty added in quadrature with the position uncertainty from Table A.1.

We also re-measured the proper motion of VLA 1623B from Sadavoy et al. (2018) using our revised errors that include adopting 15 mas pointing uncertainties. The revised proper motion for VLA 1623B is $${\mu }_{\alpha }\cos \delta =-7.8\pm 2.9\text{\hspace{0.17em}mas\hspace{0.17em}}{\text{yr}}^{-1}$$(A.3) $${\mu }_{\delta }=-29.0\pm 3.4\text{\hspace{0.17em}mas\hspace{0.17em}}{\text{yr}}^{-1}\text{, }$$(A.4)

which is nearly identical in RA and has a slight difference in Dec within errors compared to Sadavoy et al. (2018).

We combined the above proper motions for VLA 1623W and VLA 1623B to determine when the two objects had their closest approach in the plane of the sky, assuming constant velocities. Figure A.2 shows the extrapolation of both nominal proper motions with black circles showing look-back times in 500 yr intervals. The closest approach is shown with open circles at a time of ~ 1400 years, when the two sources were separated by ~ 4.6 arcsec (~ 650 au) in the plane of the sky.

Nevertheless, the proper motions for both sources have large error bars (> 10%), which make uncertainties on the closest approach difficult to constrain analytically. We used a Monte Carlo error analysis for the proper motions for VLA 1623B and VLA 1623W to estimate the error in closest separation. We assumed that the proper motion uncertainties given above correspond to the FWHM of a Gaussian distribution and then draw 10000 random values to add to the nominal proper motions. For each new set of proper motions, we re-evaluated the closest approach. Figure A.3 shows histograms from this Monte Carlo analysis for the distribution of closest approach and the time since closest approach. In both cases, we excluded roughly 200 data points for which the closest approach was given by the current epoch (e.g. in the case where the two sources are moving towards each other rather than away) as this case slightly skewed the statistics. From the error analysis, we find a median minimum separation of 4.9 arcsec (first and third quartiles are 2.7 arcsec and 6.6 arcsec) and a median time of 1150 years (first and third quartiles are 830 years and 1590 years). These results are in good agreement with the nominal case using the best-fit proper motions (without the errors). Moreover, this analysis indicates that VLA 1623W was unlikely to have a very close encounter as a plane of sky separation of < 100 au is only at the 6th percentile.

Fig. A.2 Positions of VLA 1623B and VLA 1623W extrapolated in time based on their proper motion. Arrows show the relative proper motions of both YSOs as measured here (for W) and in Sadavoy et al. (2018) (for B). Solid circles show look back times for 5000 years in steps of 500 years. The closest approach (in the plane of the sky) is represented by open circles.

Fig. A.3 Results for the minimum separation (top) and time of closest approach (bottom) for the Monte Carlo error analysis between the proper motions of VLA 1623B and VLA 1623W. The solid lines show the median values and the dashed lines show the first quartile (25th per-centile) and third quartile (75th percentile) for each distribution. Note that the minimum separation has an upper limit of 9.3 arcsec, corresponding to the current separation between the two sources.

References

All Tables

All Figures

	Fig. 1 Dust continuum map at 341.854 GHz (0.877 mm) zoomed in on the VLA 1623 sources. Protostars are labelled A (for the unresolved Aa and Ab binary), B, and W.
In the text

Fig. 2 Moment 0 maps (in Jy beam−1 km−1 s−1) of C17O (3−2), H13CO+ (4−3), SO2 (8−7), and (blended) H13CN (4−3). The velocity ranges are given in Table 1. Black contours show the Stokes I continuum at 7 and 49 mJy beam−1. The beam for each map is in the lower-left corner. The upper panels use log scaling to highlight the extended emission, whereas the bottom panels have linear scaling. The line maps are not primary-beam-corrected, for the sake of visualization. The first panel shows labels for VLA 1623W (W) and the velocity streamers (S) from Mercimek et al. (2023) for clarity.

In the text

	Fig. 3 C17O (3−2) data from Fig. 2 with contours of c-C3H2 (217.8 GHz) integrated intensity from Murillo et al. (2018b) to show the outflow. The c-C3H2 integrated intensity data are evaluated over a velocity range of 2.65−5.36 km s−1, and the contours go from 0.02 to 0.05 Jy beam−1 km s−1 in steps of 0.005 Jy beam−1 km s−1. The map resolutions are given in the bottom-left corner.
In the text

	Fig. 4 Red- and blueshifted SO2 (8−7) emission overlaid on the dust continuum data. The integrated intensities were generated using the velocity ranges indicated in the legend. Contours correspond to 10, 15, 20, 25, 30, 35, and 40 σ, where σ = 2.5 mJy beam−1 km s−1 for both the blue- and redshifted data. The SO2 map resolution is in the lower-left corner.
In the text

	Fig. 5 C17O Moment 1 maps for all three disks using restricted velocity channels (see the main text for details). Black contours show the Stokes I continuum at 7 and 49 mJy beam−1. The beam for each map is in the lower-left corner. The Moment 1 maps are masked using thresholds in Moment 0 and dust continuum to avoid noisy pixels and cleanly show the gradients in the disks.
In the text

Fig. 6 Comparison of the channel maps for VLA 1623A with the best-fit Keplerian disk model. The top two rows show (blue- and redshifted) channel maps from the observations, the middle two rows show the best-fit model, and the bottom two rows show the residuals. The central velocities of the channel maps are indicated in the top two rows. For the model and residual channel maps, we used galario to match baselines and then cleaned the synthetic data with similar imaging parameters as the data. The contours start at 5σ and continue in increments of 20σ, with dashed contours showing negative levels. We excluded channels between 3−5.75 km s−1 as these channels were not fit for VLA 1623A.

In the text

	Fig. 7 Same as Fig. 6, but for VLA 1623B. Contours start at 5Σ and continue in increments of 3Σ. Channels between −2 and 6.5 km s−1 and baselines <100 kλ are excluded.
In the text

	Fig. 8 Same as Fig. 6, but for VLA 1623W. Contours start at 5Σ and continue in increments of 3Σ. Baselines <200 kλ are excluded from the imaging, and the data are binned into 0.5 km s−1 channels to improve image signal-to-noise.
In the text

Fig. 9 Position-velocity (PV) diagrams for VLA 1623A (left), VLA 1623B (centre), and VLA 1623W (right). The background images show the observed data, and contours show the synthetic data from the best-fit model for each source. The PV diagrams were extracted with an aperture width of ~0.35″, i.e. about one beam-width, perpendicular to the extracted spatial direction for VLA 1623B and VLA 1623W as the emission for both sources is marginally resolved. As VLA 1623A is much better resolved, we instead used a width of ~0.7″ to extract its PV diagram. Velocities between 3.25 km s−1 and 5.5 km s−1 were masked for VLA 1623A to avoid channels with significant envelope emission, and velocities between −2.5 km s−1 and 6.75 km s−1 were masked for VLA 1623B to avoid confusion from emission associated with VLA 1623A. We note that these cuts lead to asymmetries in the number of channels available on either side of the line centre. We also show example Keplerian rotation curves for stars with a range of stellar masses as a simple visual comparison.

In the text

	Fig. A.1 Proper motion fits for VLA 1623W using the data from Table A.1. Offsets are shown relative to Sadavoy et al. (2019). Best-fit proper motions are given by the dashed lines. Only the 1.3 mm (230 GHz) data are used for the fit. Error bars include a 15 mas pointing uncertainty added in quadrature with the position uncertainty from Table A.1.
In the text

	Fig. A.2 Positions of VLA 1623B and VLA 1623W extrapolated in time based on their proper motion. Arrows show the relative proper motions of both YSOs as measured here (for W) and in Sadavoy et al. (2018) (for B). Solid circles show look back times for 5000 years in steps of 500 years. The closest approach (in the plane of the sky) is represented by open circles.
In the text

Fig. A.3 Results for the minimum separation (top) and time of closest approach (bottom) for the Monte Carlo error analysis between the proper motions of VLA 1623B and VLA 1623W. The solid lines show the median values and the dashed lines show the first quartile (25th per-centile) and third quartile (75th percentile) for each distribution. Note that the minimum separation has an upper limit of 9.3 arcsec, corresponding to the current separation between the two sources.

In the text