The Disk Orientations of Perseus Protostellar Multiples at $\sim$8 au Resolution

For consistency we utilized the following definitions for this work:

1. 1.

   source: A single protostar that can have a compact and/or extended disk.

2. 2.

   system: A collection of sources within a defined separation that may be interacting.

3. 3.

   aligned: A pair of sources whose dot product of the orientation vectors would correspond to a value of $<$30^∘ (Lee et al., 2019a).

4. 4.

   misaligned: A pair of sources whose dot product of the normalized orientation vectors would correspond to a value of $>$30^∘.

The definitions of alignment are for ease of qualitative referencing and is not relied on for the analysis detailed in Section 4. In the analysis, we utilized the alignment angles and the corresponding observational errors. The demarcation of alignment between two orientation vectors of 30^∘ are chosen to remain consistent with studies of simulated data (Lee et al., 2019a).

The observations were taken as part of project 2019.1.01425.S with 48 antennas included between 2019 September 11-18 at Band 6 (1.3 mm) toward 12 protostellar systems in the Perseus molecular cloud (d$\sim$300 pc). The observations were carried out in the most extended configuration C-9/10 (baselines 150 m$\sim$14.9 km) and have an effective angular resolution of 42 mas$\times$23 mas to 66 mas$\times$29 mas, with a continuum sensitivity of 13-76 $\mu$Jy beam^-1, when reconstructed with the Briggs robust 1 weighted imaging. The correlator was configured with 3 spectral windows setup for 1.875 GHz bandwidth and 3840 channels and four spectral windows used 117.19 MHz bandwidth and centered on the ¹³CO ($J=2\rightarrow 1$), C¹⁸O ($J=2\rightarrow 1$), SO ($J=6(5)\rightarrow 5(4)$), and SiO ($J=5\rightarrow 4$) transitions. However, the spectral lines were not well-detected given that the integration times were chosen for continuum sensitivity.

The first set of observations for 5 of the sources took place across a 1.5-hour block, with each time-on-source $\approx$9 minutes (scans$\sim$9.04 s). The remaining 7 sources were observed in 3 execution blocks (EBs), across 3 days, with the average time-on-source $\sim$14 minutes for the 3 EBs combined. While the absolute flux density scale is expected to be accurate to $\sim$10%, for the purpose of the results presented, all flux uncertainties only consider statistical uncertainties. A summary of the scheduled observations and execution blocks is given in Tables 1 and 2.

The raw visibility data were calibrated by the North American ARC staff using Common Astronomy Software Applications (CASA) version 6.2.1 automated pipeline. The high sensitivity of the observations (pre-self-calibration$\sim$40 $\mu$Jy) and signal-to-noise ratios (S/N) of at least 50, enabled self-calibration to be attempted for all sources except for Per-emb-2. We performed 4 to 5 rounds of phase-only self calibration, with the first round of intervals starting at the full length of the EB, then round two of intervals starting at the full length of the on-source scans, then progressing to 18.14 s, 9.07 s, and ending at the single integration timestep. Per-emb-18, Per-emb-35, NGC1333 IRAS2B, and L1448 IRS3B were unable to be phase self-calibrated down to the shortest timestep due to the S/N degrading, but were phase-only self-calibrated down to 9.07 s. The final average sensitivity resulting from the phase-only self calibrations was $\sim$30 $\mu$Jy and an average increase of S/N by a factor of 1.5. A summary of the observations is detailed in Table 1.

The data were imaged using CASA version 6.5.0-15 with the task tclean, and the images using Briggs weighting with the robust parameter 1.0 are shown in Figures 1 and 2 with image sizes of 9000$\times$9000 and 4 mas per pixel. All images are shown with square-root stretch colormaps and a common RMS value of 20 $\mu$Jy beam^-1. To restore the images, we used the Multi-(Taylor-)term Multi-Frequency Synthesis (MTMFS) with scale sizes of 0, 5, and 20 pixels and 2 Taylor terms. The scale sizes were chosen to recover dominant features in the disk and correspond to physical sizes of point source, typical size of the beam minor axis, and 2$\times$ the typical beam major axis. We utilized the “auto-multithresh” masking technique to non-interactively mask and clean the data in a reproducible manner by defining the sidelobe threshold sub-parameter to be 2.0. The final images were checked by simultaneously cleaning the data with a conservative user-defined mask.

