ANGULAR MOMENTUM EXCHANGE BY GRAVITATIONAL TORQUES AND INFALL IN THE CIRCUMBINARY DISK OF THE PROTOSTELLAR SYSTEM L1551 NE

Figures 1(a), (b), and (c) show the 0.9 mm continuum images of L1551 NE at the natural weighting, robust = +0.5 weighting, and the uniform weighting using data with projected uv distances larger than 80 kλ. These three images demonstrate the circumbinary structures at progressively higher angular resolutions without missing the primary emission components. In Figures 1(a) and (b), a compact component associated with Source A, with an emission extension to the location of Source B, is seen. This component most likely traces the circumstellar disks around the binary members. This circumstellar-disk component is embedded in an extended emission elongated along the northwest to southeast direction, which most likely traces the circumbinary disk. In the highest-resolution (072 × 036; P.A. = 91) 0.9 mm image shown in Figure 1(c), the central emission component seen in the lower-resolution images is resolved into two compact components, verifying the presence of the circumstellar disks. Gaussian fits to these two components demonstrate that they are barely if at all resolved. The southeastern component has a flux density of 0.35 Jy and a centroid location of (04^h31^m4451, 18°08'314), and the northwestern component has a flux density of 0.17 Jy and a centroid location of (04^h31^m4447, 18°08'316). The positions of these components are slightly (∼03) south of the positions of Sources A and B derived from previous observations at 3.6 cm that trace free–free emission (outflows) from the individual protostars (Reipurth et al. 2002), but consistent within ≲01 with those measured in our recent observations with the JVLA at 7 mm (J. Lim et al. 2014, in preparation). The slight positional shift over the last decade is consistent with the absolute proper motion of the L1551 region found by Jones & Herbig (1979) and Rodríguez et al. (2003) (μ_α = 0012 yr⁻¹, μ_δ = −0023 yr⁻¹). Hereafter, we regard the centroid positions of these components as derived from their Gaussian fits as the positions of the individual protostars. In addition to the circumstellar disks, the highest-resolution image reveals, for the first time, brightenings in the circumbinary disk. The most prominent internal brightenings comprise the U-shaped feature to the south of the binary system as well as protrusions to the north of Source B and northeast of Source A. A depression or gap in the dust emission also can be seen between the northern and southern parts of the circumbinary disk and the circumstellar disk of Source A.

To better separate the circumstellar disks from the circumbinary disk, we subtracted two Gaussian components with parameters derived from the fit to the two compact dust components. The resultant image is shown in Figure 1(d), where the southern U-shaped feature remains evident. The protrusions to the north of Sources A and B can now be seen to comprise another U-shaped feature, but now located to the north of the binary system. As we will show below, the northern and southern U-shaped features correspond to spiral arms (Arms A and B, respectively) produced by the gravitational torques from the binary stars.

We estimate the masses of the individual circumstellar disks and circumbinary disk (≡ M_d) from their individual continuum fluxes (≡ S_ν) using the relationship

$$\begin{equation} M_{d}=\frac{S_{\nu }d^2}{\kappa _{\nu } B_{\nu }(T_d)}, \end{equation} \tag{ 1 }$$

where ν is the frequency, d the distance, B_ν(T_d) the Planck function for dust at a temperature T_d, and κ_ν the dust opacity per unit gas + dust mass on the assumption of a gas-to-dust mass ratio of 100. We adopt the relation κ_ν = $\kappa _{\nu _{0}}$(ν/ν₀)^β, where β denotes the dust-opacity index and κ_250 μm = 0.1 cm² g⁻¹ (Hildebrand 1983), as is widely adopted in the literature for circumstellar disks around low-mass protostars. Recent multi-frequency observations of circumstellar disks in dust emission show that β is in the range of ∼0–1 (Guilloteau et al. 2011; Chiang et al. 2012), and so here we adopt β = 0.5. The dust mass opacity at 0.9 mm is then calculated to be κ_0.9 mm = 0.053 cm² g⁻¹. This value is a factor of three higher than that of Ossenkopf & Henning (1994) for grains with thin ice mantles coagulated at a density of 10⁶ cm⁻³ (κ_0.9 mm = 0.018 cm² g⁻¹). Thus, adopting the dust mass opacity by Ossenkopf & Henning (1994) provides a factor of three higher masses of the circumstellar and circumbinary disks. Regarding the dust temperature, modeling of the spectral energy distribution (SED) toward L1551 NE as measured in single-dish observations at wavelengths from 12 μm to 2 mm derived T_d ∼ 42 K (Barsony & Chandler 1993; Moriarty-Schieven et al. 1994). This temperature is, however, likely to be biased by the optically thick mid-infrared emission (∼12–100 μm), and thus the temperature is more likely the surface temperature of the disk. The dust temperature in the disk midplane of Class I sources can be as low as ∼10 K (Nakazato et al. 2003; Whitney et al. 2003). With T_d = 10–42 K and κ_0.9 mm = 0.053 cm² g⁻¹, the mass for the circumstellar disk of Source A is calculated to be ∼0.005–0.044 M_☉, that of Source B ∼0.003–0.022 M_☉. The flux density of the circumbinary disk in the 0.9 mm continuum as computed from Figure 1(d) is 0.47 Jy, corresponding to a mass for the circumbinary disk of ∼0.007–0.059 M_☉. The individual as well as combined disk masses are therefore much smaller than the inferred total binary mass of 0.8 M_☉. In our theoretical model for the distribution and dynamics of matter around L1551 NE, we therefore ignore the self-gravity of the circumstellar disks and circumbinary disk, and only consider the gravitational field of the binary protostars.

