ROTATIONALLY DRIVEN FRAGMENTATION IN THE FORMATION OF THE BINARY PROTOSTELLAR SYSTEM L1551 IRS 5

Figure 1(a) (color map, along with contours plotted at varying intervals of σ) and Figure 1(d) (the same color map, along with contours now plotted at uniformly spaced percentage levels of the peak intensity) show the map obtained for L1551 IRS 5 using natural weighting, yielding an angular resolution of 55.6 mas × 52.7 mas (7.8 au × 7.4 au). Two sources are apparent, separated by 036 and aligned almost exactly along the north–south direction. We shall henceforth refer to the more northerly source as the N source, and the more southerly source as the S source. Both sources are clearly elongated along the north–west to south–east directions. In addition, relatively weak protrusions can been seen extending away from each source toward the north–east and south–west, oriented in approximately the same directions as a pair of ionized jets imaged at 3.5 cm with the VLA by Rodríguez et al. (2003b). Rodríguez et al. (2003b) infer a position angle of 67° ± 3° for the northern and 55° ± 1° for the southern jet, as indicated by the arrows in Figure 1(a). Notice that the protrusion on the south-western side of the S source lies northwards of the inferred jet axis for this component. The same offset also is apparent for this part of the jet in its image at 3.5 cm as pointed out by Rodríguez et al. (2003b). Finally, even weaker protrusions can be seen on the northern and south-eastern sides of the S source, and perhaps (at close to marginal significance) also on the northern and southern sides of the N source.

Zoom In Zoom Out Reset image size

Figure 2(a) (color map, along with contours plotted at varying intervals of σ) and Figure 2(d) (the same color map, along with contours now plotted at uniformly spaced percentage levels of the peak intensity) show the map obtained using robust weighting, providing nearly a factor of two higher angular resolution of 33.4 mas × 31.7 as (4.7 au × 4.4 au) but a factor of about three poorer sensitivity. Apart from the strongest protrusion on the north-eastern side of the S source (related to its jet), all the other protrusions seen in the naturally weighted map are not detectable in the robust-weighted map. As before, both sources are clearly elongated along the north–west to south–east directions. At or close to their centers, both sources also appear to be elongated along the north–east to south–west directions, quite closely aligned with the inferred orientations of their ionized jets.

Zoom In Zoom Out Reset image size

