The Disk Orientations of Perseus Protostellar Multiples at ∼ 8 au Resolution

We present a statistical characterization of circumstellar disk orientations toward 12 protostellar multiple systems in the Perseus molecular cloud using the Atacama Large Millimeter / submillimeter Array at Band 6 ( 1.3 mm ) with a resolution of ∼ 25 mas ( ∼ 8 au ) . This exquisite resolution enabled us to resolve the compact inner-disk structures surrounding the components of each multiple system and to determine the projected 3D orientation of the disks ( position angle and inclination ) to high precision. We performed a statistical analysis on the relative alignment of disk pairs to determine whether the disks are preferentially aligned or randomly distributed. We considered three subsamples of the observations selected by the companion separations a < 100 au, a > 500 au, and a < 10,000 au. We found for the compact ( < 100 au ) subsample, the distribution of orientation angles is best described by an underlying distribution of preferentially aligned sources ( within 30 ° ) but does not rule out distributions with 40% misaligned sources. The wide companion ( > 500 au ) subsample appears to be consistent with a distribution of 40% – 80% preferentially aligned sources. Similarly, the full sample of systems with companions ( a < 10,000 au ) is most consistent with a fractional ratio of at most 80% preferentially aligned sources and rules out purely randomly aligned distributions. Thus, our results imply the compact sources ( < 100 au ) and the wide companions ( > 500 au ) are statistically different.


INTRODUCTION
Recent studies in the past several decades have shown nearly half of all solar-type star systems are multiples (Raghavan et al. 2010;Duchêne & Kraus 2013;Moe & Di Stefano 2017;Offner et al. 2022a).It has been discovered that stellar multiplicity is even more common for young stars (Mathieu 1994;Chen et al. 2013;Tobin et al. 2020), and protostars in the midst of the stellar assembly process have the highest multiplicity fractions (Connelley et al. 2008a;Chen et al. 2013;Tobin et al. 2016a).During the earliest stages of star formation, the deeply-embedded protostellar phase, the largest reservoir of mass is available to form multiples (Tohline 2002).This is the stage of stellar evolution that must be examined to reveal the origins of stellar multiplicity.
Multiple star systems are thought to primarily form via two processes that operate on distinct scales: massive disks undergoing disk fragmentation on 10s-100s au scales (e.g., Kratter et al. 2010a) and turbulent core fragmentation on ∼1000s au scales (e.g., Offner et al. 2010).These processes can operate simultaneously, possibly giving rise to populations of close (<500 au) and wide (>1000 au) multiple systems (Tobin et al. 2016a).However, while the scales by which these processes form multiples are distinct, the systems formed via turbulent fragmentation may migrate to ∼100s au separations (or less) in ∼10s of kyr (Offner et al. 2010;Lee et al. 2019a) depending on their relative gravitational attraction with the core and the relative velocities of the sources at the times of formation.This makes it difficult to uniquely identify the dominant formation mechanism from separation measurements alone.Studies of protostellar multiplicity within the Orion molecular cloud by Tobin et al. (2022) found that current simulations of turbulent collapse alone did not account for all of the observed multiples found between 20-500 au, and thus an additional mechanism was needed to explain the observations.Meanwhile, Murillo et al. (2016) characterized the relative evolutionary states for wide and compact YSO multiples using SED modeling and found ∼33% of multiple systems were inconsistent with "co-eval" formation mechanisms.
Distinguishing if there is a primary mechanism for close multiple star formation is important for understanding the origin of stellar multiples, their evolution, and the potential impact they might have on their circumstellar disks and the planet formation potential.A close multiple system formed within a circumbinary disk may undergo relatively smooth evolution, while a close multiple system that forms as a result of migration from large separations can greatly disturb the system, leading to misalignment of outflows (Offner et al. 2016) and possibly disrupting accretion disks.
The VLA Nascent Disk And Multiplicity (VANDAM) Survey characterized the multiplicity of the entire protostar population in the Perseus molecular cloud (Tobin et al. 2016b), finding 17 multiple systems with separations less than 2. ′′ 0 (600 au) out of 90 sources observed.This sample of close multiples was followed-up with ALMA observations of 1.3 mm continuum and molecular lines that likely trace disk kinematics ( 13 CO, C 18 O, SO, and H 2 CO), enabling the most likely formation mechanism to be inferred for 12 systems (Tobin et al. 2018).Eight out of 12 systems were found to be consistent with disk fragmentation and four were inconsistent.The systems that were consistent with disk fragmentation had apparent rotating circumbinary structures surrounding the binary/multiple system.However, even with this evidence, turbulent fragmentation could not be completely ruledout.This is because a system formed from turbulent fragmentation could migrate inward, interact, and form a close multiple system with a new circumbinary disk (e.g., Bate 2018).
However, the compact circumstellar disks around each component of the close multiple systems could provide more definitive evidence on the formation mechanism.If the close multiple system is formed via disk fragmentation, the angular momentum axis of each component and the circumbinary disks should be relatively aligned due to forming within a common disk with the same net angular momentum (Offner et al. 2016;Bate 2018).On the other hand, Lee et al. (2019b) and Offner et al. (2022b) found that in simulations of turbulent core collapse, multiples preferentially formed as randomly aligned systems whose (mis-)alignment angles would persist throughout the calculation, suggesting randomly distributed alignment angles should be a signature for this formation mechanism.
These primordial alignments can further evolve via dynamical interactions to either mis-align or realign at later stages.Highly inclined alignments with respect to the orbital plane can quickly decay, on order of 10 4 -10 5 years, while more moderate misalignments decay over 10 6 -10 7 years (typical of 1 M ⊙ companions on ∼100 au separations, e.g., Bate et al. 2000).Further studies of the more evolved protoplanetary disks are consistent with a greater likelihood of misaligned companions but these system probably experienced many orbital timescales, thus undergoing dynamical evolution (Jensen et al. 2004;Jensen & Akeson 2014;Rota et al. 2022).Furthermore, Larwood et al. (1996) found the disk inclination angle and orbital angular momentum axis evolve on timescales of order the viscous timescale, thus it is not clear if the observed relative alignment angles are inherited from the formation or if the angles are a product of the dynamical evolution.Making surveys of the circumstellar disks at such young evolutionary stages critical for breaking the degeneracy between the primordial alignments versus the more evolved protoplanetary disks.
At such a young age, protostars are deeply embedded, making direct measurements of the stellar rotation axis impossible.Bate et al. (2010) found that the stellar rotation axis (the inferred stellar angular momentum axis) would not differ significantly from the inner disk rotation axis (<5 • ).The outflow position angle can be a proxy for the angular momentum axis, can be difficult to separate for compact systems.These outflows may be entangled and/or misaligned from the rotation axis due to n-body interactions (Ohashi et al. 2022).Tobin et al. (2018) confirmed and resolved the companion separations, but given the resolution were not able to resolve the compact circumstellar disks or the individual outflows around the sources.Furthermore, Segura-Cox et al. (2018) conducted high resolution VLA observations and identified Per-emb-5 potential compact multiple, requiring highsensitivity follow-up to confirm the multiplicity.Thus high-resolution and high-sensitivity studies of compact and wide multiples are needed to accurately determine the angular momentum vectors of the disks.
We carried out a novel method of empirically testing multiple protostar formation mechanisms by observing 12 known wide and compact multiples with 25 mas angular resolution (∼8 au), with the ability to resolve small protostellar disks.This type of survey can best recover the projected disk rotation axis, the implied orientation vector, of the nascent circumstellar disk.Similar studies have been conducted in the past, observing the polarization angle of the sources to infer the circumstellar disk alignments but are highly sensitive to the amount of intervening interstellar polarization alignments which may contaminate the resulting distribution of angles Jensen et al. (2004).
We present our findings of 12 protostellar multiple systems within the Perseus molecular cloud and detail the collection of observations in Section 2. We present our analysis of the protostellar sample in Section 4. Further, we interpret our findings in the broader aspects of star formation and with the specific sources of the sample in Section 5.