3.2.1. Integrated Intensity Map

Figure 2 shows the C¹⁸O integrated intensity map at an angular resolution of 082 × 049 (P.A. = 159) (contours) superposed on the 0.9 mm dust-continuum map (colors; same as that in Figure 1(c)) of L1551 NE. The C¹⁸O map in Figure 2(a) includes an unknown contribution from dust, whereas that in Figure 2(b) has a constant level subtracted from every channel corresponding to the dust continuum emission as measured outside the velocity range spanned by the C¹⁸O line. The C¹⁸O map before continuum subtraction shows emission from the circumstellar disks, whereas that after subtracting a constant continuum component shows little if any such emission. In the circumstellar disks, both the dust and gas opacities are appreciable, resulting in partial absorption of the line emission by the dust and partial absorption of the continuum emission by the gas. In this situation, it is not straightforward to separate the line and continuum emission. The same complexity is seen, for example, in the circumstellar disk around the T-Tauri star HD 142527 (Fukagawa et al. 2013). Beyond the circumstellar disks where the dust opacity is significantly lower, the intensity distribution in C¹⁸O is similar irrespective of whether a constant continuum component is subtracted or not. The circumbinary disk can be traced out to a radius of ∼300 AU in C¹⁸O, farther out than in the dust continuum. By contrast with the dust-continuum map, the C¹⁸O map shows no depression between the circumbinary disk and the circumstellar disk of Source A, but instead increases inward in intensity toward both circumstellar disks. The observed behavior suggests that C¹⁸O is more optically thick than dust and that the region between the circumbinary disk and the circumstellar disk of Source A is not devoid of matter.

Zoom In Zoom Out Reset image size

3.2.2. Velocity Structure

Figure 3 shows the continuum-subtracted velocity channel maps of the C¹⁸O (3–2) line (contours) superposed on the 0.9 mm dust-continuum image (gray scale). In our previous observations of L1551 NE with the SMA, we established that C¹⁸O has a systemic velocity of V_LSR = 6.9 km s⁻¹, and detected this line over the velocity range ∼4.3 km s⁻¹ to ∼9.1 km s⁻¹ (Takakuwa et al. 2012; Takakuwa et al. 2013). In our observation here at a higher sensitivity using ALMA, C¹⁸O can be detected to higher blueshifted velocities of 2.8–4.3 km s⁻¹ and redshifted velocities of 9.0–10.3 km s⁻¹. As can be seen in Figure 3, the C¹⁸O emission peaks between Sources A and B at blueshifted velocities of 2.8–3.5 km s⁻¹, and to the southeast of Source A at redshifted velocities of 9.9–10.3 km s⁻¹. As we will explain later, these high-velocity components likely trace gas in the Roche lobes of the two protostars. At blueshifted velocities of 3.7–4.8 km s⁻¹ the C¹⁸O emission extends northeast from the midpoint between Sources A and B, and at redshifted velocities of 8.3–9.7 km s⁻¹ to an area spanning from south to southwest of Source A. At velocities of 5.0–6.1 km s⁻¹ and 7.4–8.1 km s⁻¹, the C¹⁸O emission originates predominantly from the northern and southern parts, respectively, of the circumbinary disk. These low-velocity blueshifted and redshifted emissions trace the outer northern and southern parts of the circumbinary disk seen in the 0.9 mm continuum emission, but extend farther out than in the continuum. Close to the systemic velocity (6.3–7.2 km s⁻¹), the C¹⁸O emission exhibits a "butterfly" pattern (particularly at 6.8 and 7.0 km s⁻¹), a characteristic signature of a Keplerian disk (Simon et al. 2000). The channel maps at velocities relatively close to the systemic velocity (5.0–8.1 km s⁻¹) trace the same regions as previously seen in our SMA observations, where the spatial-kinematic structure of the circumbinary gas out to a radius of ∼300 AU was satisfactorily modeled as a Keplerian disk. On the other hand, as we shall show, the emission extending to higher blueshifted and redshifted velocities detected here for the first time traces deviations from Keplerian rotation.