ALMA Band 3 data were taken toward three targets, a subset of the Band 6 observations: Per-emb-2, Per-emb-18, and Per-emb-5, within ALMA project 2016.1.00337.S. Data were taken in two ALMA configurations, C40-6 (2017 Aug 1) and C40-9 (2017 September 27), summarized in Table 3. For both datasets, the phase calibrator was J0336+3218, the flux calibrator was J0238+1636 and the bandpass calibrator was J0237+2848. The baselines sampled in the combined dataset ranged from 16 m to 14,500 m.

The correlator was configured with 3 spectral windows setup for 1.875 GHz bandwidth and 128 channels and two spectral windows used 58.59 MHz bandwidth and centered on the ¹³CO and C¹⁸O ($J=1\rightarrow 0$) transitions. However, the spectral lines were not well-detected given that the integration times were chosen for continuum sensitivity. The central frequency of the observations is $\sim$102 GHz.

The data were pipeline processed by the observatory using the pipeline included in CASA 4.7.2 (r39732). The data were also self-calibrated and imaged using this same version of CASA. The C40-6 data went through 3 rounds of phase-only self-calibration with solution intervals of one scan, 24.15s, and 6.05s (a single integration). The C40-9 data also went through three rounds of phase-only self-calibration but the second interval used a 12.10s solution interval, and the final solution interval of 6.05s corresponded to 3 integrations.

The final self-calibrated data were imaged together using the CASA task clean with image sizes of 2048$\times$2048 pixels and 5 mas pixels. We made use of MTMFS imaging given the wide fractional bandwidth, restoring the images with super-uniform weighting scheme to closely match the beamsize of the 1.3 mm observations, and interactively cleaned using hand-drawn masks. The images were cleaned down to $\sim$1.5$\times$ the noise in each image (Figure 3).

For the most compact companion sources (Per-emb-18, L1448 IRS3B, L1448 IRS3C, and SVS13A), the compact disk emission is only slightly larger than the size of the beam, so the deconvolved PA and i derived from image-plane analysis will be less well-constrained. Moreover, we desired an alternate fitting method to measure the source parameters independent of the images produced with the CLEAN algorithm. For completeness, we conducted the image-plane analysis in Appendix C and find the results between the uv-visibility and image plane analysis are consistent within 3 $\sigma$.

In order to utilize the full spatial constraints afforded by the observations, we constructed a number of Gaussians equal to the number of sources in the uv-visibilities. Similar techniques were applied to protoplanetary disks (Jennings et al., 2022a, b) to recover substructure with $>2\times$ longer effective baselines than the reconstruction from CLEAN by fitting the visibilities directly. We constrained the sources using Bayesian inference, fitting multi-component 2-D Gaussians to the visibilities of each individual source using (emcee, dynesty, pdspy; Foreman-Mackey et al., 2013; Speagle, 2020; Sheehan, 2022). We restricted the uv-visibilities fitting to scales smaller than 0.^′′5 (by restricting the uv-distance $>$400 k$\lambda$) to ensure we fit the compact disk and not the extended emission of the envelope or circum-multiple material that was not resolved out in the observations. We also limited the fit phase center to be within 0.^′′5 of previously published results, and in the cases of new detections, we utilized the centering from the CASA task imfit to form the prior (a summary of the imfit results and the comparison with the uv-visibility results is described in Appendix C).

A summary of the fitted parameters is provided in Table 4 , and a summary of the projected 3D orientation vectors solved from fitting the uv-visibilities is shown in Figure 4. The errors reported are derived as the 1 $\sigma$ uncertainty from the median of the sampled posterior.

While all sources could be described by a Gaussian, the source L1448 IRS3A was best described with a ring (see Appendix A). The detailed analysis of L1448 IRS3A falls outside of the scope of this paper and we leave analysis of its disk structure for a future paper. While the uv-visibilities enable a more complete picture of the system and accurate representation of the sources, without bias from the inherent beam geometries and subsequent clean procedure.