In our previous observation of L1551 IRS 5 at 7 mm with the VLA (before its upgrade), we attained a far inferior sensitivity along with a somewhat poorer angular resolution, despite utilizing the Pie Town antenna as part of the array (Lim & Takakuwa 2006). In that observation, we detected a relatively compact feature at a 5σ significance level southeast of the N source. In the observation reported here, apart from the N and S sources, we detect no other source in our field of view (extending far beyond the binary system to a full-width half-maximum of 1'). Although we cannot rule out a transient source, it seems more likely that the weak and relatively compact feature detected in our previous observation is an artifact produced by the poorer quality of the calibration possible at the time (compared with the case here where, because of the much higher instantaneous sensitivity, corrections for variations in atmospheric phase can be made on much shorter timescales and at a much higher precision).

As pointed out by Rodríguez et al. (1998), at 7 mm the emission from the binary components in L1551 IRS 5 is dominated by dust from their circumstellar disks. When observed at a sufficiently high sensitivity, however, free–free emission from their ionized jets also is detectable at 7 mm, as found by Lim & Takakuwa (2006) and as is clearly apparent in Figure 1(a) and to a lesser extent also Figure 2(a). To determine the geometry of the circumstellar disks, one of the main goals of this work, we therefore have to remove the contributions from the ionized jets.

As can be seen in the naturally weighted map of Figure 1(a) in particular, the ionized jets from both protostars have a complicated structure. Fortunately, to separate the ionized jets from the circumstellar disks, we only have to subtract the portion of the jet projected against each circumstellar disk rather than the entire jet structure. As we shall show, a two-component, two-dimensional Gaussian function fitted simultaneously to each source provides an adequate decomposition of the projected portions of the jets from the respective circumstellar disks. The parameters of the dual Gaussian components that provide the best fit to each source in the natural- and robust-weighted maps are listed in Table 2. These components, convolved with the corresponding synthesized beams (also two-dimensional Gaussian functions with parameters as listed in Table 1), are plotted in contours in Figures 1(b) and 2(b), overlaid on the color maps of Figures 1(a) and 2(a), respectively. To provide a simple visual comparison of the relative sizes of the Gaussian components fitted to each source, the contour levels for each Gaussian component are plotted at the same percentage levels of the peak intensity of that component. As can be seen, for each source, one of the fitted Gaussian components has its major axis closely aligned with the inferred orientation of the jet from that source (indicated by the arrows in Figures 1(a) and 2(a)), and must therefore correspond to the inner portion of the jet. The other component has its major axis closely orthogonal to the inferred jet orientation, and must therefore correspond to the circumstellar disk.

The residuals after subtracting the two-component, two-dimensional Gaussian components fitted to each source are shown in contours in Figure 1(c) for the natural- and Figure 2(c) for the robust-weighted map, again overlaid on the color maps of Figures 1(a) and 2(a), respectively. Although significant residuals (≤−3σ and/or ≥3σ) are apparent in Figure 1(c), these residuals coincide primarily with protrusions outside the main bodies—and therefore also circumstellar disks—of the N and S sources. Within the main bodies of the two sources, all the residuals are <4σ. Thus, in the naturally weighted maps, the fitted Gaussian components provide a reasonably good decomposition of the projected portions of the jets from the circumstellar disks. In Figure 2(c), the only significant residual (3.5σ) is that at the centroid of the S source. Thus, in the robust weighted maps, the fitted Gaussian components provide an excellent decomposition of the projected portions of the jets from the circumstellar disks. Subtracting the Gaussian components fitted to the jets leave the jet-subtracted images of the N and S sources shown in Figure 3 (first column), used in Section 4 for determining the physical parameters of the circumstellar disks.

Zoom In Zoom Out Reset image size

As a simple demonstration of the need to subtract the jets so as to obtain a good fit for the circumstellar disks, in Figures 1(e) and 2(e) we plot in contours just a single two-dimensional Gaussian function fitted to each of the two sources. Once again, these Gaussian components are plotted at the same percentage levels of their individual peak intensities. To provide a simple visual comparison of how well these components actually match the main bodies of the N and S sources, in Figures 1(d) and 2(d) we plot contour levels for the two sources at the same percentage levels of their individual peak intensities. As can be seen, there are clear differences between the fitted and observed structures in the central region of the N source as well as the north-eastern side of the S source. These differences correspond to the projected portions of the jets captured in the two-component fit to each source as described above. The same differences appear as relatively strong residuals in Figures 1(f) and 2(f) after subtracting the single Gaussian component fitted to each source.

As listed in Table 2, both the peak and integrated intensities of the Gaussian components fitted to the circumstellar disks are higher (by up to an order of magnitude) than those fitted to their corresponding ionized jets. Thus, at 7 mm, emission from the circumstellar disks dominate over that from the jets projected closely against these disks, consistent with that found previously. Furthermore, the centroids of the Gaussian components fitted to the circumstellar disks coincide with their respective source centroids. In contrast, the centroids of the Gaussian components fitted to the jets are clearly displaced from their respective source centroids in the natural- and, for the S source, also the robust-weighted maps. Finally, the circumstellar disks are resolved along their minor as well as major axes in both the natural- and robust-weighted maps. In contrast, only the southern jet is resolved along both its major and minor axes, and then only in the natural-weighted map.

Leaving all parameters in the NUKER function to be freely fitted, the GALFIT program was not able to find satisfactory fits for either circumstellar disk. On the other hand, by fixing α (which determines the sharpness of the transition between the inner and outer power laws) at different values, we found that the position angles of the major axis and inclination (as determined from the ratio of the minor to major axes) of both circumstellar disks were essentially constant independent of α in both the natural- and robust-weighted maps. Furthermore, these parameters are similar to those derived by fitting a two-dimensional Gaussian function to each circumstellar disk (to within the uncertainties of the Gaussian fits, as listed in Table 2). Evidently, although the radial intensity profiles of the circumstellar disks cannot be uniquely constrained, their position angles and inclinations are well determined. We therefore fixed these parameters at their formal values derived from the Gaussian fits.

When deriving their inclinations, we assumed implicitly that both circumstellar disks have circularly symmetric radial intensity profiles and are geometrically thin (i.e., a vertical thickness much smaller than the disk diameter). Pichardo et al. (2005) show that the radius at which orbits in the circumstellar disks significantly depart from near circular only occurs beyond about 80% of their tidally truncated radius. The Gaussian functions fitted to the circumstellar disks have radii at FWHM (the dimensions used to infer the disk inclinations) that are only about 60%–70% of their break radii found by fitting NUKER functions (Section 4.3). As a consequence, tidal distortions should not significantly affect the inclinations we infer for the circumstellar disks. More problematically, if the circumstellar disks are strongly flared, then their inclinations may be underestimated. Vertical settling of dust to the mid-plane, if it occurs, can mitigate this problem for disks traced in dust emission, as is the case here. As we will show in Section 6.2, the inferred inclinations of both circumstellar disks are systematically lower than the inclination of the orbital plane for many of the best-fit acceptable orbits, which span a wide range of eccentricities. The possibility that the circumstellar disks are strongly flared, and hence their inclinations are underestimated, may explain their apparent misalignment with the orbital plane. On the other hand, as we show in Section 5.1, both circumstellar disks are closely parallel to their surrounding flattened envelope, arguing against a large systematic error in their inferred inclinations.

As listed in Table 2, whether derived from the natural- or robust-weighted maps, the position angles of both circumstellar disks are identical to within the measurement uncertainties. From the natural-weighted map, the position angle is 1482 ± 07 for the circumstellar disk of the N source and 1477 ± 12 for the circumstellar disk of the S source, similar to within 05 ± 14. The uncertainties in the position angles are larger in the robust-weighted maps, where the position angles of the two circumstellar disks are similar to within 23 ± 74. The inclination of each circumstellar disk as derived from the natural- and robust-weighted maps is also similar to within the measurement uncertainties, although there is a small difference between the inclinations of the two disks in both these maps. In the naturally weighted map, the circumstellar disk of the N source has an inclination of 481 ± 05 and that for the S source an inclination of 443 ± 08, a difference of just 38 ± 09. The uncertainties in the inclinations are larger in the robust-weighted maps, where the inclinations of the two circumstellar disks differ by 91 ± 50. Henceforth, we shall adopt the position angles and inclinations obtained from the natural-weighted map, where the uncertainties in both these parameters are significantly smaller than those obtained from the robust-weighted maps; this map has a higher sensitivity than the robust-weighted map, and hence samples a larger extent of the circumstellar disks. The circumstellar disks of the binary protostars in L1551 IRS 5 are therefore very closely parallel, providing the first of several constraints that we will impose on models for the formation of this binary system.

In binary systems such as L1551 IRS 5 where ionized jets can be traced far beyond the circumstellar disks, the alignment between the jets as projected onto the sky plane can potentially serve as a sensitive probe of the alignment between the circumstellar disks. As mentioned earlier, from a high angular-resolution map at 3.6 cm where there is no apparent emission from the circumstellar disks, Rodríguez et al. (2003b) infer a position angle of 67° ± 3° for the northern and 55° ± 1° for the southern jet. If the jets emerge perpendicular to their respective circumstellar disks, then the major axes of the two circumstellar disks should have positions angles that differ by 12° ± 3°. Instead, we find that the major axes of the two circumstellar disks are aligned to within 05 ± 14.

As we will now show, the orientations that Rodríguez et al. (2003b) infer for the two jets may differ from the true orientations at their bases. At 3.6 cm, both jets in L1551 IRS5 are detected as a series of knots. The knots are not aligned along a single axis for either the northern or southern jets. As a consequence, to infer the position angles of the individual jets, Rodríguez et al. (2003b) fitted a two-dimensional Gaussian function to the knot closest to the individual protostellar components (the latter defined by the centroids of their circumstellar disks as measured by Rodríguez et al. 1998). The knots that Rodríguez et al. (2003b) selected are the same as those fitted by the Gaussian components corresponding to the jet of the N source in Figure 1(b) and the jet of the S source in Figures 1(b) and 2(b). As can be seen even in the map at 3.6 cm presented by Rodríguez et al. (2003b), the centroid of the knot closest to the S source is visibly displaced to the north–east of this source. We see the same displacement at 7 mm in both Figures 1(b) and 2(b), confirming that this knot does not trace the base of the southern jet. In the map at 3.6 cm presented by Rodríguez et al. (2003b), the centroid of the knot closest to the N source appears to coincide with this source. At the higher angular resolution of our natural-weighted map at 7 mm (Figure 1(b)), however, the centroid of this knot can be seen to be displaced to the south–west of the N source. Thus, this knot does not actually trace the base of the northern jet.

From our naturally weighted map at 7 mm (Figure 1(b)), we infer a position angle of 423 ± 89 for the northern and 272 ± 62 for the southern jet (Table 2) based on the same knots as those used by Rodríguez et al. (2003b). These values are significantly different (by 247 ± 94 for the northern and 278 ± 63 for the southern jet) from those determined by Rodríguez et al. (2003b) at a lower angular resolution at 3.6 cm. Figure 4 shows the structure of the knots closest to the N and S sources after the Gaussian components fitted to their circumstellar disks have been subtracted from the natural-weighted map of Figure 1. Evidently, the knot closest to the S source has a complicated structure. Indeed, this knot is resolved not only along its major but also its minor axes; thus, the two-dimensional intensity distribution of this knot affects the position angle inferred for its major axis. The same also is true for this knot at 3.6 cm, where it is resolved along both its major and minor axes (see Table 1 of Rodríguez et al. 2003b). Thus, in L1551 IRS 5, where the knots defining either of the two jets do not lie along a single axis and the knot closest to each protostar does not actually trace the base of the jet, we caution against using the apparent misalignment between the jets to infer a corresponding misalignment between their associated circumstellar disks.

Zoom In Zoom Out Reset image size

Fixing the position angles and inclinations of the circumstellar disks to the values listed in Table 3 (their formal values as derived from  Gaussian fits to the N and S sources and listed in Table 2), we then varied the value of α to investigate how the remaining free parameters changed with α. For the jet-subtracted image of the N source, we found that the solutions converge to a common set of values for γ, β, r_b, and I_b as α increases (i.e., a sharper transition between the inner and outer power laws). For the circumstellar disk of the S source, we found two families of acceptable fits for both the natural- and robust-weighted images. Only one of the two families in the natural-weighted image (but neither family in the robust-weighted image), converges to a common set of values for γ, β, r_b, and I_b as α increases. Table 3 lists the parameters from the NUKER fits to the jet-subtracted images the N source (both natural- and robust-weighted maps) and S source (natural-weighted map only, and for the fit that converges to a common set of parameter values) at the largest value of α for which a solution exists. The second column of Figure 3 shows images of the fitted NUKER profiles convolved by the synthesized beam. The third column shows the residuals, plotted in contours, after subtracting the fitted NUKER profile from the source. The last column shows the actual intensity profile of the fitted NUKER function along the major axis.

Although the fitted NUKER profiles have different values of γ and β in the natural- and robust-weighted maps for the circumstellar disk of the N source, the inferred break radius of this disk (where its brightness tapers off sharply) is quite closely comparable in both maps (a difference of less than 10%). As the natural-weighted image has a higher sensitivity and is therefore able to trace the circumstellar disk further out, we take the derived break radius of 87.4 mas (12.2 au) as a measure of the outer radius of the circumstellar disk of the N source. For the circumstellar disk of the S source, we take the derived break radius of 73.9 mas (10.3 au) as a measure of its outer radius. Consistent with the relative dimensions of the two disks as inferred from the Gaussian fits, the circumstellar disk of the N source is larger than that of the S source. The ratio of the outer disk sizes based on our power-law fits is 0.85, whereas that based on the FWHMs of the Gaussian fits is 0.75.

In Figure 5(a), we show a continuum map at 0.8 mm tracing dust in the envelope around L1551 IRS 5. This map was made from data obtained by Chou et al. (2014) with the SubMillimeter Array (SMA), but using only those data taken in the extended configuration so as to trace the inner region of the dust envelope. The data used were taken in 2013 January 14, close in time to our observations of L1551 IRS 5 with the VLA. The positions of the two protostars (as determined from our VLA observations and indicated by the two crosses in Figure 5) relative to the dust envelope can therefore be accurately located without having to correct for absolute proper motion. The emitting region traced in this map is about an order of magnitude smaller than the centrally condensed portion of the envelope traced in C¹⁸O by Momose et al. (1998).

Zoom In Zoom Out Reset image size

In Figure 5(b), we show in contours a two-dimensional Gaussian function with parameters listed in Table 4, fitted to the continuum source. The residuals after subtracting this Gaussian function from the map are shown in Figure 5(c). As can be seen, the fitted Gaussian function provides a good representation of the continuum source except around its central and outer eastern regions. With a radius at half intensity along its major axis of 072 ± 0.01, the continuum source extends well beyond the putative circumbinary disk (radius of ∼046). Based on its deconvolved dimensions at FWHM, we infer an inclination of 541 ± 10 for the flattened envelope, closely comparable to the inclinations of the circumstellar disks in L1551 IRS 5 (a difference of 60 ± 11 for the circumstellar disk of the N source, and 98 ± 13 for the circumstellar disk of the S source).

Just as we did for the circumstellar disks, when deriving the inclination of the envelope, we assumed implicitly that the envelope has a circularly symmetric radial intensity profile and is geometrically thin. Momose et al. (1998) investigated a number of simple analytical models to reproduce the observed position–velocity (PV-)diagram measured in C¹⁸O for the envelope around L1551 IRS5. They found that the observed PV-diagram is better fit by an envelope with a finite thickness rather than an infinitesimally thin disk. For a vertical thickness to radius of ∼1/5 (at a given isodensity contour; their preferred model for a sheet-like envelope), the inferred inclination of the envelope is only slightly changed to ∼60°.

From the same two-dimensional Gaussian function fitted to the continuum source, we determined a position angle for the major axis of the flattened envelope of 1566 ± 13. The position angle of the envelope is therefore closely comparable to those of the two circumstellar disks (a difference of 83 ± 12 for the circumstellar disk of the N source, and 89 ± 16 for the circumstellar disk of the S source). Given also their comparable inclinations, the circumstellar disks of the binary protostars in L1551 IRS 5 are closely parallel, not just to each other but also to their surrounding flattened envelope, providing a second constraint that we will impose on models for the formation of this system.

As described by Pichardo et al. (2005) and Pichardo et al. (2008), tidal torques exerted by both protostars cause orbits in their circumbinary disk to become eccentric. In a manner analogous to the circumstellar disks, the smallest non-intersecting orbit in the circumbinary disk defines the innermost stable orbit and hence the minimum size of a central cavity. In Figure 5(a), no central depression, let alone cavity, is evident in the image of the dust emission (or indeed in images of molecular-gas emission as presented by Chou et al. 2014). At an angular resolution of 080 × 052 in Figure 5(a), however, emission from the circumstellar disks cannot be separated from that of the circumbinary disk. Thus, even if the circumbinary disk has a central cavity, no central depression or cavity would be evident in the image shown in Figure 5(a).

The residual map of Figure 5(c) shows a positive residual close to the centroid of the continuum source. This positive residual is point-like, consistent perhaps with a contribution from the more intense circumstellar dust disk of the N source. Surrounding the central positive residual is a ring-like distribution of negative residuals. The major axis of this ring lies at the same position angle as the fitted Gaussian function. Evidently, away from the center, the dust emission has a profile that is not as strongly centrally peaked as a Gaussian function, which otherwise provides a good fit to the bulk of the dust emission. The presence of the ring-like negative residuals may indicate a central depression or cavity in the dust envelope. On the other hand, the intensity profile of the envelope may simply deviate from a Gaussian function in its inner regions.

Figure 6 shows the separation and position angle of the S source relative to the N source over the last ∼30 years. The most recent measurement, in 2012, is that reported here, based on the centroids of the Gaussian components fitted to the circumstellar disks of the N and S sources in the robust-weighted map. Emission from the jets is much weaker in this map than in the natural-weighted map, and is less likely to affect the determination of the centroids of the circumstellar disks. The second last measurement, in 2002, also was made at 7 mm with the VLA (Lim & Takakuwa 2006). Similarly, the third measurement, in 1997, was made at 7 mm with the VLA, although at a time when only about half of the array was equipped with the necessary receivers (Rodríguez et al. 1998). The remaining measurements, in 1983, 1985, 1995, and 1998, were all made at 2 cm with the VLA (Rodríguez et al. 2003a). As can be seen, the two most recent measurements provide the strongest constraints on the relative proper motion of the binary protostars, and largely motivate the full exploration of possible orbital solutions as described next in Section 6.2. These two measurements confirm the trend for an increasing separation between the binary protostars over time, and that the position angle of the S source as measured from the N source (anticlockwise from north) is decreasing over time (indicating clockwise orbital motion).

Zoom In Zoom Out Reset image size

The flattened envelope around L1551 IRS 5 has its major axis (at a position angle of ∼157°; Table 4) closely orthogonal to that of the ionized jets (at position angles of about 60°; Rodríguez et al. 2003b) and bipolar molecular outflow (at a position angle of about 45°; Moriarty-Schieven & Snell 1998) from the binary protostars as projected onto the sky plane. If the equatorial plane of the flattened envelope is indeed orthogonal to the outflow axis (i.e., so that its eastern side is the nearer side), then the kinematics measured in molecular lines imply that the envelope is rotating in the clockwise direction. In this case, as noted by Lim & Takakuwa (2006), the orbital motion of L1551 IRS 5 is in the same direction as the rotational motion of its surrounding envelope, providing a third constraint that we will impose on models for the formation of this system.

We fitted orbital solutions, at selected orbital periods as indicated by the data points in Figure 7, to the measurements of relative proper motion shown in Figure 6 for eccentricities, e, of e = 0 (circular orbit) and e = 0.1, 0.2, 0.4, 0.6, and 0.8 (elliptical orbits). Figures 7(a)–(b) shows the reduced χ² of the best-fit orbit(s) at the selected orbital periods for the different eccentricities considered, the latter as indicated by the different colors with corresponding eccentricities defined in the upper right corner of Figure 7(a). For circular orbits, there is a unique acceptable (i.e., reduced χ² ≤ 1) best-fit orbit at a given orbital period. For elliptical orbits, we find two families of acceptable best-fit orbits at a given orbital period: (i) where the system is leaving periastron or approaching apastron (family A), plotted in the left column of Figure 7; and (ii) where the system is leaving apastron or approaching periastron (family B), plotted in the right column of Figure 7. In family A (B), the system is closer to apastron (periastron) as the period decreases (increases). As can be seen, acceptable orbital solutions exist at all the different eccentricities considered. The orbital periods of acceptable circular orbits have a lower bound of ∼250 years, but no upper bound. The orbital periods of acceptable elliptical orbits have both a lower and an upper bound, with a range that increases with increasing eccentricity from just ∼40–520 years at e = 0.1 to ∼20–5000 years at e = 0.8.

Zoom In Zoom Out Reset image size

Figures 7(c)–(d) shows the semimajor axis of the best-fit orbits at  the selected orbital periods, with the different colors again corresponding to orbits with the different eccentricities considered. Over a broad range of orbital periods (∼20–5000 years), the orbital semimajor axis spans the relatively narrow range ∼03–15 (∼40–210 au). In family A, the semimajor axis is nearly independent of eccentricity at a given orbital period. In family B, the semimajor axis increases with eccentricity at a given orbital period. Figures 7(e)–(f) shows the corresponding total mass of the binary system as a function of orbital period, with the different colors corresponding to the different orbital eccentricities considered. For circular obits, the total binary mass attains a minimum value of ∼0.5 M_☉ at orbital periods in the range ∼600–1000 years. For elliptical orbits in family A, the total binary mass is confined to a relatively narrow range spanning ∼0.5–1.0 M_☉ over a broad range of orbital periods spanning ∼30–2000 years. For elliptical orbits in family B, the total binary mass increases with eccentricity, spanning the range from ∼1 M_☉ at e = 0.1 to ∼9 M_☉ at e = 0.8.

Figures 7(g)–(h) show the inclination of the orbital plane for the best-fit orbits at the selected orbital periods and eccentricities. We note that the inclination of the orbital plane is degenerate insofar as, without knowledge of whether the N or S source is approaching us, we do know which of the nodes correspond to the ascending and which to the descending node. (In the same way, we do not know on which sides of the circumstellar disk, along its minor axis, the near and far sides are.) In the same panels, the inclinations of the circumstellar disks derived from the natural-weighted maps (Section 4.1) are indicated by shaded strips (the widths of the strips reflecting the ±1σ uncertainty in the derived inclinations; Table 2). Locations where the shaded strips intersect solutions for the orbital inclination indicate orbits that have the same inclination as one of the circumstellar disks. For circular orbits, the orbital inclination is systematically higher than the disk inclinations by ∼15°–25° (based on their formal values) over the range in orbital periods spanning ∼300–1000 years. For elliptical orbits in family A, the orbital inclination differs from the disk inclinations by no more than about ±20° (based on their formal values), a difference that is usually much smaller than 3σ. For elliptical orbits in family B, the orbital inclination is systematically higher than the disks inclinations, increasing with eccentricity from ∼20° at e = 0.1 to ∼35° at e = 0.8.

Figures 7(i)–(j) show the position angle of the line of nodes  for the best-fit orbits at the selected orbital periods and eccentricities. The shaded strips indicate the position angles of the major axes of the circumstellar disks for the N and S sources, with widths reflecting the ±1σ uncertainties, as inferred from the natural-weighted maps (Table 2). Locations where the shaded strips intersect solutions for the position angle of the line of nodes indicate orbits that cut through the plane of the sky at the same position angle as the circumstellar disks. For circular orbits, the position angle of the line of nodes  differs from the position angle of the disk major axes by no more than ∼20° (based on their formal values). For elliptical orbits in family A, the position angle of the line of nodes differs by no more than ∼30° from the position angle of the disk major axes over the eccentricity range 0 ≤ e ≤ 0.4. This difference can be as large as ∼50° at e = 0.6, and as large as ∼70° at e = 0.8. As we will show in Section 6.3.1, however, constraints on the circumstellar disk sizes as imposed by tidal truncation limit acceptable orbits with e = 0.6 in family A to those where the position angle of the line of nodes axis is quite closely aligned (within ∼20°) with the disk major axes. For elliptical orbits in family B, the position angles of the line of nodes differs by no more than ∼20° from the position angle of the disk major axes.

In summary, although acceptable orbital solutions can be found over a broad range of eccentricities, all (omitting those in family A with eccentricities e ≥ 0.6 for the reason mentioned above) imply a relatively close alignment (mostly ≲25°) between the circumstellar disks and the orbital plane. The circumstellar disks of the binary protostars in L1551 IRS 5 are therefore quite closely aligned with the orbital plane, possibly essentially coplanar, although more probably tilted by up to typically ∼25° from the orbital plane. This condition provides the fourth constraint that we will impose on models for the formation of L1551 IRS 5.

Figure 8 shows examples of best-fit orbits all with a common orbital period of 500 years. The first and third columns show the projected orbits for families A and B, respectively. The second and fourth columns show the deprojected orbits (as seen along the orbital axis) for families A and B, respectively. At this orbital period, the orbit for e = 0.8 (not shown) is very similar to that for e = 0.6 in a given family. As can be seen, for all the orbits shown, the available measurements of relative proper motion span only a small fraction (less than one-tenth) of a complete orbit, explaining why the orbital eccentricity is poorly constrained. Despite this obvious limitation, the available measurements already constrain the best-fit orbits to being quite closely if not essentially coplanar with the circumstellar disks.

Zoom In Zoom Out Reset image size

6.3.1. Constraints from Circumstellar Disks

In a coplanar binary system, the radii of the two circumstellar disks combined obviously cannot be larger (and, indeed, must be significantly smaller) than the binary separation. Furthermore, for elliptical orbits, Pichardo et al. (2005) find that the radius of the outermost non-intersecting orbit in a circumstellar disk varies little with orbital phase, generally by only about 10%. Thus, in a coplanar binary system, the orbital separation at periastron (closest approach) sets a limit on the maximal possible sizes (as imposed by tidal truncation) for the circumstellar disks of the binary components. The tidally truncated radius of each circumstellar disk also depends on the mass ratio of the binary companions, with the more massive companion having a larger tidally truncated radius for its circumstellar disk. For a given mass ratio, the inferred sizes of the circumstellar disks in L1551 IRS 5 can therefore be used to rule out orbits with too small a periastron separation, a point we previously touched upon in Lim & Takakuwa (2006) and fully explore below.

To begin, we assume equal protostellar masses (i.e., q = 1) and hence equal maximum possible disk sizes. Note that Pichardo et al. (2005) define q = m_s/(m_p + m_s), where m_p is the mass of the primary and m_s the mass of the secondary, whereas we define q = m_s/m_p as is commonly used. The data points in Figures 9(a)–(d) indicate the predicted tidally truncated radii (averaged over an orbit) of the circumstellar disks (identical for the two protostars) for the best-fit orbits with the different eccentricities considered, the latter as indicated in the different colors. Like in Figure 7, the left column corresponds to the best-fit orbits in family A and the right column those in family B. At a given orbital period, the predicted tidally truncated radius decreases as the orbital eccentricity increases, as is expected. The break radius for the circumstellar disk of the N source is indicated by the horizontal line in Figures 9(a)–(b), and that for the circumstellar disk of the S source by the horizontal line in Figures 9(c)–(d). The data points lying above the horizontal line in each panel therefore correspond to orbits where the predicted tidally truncated radius of the circumstellar disk is larger than its observed break radius, whereas the data points lying below the horizontal line to orbits where the predicted tidally truncated radius of the circumstellar disk is smaller than its observed break radius. The curves joining the data points indicate the best-fit orbits for which  the predicted tidally truncated radius is, to within an uncertainty of 3σ, no smaller than the observed break radius. As can be seen,  the condition that the predicted tidally truncated radius cannot be smaller than the observed break radius rules out an increasingly larger range in orbital periods as the orbital eccentricity increases. For q = 1, the larger circumstellar disk of the N source imposes a more stringent constraint than the smaller circumstellar disk of the S source. For the former disk, the aforementioned condition can be met over nearly the entire range of orbital periods spanned by the best-fit orbits in family A only for e ≤ 0.1. This condition can be met over the entire range of orbital periods spanned by the best-fit orbits in family B only for e ≤ 0.2. For both families, this condition is satisfied over only a limited range in orbital periods spanned by the best-fit obits for 0.4 ≤ e ≤ 0.6, and cannot be at all satisfied for e = 0.8.

Zoom In Zoom Out Reset image size

Whether fitted by a Gaussian or a NUKER profile, the circumstellar disks of the N and S sources have different widths as measured at their FWHM (Table 2) or their break radii (Table 3). This difference raises the possibility that the two protostars in L1551 IRS 5 have different masses, resulting in different tidally truncated radii for their circumstellar disks. Based on the ratio in their break radii of about 0.85, the mass ratio of the two protostars according to the calculations of Pichardo et al. (2005) is about 0.72. For q = 0.72, the data points in Figures 9(e)–(f) indicate the predicted tidally truncated radius for the circumstellar disk of the N source, and those in Figures 9(g)–(h) the predicted tidally truncated radius for the circumstellar disk of the S source, for the best-fit orbits with the different eccentricities considered. Like before, the left column corresponds to the best-fit orbits in family A and the right column those in family B. As can be seen, at a given orbital period, the tidally truncated radius of the circumstellar disk for the more massive northern protostar is larger than that for the less massive southern protostar. Compared to q = 1, the tidally truncated radius for the circumstellar disk of the N source is larger (except for e = 0.6, due to the complex manner in which the tidally truncated radius changes with orbital phase), whereas that for the circumstellar disk of the S source is smaller. For both circumstellar disks, the changes can be quite small at high eccentricities so as to be unnoticeable. The break radius for the circumstellar disk of the N source is indicated by a horizontal line in Figures 9(e)–(f), and that for the circumstellar disk of the S source by a horizontal line in Figures 9(g)–(h). As before, the curves joining the data points indicate the best-fit orbits for which  the predicted tidally truncated radius is, to within an uncertainty of 3σ, no smaller than the observed break radius. Coincidentally, this condition rules out nearly the same range in orbital periods at a given eccentricity for both protostars as the circumstellar disk of the N source in the case where q = 1.

In making the aforementioned comparisons, we have assumed that the break radii of the circumstellar disks  provide a good measure of the lower bound on their tidally truncated radii. We now consider situations where the break radii of the circumstellar disks may be driven to either extend beyond or lie within their theoretically predicted tidally truncated radii. First, to continue to grow, both protostars must accrete through their circumstellar disks, which in turn must be replenished by material from the surrounding envelope. Accretion of infalling matter onto the circumstellar disks allows the disks to extend, albeit temporarily (within a disk dynamical timescale) if the accretion is intermittent, beyond their theoretically predicted tidally truncated radii. Theoretical simulations that assume circular coplanar orbits find that accretion onto the circumstellar disks of binary protostars is highly unequal, with some models favoring accretion onto the circumstellar disk of the secondary component (Artymowicz 1983; Bate 1997; Bate & Bonnell 1997) and others favoring the opposite (Ochi et al. 2005). For eccentric coplanar systems, Artymowicz & Lubow (1996) find that accretion is stronger onto the circumstellar disk of the secondary component. A recent study by Dunhill et al. (2015), however, paints a more complicated picture for accretion in eccentric coplanar systems. They performed hybrid smoothed-particle hydrodynamics (SPH)/N-body simulations of a binary system with a moderately high eccentricity (e = 0.6) and a moderately small mass ratio (q = 0.64) accreting from its circumbinary disk. The binary system clears a large cavity in its circumbinary disk and drives elliptical orbits at the cavity edge. Inflow from the circumbinary disk onto a given circumstellar disk is strongly modulated with orbital phase, such that accretion is strong only during a relatively short interval centered just before (∼0.1 orbital phase from) periastron. The binary system also drives precession in the elliptical orbits at the cavity edge. Accretion onto the circumprimary disk is much stronger than that onto the circumsecondary disk for half of this precession period (which is more than an order of magnitude longer than the binary orbital period), and vice versa. Averaged over time, accretion onto the circumprimary disk is similar to or slightly higher than that onto the circumsecondary disk, opposite to that found by Artymowicz & Lubow (1996). Despite their different predictions, all these studies find that accretion, when ongoing, is stronger onto one than the other circumstellar disk. Thus at least one of the circumstellar disks in L1551 IRS 5 is likely to have a size that is not affected by recent accretion. (We note that rotational periods at the predicted tidally truncated radii in the circumstellar disks of L1551 IRS5 are, for q ∼ 1, considerably shorter than the orbital periods of the acceptable best-fit orbits. In other words, dynamical times in both circumstellar disks are considerably shorter than the binary orbital period.)

Elliptical orbits in the circumstellar disks may be damped by gas pressure, providing an alternative means by which circumstellar disks can extend beyond  their largest non-intersecting orbit. If the dust-to-gas ratio is constant and the gas well coupled to the dust at all radii in the circumstellar disks, then the gas pressure ρ_gas ∝ r^−Γ. Equating Γ = γ, then given the inferred power-law index of γ ∼ 0.60 for the circumstellar disk of the N source and γ ∼ 0.81 for the circumstellar disk of the S source (Table 3), the gas pressure only decreases modestly with increasing radius. Although on this basis we cannot rule out that gas pressure may damp elliptical orbits in the outer regions of the circumstellar disks, we note that the break radii of both circumstellar disks are closely comparable with their predicted tidally truncated radii for the best-fit circular (for orbital periods <1000 years) or low-eccentricity (e ≤ 0.2) orbits. It seems to us most unlikely  that, if L1551 IRS 5 has has a circular or low-eccentricity orbit, gas pressure conspires to extend both its circumstellar disks so as to closely match their predicted tidally truncated radii.

Rather than extending beyond their predicted tidally truncated radii, it is more likely that one or both circumstellar disks do not extend to their tidally truncated radii because of depletion by accretion onto their respective protostars. In addition, or alternatively, magnetic braking can effectively shrink the circumstellar disks through the efficient outwards transfer of angular momentum (Machida et al. 2011). The possibility that the inferred break radii of the circumstellar disks are significantly smaller than their predicted tidally truncated radii rules out even more of the best-fit orbits derived for L1551 IRS 5 (i.e., the actual tidally truncated radii of the circumstellar disks lie above the horizontal lines drawn in Figure 9), thus rejecting most if not all orbits with relatively high eccentricities (e ≥ 0.4). In short, the break radii inferred for the circumstellar disks of the binary protostars in L1551 IRS 5 strongly favor orbits with relatively low eccentricities (e ≲ 0.2).

6.3.2. Constraints from the Envelope

We begin by summarizing the constraints imposed on models for the formation of L1551 IRS 5 as derived in Sections 4–6. From the measured geometry of its circumstellar disks and its surrounding envelope:

• The circumstellar disks of the binary protostars are very closely parallel and
• The circumstellar disks also are very closely parallel to their surrounding flattened envelope.

Based on the measured relative proper motion of the binary protostars and published measurements of the kinematics of its surrounding envelope:

• The orbital motion is in the same direction as the rotational motion of the surrounding envelope.

The close geometrical and dynamical relationship between the circumstellar disks and surrounding flattened envelope suggests that the angular momentum of material that gave rise to the protostellar system and that comprising its parental core are aligned.

Finally, although orbits spanning a wide range of eccentricities provide acceptable fits to the available measurements of the relative proper motion for the binary protostars in L1551 IRS5, theoretical predictions for limits to the sizes of the circumstellar disks as imposed by tidal truncation rule out the vast majority of best-fit orbits with moderate eccentricities (0.4 ≲ e ≲ 0.6) and all best-fit orbits with high eccentricities (e ≳ 0.8). For the remaining best-fit orbits:

• In the majority of cases, there is a systematic tilt (by up to typically ∼25°) between the circumstellar disks and the binary orbital plane; nevertheless, there is a clear preference for a relatively close alignment if not coplanarity.

As explained in Section 1, current models invoke either local (small-scale) turbulence in cores or the bulk rotation of cores to drive fragmentation. In cores that have little or no bulk rotation, turbulence introduces velocity and density inhomogeneities that can seed and drive the growth of multiple density perturbations to become self-gravitating (e.g., Bate et al. 2002, 2003; Delgado-Donate et al. 2004b, 2004a; Goodwin et al. 2004a, 2004b, 2006; Bate & Bonnell 2005). Thereafter, dynamical interactions between the fragments can play an important role in determining the final properties of the binary system (e.g., the number and masses of component stars, together with their orbital separation and eccentricity). Multiple fragments produced in different turbulent cells are predicted to exhibit random orientations between the circumstellar disks or spin axes of the binary components, and no particular relationship between the circumstellar disks and orbital plane, or the spin and orbital axes. If multiple fragments are produced in a common region where turbulence conspires to create local angular momentum, however, the binary system thus assembled can exhibit quite well aligned circumstellar disks and orbital planes, or spin and orbital axes.

Rather than through turbulence, the large-scale ordered rotation of the core can also drive dynamical instabilities to induce fragmentation during collapse. In such models, conservation of angular momentum forces cores to become increasingly flattened as they collapse. As a result, a disequilibrium disk-like (i.e., flattened and rotating) structure forms at the center of the core. The central region of the core can become particularly flattened if magnetic fields are invoked to direct infalling matter onto the mid-plane of the disk-like structure; the resulting structures closely resemble, at least morphologically, rotationally supported disks, and are therefore referred to as pseudodisks (Galli & Shu 1993a, 1993b). By introducing an initial density or velocity perturbation, the large-scale ordered rotation of the core can drive dynamical instabilities in the form of a spiral, bar, or ring in its central flattened region (Matsumoto & Hanawa 2003; Cha & Whitworth 2003; Machida et al. 2008). The pattern of the dynamical instability does not depend on the nature of the initial velocity or density perturbations. Furthermore, because the perturbations introduced can either promote or hinder fragmentation (in the latter case through the effective removal of angular momentum by gravitational torques from the resulting dynamical instability), increasing the amplitude of the initial perturbation does not necessarily increase the likelihood of fragmentation. Fragments form in localized regions of the resulting dynamical instabilities that are gravitationally unstable (according to the Toomre criterion) and have masses exceeding the local Jeans mass. Theoretical considerations suggest that, when driven by the large-scale ordered rotation of the core, fragmentation most likely occurs during the adiabatic phase (when a first hydrostatic core forms and grows) or the protostellar phase (see arguments in Machida et al. 2008).

Binaries that form through the rotation fragmentation of disk-like structures should naturally exhibit a close alignment between the spin axes of the component stars. Because fragments can be produced at different heights from and perhaps even on opposite sides of the mid-plane, however, the circumstellar disks and orbital planes of the resulting binary system need not be closely aligned but can span a  range of angles. The mass ratio and orbital parameters (i.e., orbital separation and eccentricity) of such binary systems depend not only on the initial mass ratios and orbits of their parental fragments, but also on how these fragments interact with each other (particularly in systems initially comprising three or more fragments) as well as with the surrounding material from which they accrete.

As explained above, both turbulent and rotationally driven fragmentation can produce binary systems in which the circumstellar disks are parallel to each other, as well as quite closely parallel to if not coplanar with the binary orbit, satisfying some of the observed physical properties of L1551 IRS 5. The crucial distinction between the two fragmentation models, however, is the physical relationship between the binary system and its surrounding bulk envelope, constituting the remnant parental core. In models that invoke turbulent fragmentation, even in the situation where multiple fragments are produced in a common region where turbulence conspires to create local angular momentum (and thus produce an aligned binary system), the binary system need not orbit in the same direction as the rotation of the bulk envelope, let alone have circumstellar disks aligned with any flattening of the bulk envelope. On the other hand, in models that invoke rotational fragmentation, the binary system will of course orbit in the same direction as the rotation of its surrounding envelope, and its circumstellar disks will be aligned with the surrounding flattened envelope. The observed physical properties of L1551 IRS 5 therefore directly point to its formation through rotational fragmentation. In support of this argument, both the measured rotational velocity and the high degree of flattening exhibited by the envelope around L1551 IRS 5 suggest that its parental core possessed considerable angular momentum.

As mentioned above, even when formed by rotational fragmentation, the circumstellar disks of the resulting binary protostars, although parallel to each other, need not be accurately aligned with the orbital plane. Bate et al. (2000) argue that tidal interactions can quite rapidly align the circumstellar disks of binary protostars with their orbital planes, such that low-mass binary protostellar systems with orbital separations ≲100 au should have coplanar disks and orbital planes. If this process had indeed driven the circumstellar disks of L1551 IRS 5 toward close alignment with the binary orbital plane, the same process would have led to a misalignment between the circumstellar disks/orbital plane and the surrounding flattened envelope. Consider the situation where the two protostars formed at different heights from, and perhaps even on opposite sides of, the mid-plane of a flattened core, as shown in the schematic diagram of Figure 10(a). At birth, the circumstellar disks of the two protostars are therefore closely parallel to their surrounding flattened envelope, but not coplanar with each other or the orbital plane. If the circumstellar disks are then brought into alignment with the orbital plane by tidal interactions, the mutual plane of the binary system now becomes misaligned with the surrounding flattened envelope, as shown in Figure 10(b). In contrast, we find that the circumstellar disks of L1551 IRS 5 are very closely parallel to their surrounding flattened envelope, but probably misaligned by up to typically ∼25° from the orbital plane, a geometry more akin to that shown in Figure 10(a). The probable misalignment between the circumstellar disks/flattened envelope and the binary orbital plane may therefore reflect the formation of the binary fragments in different planes sharing the same rotational axis.

Zoom In Zoom Out Reset image size

By imposing the condition that the inferred break radii of the circumstellar disks cannot be larger than their theoretically predicted tidally truncated radii, only a very narrow range of orbital periods, if any, are permitted for the best-fit orbits with eccentricities e ≳ 0.4, as explained in Section 6.3.1. Thus, most probably, L1551 IRS 5 has a circular or low-eccentricity orbit. Early in the protostellar and perhaps even the first core phase when fragmentation occurs, the fragments have a very low mass compared to their surrounding envelope (parental core). In the situation where, at apastron, the orbital velocity of the fragments is smaller than the rotational velocity of their surrounding envelope (in the immediate vicinity), tidal interactions between the envelope and the fragments transfer both angular momentum and energy from the envelope to the binary. This action results in an increase in the orbital velocity and hence separation during the subsequent passage through periastron, and consequently a reduction in the orbital eccentricity. Similarly, at periastron, if the orbital velocity of the fragments is larger than the rotational velocity of their surrounding envelope (in the immediate vicinity), tidal interactions between the fragments and the envelope transfer both angular momentum and energy from the binary to the envelope. This action results in a decrease in the orbital velocity and hence separation during the subsequent passage through apastron, and once again a reduction in the orbital eccentricity. Thus, the probable low orbital eccentricity of L1551 IRS 5 may have been imprinted at birth.

The total mass of L1551 IRS 5 is at least ∼0.5 M_☉ (based on the best-fit orbits), compared with a mass for its envelope of about 0.1 M_☉ (Momose et al. 1998). Thus, during its remaining protostellar phase, L1551 IRS 5 is not likely to accrete more than ∼20% of its present total mass and perhaps much less. Thus, both the total mass and the mass ratio of the binary system should be largely preserved as it evolves onto the pre-main-sequence and then the main-sequence. Whether or not the orbital parameters of L1551 IRS 5 will evolve significantly during its remaining protostellar phase is not so clear. In the situation where the binary system dominates in mass over its circumbinary disk (so that, even at apastron, the binary orbits faster than the circumbinary disk rotates), as must be the case here, theoretical simulations show that tidal interactions between the binary system and its (coplanar) circumbinary disk result in a transfer of both angular momentum and energy from the binary to its circumbinary disk (Artymowicz et al. 1991; Bate 2000). The effect is to drive up the orbital eccentricity, opposite to the situation described above where the envelope dominates in mass over the binary fragments. By preserving the apastron distance, a(1 + e), the effect also is to decrease the orbital semimajor axis.

Thus, L1551 IRS 5 will likely evolve into a solar-mass binary system that has roughly equal companion masses, quite closely coplanar circumstellar disks and orbital planes (in particular if tidal interactions between one protostar and the circumstellar disk of the other protostar bring their circumstellar disks further into alignment with the orbital plane), and hence presumably also closely aligned stellar spin and orbital axes. Unless its orbital parameters are dramatically altered during its subsequent evolution, L1551 IRS 5 will retain a likely circular or low-eccentricity orbit with a semimajor axis of many tens to perhaps a few hundred au, and a corresponding orbital period of many hundreds to perhaps a few thousand years. The physical properties that we project L1551 IRS 5 will have in its final state closely resemble those observed for binary systems on the pre-main-sequence and main-sequence with similar masses and orbital separations (see, e.g., the review by Duchêne & Kraus 2013). Such systems do not show a preferred mass ratio or orbital eccentricity, but do exhibit preferentially aligned circumstellar disks on the pre-main-sequence and preferentially aligned stellar spin and orbital axes on the main-sequence. For example, the circumstellar disks of T Tauri star binaries with projected orbital separations spanning the range ∼100–1000 au are preferentially aligned to within a few tens of degrees (e.g., Wolf et al. 2001; Jensen et al. 2004; Monin et al. 2006). Main-sequence solar-type binaries with orbital separations smaller than about 30–40 au (close to the median orbital separation) have approximately aligned spin and orbital axes, whereas those with larger orbital separations commonly show misaligned spin and orbital axes (Hale 1994). Many of the systems that are closely aligned may have formed in a manner similar to L1551 IRS 5, rather than having being brought into alignment by internal or external interactions.