Zoom In Zoom Out Reset image size

To help in visualizing the spatial distribution of the newly detected high-velocity components, in Figure 4 (upper panels) we show maps of the high-velocity (≳2 km s⁻¹) blueshifted (blue contours) and redshifted (red contours) C¹⁸O (3–2) emission superposed on the 0.9 mm dust-continuum image (gray scale). These high-velocity components also are detected in ¹³CO (3–2), CS (7–6), and SO (7₈–6₇) (no emission was detected in HC¹⁸O⁺). For example, in Figure 4 (lower panels), we show the corresponding channel maps in ¹³CO (see the full velocity channel maps of the ¹³CO emission in Appendix A). At the highest blueshifted velocities (left panels), both the C¹⁸O and ¹³CO emissions peak between Sources A and B. On the other hand, at the highest redshifted velocities (left panels), both the C¹⁸O and ¹³CO emissions peak to the south-east of Source A. At slightly lower velocities (middle panels), the blueshifted C¹⁸O and ¹³CO emissions peak to the north from the midpoint between Sources A and B and curl northeast, whereas the redshifted C¹⁸O and ¹³CO emissions peak to the southeast of Source A and curl to the southwest. The trajectory of the emission peaks at blueshifted and redshifted velocities follows an S-shaped curve as highlighted by the yellow curves. At even lower velocities (right panels), the emission in both C¹⁸O and ¹³CO now encompasses the circumbinary disk as traced in the continuum.

Zoom In Zoom Out Reset image size

Figure 5(a) shows the intensity-weighted mean velocity map in C¹⁸O. Whereas the velocity field measured in our previous SMA observations at a lower angular resolution and sensitivity can be satisfactorily modeled as pure Keplerian rotation (Figure 5(b)), the velocity field measured in our ALMA observation shows clear deviations from Keplerian rotation. If the rotation was purely circular, the line-of-sight (LOS) velocity would be symmetric about the major axis of the disk (tilted vertical dashed lines in Figure 5) and hence the ridge of peak LOS velocities aligned with the major axis. Even in the case of aligned Keplerian elliptical orbits, the ridge of peak LOS velocities should lie along the straight line for any projection in the sky. The ALMA mean-velocity map shows, however, that this ridge deviates from the major axis in a systematic way that cannot be explained as either circular or aligned elliptical orbits in a Keplerian disk. In addition, a velocity gradient is apparent along the minor axis of the disk, such that the eastern side of the disk is redshifted and the western side blueshifted.

Zoom In Zoom Out Reset image size

In Figure 5(c), we show the residual mean-velocity map after subtracting the best-fit Keplerian model derived from our SMA observations (i.e., Figure 5(b)) from the mean velocity map observed with ALMA (Figure 5(a)). Details of the subtraction process are described in Appendix B. For comparison, we overlay the ALMA 0.9 mm dust-continuum image after subtracting the circumstellar disks so as to highlight features in the circumbinary disk (white contours). The subtracted mean-velocity map shows that, away from the major axis of the circumbinary disk, the LOS velocities are, in general (except relatively close to the protostellar system), more blueshifted to the north and more redshifted to the south than expected if the rotation was purely Keplerian (and circular). Along the purple and red curves drawn in Figure 5(c), coinciding with a portion of the individual U-shaped features seen in the continuum, the LOS velocities are ≳0.5 km s⁻¹ more blueshifted and redshifted, respectively, than expected for pure (circular) Keplerian rotation. Just beyond the protostellar components, the subtracted mean-velocity map shows a patch of residual redshifted velocities to the north and blueshifted velocities to the south. These patches indicate much lower LOS velocities than Keplerian in the close vicinity of the binary system, implying significant changes in orbital velocities close to and within the Roche lobes (≲10; see Section 4.1).