This source was previously reported to be a close multiple $<$50 au with the VLA at 9 mm (Tobin et al., 2016a) and was further observed in Tobin et al. (2015a, 2018) as a smooth continuum surface brightness distribution at 1.3 mm (Figure 5). We resolved the compact binary (a$\sim$80 mas$\approx 24$ au) in both the 1.3 mm and 3 mm observations. Our observations resolved much of the large scale emission, but further revealed a possible additional 3 compact sources with separations of $\sim$0.^′′431 ($\sim$129 au), $\sim$1.^′′432 ($\sim$430 au), and $\sim$0.^′′5 ($\sim$150 au) for the southern (SNR$\approx$25), northern-most (SNR$\approx$22), and southern-most (SNR$\approx$15) sources respectively, relative to the Per-emb-2-A source of the compact binary. All of these compact sources (except the southern-most companion) are found within regions of relatively enhanced surface brightness within the Tobin et al. (2018) observations, presumed to be a massive extended disk. The sources are present at least the 5 $\sigma$ level in the 1.3 mm and the 3 mm observations. It should be noted that only the compact sources -A,-B and a marginal detection of the brighter southern source -C appear in 9 mm VLA observations (Tobin et al., 2016a) whereas more diffuse sources, the northern-most -D and and southern-most -E sources do not.

Per-emb-5 was also previously reported to be a close multiple $<$50 au with the VLA at 9 mm (Tobin et al., 2016a). We instead found continuum emission that appears consistent with that of a disk surrounding what appeared as two peaks in the VLA data (Figure 6). The disk appeared to have a central cavity (centered between the two VLA peaks), and maybe a single spiral arm to the west. There is an asymmetry across the minor axis of the disk and a flux enhancement in the south-east portion of the disk.

L1448 IRS3B is certainly a exceptional source (Figure 7). The system is home to at least 4 compact continuum sources within 8^′′, 3 of which are within 1^′′ of each other. The brightest feature, the tertiary companion commonly known as L1448 IRS3B-C is likely optically thick. The tertiary companion is embedded within one of the large spiral arms that stems from the inner disk to the outer disk. Zooming into the center of the system, two bright continuum sources are obvious, one just inside of an inner spiral arm/disk structure and one just outside. We apparently resolved the “clump” as reported in Reynolds et al. (2021) as the north-east portion of the inner disk and now report an additional faint compact source near the geometric center of the inner disk. The bright point L1448 IRS3B-A is now resolved as the south-west portion of the ring and a bright source just inside of the “ring” of the inner disk. L1448 IRS3B-B is just outside of the “ring”.

We reproduced Figure 16 from Reynolds et al. (2021) in Figure 7, denoting the locations of the kinematic centers for the L1448 IRS3B system using various disk tracing molecular line observations and techniques. We visually depict the new geometric center of the inner ring, which coincides with the center of the “deficit” previously reported and overlaps, within observation uncertainties, with the center-point of the kinematic centers. It is possible this newly resolved source is another deeply embedded protostellar source and the disk could harbor as many as 4 protostellar sources. For the purpose of the analysis conducted later in the paper, we do not consider this small point source at the center as another independent source; instead, we only consider the two confirmed compact continuum sources as protostar sources. The new source, while confidently detected, is too faint to have its geometric parameters well-constrained from Gaussian fitting. We refer to this new source as L1448 IRS3B-D, centered at 03h25m36.326s 30:45:14.93. The designation of sources in L1448 IRS3B may need reassignment in the future once the nature of the source and the inner disk are better characterized with additional observations.

Within the same field-of-view as the L1448 IRS3B observations, we resolved the disk around L1448 IRS3A (Figure 8). Previous ALMA 870 $\mu$m observations hinted of substructure in the surrounding circumstellar disk of L1448 IRS3A, but were unable to determine the nature of the features (Reynolds et al., 2021). Additionally, Reynolds et al. (2021) determined the structure was likely not spiral arms due to the disk being relatively stable against collapse. With the higher resolutions afforded by these observations, we resolve the circumstellar material to be a ring surrounding a compact source which remains unresolved in these observations. The compact continuum source appears slightly off the geometric center of the ring, but this could be explain by projection effects of an inclined disk on the observations.

In order to facilitate the analysis of the relative disk orientations in each multiple system, we must first determine which systems are associated with each other. While companions $<$100 au are mostly trivial to assign, systems with $>$3 sources and separations that range out to $\sim$10,000 au require a more automated approach. We made use of a method similar to the one implemented by Tobin et al. (2022) which will automatically create the companion associations by calculating the line-of-sight (L.O.S.) distances given an input catalog of positions. The systems constructed by this algorithm are verified to be consistent with prior studies (Tobin et al., 2016a, 2018) and the algorithm is discussed in detail in Appendix B.