DEFINITIONS
For consistency we utilized the following definitions for this work: 1. source: A single protostar that can have a compact and/or extended disk.
2. system: A collection of sources within a defined separation that may be interacting.

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

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).

OBSERVATIONS AND DATA ANALYSIS
The Atacama Large Millimeter/submillimeter Array (ALMA) is a state-of-the-art interferometer located on Llano de Chajnantor in the Atacama region of Chile at an elevation of ∼5000 meters).We conducted observations of protostellar multiple systems in Perseus primarily using Band 6 (1.3 mm) with some supplementary observations in Band 3 (3 mm).

Band 6 (1.3 mm) Observations
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∼300 pc).The observations were carried out in the most extended configuration C-9/10 (baselines 150 m∼14.9 km) and have an effective angular resolution of 42 mas×23 mas to 66 mas×29 mas, with a continuum sensitivity of 13-76 µ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.19MHz bandwidth and centered on the 13 CO (J = 2 → 1), C 18 O (J = 2 → 1), SO (J = 6(5) → 5(4)), and SiO (J = 5 → 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 timeon-source ≈9 minutes (scans∼9.04s).The remaining 7 sources were observed in 3 execution blocks (EBs), across 3 days, with the average time-on-source ∼14 minutes for the 3 EBs combined.While the absolute flux density scale is expected to be accurate to ∼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∼40 µ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, Note-L1448IRS3B also contains L1448IRS3A within the field of view.SVS13A also contains SVS13B within the field of view.Per-emb-18 also contains Per-emb-21 within the field of view.
Note-RMS is specified as the root-mean-squared of the Briggs robust 1 weighted image with the upper 95% of emission clipped.The beam specified is the self-calibrated, multi-frequency synthesis clean beam using Briggs robust 1 weighting.S/N is the signal-to-noise ratio defined as the emission peak divided by the respective RMS.Class, T bol , and L bol are given in Tobin et al. (2018); Connelley et al. (2008b).

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 ∼30 µ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×9000 and 4 mas per pixel.All images are shown with square-root stretch colormaps and a common RMS value of 20 µ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× 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.
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 13 CO and C 18 O (J = 1 → 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 ∼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 selfcalibration 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×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 ∼1.5× the noise in each image (Figure 3).

Gaussian Fitting the uv-visibilities
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 Note-Per-emb-18 also contains Per-emb-21 within the field of view.
Note-RMS is specified as the root-mean-squared of the superuniform weighted image with the upper 95% of emission clipped.The beam specified is the self-calibrated, multi-frequency synthesis clean beam using superuniform weighting.S/N is the signal-to-noise ratio defined as the emission peak divided by the respective RMS.
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 σ.
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× 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λ) 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 σ 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.Note-The separations are given in units of au, which is derived based on the average distance to the Perseus molecular cloud of 300 pc.Companion separations are defined as the distance from the first target listed.
Note-The uv -plane nested sampling results.The models fit multiple Gaussians to the uv data rather than the image-plane, which is sensitive to image reconstruction efforts and the specific beam shape, particularly towards marginally-resolved or un-resolved sources.

RESULTS
With these observations, we detected all the circumstellar disks toward each multiple system within the survey at ∼8 au resolution for the 1.3 mm observations and ∼26 au resolution for the 3 mm observations, and most circumstellar disks are at least marginally-resolved.The ALMA images are shown in Figure 1 with the respective beam sizes in the lower-right corner.While we resolved out much of the >100 au (>0.′′ 3) scale disk structures previously resolved (Tobin et al. 2018), we did recover a large variety of disk sub-structures never previously resolved toward these sources.We briefly detail some key observations here and further discuss the morphologies of the individual sources in Appendix A.

PER-EMB-2
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. (2015aTobin et al. ( , 2018) ) as a smooth continuum surface brightness distribution at 1.3 mm (Figure 5).We resolved the compact binary (a∼80 mas≈ 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 ∼0.′′ 431 (∼129 au), ∼1.′′ 432 (∼430 au), and ∼0.′′ 5 (∼150 au) for the southern (SNR≈25), northern-most (SNR≈22), and southern-most (SNR≈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 σ 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
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
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 Per-emb-2 Figure 5. Left side is a 2. ′′ 5 image of the Per-emb-2 system whose inner binary has been resolved.The right side is a 2× zoom in on the source.The top row is the ALMA 1.3 mm observations and the bottom row is the ALMA 3 mm observations.In both observations, we resolve the compact inner binary but also report potentially 3 additional companions, 2 more within the disk and the southern-most companion falls just on the edge of the image.This would make Per-emb-2 a possible 5 companion protostellar system.A 0. ′′ 5 (150 au) scalebar is shown in the left panels and a 0. ′′ 25 scalebar in the right panels.The beam is located in the lower right of both images.The colormap is square-root scaled.
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.326s30: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.

L1448 IRS3A
Within the same field-of-view as the L1448 IRS3B observations, we resolved the disk around L1448 IRS3A (Figure 8).Previous ALMA 870 µm observations hinted of substructure in the surrounding circumstellar disk of L1448 IRS3A, but were unable to determine the nature of the features (Reynolds We now resolve the inner disk to be a ring and the so called "clump" (see Figure 2; Reynolds et al. 2021), is the north-east side of the ring.We report the L1448 IRS3B-A source is now resolved as just inside of the ring and L1448 IRS3B-B is now resolved as just outside of the ring.The separation between the two sources is well resolved and the L1448 IRS3B-B source is resolved.Towards the geometric center of the ring there is compact emission.We also reproduce the kinematic centers as the numbers "1" and "2" in the plot, as reported in 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.

Companion Finding
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 ∼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(Tobin et al. , 2018) ) and the algorithm is discussed in detail in Appendix B.

Geometric Orientations
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( σ minor σ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.

Models of Companion Orientations
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 (±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.

Statistical Tests
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.75 UC 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 σ), such that if the probability test can reject the null hypothesis at the 3 σ 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(Tobin et al. , 2020;;Encalada et al. 2021;Tobin et al. 2022) found the average separation of companions to be ∼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 ∼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 uvvisibility analysis, but the results from the uv -plane analysis tend to be favored as the uncertainties in the fits are more constrained.

Formation Pathways
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 ∼100s au scales in massive disks ( M d M * ≳ 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 ≈ 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.

Formation Mechanism of Individual Systems
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 ≈100s au) and the kinematics of the system to be organized.Whereas, formation of system via turbulent fragmentation happens on much larger scales (a ≈1,000-10,000s au).The sources formed at larger scales could further migrate to ∼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.

Disk Fragmentation Candidates
We detected circum-multiple material in nine of the 11 observed sources (Per-emb-5 is now classified as a single source) in these observations (Figure 2); however our array configuration is less sensitive to extended emission.In all detected cases, the circum-multiple material was observed around Class 0 sources, while sources that show significant misalignment do not show much circum-multiple material within the detection limits of our observations.In selecting a subsample consisting of only Class 0 sources from our full sample of sources, we found the median relative orientation angle between the sources in each system is 24 • for the Class 0 sources.
Seven of the observed multiple systems are consistent with having their <500 au companions formed via disk fragmentation.The sources Per-emb-2, Per-emb-17, Per-emb-18, Per-emb-22, L1448 IRS3B, L1448 IRS3C, and SVS 13A appear most consistent with the disk fragmentation formation pathway.
In particular for Per-emb-2, Pineda et al. (2020) found the mass-infall rate exceeds that of the accretion rate derived from the bolometric luminosity.This scenario provides favorable conditions to trigger gravitational instabilities that lead to the fragmentation of the disk (Kratter et al. 2010a).
These compact circumstellar disks appear to be strongly aligned and deeply embedded Class 0 systems, with bright circum-multiple disks (in some cases).

Turbulent Fragmentation Candidates
There are several other systems with separations a <1,000 au which are not consistent with disk fragmentation and likely formed via a combination of other methods.The sources L1448 IRS1 and NGC1333 IRAS2A are relatively compact proto-multiple systems, with strongly misaligned disks.
The companion (a ≈ 418 au) of L1448 IRS1 is nearly orthogonal to the brighter primary source.This source could have formed via turbulent fragmentation and then evolved via dynamical interactions to more compact scales (a∼1.′′ 4≈420 au).
The companion of NGC1333 IRAS2A is nearly orthogonal to the brighter primary source with a much higher inclination within the limits of our observations.Our observations for the disk position angle are consistent with an angle orthogonal to prior observations of the outflow position angles (Tobin et al. 2015b).Similar to L1448 IRS1, this source could have formed via turbulent fragmentation and then evolved via dynamical interactions to more compact scales (a∼617 mas≈185 au).
Per-emb-35 appears to be two similar sources, with aligned relative disk orientations, at a wide separation of 1. ′′ 9 (∼570 au).While the orientations are aligned with an average relative orientation angle of ∼19 • , the source is not likely to have formed via gravitational instability.We report no detection of circumbinary material or any dust continuum between the two sources.This coupled with the wide separation means Per-emb-35 is likely to have formed via large-scale (1000s of au) mechanisms, likely turbulent fragmentation or dynamical capture, and migrated to more compact scales (100s of au).
The sources Per-emb-21, L1448 IRS3A, and SVS13 B are wide (a > 1000 au) companion sources and additionally appear to not be preferentially aligned towards the outer sources in the respective systems, and are thus likely formed via turbulent fragmentation.Figure 9. Left panel: corresponds to the subsample of compact companion separations (< 100 au); Middle panel: the full companion sample with separations (a<10,000 au); Right panel: subsample of extended companions (a>500 au).Median CDFs of the dot-product of the orientation vectors the two fitting techniques and of the various distributions of fractional alignment ratio.Each of the colored lines corresponds to a particular sample of fractional ratio distributions between a fully preferentially aligned distribution to fully randomly aligned distribution, with the given ratio in the colorbar, as detailed in Section 4.4.We resample the data with Gaussian errors considering the observation uncertainties and construct empirical CDFs.The solid black line is a median CDF for the CASA imfit results.The dashed black line is a median CDF for the uv -plane fit results.The red, horizontal error bars associated with each ECDF, represent the 1 σ uncertainty for the lower and upper bound, respectively.The compact (<100 au) subsample is consistent with preferentially aligned distributions, with >80% of all sources having a dot-product within 30 • .The full sample and extended companion subsample are consistent with a ratio of aligned orientations, with 60% of all sources having a dot-product within 30 • .The full statistical analysis is given in Table 5 and further discussed in Appendix E.

Ambiguous Formation
It is not readily apparent via which formation pathway NGC 1333 IRAS2B (a ≈ 95 au) formed.The orientation vectors are relatively aligned within 40 • but this could be happenstance as there is no circum-multiple material found in these or prior observations.The system could have formed at long separations and migrated to compact scales.

CONCLUSIONS
We have presented very high resolution observations (∼8 au) towards 12 Perseus protostellar multiple systems, resolving and detecting 32 sources within the field of view of the targeted sources.We observed the dust continuum at 1.3 mm and provide spatial resolution at 3 mm data as well.Our results can be summarized as follows: 1. We detected and confirmed all previously detected multiple systems from the VANDAM observations (Tobin et al. 2018(Tobin et al. , 2016a)), except one, and we detected a total of 32 sources.Per-emb-5 is the only source that we reclassified as a single source.Note-The fractional alignment column details the fractional ratio of aligned vs randomly aligned as indicated by the subscripts.% KS , % AD , and % ES correspond to the % of resampled observation dot products that cannot reject the null hypothesis in favour of the alternative for Kolmogorov-Smirnov, Anderson-Darling, and Epps-Singleton statistical tests respectively.The rejection criteria is evaluated at 0.3%.The total number of resampled observation dot products are 10,000.The null hypothesis tested is the empirical distributions and the corresponding constructed fractional distributions are drawn from the same underlying distribution.
2. We detected circum-multiple continuum emission in seven of the 12 systems consistent with other observations at lower resolution.
3. We statistically characterized our full sample of 11 Perseus protostellar multiples (N pairs =21), with separations <10,000 au, to be consistent with forming from a combination of gravitational instability and turbulent fragmentation pathways, with the underlying distribution described as 40%-80% preferentially aligned systems (or conversely 60%-20% randomly aligned systems).
4. If we select a compact subsample of the sources, with separations <100 au (N pairs =7), we found the underlying distribution of at least 40% preferentially aligned and we ruled out distributions that are randomly aligned dominated.
5. Futhermore, if we select an extended subsample of the sources, with separations >500 au (N pairs =15), we found the underlying distributions of 40%-80% preferentially aligned companions are equally likely.
6. Combining our statistical approach with prior observations, we determined 7 of 12 of our systems are more consistent with disk fragmentation; while 3 systems (and 3 wide companions) are more consistent with turbulent fragmentation.One system, NGC1333 IRAS2B has an ambiguous formation pathway.Determining the formation mechanism via the statistical approach gives the likelihood of the sample being consistent with some underlying distribution of aligned and misaligned disk; whereas, combining this approach with multi-band observations detailing larger scale structures and molecular line features will give the most holistic determination.
7. Towards Per-emb-2, we detected the previously reported compact (a =24 au) binary Per-emb-2-A/-B, 2 additional possible companions Per-emb-2-C/-D, and additionally a potential fifth companion Per-emb-2-E, embedded within the northern part of the disk.
8. Per-emb-5 is now resolved and appears as a ringed disk with a single spiral arm extending to the west.The double peaks as appeared in the VLA data were likely associated with bright peaks along the ring and the underlying disk structure was below the surface brightness limit of the VLA. 9. Toward L1448 IRS3B, we resolved the compact ring of the inner disk and resolved the previously detected IRS3B-A, -B as compact sources just inside and outside of the ring, respectively.Additionally, toward L1448 IRS3B, while faint, we confidently detected an additional continuum source at the geometric center of the inner disk ring.This might be emission surrounding a more central protostellar source that could dominate the gravitational potential given that it is near the two kinematic centers as described with the molecular lines C 17 O and C 18 O.
10. Towards L1448 IRS3A we resolved the disk to be a ring and an unresolved compact continuum source surrounding the central protostar at the geometric center of the ring.
These results, while unable to determine with certainty the relative contribution of the different formation pathways due to the low number of systems, suggest a path for characterizing multiple star formation with future surveys.To best utilize the methods described, the highest resolution surveys towards larger samples of multiples are required to confidently resolve the most compact scales of the circumstellar disks around the protostars.NKR, JJT, and NAK gratefully acknowledge support from NSF AST-1814762.LWL acknowledges support from NSF AST-2108794.ZYL is supported in part by NASA 80NSSC20K0533 and NSF AST-2307199.DMSC is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2102405.DMSC is grateful for support from the Max Planck Society.This paper makes use of the following ALMA data: ADS/JAO.ALMA#2019.1.01425.S and 2017.1.00337.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.We continue our discussion of Per-emb-2, a previously reported close multiple (a<50 au).We resolved the compact binary (a∼80 mas≈ 24 au) in both the 1.3 mm and 3 mm observations and from our observations further revealed a possible additional 3 compact sources.
Tobin et al. ( 2018) observed a velocity gradient in 13 CO and C 18 O centered on the compact binary of Per-emb-2-AB.Several other studies have also detected the outflow from Per-emb-2-AB (e.g., Tobin et al. 2015b;Yen et al. 2015;Stephens et al. 2018), but none have had the resolution or sensitivity to detect which of the two compact components is driving the outflow.Also, there do not appear to be detectable outflows originating from any of the companion sources.

PER-EMB-5
Previous Submillimeter Array (SMA) observations towards Per-emb-5 at 1.3 mm targeting C 18 O (2-1) and CO (2-1) molecular lines indicate a rich gas presence in the disk and envelope (Stephens et al. 2018;Heimsoth et al. 2022).The source also exhibits bi-polar outflows as evident in the SMA data, which are orthogonal to the major axis of the continuum disk in these observations.The C 18 O observations shown in Heimsoth et al. (2022) exhibit a velocity gradient that is in the same direction as the outflow with ∼3 ′′ (900 au) resolution, along the minor axis of the disk that we detect.Thus, the envelope may be influenced by the outflow (Arce & Sargent 2006), or there is infall from a flattened envelope (Yen et al. 2011).Higher resolution kinematic observations are required to determine at what scale the disk rotation is detectable.

L1448 IRS3B
Reynolds et al. ( 2021) observed C 17 O in the disk of IRS3B.C 17 O appears to trace well-ordered Keplerian rotation on the scale of the continuum disk (see Figure 4; Reynolds et al. 2021).Considering only the high velocity Doppler-shifted channels of the C 17 O emission, which should trace scales closest to the center of mass, the kinematic center of C 17 O coincides with the position of L1448 IRS3B-D.In Figure 7, we compile our observations at 1.3 mm and observations by Reynolds et al. (2021).We plot the kinematic centers of the C 17 O molecular line observations using the position-velocity diagram fitting technique (i.e."1") and the Bayesian analysis of kinematic flared disk uv -visibility modeling ("2") with the pdspy software (Sheehan 2022).We also indicate the position of L1448 IRS3B-A,-B, derived from previous observations, which we now resolved as two sources just inside and outside of the ring, respectively.Finally we draw a visual aid indicating the proposed geometric configuration of the inner disk, which now appears as a ring with a faint, but confidently detected, continuum source located at the geometric center of the ring.
NGC1333 IRAS2A -We resolved the binary of NGC1333 IRAS2A and resolved the brighter compact source NGC1333 IRAS2A-A.We also detected faint extended emission around NGC1333 IRAS2A-A at low surface brightness.The compact faint component NGC1333 IRAS2A-B does not appear to have much extended emission and is marginally resolved in these observations.The sources appeared oriented nearly orthogonal to each other in continuum.Bipolar outflows have been observed toward IRAS2A in several studies (Jørgensen et al. 2005;Plunkett et al. 2013;Codella et al. 2014;Tobin et al. 2015b;Jørgensen et al. 2022).The two outflows are nearly orthogonal to each other and appear to originate from the two components of the system.The outflow position angles are approximately orthogonal to the major axis of each presumed driving source.Other molecular gas emission has been characterized toward IRAS2A (Tobin et al. 2018).There is a possible velocity gradient in C 18 O, but the kinematics do not appear highly organized on scales >0.′′ 3.
Per-emb-17 -We resolved the compact disks toward each component of the binary (a∼86 au) and detected the faint circumbinary material; circumbinary gas emission was also previously detected by Tobin et al. (2018).The northern source, Per-emb-17-A, is brighter and less extended than the southern source, Per-emb-17-B, which appears appears more edge-on.The circumbinary gas detected in C 18 O and 13 CO show a clear velocity gradient in the same plane as the companions.A bipolar outflow is also seen by (Tobin et al. 2018), but the extended emission of the outflow is better recovered by observations from Stephens et al. (2018).The outflow position angle is orthogonal to the orientation of the binary system and the orientation of Per-emb-17-A.
Per-emb-18 -We detected both the extended circumbinary structure reported in Tobin et al. (2018) and resolved the separation between the compact binary (a∼32 au) as reported in Tobin et al. (2016a).The binaries are situated in the near geometric center of the circumbinary disk (Figure 10).There is an apparent surface brightness asymmetry that could appear as an azimuthal asymmetry if viewed more face-on (e.g., van der Marel 2015).When viewed at 9 mm the circumbinary disk appeared one-sided, which could be due to a combination of lower surface brightness sensitivity from the VLA and dust trapping due to a vortex created by the inner binary pair (van der Marel 2015).The binaries themselves are similar in physical structure and brightness between the 1.3 mm and the 3 mm observations.The circumbinary material is shown to have a clear velocity gradient in 13 CO, C 18 O, and H 2 CO in Tobin et al. (2018).These velocity gradients are orthogonal to the outflow traced by 12 CO (Stephens et al. 2018;Tobin et al. 2018;Heimsoth et al. 2022).2015a).The velocity gradient was more evident at ∼1 ′′ (300 au) resolution than at ∼0. ′′ 3 (90 au) resolution in Tobin et al. (2018).There is a clear outflow from Per-emb-22 (Stephens et al. 2018) that is at a ∼45 • angle with respect to the plane of the binaries.Higher resolution maps of CO from Tobin et al. (2018) seem to show that the more prominent outflow originates from Per-emb-22-B while there might be some evidence for a second outflow from Per-emb-22-A.
L1448 IRS3C -We resolved both disks of the compact (a∼72 au) binary L1448 IRS3C.We detected faint circum-binary continuum emission at the limit of our observations.The disks appear nearly aligned and the southern source is much brighter than the northern source.This system is also a wide companion to the L1448 IRS3A/IRS3B system, located ∼5300 au away.Tobin et al. (2018) reported observations of L1448 IRS3C with disk tracing molecules C 18 O and 13 CO, and found position angles of the velocity gradients corresponding to 220 • , consistent with the major axes of the disk we observed here.There is also a velocity gradient in the circumbinary emission, consistent with orientation of the binary, in both 13 CO and C 18 O.Outflows have been observed from L1448 IRS3C (Lee et al. 2015;Tobin et al. 2015a), but even the highest resolution 12 CO maps thus far from Tobin et al. (2018) only detect a single outflow.
SVS13 A -is a compact (a∼92 au) binary with prominent spiral arms that were also seen in Tobin et al. (2018) and Diaz-Rodriguez et al. (2022).We detected the circumbinary one-armed spiral continuum and resolved each compact circumstellar disk.We also detected the unresolved wide (a∼1830 au) companion SVS13 A2.Additionally, we detected the wide (a∼4770 au) companion SVS13 B. We resolved SVS13 B, but due to limited sensitivity far out in the primary beam, the image is very noisy.Tobin et al. (2018) found apparent rotation in the circumbinary spiral structure, while higher resolution data from Diaz-Rodriguez et al. (2022) found rotation across the two binary sources and circumstellar rotation from one.Stephens et al. (2018) found the outflow position angle to be 150 • ±10 • , which is approximately orthogonal to our reported disk position angle of ∼70 • .Tobin et al. (2018) found the disk position angle to be ∼220 • (180 • symmetry yields angle ∼50 • ) by evaluating the angle of the velocity gradient in the disk tracing molecules C 18 O and 13 CO.This is consistent with the average of the position angles for the two compact sources SVS13A-AB.
RAC1999 VLA20 -We detected the continuum source RAC1999 VLA20.However it is believed to be an extra-galactic source (Rodríguez et al. 1999;Tobin et al. 2016a) so the source is not considered for the purpose of the analysis conducted in this paper.A summary of the observed source is provided in Table 6.
A.2. Class I/II Protostars L1448 IRS1 -The bright disk of L1448 IRS1-A, the northern source, is well-resolved and is likely the most evolved source in our sample.This source has a bright, disk (25 ×12 mas≈75×36 au)..The fainter, southwest companion is marginally resolved orthogonal to the beam.Previous observations of L1448 IRS1 reported a disk position angle of 24 • based on continuum observations (Tobin et al. 2018), and we report a consistent position angle of 25 • .Analyzing the CO emission, which is tracing the disk toward L1448 IRS1-A, coincides with the major axis of the disk and is likely tracing only the gas component of the disk and not the outflows.
L1448 IRS3A -Within the pointing of L1448 IRS3B, we detected this wide (a∼2280 au), more evolved companion.Furthermore, we resolved the substructure of the continuum disk to be a clear circumstellar ring.The compact inner disk centered at the geometric center of the ring remains unresolved in these observations, and the relative orientation cannot be determined to be different with respect to the outer ring.
Previous observations of L1448 IRS3A reported the disk position angle to be ∼130 • (Tobin et al. 2018;Reynolds et al. 2021) based on the PA of the ring and apparent orthogonality to the outflows observed by Lee et al. (2015); Tobin et al. (2018).As such, for the purpose of the statistical analysis later in the paper, we assumed the inner compact disk orientation would be comparable to the orientation of the ring.Rotation of the outer disk has been traced in C 18 O, 13 CO, and C 17 O, and the inner disk, close to the scales of the inner compact continuum source, has been traced in SO 2 , consistent with the outer disk (Tobin et al. 2018;Reynolds et al. 2021).The origin of the ring is Note-The CASA imfit fitting results for the source RAC1999 VLA20, which was detected in the field of SVS13 A. This source is believed to be extra-galactic and not considered part of the SVS13 A system.
beyond the scope of this paper, but we can report that we do not detect a continuum source within the gap between the inner continuum source and the ring.
Per-emb-35 -We marginally resolved the continuum emission toward each source, which are separated by ∼644 au and report no obvious detection of any continuum emission between the two compact sources.
Observations of the the CO outflows in Stephens et al. (2018) show an outflow position angle of 123 • , and Tobin et al. (2018) detected what appear to be parallel outflows from each protostar.Due to the compact nature of the circumstellar disks, rotation within them is not clear, though Tobin et al. (2018) did detect compact blue and red-shifted H 2 CO and SO that could be from the disks.With our higher resolution 1.3 mm observations, we report a disk position angle of 31 • for Per-emb-35-A, which is consistent with being orthogonal to the outflow as reported in Stephens et al. (2018).
NGC1333 IRAS2B -The disk around the brighter southern source, NGC1333 IRAS2B-A is wellresolved.The companion source ∼95 au away, is marginally resolved and visually appears mis-aligned to the brighter source.Stephens et al. (2018) reported the outflow position angle of 24 • ±10 • , which is consistent with an earlier map by Plunkett et al. (2013).We reported the disk position angle of IRAS2B-A to be 104 • , which is nearly orthogonal with the reported outflow position angles.While the companion IRAS2B-B appears relatively aligned (θ ≈30 • ) with the brighter IRAS2B-A source.Clear rotation across the binary or toward either disk has not been observed, although Tobin et al. (2018) found emission to generally be concentrated northeast of IRAS2B-A, towards IRAS2B-B, with SO having the greatest likelihood of exhibiting a velocity gradient. 13CO and H 2 CO appear on the north-east side of IRAS2B-A while C 18 O appears almost entirely on the west side of IRAS2B-B.Step 6 Single: 1 Binary: 0 Triple: 1 Quad: 1 # Pairs: 9 Figure 11.Demonstration of the companion finding algorithm used to associated the observed sources within each system.In this case, the companion maximum distance is 500 au.The algorithm first attempts find the most compact pair within the maximal defined distance and groups them together, re-inserting the geometric center of the group into the source list.The algorithm then finds the most compact pair within the maximum distance allowed and groups these two.In the plot, all of the example sources are shown in red crosses.At "Step 1", the algorithm finds the most compact pair.In "Step 2" the next most compact pair is a companion to the first binary, thus making a triple source system.The algorithm continues until no more sources meet the distance criteria, arriving at a quad-source system, a triple-source system, and a single source system for this example.-90 • and in the red "+", the position angle orientations are shown between 0 • -180 • .The respective errorbars for each fit are also shown, where the uv -visibility errors correspond to the isometric 1 σ average uncertainty in the posterior.The one-to-one line demonstrating the two fit methods return the same result is plot in black.There appears to be no clear systematic difference for the two fitting schema, as all sources are within 3 σ of the one-to-one line..The uv -visibility fit would be able to achieve higher resolution fits as compared to the image-plane fit, thus the uv -visibility fit is likely more consistent with the true compact disk orientation.

Comparison of Orientation Angles
Figure 13.Comparison of the uv-visibility fit and imfit derived dot-product orientations for each of the sources.The red errorbars represent the 1σ confidence intervals while the black points represent the median of the resampled distributions.The one-to-one line demonstrating the two fit methods give comparable results is plot in black, with the upper left half corresponding to targets that appear more aligned with the uv-visibility fitting while the lower right half corresponding to targets that appear more aligned with the image-plane fitting.There appears to be no clear systematic difference for the two fits.The only outlier system L1448 IRS3A.Note-Same as Table 5, except using the image-plane results from the deconvolved imfit task from CASA.

X
Note-The separations are given in units of au, which is derived based on the average distance to the Perseus molecular cloud of 300 pc.
Note-Compiled comparison table between the uv -plane derived results and the image-plane derived results for each of the sources.The "source A" and "source B" columns correspond to each pair where the comparison was performed to derive the orientation vector correlations.Where "source A" is not listed, the previously listed "source A" is to be assumed.All separations are defined with respect to the "source A".The dot-product is given in cosine coordinates.The final column, "Aligned (<30 • )" compares the derived dot-product values by averaging the two values and if the value corresponds to <30 • , the sources are considered aligned, while values >30 • are considered misaligned.
E. STATISTICAL TESTS Across all scientific fields, many statistical tests have been used to test whether two independent empirical samples were drawn from the same underlying distribution.This procedure of testing the the distributions usually involves solving the inverted hypothesis, whether the "null hypothesis" can be rejected.The null hypothesis being that the two distributions are drawn from the same underlying distribution against the alternative hypothesis, that the two distribution are not drawn from the same underling distribution.If we reject the null hypothesis at some significance level, then we can confidently say the underlying distributions are not identical at that significance level.However, careful consideration is needed as proving the inverse is not necessarily proving the distributions are drawn from the same underlying population.That is to say if we cannot reject the null hypothesis in favor of the alternative, that simply means the underlying distributions are consistent with being drawn from the same distribution at a particular significance level, and does not state the two distributions themselves are identical.
The most well-known statistical test is the Kolmogorov-Smirnov test (KS test), in the case of k-samples, tests the empirical cumulative distribution functions (ECDFs) of each sample against the null hypothesis that they were drawn from the same distribution (under the assumption one distribution is continuous).The KS test is sensitive to changes in the mass center and shape of the ECDF by computing the maximum difference of the distributions (higher sensitivity towards CDF center).The Anderson-Darling test (AD test) belongs to a family of quadratic rank tests and thus places more weight on the differences in the tails of the distributions.However, both the KS and AD tests assume the samples are drawn from continuous distributions and that the underlying distribution is stochastically larger than the drawn distribution.This is the so-called directional alternative hypothesis.
The Epps-Singleton test (ES test), is an empirical rank test that compares the characteristic functions (CFs) of the distributions, not the distributions themselves as compared to the KS and AD test.The test performs a quadratic form of differences between the CFs.The ES test is more powerful in that it does not assume the directional-alternative hypothesis or the continuous criteria, as it can be applied to discrete distributions.The primary prerequisite for the ES test, however, is the samples are fully independent, both across samples and within the samples (Epps & Singleton 1986;Goerg & Kaiser 2009), but otherwise is shown to provide a more robust statistical test.
For the purpose of our analysis, we apply all three tests, KS, AD, and ES test (as implemented in scipy), with k-samples, against the empirical distributions of the Perseus samples and the suites of mixture models.To satisfy the conditions of the KS and AD tests' directional alternative hypothesis and continuous prerequisite, we construct well-sampled models of 10,000 systems (each system comprised of at least 2 sources).To statistically characterize the observations, we resample all of the observation parameters (α, δ, σ major , σ minor , i, and p.a.) with Gaussian errors representative of the uncertainties of the observations 10,000 times.We then compute the statistic and probability for each of the KS, AD, and ES tests for every resample of the observations against every constructed distribution of fractional ratios.We finally compile the number of resampled distributions that cannot reject the null hypothesis at the 0.3% (3 σ) significance level and provide these numbers in Table 5.This yields six numbers per mixture ratio suite that describes the number of distributions that cannot reject the null hypothesis.

Figure 1 .Figure 2 .
Figure 1.ALMA 1.3 mm continuum images of the Perseus multiples, constructed with a Briggs robust weighting parameter of 1.The sources are contained within 12 pointings detailed in Table1.The distance given underneath the source name indicates the distances of that companion to the primary pointing source (the prior source in the list).If no distance is given, the center of the image is near to the center of the primary pointing.A 0. ′′ 5 (150 au) scalebar is shown in the lower left, and the respective restoring beam is shown in the lower right.The colorscale is square-root scaled, with the lower bound set by a common RMS value of 20 µJy beam −1

Figure 3 .Figure 4 .
Figure3.Continuum images at 3 mm of the Perseus multiples.Per-emb-2, Per-emb-5, and Per-emb-18 are constructed with superuniform weighting with an average restoring beam of 0. ′′ 09×0.′′ 04; while Per-emb-21 is constructed with Briggs robust 0.5, having a restoring beam of 0. ′′ 11×0.′′ 05, in order to recover the source from the noise.The sources are contained within 3 pointings detailed in Table3.The distance given underneath the source name indicates the distances of that companion to the primary pointing source (the prior source in the list).If no distance is given, the center of the image is near to the center of the primary pointing.A 0. ′′ 5 (150 au) scalebar is shown in the lower left and the respective restoring beam is shown in the lower right.

Figure 6 .Figure 7 .
Figure 6.Multiwavelength observations of the Per-emb-5 system.The upper panels are naturally weighted clean images from the 1.3 mm ALMA data.The middle panels are from Tobin et al. (2016a) robust 0.25 VLA 9 mm.The bottom panels are super-uniform clean images from the 3 mm ALMA data set reported here.The disk is asymmetric and appears as a ring and a single-arm spiral structure.We overlay an ellipse and line to assist in defining the substructure of Per-emb-5.A 0. ′′ 1 (30 au) scalebar is shown in the lower left and the representative beam is in the lower right.The white contours start at the 5 σ level and iterate by 20 σ, where the 1 σ = 50 µJy beam −1 .The colormap is square-root scaled.

Figure 8 .
Figure7.The left side is the 2 ′′ view of the system L1448 IRS3B, in which you can clearly see the triple system and the spiral arm substructures.Right side is a ∼4× zoom in on the inner disk of the binary L1448 IRS3B-AB.The top row is reproduced fromReynolds et al. (2021) using ALMA 870 µm continuum emission.The upper-middle row is the current 1.3 mm, 8 au resolution observations.The lower-middle row has the 870 µm data as the background and the 1.3 mm data is overlayed as white contours.The lower row is the same 1.3 mm data but with diagrams to aid in visualizing the likely configuration of the disk.We now resolve the inner disk to be a ring and the so called "clump" (see Figure2;Reynolds et al. 2021), is the north-east side of the ring.We report the L1448 IRS3B-A source is now resolved as just inside of the ring and L1448 IRS3B-B is now resolved as just outside of the ring.The separation between the two sources is well resolved and the L1448 IRS3B-B source is resolved.Towards the geometric center of the ring there is compact emission.We also reproduce the kinematic centers as the numbers "1" and "2" in the plot, as reported inReynolds et al. (2021).The 0. ′′ 25 (75 au) scalebar is shown in the lower left and the representative beam is in the lower right.The white contours start at the 3 σ level and iterate by 10 σ, where the 1 σ = 20 µJy beam −1 .The colormap is square-root scaled.The contour representative beam is overlayed in black at the lower right.

Figure 10 .
Figure10.Left side is an 1. ′′ 0 image of the Per-emb-18 system, and the right column is a ∼4× zoom-in.The upper panels are Briggs weighted with robust parameter 1, clean images from the 1.3 mm ALMA data.The bottom panels are super-uniform clean images from the 3 mm ALMA data.We resolve the compact inner binary of Per-emb-18 and resolve the circumbinary disk.We report in both 1.3 mm and 3 mm observations the circumbinary disk is asymmetric in flux, with the east side being enhanced as compared to the west side.The compact inner binary is located at the geometric center of the circumbinary disk.The compact disks of the individual sources remain unresolved in both observations.A 0. ′′ 1 (30 au) scalebar is shown in the lower left and the representative beam is in the lower right.The colormap is square-root scaled.

Figure 12 .
Figure12.Comparison of the uv -visibility fit and imfit derived inclination and position angle orientations for each of the sources.In the cyan, "Y" the inclination orientations are shown between 0 • -90 • and in the red "+", the position angle orientations are shown between 0 • -180 • .The respective errorbars for each fit are also shown, where the uv -visibility errors correspond to the isometric 1 σ average uncertainty in the posterior.The one-to-one line demonstrating the two fit methods return the same result is plot in black.There appears to be no clear systematic difference for the two fitting schema, as all sources are within 3 σ of the one-to-one line..The uv -visibility fit would be able to achieve higher resolution fits as compared to the image-plane fit, thus the uv -visibility fit is likely more consistent with the true compact disk orientation.

Table 2
Note-# E.B. is the number of execution blocks.

Table 7
Note-The Gaussian image-plane fit results utilizing the CASA imfit routine.The sources were given estimates to begin the fitting routine but were otherwise not restricted by any bounding values.Note-The separations are given in units of au, which is derived based on the average distance to the Perseus molecular cloud of 300 pc.Companion separations are defined as the distance from the first target listed.aThevalues for PA are defined as the angle of the major axis, oriented east-of-north.The values for σ