Figure 6 (Left) shows the observed position–velocity (P–V) diagrams of the C¹⁸O emission along the major and minor axes of the circumbinary disk (N–S and E–W in Figure 4(a)), as well as along the northern and southern sides of the circumbinary disk parallel to the minor axis (EN–WN and ES–WS in Figure 4(a)). Along the major axis of the circumbinary disk (N–S), the observed P–V diagram shows primarily Keplerian rotation (green curves from the Keplerian model derived in our SMA observations). Along the minor axis passing through the disk center (E–W), however, the P–V diagram is not symmetric as would be expected for Keplerian rotation. Instead, the P–V diagram along this axis shows blueshifted emission extending to the highest velocities (<−2 km s⁻¹ from the systemic velocity) in the west and redshifted emission extending to the highest velocities (>+2 km s⁻¹ from the systemic velocity) in the east. Although the P–V diagram therefore indicates a blueshift to redshift velocity gradient from west to east along the minor axis passing through the disk center, the P–V diagram parallel to the minor axis but to the south of the disk center (ES–WS) shows that the emission increases in redshifted velocity to the west (red dashed line) rather than to the east. In addition, the highest redshifted velocity is located to the west (arrow) rather than the east. In the case of circular Keplerian rotation, the P–V diagrams along the on- and off-center minor axes are symmetric and do not show any velocity gradients (Takakuwa et al. 2013). In the case of a rotating disk exhibiting a radial velocity component due to infall, the P–V diagram exhibits a velocity gradient along the minor axis reflecting the infalling motion (Yen et al. 2010; Takakuwa et al. 2013), but the sense of the velocity gradient is the same for all P–V diagrams parallel to the minor axis. Thus, the different signs of the velocity gradients along and parallel to the minor axis cannot be simply explained by combination of rotation and infall.

Zoom In Zoom Out Reset image size

Thanks to the higher angular resolution and sensitivity afforded by ALMA, we have been able to detect internal structures in the circumbinary dust disk and deviations from Keplerian motion in the circumbinary gas disk of L1551 NE not previously seen in our observations with the SMA. In the following, we discuss the physical causes of these newly detected features. As will become apparent, the explanation for these features provides insights into the mechanisms that drive the exchange of angular momentum in the circumbinary disk of L1551 NE so as to permit infall through this disk.

In the Roche potential of binary systems, circumbinary orbits close to the binary system experience strong gravitational torques and cannot be described by circular or elliptical orbits as those around single stars. Instead, torques that impart angular momentum to the circumbinary disk create shocks where matter that gained angular momentum collides supersonically with matter located downstream having less angular momentum. Theoretical simulations show that these shocks form a two-arm spiral pattern that co-rotates with the binary system. One arm extends from the primary through the L3 Lagrangian point in the Roche potential, and the other arm from the secondary through the L2 point (see Figure 9(a)). These Lagrangian points correspond to locations where circumbinary material flows into the Roche lobes of the respective binary components (Artymowicz & Lubow 1996); note, however, that infall does not occur along the spiral arms, but rather in regions between the spiral arms as we will explain later. Although the degree of density enhancement and the opening angle of the spiral arms depend on the binary parameters and physical conditions (density and temperature) in the circumbinary disk, the two-arm spiral density pattern is a generic prediction of theoretical models of protobinary systems surrounded by a circumbinary disk (Artymowicz & Lubow 1996; Bate & Bonnell 1997; Bate 2000; Günther & Kley 2002; Ochi et al. 2005; Hanawa et al. 2010).

To model the internal structure and dynamics of the circumbinary disk in the specific case of L1551 NE, we performed a hydrodynamic simulation using the adaptive mesh refinement (AMR) code SFUMATO (Matsumoto 2007). The relevant parameters adopted in our simulation are summarized in Table 2. We assumed that the circumbinary disk (indeed, all the material beyond the two protostars) is isothermal at a temperature of 13.5 K, and that its equatorial plane is aligned with the orbital plane of the binary system. Material from the circumbinary disk therefore accretes onto the individual circumstellar disks of the binary protostars along the same plane. We assumed that the binary system has a circular orbit, and adopted a total mass for the binary of 0.8 M_☉, a mass ratio of 0.19, and a binary orbital separation 145 AU as determined from our previous SMA observations (Table 2). The sizes of the Roche lobes of the primary and secondary are thus ∼150 AU and ∼70 AU, respectively (see Figure 9(a)). Gas is injected at the cylindrical computational boundary located at 1740 AU from the center of the mass of the binary system. The specific angular momentum of the injected gas (≡ j_inj) is chosen to have a centrifugal radius of 300 AU (i.e., j_inj = 461.5 AU km s⁻¹). This radius corresponds to the outer radius of the circumbinary disk, where the gas motion changes from infall in the surrounding envelope to the Keplerian rotation in the circumbinary disk as found in our previous SMA observations (Takakuwa et al. 2013). The choice of specific angular momentum is guided also by the need to reproduce the inner depression in dust emission in the circumbinary disk. As previously shown by Bate & Bonnell (1997), the size of the inner depression in the circumbinary disk is dependent on the specific angular momentum of the injected gas. Our previous SMA observations show that the mass-infall rate from the surrounding envelope to the circumbinary disk is of the order of ∼10⁻⁶ M_☉ yr⁻¹ (Takakuwa et al. 2013). In our model, we assume a density for the injected gas of 1.5 × 10⁵ cm⁻³ (Table 2), which gives a comparable mass-infall rate from the envelope to the circumbinary disk of 3.3 × 10⁻⁶ M_☉ yr⁻¹. We adopted the simulation result after the ∼20 orbital periods of the binary (∼4 × 10⁴ yr). By that stage of the simulation, the gas distribution and kinematics have settled into a quasi-steady state and the circumbinary disk developed spiral arms, but the accreted masses onto the protostellar binary have not yet changed the individual protostellar masses significantly.