To investigate the proto-multiple disk orientation alignments, we measured the projected angular difference of the inclination and position angle of the disks. With the full 2-D Gaussian fit parameters provided from the CASA imfit and uv-visibility fitting, we constructed the orientation vectors for each of the sources. We assumed that the disks are axisymmetric and geometrically thin and solved for the inclination of the disk from the major and minor axis lengths, deconvolved from the synthesized beam, via arccos($\frac{\sigma_{minor}}{\sigma{major}}$). The position angle of the disk is directly output from the 2-D Gaussian fits as the angle of the major axis, deconvolved from the synthesized beam, oriented to the standard east-of-north.

From a single continuum observation alone, we are insensitive to the orbital motions of the disk, thus we were unable to differentiate between aligned and anti-aligned disks. So for our analysis, we restricted i to range from 0-90^∘ (such that face-on is 0^∘) and the PA to range from 0-180^∘. The resulting dot product of the angular momentum vectors is normalized to fall between 0-1, and we constructed the summary plot in Figure 4 by overlaying resulting orientation vectors on the continuum image from Figure 1. In most sources the image-plane and the uv-visibility derived orientation vectors are nearly aligned, however for the sources where the compact disks are either unresolved or marginally resolved, the methods yield slightly different vectors. In these cases, for the purpose of interpretation, we favored the uv-visibility derived results as these results will take full advantage of the spatial-dynamic range in the data. We present the image-plane fits in Appendix C and the results of the statistical tests in Appendix D.

Our goal is to statistically analyze the companion orientations within the sample and provide a robust way to characterize the ensemble sample of orientations. To do so, we need to generate a sample of model protostellar configurations to determine the most probable formation pathways for the observed sample. The protostellar configurations we considered for the analysis are disks that are preferentially aligned with each other, representing the expectation of companions formed via disk fragmentation, and randomly aligned disks as would be expected to result from turbulent fragmentation.

To generate these configurations, the parameters we considered to generate a single system are: stellar multiplicity, the position angle, and inclination of the compact circumstellar disks (ignoring circum-multiple disks and under the assumption the stellar angular momentum axis and disk angular momentum axis are aligned). From these parameters alone, we constructed our model distributions. To model an empirically-driven distribution of proto-multiple systems, we sampled multiplicity for separations between 20-10,000 au, following the distribution of Perseus protostellar multiples of Tobin et al. (2022).

From these constructed systems, we generated a range of possible protostellar configurations with various fractions of aligned and randomly-aligned orientations. That is to say, a system with a multiplicity of 4 sources and a fractional orientation of 75% preferentially aligned sources and 25% randomly aligned sources, will have on average, 3 preferentially aligned sources and 1 randomly aligned source. We inserted no bias or preferential weighting corresponding to the separations of the sources with regards to alignment.

Preferentially aligned systems are constructed by drawing the position angle and inclination from a normal distribution described by a full-width-half-maximum value of 30^∘ (Lee et al., 2019a). The choice of 30^∘ is chosen to be consistent with similar studies carried out but the results presented are not significantly changed if moderate deviations ($\pm$5^∘) from 30^∘ is chosen. The resulting distributions of predominantly preferentially aligned systems are similar to population synthesis simulations conducted by Bate (2018). Randomly aligned systems are constructed by drawing the position angle and inclination from a uniform distribution of inverse cosine from 0 to 1.

For a rigorous statistical analysis, the observations need to be compared against a continuous distribution of models. Thus we generated 10,000 model systems (each system will have at least two sources) for each of the different fractions of systems with fractional orientations, ranging from fully preferentially aligned to completely random alignments. When extracting the orientation vectors of the constructed systems, we applied the same 2D projection bias the empirical continuum observations are subject to (i.e. we normalized the p.a. to $<$180^∘ and the inclination to $<$90^∘). This provides an observation-like set of models to compare directly with the observations. A visual representation of the median empirical cumulative distribution function (ECDF) for each of the constructed fractional alignment distributions is shown in Figure 9, with the median ECDF of the observations with uncertainties is overlaid in black.

For the observational and model data, we evaluated the dot product of orientation vectors for unique disk pairs in every system (counting each disk pair only once; Table 4). We utilized an algorithm similar to the companion finding approach, where each compact binary within a system is compared first. We then randomly selected a source from each compact binary to use for further comparison. This gives us N${}_{s}-1$ number of pairs per system where N_s is the number of sources within a given system (more detailed explanation in Appendix D). This resulting dot product is a derived from the results of the uv-visibility fitting (imfit analysis remains entirely consistent with the uv-visibility results and is detailed in Appendix C). To account for the uncertainties in the observations and to resample the particular sources chosen for the pairs, we recalculated each of the dot products 10,000 times sampling the fit errors assuming Gaussian uncertainties. This gave a suite of 10,000 realizations of the empirical distribution which are all consistent with the observations within the uncertainties.