Having computed the gas and also dust (by simply assuming a gas-to-dust mass ratio of 100) distribution from our hydrodynamic simulation, we then performed radiative transfer calculations to produce the theoretically predicted 0.9 mm dust-continuum and C¹⁸O images. To better capture the intensity distribution of both the gas and dust, when computing the radiation transfer we assumed a combination of two power-law temperature distributions centered, respectively, on Sources A and B. The parameters of the individual temperature profiles are listed in Table 2. We first fixed the temperature at a radius of 300 AU, the outer radius of the circumbinary disk, to approximately reproduce the observed 0.9 mm continuum and C¹⁸O line intensity peaks in the circumbinary disk. Then, the power-law index of the temperature profile was determined to closely reproduce the observed flux density of the compact continuum component (circumstellar disk) associated with each protostar. The choice of the power-law index is controlled by the circumstellar-disk components, and we confirmed that changing the power-law index within the reasonable value (0–0.5) does not affect the appearance of the spiral-arm features in the circumbinary disk much, and thus the main points of the present paper. We also note that the use of power-law temperature profiles is in conflict with our model computations which assume that the matter around both protostars is isothermal. However, the adopted power-law index results in a change in the temperature of only a factor of 1.4 through the circumbinary disk (from 33 K at r ∼ 50 AU to 23 K at r ∼ 300 AU). The peak brightness temperatures of the ¹³CO and C¹⁸O image cubes on the northern blueshifted side of the circumbinary disk are 28.9 K and 19.8 K, and those in the southern redshifted side 25.9 K and 24.7 K, respectively. On the assumption of local thermodynamic equilibrium (LTE) and X (¹³CO)/X (C¹⁸O) = 7.7 (Wilson & Rood 1994), the excitation temperatures and optical depths can be estimated from these isotopic lines. We find, for C¹⁸O, an excitation temperature of ∼36 K and optical depth of ∼1.2 at the peak brightness on the northern side, and an excitation temperature of ∼33 K and optical depth of ∼3.0 at the peak brightness temperature on the southern side. Thus, the C¹⁸O emission is close to or completely optically thick at the emission peaks, and the derived excitation temperatures at these locations correspond to the gas temperatures at the surface of the circumbinary disk. More typically, the peak brightness temperature of C¹⁸O in the channel maps ranges from ∼15 K to ∼20 K, suggesting that the C¹⁸O emission is mostly optically thin. The adopted temperature of 13.5 K in our hydrodynamic simulations is lower than the typical peak brightness temperature measured in C¹⁸O, but may be more representative of the bulk of both the gas and dust in the circumbinary disk that is concentrated close to the midplane as discussed in Section 3.1. Thus, the different temperatures used in our hydrodynamic simulation and in the radiative transfer calculations may not present as severe an internal inconsistency as might appear on first sight. As a check, we confirmed that using a temperature of 25 K does not produce any significant differences in the gas distributions in our hydrodynamic simulations. Following the radiative transfer calculations, we used the CASA task "simobserve" to create the simulated visibility data for the model images with the same antenna configuration, hour angle coverage, bandwidth and frequency resolution, integration time, and noise level as those of the real ALMA observation. We further performed flagging of the simulated data to match the simulated data with the real processed data, and then made simulated theoretical images with the same imaging methods as described in Section 2. More details of our theoretical simulations are described in Appendix C.

Figure 7(a) shows the model continuum image (i.e., before passing through the ALMA simulator). This image clearly shows two spiral arms, along with a depression or gap between the spiral arms and the two circumstellar disks. The northern arm labeled A connects to the circumstellar disk around Source A, and that labeled B to the circumstellar disk around Source B. In Figures 7(b) and (c), the simulated theoretical continuum image (i.e., after the model image has been passed through the ALMA simulator) and the observed ALMA 0.9 mm continuum image (same as that in Figure 1(c)) are shown, respectively. The observed U-shaped features to the north and south of the binary system coincide nicely with the spiral arms. Furthermore, the simulated theoretical continuum image closely reproduces the inner depression observed in the ALMA continuum image.

Zoom In Zoom Out Reset image size

Figure 8 shows the simulated theoretical velocity channel maps of the C¹⁸O (3–2) line (contours) superposed on the simulated theoretical 0.9 mm dust-continuum image (gray scale). These channel maps closely capture all the important features seen in the ALMA image cube as described above. At the highest blueshifted velocities of 2.8–4.3 km s⁻¹, the C¹⁸O emission is located between Sources A and B. At lower blueshifted velocities of 4.6 to 5.2 km s⁻¹, the C¹⁸O emission extends toward the north before curling northeast from the midpoint between Sources A and B. By comparison, at the highest redshifted velocities of 9.9–10.3 km s⁻¹, the C¹⁸O emission is located to the southeast of Source A. At slightly lower redshifted velocities of 8.6–9.7 km s⁻¹, the C¹⁸O emission is located to the southeast of Source A and curls to the southwest. At lower blueshifted and redshifted velocities closer to the systemic velocity, the channel maps primarily exhibit Keplerian rotation. The P–V diagrams constructed from the simulated theoretical images are shown in Figure 6 (right), and closely reproduce the velocity features in the observed P–V diagrams as described above, i.e., the velocity gradient along the minor axis, and the different signs of the velocity gradient along the minor axis compared to the velocity gradients along cuts parallel to the minor axis on either side of the circumbinary disk.

Zoom In Zoom Out Reset image size

The simulated theoretical C¹⁸O mean-velocity map shown in Figure 5(d) shows a systematic deviation in the peak LOS velocity away from the major axis of the circumbinary disk. The ridge in peak LOS velocities traces an S-shaped pattern as highlighted by the yellow line in Figure 5(d) that is seen also in the velocity field of our ALMA C¹⁸O map (Figure 5(a)). This deviation results in higher blueshifted LOS velocities at the northeastern (purple curve in Figure 5(c)) and higher redshifted LOS velocities at the southwestern (red curve) portions of the circumbinary disk than predicted for Keplerian rotation. These portions of the circumbinary disk correspond, respectively, to the eastern part of Arm A and the western part of Arm B. Because the spiral arms are produced by gravitational torques from the binary protostars that impart angular momentum onto matter in the circumbinary disk, the matter along the arms orbits faster than Keplerian and expands outward, thus producing a peak in the LOS velocity (see Figure 9). The east (blueshifted) to west (redshifted) velocity gradients along the offset minor axes seen in the observed and theoretical P–V diagrams (Figure 6) reflect the expanding motion of the arms, since the disk near side is on the eastern side and the far side is on the western side. On the other hand, the minor axis of the circumbinary disk does not cross the spiral arms, but instead crosses regions where our theoretical model predicts torques largely extract angular momentum from the circumbinary disk (Figure 9). Here, our theoretical model predicts infall occurs, thus explaining the observed velocity gradient along the minor axis of the circumbinary disk. As the observed velocities on both sides of the center along the minor axis reach ±0.5 km s⁻¹, the infall velocity corresponds to ∼0.5 km s⁻¹/sin i ∼ 0.6 km s⁻¹ (where i is the inclination of the disk). The inferred infall velocity closely matches that predicted in our model.

Zoom In Zoom Out Reset image size

In summary, both the U-shaped dust brightenings in the circumbinary disk of L1551 NE and its velocity field as measured in C¹⁸O can be explained solely by the effects of gravitational torques from the binary components. The gravitational torques, rather than magnetic fields or other mechanisms for generating turbulence, constitute the primary driver for exchanging angular momentum in the circumbinary disk so as to permit infall through this disk and hence growth of L1551 NE. In Figure 10, we present a schematic picture that shows the major features so far discovered in our study of L1551 NE.

Zoom In Zoom Out Reset image size