We need to construct a distribution of alignments that could form the basis for the observational underlying distribution. Since this underlying distribution is not known, we constructed a grid of distributions that would cover the range of possible alignment distributions. Each particular constructed distribution is an aggregate sample of preferentially aligned and randomly aligned distributions, sampled via some fractional ratio (e.g. C_0.75UC${}_{0.25}=$3-to-1, aligned to randomly-aligned, etc).

We then calculated the probability of each of the 10,000 resampled observed distributions being drawn with each of the fractional ratios (i.e. the null hypothesis, see Appendix E) by utilizing the 2-sample Kolmogorov-Smirnov, Anderson-Darling, and Epps-Singleton probability tests (Scholz & Stephens, 1987; Hodges, 1958; Epps & Singleton, 1986; Goerg & Kaiser, 2009, respectively). We set a null-hypothesis rejection threshold of 0.3% (3 $\sigma$), such that if the probability test can reject the null hypothesis at the 3 $\sigma$ threshold, we discarded that particular empirical distribution. We then counted up all of the distributions that pass this threshold and summarize the full statistics in Table 5. A full description of the tests and methodologies is given in Appendix E.

To analyze the potential signatures of formation mechanism pathways, we divide the full set of observations into three subsamples and perform analysis on these subsamples. The subsamples chosen correspond to the maximum separation of companions, derived from the fit parameters: $<$100 au, $<$10,000 au, and $>$500 au. Surveys of protostellar multiplicity (Tobin et al., 2016a, 2020; Encalada et al., 2021; Tobin et al., 2022) found the average separation of companions to be $\sim$75 au, thus we chose 100 au as the compact subsample. Since turbulent fragmentation can form on 1000s au scales and then migrate down to 100s au scales, we select scales $>$500 au to select the subsample with a different underlying distribution. However, our selection of subsample noticeably foregoes the sources that fall within the range of separations $100<a<500$ au. This selects out $\sim$13 source pairs. While performing such a selection cut reduces the overall number of sources included in the subsamples, the authors do not detect any major difference that would change the findings. This particular cut of $100<a<500$ au was chosen to ensure sources selected by the two subsamples $<100$ au and $>$500 au would probe different underlying distributions, minimizing the overlap in these particular distributions, making the resulting statistical test more sensitive to differences in the underlying distributions.

When the sample is limited to $<$100 au separations (Left panel: Figure 9), the statistical tests imply the fractional ratio of orientations is comprised primarily of preferentially aligned sources with distributions at least 40% preferentially aligned.

When the sample is limited to only extended companions with separations $>$500 au separations (Right panel: Figure 9), 15 source pairs are analyzed and we rule out distributions of fully preferentially aligned companions down to 80% preferentially aligned. The results of the statistical results for each of the distributions appears equally likely for fractional alignments between 50%-70% preferentially aligned, thus we are not able to conclusively determine the exact alignment ratio.

Our analysis found the full observed sample (Middle panel: Figure 9), separations out to 10,000 au, is most consistent with a hybrid population of multiples with some contribution from multiple systems with preferentially aligned orientation vectors (up to 80% aligned sources) but is consistent with distributions down to 40% fractional alignments. We strongly ruled out fully aligned and fully randomly aligned distributions of disk orientation pairs. While it would be expected for distributions of highly separated sources to consist of a higher fractional ratio of randomly aligned systems, we found this distribution is consistent with roughly an equal fractional ratio.

The results are not significantly different whether we used the result from the image-plane or uv-visibility analysis, but the results from the uv-plane analysis tend to be favored as the uncertainties in the fits are more constrained.