An important and as yet unresolved question in binary star formation is once one protostar gains a significantly larger mass than the other, which protostar accretes at a higher rate and what factors determine the mass ratio of the final binary product? At approximately solar masses, observations suggest that binary stars do not exhibit a preferred mass ratio. Duquennoy & Mayor (1991) showed that the mean mass ratio of the solar-type stars is ∼0.4, and Reid & Gizis (1997) found that the mass ratio distribution is approximately flat from ∼0.1 to 1. Theoretical simulations have produced conflicting results on whether the primary or secondary protostar accretes at a higher rate, and thus whether the mass ratio of binary systems is driven away from or toward unity. Bate (1997), Bate & Bonnell (1997), and Bate (2000) suggest that the secondary protostar accretes at a higher rate because it is located further away from the center of the mass of the binary system, and thus sweeps up more material in its orbit than the primary protostar. On the other hand, numerical simulations by Ochi et al. (2005) and Hanawa et al. (2010) employing higher resolutions show that the accretion rate of the primary is larger than that of the secondary. They argue that high-resolution simulations are necessary to properly follow the complex gas motions around the L2 and L3 Lagrangian points where material from the circumbinary disk flows into the Roche lobes. Similar results are also obtained by Fateeva et al. (2011).

We anticipate that observations that can properly resolve the gas kinematics around the L2 and L3 Lagrangian points and within the Roche lobes will be crucial for addressing the factors that determine the relative accretion rates onto the individual components of a binary system. Future observations with ALMA will be able to resolve the gas kinematics around the Lagrangian points and within the Roche lobes, and to provide important observational insights on the mass accretion onto the primary and secondary of the protostellar binaries.

To help interpret the two U-shaped brightenings and velocity field of the circumbinary disk in L1551 NE as observed with ALMA, we performed a hydrodynamic simulation of an accreting binary system tailored to the specific properties of L1551 NE. In our model, we assume that the binary stars have a circular orbit. The masses and separation of the protostars are as inferred from our previous observations of L1551 NE with the SMA, corresponding to a mass of 0.675 M_☉ for Source A, mass of 0.125 M_☉ for Source B, and a separation of D = 145 AU. Using cylindrical coordinates (R, φ, z), we placed the protostars at the midplane (z = 0 plane) such that their center of mass coincides with the origin (R = 0). The computational domain is defined as 0 ⩽ R ⩽ R_bound, 0 ⩽ φ ⩽ 2π, and −Z_bound/2 ⩽ z ⩽ Z_bound/2, where R_bound = Z_bound = 12D = 1740 AU. Gas having a constant density is injected at the cylindrical boundary at R = R_bound with a velocity

$$\begin{equation} \left(\begin{array}{@{} c @{}}v_R \\ v_\varphi \\ v_z \end{array} \right) = \left(\begin{array}{@{} c @{}}\displaystyle - \sqrt{\frac{2GM_\star }{R_\mathrm{bound}} - \left(\frac{j_\mathrm{inj}}{R_\mathrm{bound}} \right)^2 } \\ \displaystyle \frac{j_\mathrm{inj}}{R_\mathrm{bound}}\\ 0 \end{array} \right). \end{equation} \tag{ C1 }$$

The radial velocity corresponds to the infall velocity for material experiencing free-fall from infinity in the mid-plane. The specific angular momentum of the injected gas is given by

$$\begin{equation} j_\mathrm{inj} = \sqrt{ GM_\star R_\mathrm{cent} }, \end{equation} \tag{ C2 }$$

where R_cent (=300 AU) is the centrifugal radius of the gas as determined from our observations with the SMA (Takakuwa et al. 2012; Takakuwa et al. 2013). The mass-injection rate into the computational domain depends on the density of the injected gas, and is given by

$$\begin{eqnarray} \dot{M} &=& {-} 2 \pi R_\mathrm{bound} Z_\mathrm{bound} \rho _\mathrm{inj} v_R \nonumber\\ &=& 2.22 \times 10^{-6} \left(\frac{n_\mathrm{inj}}{10^5 \,\mathrm{cm}^{-3}} \right) \, M_\odot\, \mathrm{yr}^{-1}, \end{eqnarray} \tag{ C3 }$$

where ρ_inj and n_inj denote the mass density and the number density, respectively, of the injected gas. For simplicity, we assume that the gas is isothermal throughout the simulation with a corresponding sound speed of $c_s = 0.1 \sqrt{{GM}_\star /D} = 0.221$ km s⁻¹. The gas temperature is therefore 13.5 K, a value that is typical for dense molecular cloud cores and infalling envelopes. The self-gravity of the gas can be safely ignored because, based on both our SMA and ALMA observations, the mass of the circumbinary disk is considerably smaller than the total mass of the binary stars (see Section 3.1). The gas is therefore attracted inward solely by the gravity of the protostars.

The calculation was performed using SFUMATO, a three-dimensional AMR code (Matsumoto 2007). The hydrodynamic scheme was modified to have a third order of accuracy in space and second order in time. The computational domain of (2R_bound)² × Z_bound is resolved by the base grid (l = 0) with 128² × 64 cells. In the AMR hierarchy, the maximum grid level is set at l = 5. The cell width is Δx = 0.85 AU on the finest grid of l = 5, compared with Δx = 27 AU on the base grid of l = 0. The effective resolution therefore corresponds to 4096² × 2048. The initial hierarchical grid was fixed during the calculation, and no further grid refinement during the time integration is performed. The two protostars are represented by two point masses that comprise sink particles. The sink particles accrete gas within a given radius referred to as the sink radius. The sink radius is set at 3.4 AU, which is considerably smaller than the binary separation of D = 145 AU, and we adopted the accretion method described by Krumholz et al. (2004). The simulation is performed in a rotating frame where the sink particles are at rest, and takes into account both centrifugal and Coriolis forces.

In the simulation, the injected gas falls toward the center, and initially forms a circumbinary disk having an inner radius of ∼200 AU and an outer radius of ∼400 AU. At this early stage, the circumbinary disk, which is surrounded by an infalling envelope, therefore exhibits a central gap. Matter at the inner edge of the circumbinary disk is prevented from falling into the central gap by its angular momentum. After ∼20 orbital periods, by which time the gas has settled into a quasi-steady state, the circumbinary disk has developed spiral arms as shown in Figure 9(a). The spiral arms in the circumbinary disk are excited by the orbiting binary protostars; the rotation of their non-axisymmetric potential generates torques that act on the gas in the circumbinary disk. This non-axisymmetric potential is expressed as the Roche potential drawn in Figure 9(a). The spiral arms co-rotate with the protostellar binary, and sweep through the circumbinary disk and generate the shock waves, as their angular velocity is higher than that of the circumbinary disk (see Figure 9(b)). The shock wave pushes gas downstream in the circumbinary disk. In other words, the gas in the spiral arms obtains angular momentum which is transferred from the orbital angular momentum of the protostellar binary via the non-axisymmetric potential and shock waves. The higher orbital velocity in the spiral arms causes the positive radial velocity there as shown in Figure 9(c). Where the angular momentum is lowered by the non-axisymmetric potential, as indicated by regions having negative radial velocities in the circumbinary disk shown in Figure 9(c), infall occurs. The positive and negative radial velocities coincide with the downstream and upstream sides, respectively, of the shock waves.

Infall through the circumbinary disk promotes accretion onto the circumstellar disks through the gap in the circumbinary disk. Gas falls into this gap from the upstream side of the spiral arms, and then onto the circumstellar disks orbiting inside the gap. In our simulation, we find that the mass-accretion rate onto the secondary star (Source B) is about one order of magnitude higher than that onto the primary star (Source A). Thus, in our simulation, the binary mass ratio increases as the protostellar binary evolves. Similar results have also been reported in previous works (e.g., Bate & Bonnell 1997).

The simulation shows that the accretion rate onto Source B is roughly equal to that onto the circumbinary disk in the quasi-steady state. In our simulation we assumed the number density of the injected gas as 1.5 × 10⁵ cm⁻³ (see Table 2 and Equation (C3)). This assumed number density provides an accretion rate from the infalling envelope to the circumbinary disk of 3.3 × 10⁻⁶M_☉ yr⁻¹, similar to that estimated from our previous SMA observations (Takakuwa et al. 2013). By using this number density, the radiative transfer calculation reproduces the dust emission as shown in Figure 7(a). This indicates that the mass accretion rates onto Sources A and B are ∼10⁻⁷M_☉ yr⁻¹ and ∼10⁻⁶M_☉ yr⁻¹, respectively.

From the model distribution of matter in the circumstellar and circumbinary disks, we performed radiative transfer calculation so as to permit a direct comparison between the theoretical images and those obtained in our ALMA observation. To derive the emission distributions of the model images, we simply integrated each specific intensity in each cell along the light of sight using the analytic formula of the radiative transfer and assuming the local thermodynamic equilibrium (LTE) condition (i.e., T_ex = T_k). For the purpose of producing more realistic theoretical images, we assumed a power-law temperature distribution of T∝r^−0.2, whereas for simplicity we assumed the gas to be isothermal in the hydrodynamic simulation (Detailed discussion on the assumed temperature is given in Section 4.1.). The assumed dust opacity and C¹⁸O abundance are listed in Table 2. Finally, the theoretical images are passed through the CASA task "simobserve" to create simulated visibility data, and the simulated visibility data are Fourier-transformed and de-convolved to create the simulated images, shown in Figures 5–8. In this step, we replicated the same antenna configuration, observing time and hour angle convergence, data flagging, thermal noise, and imaging parameters as in the actual observation and data processing.