There are two primary mechanisms for multiple star formation which operate at various scales, disk fragmentation and turbulent fragmentation during core collapse. Disk fragmentation operates on $\sim$100s au scales in massive disks ($\frac{M_{d}}{M_{*}}\gtrsim 0.1$) that are gravitationally unstable (i.e. Toomre Q $<$ 1), if the disk cools sufficiently fast (Gammie, 2001). The outcomes of gravitational instability can be observed directly (such as; Tobin et al., 2016b; Reynolds et al., 2021) but also has been extensively modeled (Kratter & Matzner, 2006; Boley, 2009; Kratter et al., 2010b; Vorobyov et al., 2013; Vorobyov & Elbakyan, 2018). However clear cases of ongoing disk fragmentation are somewhat elusive aside from L1448 IRS3B, perhaps due to the short timescales before the disk self-stabilizes by redistributing the angular-momentum and/or fragments. Turbulent fragmentation typically operates on 1000s of au scales in turbulent cloud cores. Moreover, simulations show that systems may migrate significantly from their nascent locations (Ostriker, 1999; Lee et al., 2019b) due to gravitational attraction and initial velocities with respect to the cloud core. Stars in these environments can become bounded and unbounded with empirical evidence that suggests multiples frequently form via turbulent fragmentation (Murillo et al., 2016). Therefore, it is likely that gravitational instability produces companions with compact separations ($<100$ au), both pathways can produce companions with moderate separations ($<$500 au), and a single mechanism may primarily populate the more extended configurations at scales of $>$500-10,000s of au.

While we cannot directly infer the evolution of any one particular system, we can use statistical approaches to modeling and determine the most probable formation pathway for an ensemble of multiple systems. Particularly, Bate (2018) conducted hydrodynamic simulations of proto-multiples undergoing gravitational collapse in a viscous medium, and found for a set of unbiased sources, the relative alignment angles are not correlated with hierarchy number, separation distributions, or age (see Figures 19 and 24 in Bate, 2018). Additionally Offner et al. (2016) showed that systems formed through turbulent fragmentation are randomly aligned and found that partial misalignment persists even after inward orbital migration. This was further supported by Lee et al. (2016) in observing companions with separations $a>$1,000 au who found a nearly complete random distribution of alignment angles. Analyzing our results of statistical tests (Table 5), we can statistically infer the fractions of preferentially aligned and randomly aligned orientation pairs that are in the given sample.

We found for our particular Perseus sample, for the distribution of companions with separations $<$10,000 au, distributions of a majority preferentially aligned orientation pairs are ruled out with at most 80% preferentially aligned. The $a<$10,000 au sample also disfavors distributions of less than 40% preferentially aligned sources. With the relatively low number of statistics (at most 21 pairs of sources), we were not able to determine the exact underlying fractional alignment ratio that best describes the sample of sources. This characterization of the subsamples and full sample is consistent with synthetic radiation hydrodynamical simulations of protostellar clusters (Bate, 2018).

It is likely the more compact companions of our sample ($a<$100 au, left panel of Figure 9 and Table 5), which do not rule out the null hypothesis in fractional ratio tests from 100% preferentially aligned distributions down to 40% preferentially aligned distributions, are formed primarily via gravitational instability, forming preferentially aligned pairs. This is similar to surveys of compact ($a\approx 50$ au) star-planet binaries (Dupuy et al., 2022) which were observed to have mutual orbital inclinations $<$30^∘. The large spread in the fractional ratio tests is likely due to the low number of sources (n$=7$) that have separations $a<$100 au). Bate (2018) shows an EDCF of compact companions that appears to be most consistent with 80% preferentially aligned sources.

The extended samples ($a>$500 au, right panel of Figure 9 and Table 5) appear to reject distributions with more than 80% preferentially aligned sources and distributions with less than 40% preferentially aligned sources.

The full sample ($a<$10,000 au, middle panel of Figure 9 and Table 5) are likely more constrained due to the higher number of sources (n$=$21) and appear consistent with distributions of 80% to 40% preferentially aligned sources.

The statistical tests alone are not able to tell a complete picture of an individual system, but it can provide a statistical way to characterize a sample of systems. To best determine the formation mechanism of any particular system, a combination of multi-wavelength observations of continuum and molecular lines are needed. We would expect sources that have formed via gravitational instability to have separations on the order of the continuum disk size ($a\approx$100s au) and the kinematics of the system to be organized. Whereas, formation of system via turbulent fragmentation happens on much larger scales ($a\approx$1,000-10,000s au). The sources formed at larger scales could further migrate to $\sim$100s of au, preserving no preferential alignment of the kinematics. We discuss the likelihood of the observed systems deriving from these two formation mechanisms based on our 1.3 mm and 3 mm ALMA continuum observations combined with prior observations.