Structure of a Protobinary System: An Asymmetric Circumbinary Disk and Spiral Arms

A fiducial model with (q, j_inf, c) = (0.2, 1.2, 0.1) is examined here because the model parameters are similar to those for the Class I binary L1551 NE (Takakuwa et al. 2014, 2017). The gas is injected at the boundary surfaces and falls toward the stars. At ∼5 revolutions of the binary, a circumbinary disk forms with a radius that is approximately equal to the centrifugal radius of the infalling gas, R_cent = 1.44. Since the primary and secondary stars each have circumstellar disks, they are hereafter referred to as the circumprimary and circumsecondary disks.

Figures 1 and 2 show the volume density and surface density distributions, respectively, at ∼100 revolutions of the binary. At this stage, the circumbinary disk extends up to a radius of R ∼ 6. The circumbinary disk has a gap with a radius of R ∼ 1.5, in which both the stars rotate. As shown in Figure 1, the circumstellar disks are geometrically thin, while the circumbinary disk exhibits a thick structure and extends up to z ∼ ±3 in the z-direction.

Zoom In Zoom Out Reset image size

The circumstellar disks have clear outer edges, the sizes of which are consistent with Pichardo et al. (2005), who estimated the disk sizes by exploring test particle orbits. In contrast, the gap radius (R ∼ 1.5) is smaller than their estimation of R ∼ 1.8. The gas inflow along the spiral arms may have reduced the gap.

Two spiral arms exist in the gap of the circumbinary disk, as shown in Figure 2. Hereafter, we call the spiral arms on the sides of the primary and secondary stars the primary and secondary arms, respectively. The secondary arm is connected to the circumsecondary disk, while the primary arm runs around the circumprimary disk and is connected to the bridge between the circumstellar disks. The spiral arms are caused by the gravitational torque of the stars. The bridge is a shock wave that is caused by the convergence of the gas flows inside the gap. These spiral arms and bridge were commonly seen in previous simulations (e.g., Bate & Bonnell 1997). While Kaigorodov et al. (2010) and Fateeva et al. (2011) have argued that two bow shocks associated with the circumstellar disks play an important role in opening the gap of the circumbinary disk, no bow shock is associated with the circumsecondary disk (Figure 2). Moreover, the gas flow along the primary arm suggests that the primary arm is a stream rather than a bow shock (see also the streamlines in Figure 4).

The mass and radius of the circumbinary disk increase as time proceeds (Figure 3). The circumsecondary disk reaches a constant mass at ∼30 revolutions of the binary, while the circumprimary disk continues to show an increasing trend at 100 revolutions. At the last stage, the ratio of the masses of the circumbinary disk, circumprimary disk, and circumsecondary disk is M_CBD : M_CPD : M_CSD = 1.00 : 0.043 : 0.045.

Zoom In Zoom Out Reset image size

Of the gas injected at the boundary surfaces, $\sim 40 \%$ is eventually accreted onto either of the sink particles (see dashed lines in Figure 3), and the remaining ∼60% accumulates in the circumbinary disk, which grows in mass and radius. Because the growth of the circumbinary disk is slow compared to the rotation timescale, ${\dot{M}}_{\mathrm{CBD}}/{M}_{\mathrm{CBD}}\sim {10}^{-3}$ at ∼100 revolutions, the structure of the circumbinary disk exhibits a quasi-steady state during the simulation. In the measurement of the masses of the disks, each circumstellar disk is defined by gas with a density higher than 10⁵ that is located inside each Roche lobe, and the circumbinary disk is defined by gas with a density higher than 10 that is located outside the two Roche lobes (see Figure 1).

Note that the accreted masses predominate over the masses of the circumstellar disk for both the primary and secondary in the long-term simulation here. This suggests that models without sink particles are suitable for short-term simulations (e.g., Ochi et al. 2005; Hanawa et al. 2010).

The circumbinary structures exhibit small-scale fluctuations. In order to erase such fluctuations, the temporal average is taken over the period of 80–100 revolutions, as shown in Figure 4. The purpose of the temporal average is to investigate the dynamics of the gas, not for comparison with observations. With the temporally averaged distributions, we can present long-lived, essential physical phenomena in the circumbinary disk. While the temporally averaged surface density exhibits a smooth distribution compared to the snapshot, the spiral arms are clearly seen. The bridge between the circumstellar disks is blurred in the temporally averaged surface density because the bridge changes as it wavers in the vicinity of the L1 point.

Zoom In Zoom Out Reset image size

The streamlines in Figure 4 show the temporally averaged density-weighted velocity distribution. The gas falls along the spiral arms from the circumbinary disk to the circumstellar disks. The gas stream along the secondary arms falls directly onto the circumsecondary disk. The gas stream along the primary arm is connected to the bridge between the two circumstellar disks. These streams along the spiral arms cause the accretion from the circumbinary disk to the circumstellar disks.

The streamlines show eddy-like streams around the L4 and L5 points in the gap of the circumbinary disk. These streams are similar to a tadpole orbit, which is one of the closed orbits in the restricted three-body problem. The eddy-like streams are smoothly connected to those in the spiral arms and the circumbinary disk. Note that there are two stagnation points near the L2 and L3 points, at (x, y) ≃ (−1.2, 0.5) and (0.9, −0.6), respectively.

Figure 5 shows the gas flow more quantitatively. The two spiral arms have a rotation velocity faster than the Keplerian velocity (Figure 5, left panel). The gas of the spiral arms is accelerated by the gravitational torque from the orbiting stars, thereby showing fast rotation.

Zoom In Zoom Out Reset image size

Along the spiral arms, the radial velocity indicates both infall and outflow (Figure 5, middle panel). The inner portion of each spiral arm exhibits inflow, while the outer portion exhibits outflow. The boundary between the inflow and outflow coincides with the location of the stagnation point near the L2 or L3 point for each spiral arm. Note that the circumbinary disk (outside the outermost contour of the Roche potential) has both inflow and outflow regions. In the circumbinary disk, outflow regions coincide with the spiral arms, while inflow regions coincide with the inter-arm regions. Such inflow and outflow are caused by gravitational torque from the orbiting stars and have definitely been detected by ALMA observations (Takakuwa et al. 2014).

The radial mass flux is examined using the distribution of ${\rm{\Sigma }}{v}_{r}$ in the x–y plane (Figure 5, right panel). Inside the gap, the secondary arm exhibits higher inflow (deeper blue) than the primary arm does. This is consistent with the fact that the secondary has a higher accretion rate than the primary (Figure 6). The bridge between the circumprimary and circumsecondary disks shows a positive radial flux (pale red), indicating that the gas is transferred from the primary side to the secondary side over the bridge. However, this flux is considerably lower than that through the secondary arm, thereby indicating that accretion onto the circumsecondary disk occurs mainly through the secondary arm.

Zoom In Zoom Out Reset image size

Figure 6 shows the accretion rates for the primary and secondary for the fiducial model. To reduce the short-term variability of the accretion rates, we measure the accretion rates as average rates, ${\dot{M}}_{i}(t)=({M}_{i}(t-{\rm{\Delta }}t/2)-{M}_{i}(t+{\rm{\Delta }}t/2))/{\rm{\Delta }}t$ for i = 1, 2 (the primary and secondary), where M_i denotes the accumulated mass of the gas that accretes onto each star and Δt ≃ 1.5 (≃0.25 rotation period). The accretion rate for the secondary (${\dot{M}}_{2}$) reaches a constant value by the early phase, while that for the primary (${\dot{M}}_{1}$) gradually increases with significant modulation. The secondary has a higher accretion rate than the primary by a factor of ${\dot{M}}_{2}/{\dot{M}}_{1}\sim 4\mbox{--}6$ after t ∼ 300.

Figure 6 also shows the dimensionless variable of

$$\begin{eqnarray}&&{\rm{\Gamma }}=\displaystyle \frac{\dot{q}/q}{{\dot{M}}_{\mathrm{tot}}/{M}_{\mathrm{tot}}}=\displaystyle \frac{(1+q)\left({\dot{M}}_{2}-q{\dot{M}}_{1}\right)}{q\left({\dot{M}}_{1}+{\dot{M}}_{2}\right)}.\end{eqnarray} \tag{ 1 }$$

The variable Γ is the change rate of the mass ratio q normalized by the accretion rate onto both of the stars. The fiducial model exhibits Γ ∼ 4.5–5, which is in agreement with corresponding models in Bate & Bonnell (1997) and Young et al. (2015). The temporal average of the last 20 revolutions is Γ = 4.76 (Table 1). The positive Γ indicates an increase in the mass ratio.

For many models, a circumbinary disk shows asymmetry in the density and surface density distributions. For the fiducial model with (q, j_inf, c) = (0.2, 1.2, 0.1), as shown in Figure 1(b), the circumbinary disk has higher density on the lower side at the stage shown there. In Figure 2, the surface density of the circumbinary disk is higher on the lower left side. The asymmetric pattern rotates in the azimuthal direction at an angular velocity lower than that of the binary Ω_⋆. Therefore, the relative position angle between the asymmetry and the binary stars changes as time proceeds.

Note that the spiral arms rotate at the binary angular velocity Ω_⋆ in the gap and inner region (R ≲ 1.5 for the fiducial model) of the circumbinary disk. The waves of the spiral arms accumulate to form a local density bump in the circumbinary disk, which causes the asymmetry there, as can be clearly observed in Figure 2.

In order to evaluate the asymmetry of the circumbinary disk, we measure the surface density in thin cylindrical shells with radii of R = R_ring ± 0.1 for R_ring = 2, 4, 6 for each stage, as shown in Figure 7. The innermost shell (R_ring = 2) has the largest asymmetry amplitude among the three shells, thus indicating that the asymmetry is larger in the inner region of the circumbinary disk. The inclined stripes in the figure indicate that the asymmetry rotates at an almost constant angular velocity. Note that the slope of the stripes coincides with the angular velocity of the asymmetry.

Zoom In Zoom Out Reset image size

In order to more quantitatively analyze the asymmetry of the circumbinary disk, the Fourier transform is performed on the surface density in the azimuthal distribution for each stage, yielding the Fourier components of the mth mode, a_m. The amplitudes of the mth modes, $| {a}_{m}|$, and their phases (or position angles) $\arg ({a}_{m})$ are therefore evaluated. The time differential of the phase coincides with the angular velocity of the mth mode.

The amplitude of the m = 0 mode, $| {a}_{0}|$, denotes the azimuthal mean surface density (Figure 8, bottom panel). The mean surface density at R_ring = 2 increases almost linearly. Note that the ordinates are on a logarithmic scale. The radial growth of the circumbinary disk is also seen in the mean values of the surface density at R_ring = 4 and 6. These values also show an oscillation with a period of ∼10 revolutions, indicating that the circumbinary disk oscillates radially with that period. The oscillation at R_ring = 6 is reflected by the oscillation of the outer edge of the circumbinary disk.

Zoom In Zoom Out Reset image size

The ratio $| {a}_{1}| /| {a}_{0}|$ denotes the relative amplitude of the asymmetry (Figure 8, middle panel). At t ≳ 200, the relative amplitude reaches ∼50% for the surface density at R_ring = 2 and 4 with considerable oscillation, thereby indicating that the azimuthal contrast is ∼3 (=150%/50%) at most. The m = 1 mode has the largest amplitude among the modes with m ≥ 1.

The angular velocity of the asymmetry is measured by ${{\rm{\Omega }}}_{p}=\left[{\phi }_{1}(t+5{T}_{\star })-{\phi }_{1}(t-5{T}_{\star })\right]/(10{T}_{\star })$, which gives a temporal mean of the angular velocity over 10 revolutions (10T_⋆), where ${\phi }_{1}(t)=\arg \left[{a}_{1}(t)\right]$ (Figure 8, top panel). Note that the angular velocity Ω_p is for the asymmetric pattern and not for the gas velocity. The angular velocity of the asymmetry exhibits a constant value when the relative amplitude has a significant value: Ω_p ∼ 0.25Ω_⋆ at t ≳ 200 for R_ring = 2 and 4. This angular velocity suggests that the rotation of the asymmetric pattern is caused by resonance with the binary orbit, showing Ω_p : Ω_⋆ ≃ 1 : 4. More precisely, the mean angular velocity of the asymmetry is measured as Ω_p = 0.266Ω_⋆ in the period of 80–100 revolutions for R_ring = 2, as shown in the fourth column of Table 1. Thus, the asymmetry has a corotation radius of R = 2.42 when assuming Keplerian rotation. The dependence of Ω_p on the model parameters is shown in Sections 4.6 and 4.7.

The rotating asymmetric pattern causes the vortex in the circumbinary disk. Figure 9 shows the temporal averages of the velocity, surface density, and density distributions. The temporal averages are calculated in the rotating frame with the angular velocity Ω_p, which is the angular velocity of the asymmetric pattern. As shown in the face-on view (Figure 9, left and middle panels), the asymmetry is clearly seen in both the surface density and density distributions. The peak surface density and density are, respectively, ∼4 and ∼6 times as high as those on the opposite sides, at the stages shown here. In the figure, the position of the density peak nearly coincides with that of the surface density, (x, y) ≃ (−2, 0).

Zoom In Zoom Out Reset image size

The streamlines show the vortex, the center of which nearly coincides with the surface density and density peaks, but it is slightly shifted in the rotational direction. The aspect ratio of the vortex (the aspect ratio of the closed streamlines) is ∼2. Although it has been suggested that such a thick vortex (a small aspect ratio) would be unstable (e.g., Lesur & Papaloizou 2009; Richard et al. 2013), the vortex here remains for a long time. This suggests that the vortex continues to be excited.

The thickness of the circumbinary disk is roughly the same between the dense and less dense sides (Figure 9, right panel). This indicates that asymmetry in the surface density can be attributed to the density rather than the thickness of the disk.

The left panel of Figure 10 shows the specific angular momentum distribution on a color scale. The stage and spatial scale for this snapshot are the same as those in the upper right panel of Figure 2. These figures show that the spiral arms have a high specific angular momentum. The outermost spiral arms have an especially high value of j ≳ 1.8. This high specific angular momentum is attributed to the gravitational torque of the binary. Just inside this spiral arm, a vortex is excited, as shown in the streamlines around (x, y) ∼ (0, −2). The surface density bump is associated with the vortex at roughly the same position (see Figure 2).

Zoom In Zoom Out Reset image size

This structure is observed clearly when we take a temporal average of the disk (Figure 10, right panel). The center of the vortex (black streamlines) nearly coincides with the density bump (red contours). It is located at the corotation radius with the angular velocity of the asymmetry ${{\rm{\Omega }}}_{p}$ for the Keplerian rotation (see white circle). We also find that the peak of the specific angular momentum j (blue contours) is located in almost the same orientation as the density peak and the center of the vortex with respect to the origin, and confirm that this configuration is maintained in the other periods, i.e., the 30–50 and 50–70 revolutions periods. This indicates that the density bump and vortex are caused by the locally enhanced angular momentum, which results from the gravitational torque. Further discussion of the asymmetry is made in Section 5.1.

Figure 11 shows the dependence of the surface density distribution on the temperature of the gas or the sound speed. The cold models (small c) have circumbinary disks with narrow ring shapes, while the hot models (large c) have circumbinary disks with large radii. Moreover, a colder disk exhibits sharper contrast in the surface density distribution. This temperature dependence is due to changes in the thermal pressure.

Zoom In Zoom Out Reset image size

The cold models also have asymmetric surface density distributions in the circumbinary disks. The model with (q, j_inf, c) = (0.2, 1.2, 0.01) has a narrow circumbinary disk that is off-center with respect to the barycenter of the binary stars (upper left panel of Figure 11). While this circumbinary disk is almost concentric in the early phase, it becomes eccentric at ∼50 revolutions of the binary. The narrow circumbinary disk is fitted by an elliptic orbit, the focus of which coincides with the barycenter of the binary stars (white ellipse in the upper left panel of Figure 11). The elliptic orbit is parameterized by the angular momentum, eccentricity, and periapsis. The angular momentum is obtained by the density-weighted average in the circumbinary disk, and the eccentricity and periapsis are obtained by the least mean squares method. The angular momentum and the eccentricity gradually increase as time proceeds because of the gravitational torque of the stars. The periapsis of the circumbinary disk rotates slowly. The angular velocity is Ω_p = 0.00504Ω_⋆ in the period of 80–100 revolutions, as shown in Table 1. This low angular velocity is also shown in Figure 12, in which the peak surface density diagrams show a shallow slope. Note that the two peaks of the surface density in the panel of R_ring = 2 correspond to the two intersections between the eccentric circumbinary disk and the cylindrical ring with R_ring = 2. Figure 13 shows low Ω_p in the stages where the asymmetry $| {a}_{1}| /| {a}_{0}|$ has a large amplitude.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The slow rotation of the asymmetry in the cold models is qualitatively different from that in the fiducial model. Moriwaki & Nakagawa (2004) derived a formula for the precession rate of a test particle orbiting a binary using a secular perturbation theory, and Thun et al. (2017) modified it as follows (see also Leung & Lee 2013; Dunhill et al. 2015):

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{p}=\displaystyle \frac{3}{4}\displaystyle \frac{q}{{\left(1+q\right)}^{2}}{\left(\displaystyle \frac{D}{{a}_{p}}\right)}^{7/2}{\left(1-{e}_{p}^{2}\right)}^{-2}{{\rm{\Omega }}}_{\star },\end{eqnarray} \tag{ 2 }$$

where a_p and e_p are the semimajor axis and eccentricity of the test particle, respectively. Equation (2) is an expression for the case of a circular orbiting binary. For the model with c = 0.01, the surface density of the circumbinary disk is fitted to an elliptic orbit using the method described above, yielding a_p = 3.0 and e_p = 0.42 as the temporal averages in the period of 80–100 revolutions. Using these orbital parameters, Equation (2) gives Ω_p = 0.0032Ω_⋆. For the model with c = 0.02, we obtain a_p = 2.9, e_p = 0.40, and Ω_p = 0.0033Ω_⋆. These values are consistent with the value listed in Table 1, indicating that the secular motion of the ballistic particle causes the rotation of the asymmetry of the circumbinary disk and hydrodynamical effects are minor in the cold models.

The hot models show the circumbinary disks extended up to large radii as shown in the lower panels of Figure 11. A circumbinary disk extends faster with a higher sound speed. Half of the outer edge of the circumbinary disk touches the computational boundary at $R=12D$ at the stages of 80, 50, and 31 revolutions for the models with c = 0.15, 0.175, and 0.2, respectively. For models with higher sound speed, both the surface density contrast and the amplitude of the asymmetry are lower. The extension of the circumbinary disks and the low contrast are responsible for the high pressure seen in the hot models.

Figures 14 and 15 are the same as Figures 7 and 8, respectively, but for the model with the highest sound speed (c = 0.2). For this model, the amplitude of the asymmetry is as low as $| {a}_{1}| /| {a}_{0}| \sim 0.2$ (Figure 15). The angular velocity of the asymmetric pattern is a very low value of Ω_p = 0.0335Ω_⋆ in the period of 80–100 revolutions (see also Table 1). The corresponding corotation radius is R = 9.62. This slow rotation is observed in the slope of the diagrams in Figure 14.

Zoom In Zoom Out Reset image size

The accretion rates are affected by the temperature. All the models show a tendency to decrease in Γ with time, while the cold models with c ≤ 0.02 show a rapid decrease in Γ, resulting in low Γ at the end of the calculations (Table 1). The decrease in Γ is mainly attributed to an increase in ${\dot{M}}_{1}$, as observed in the fiducial model (Figure 6). Young et al. (2015) reported a model with hot gas (c = 0.25) showing variability in accretion rates, the timescale of which was t ∼ 90. Our model with c = 0.2 also exhibits variability mainly in the accretion rate of the secondary, and the timescale is longer (t ∼ 100–200) than that of Young et al. (2015).

Zoom In Zoom Out Reset image size

Figure 16 shows representative models with different specific angular momenta for the infalling gas. We observe that the models with ${j}_{\inf }\,\leqslant \,0.8$ do not have circumbinary disks, while the models with j_inf ≥ 0.9 have circumbinary disks. This is in agreement with the results of the SPH simulations by Bate & Bonnell (1997). They implicitly used the condition for circumbinary disk formation, in which the centrifugal radius of the gas, ${R}_{\mathrm{cent}}={j}_{\inf }^{2}/({{GM}}_{\mathrm{tot}})$, is larger than the orbital radius of the secondary star, a₂ = D/(1 + q), yielding

$$\begin{eqnarray}&&{j}_{\inf }\,\geqslant \,{\left(\displaystyle \frac{{{GM}}_{\mathrm{tot}}D}{1+q}\right)}^{1/2}.\end{eqnarray} \tag{ 3 }$$

For q = 0.2, Equation (3) gives j_inf ≥ 0.91, which is in good agreement with our simulations.

Zoom In Zoom Out Reset image size

The models with circumbinary disks (j_inf ≥ 0.9) show ${\dot{M}}_{2}\gt {\dot{M}}_{1}$ throughout the evolution. On the other hand, the model with j_inf = 0.8 shows ${\dot{M}}_{1}\gt {\dot{M}}_{2}$ after ∼50 revolutions. The model with j_inf = 0.5 shows ${\dot{M}}_{1}\gt {\dot{M}}_{2}$ throughout the evolution and has a negative Γ (Table 1). The dependence of j_inf on the accretion rates is consistent with Bate & Bonnell (1997).

The models with large j_inf have extended circumbinary disks, as shown in the lower panel of Figure 16. The surface density distribution for the circumbinary disk has a peak around the centrifugal radius specified by j_inf (e.g., R_cent = 2.56, 3.24, 4.0 for j_inf = 1.6, 1.8, 2.0, respectively). For larger j_inf, the circumbinary disk gap and the circumprimary and circumsecondary disks all show lower surface densities. The circumprimary disk is influenced more by the change in j_inf than the circumsecondary disk. The mass ratio for the circumstellar disks is M_CSD/M_CPD = 1.26 at 100 revolutions for the model with j_inf = 2.0, and the circumsecondary disk is more massive than the circumprimary disk.

The asymmetry of the circumbinary disk is qualitatively different when j_inf ≥ 1.8. In these models, there are several mth modes with considerable amplitudes, and the m = 1 mode does not have the largest amplitude (see also Figure 17 for the model with j_inf = 1.8). Fourier analysis shows that the amplitudes $| {a}_{m}|$ oscillate and that the largest mode with m ≥ 1 changes with time.

Zoom In Zoom Out Reset image size

Figure 18 compares models with different mass ratios q. The models with q = 0 and 1 correspond to a single star and an equal-mass binary, respectively. In the binary cases ($q\ne 0)$, the structures of the circumbinary disks are not very sensitive to the mass ratio q. All the binary models examined here show asymmetry in their circumbinary disks, with almost the same amplitude and angular velocity Ω_p irrespective of q (see also Table 1).

Zoom In Zoom Out Reset image size

For all the binary models, the radial extent of each circumbinary disk oscillates as seen in the fiducial model. The radial extent is more sensitive to the stage than to the mass ratio because of this oscillation.

The model with q = 0 corresponds to the case of a single star (the leftmost panel of Figure 18). The disk has a narrow ring shape with almost constant specific angular momentum j_inf. The central star does not cause gravitational torque, and the static gravity of the point mass acts on the circumstellar disk. We observe asymmetry of the m = 1 mode in the ring, the amplitude of which is small compared to those in the binary models: $| {a}_{1}| /| {a}_{0}| \lesssim 0.2$ at R_ring = 2. The angular velocity of the asymmetric pattern is Ω_p = 0.589 (see also Table 1), which coincides with that of Keplerian rotation with an angular momentum of j_inf = 1.2. This angular velocity is considerably different from those of the binary models (the models with $q\ne 0$). Moreover, there is no significant accretion onto the star ${\dot{M}}_{\star }/{\dot{M}}_{\mathrm{env}}\lesssim 0.1 \%$, where ${\dot{M}}_{\star }$ is the accretion rate and ${\dot{M}}_{\mathrm{env}}$ is the gas injection rate into the computational domain. In summary, the disk in the single-star model shows different features from the circumbinary disks, indicating that the star's gravitational torque has considerable impact on the circumbinary disks.

In all the models examined above, the gas is injected from the computational boundaries, mimicking the infalling envelope. We now consider an additional model in which the gas injection at the boundaries is terminated at 50 revolutions of the binary.

Figure 19 shows the surface density distribution at 150 revolutions. This stage occurs 100 revolutions after the gas injection has ended. This model exhibits a smooth surface density distribution compared to the models with gas infall (compare with the top right panel of Figure 2 and the top right panel of Figure 11). This indicates that the accretion onto the circumbinary disk excites a turbulent feature there. We also observe that the circumbinary disk extends after gas injection ends, indicating that the ram pressure of the infalling gas confines the circumbinary disk to relatively compact sizes.

Zoom In Zoom Out Reset image size

The asymmetry of the circumbinary disk appears regardless of gas injection, as shown in Figure 19. The angular velocity of the asymmetry is measured as Ω_p = 0.254 and 0.235 for 80–100 and 130–150 revolutions, respectively. The angular velocity is insensitive to gas injection.

The accretion rate for each star is influenced by gas injection. As shown in Figure 20, ${\dot{M}}_{2}$ decreases down to the level of ${\dot{M}}_{1}$ after gas injection ends.

Zoom In Zoom Out Reset image size

The masses of the circumbinary disk and circumstellar disks are also influenced by the termination of the gas injection. As shown in Figure 21, the mass of the circumbinary disk decreases slowly after the gas injection ends, in contrast to the model shown in Figure 3, in which the gas accretion continues. This slow decrease is due to the gas infall onto the circumstellar disks. The masses of the circumstellar disks also decrease after the gas injection ends. In particular, the mass of the circumsecondary disk decreases by an order of magnitude in 60 revolutions after the gas injection ends, and then it remains almost constant. At 150 revolutions, the ratio of the circumstellar disk is M_CSD/M_CPD = 0.37, while both the disks have almost the same mass in the model with the same parameters but with uninterrupted gas injection.

Zoom In Zoom Out Reset image size

In addition, we consider two rigidly rotating envelope models in which the infalling gas has a constant angular velocity at computational boundaries. In these models, the gas envelope has a range of specific angular momenta 0 ≤ j ≤ j_inf. Figure 22 shows the surface density distribution for the models with j_inf = 1.2 (see Section 2). The surface density distribution resembles that for the model with constant specific angular momentum (the upper right panel of Figure 2), showing asymmetry in the circumbinary disk. The angular velocity of the asymmetric pattern Ω_p has almost the same value as in the fiducial model (see Table 1).

Zoom In Zoom Out Reset image size

The spiral arms remain in this model, but their density contrast is lower than that in the fiducial model because the gap has a higher density. This is because gas with low angular momentum falls along the polar directions and fills the gap. Such gas increases the accretion rate for the primary. The primary has an accretion rate that is two times as high as that of the secondary (Figure 23). The change rate for the mass ratio is low, Γ ∼ 1, but it is still positive.

Zoom In Zoom Out Reset image size

For the envelope model with a higher rotation of j_inf = 2.0, the circumbinary disk evolves in a way similar to that for the corresponding model with constant specific angular momentum (the lower right panel of Figure 16) in the early stages. After ∼60 revolutions, significant asymmetry of the m = 1 mode appears at R ∼ 2–3. As in the previous model, the gap has higher density than in the model with constant specific angular momentum. The primary star has a higher accretion rate than the secondary star by a factor of 2–3 throughout the evolution. These differences mainly come from gas with low angular momentum falling from the polar regions.

The three models in this subsection demonstrate that the accretion rates for the binary stars are sensitive to the angular momentum of the infalling envelope.

As shown in the fiducial model (Figure 10), the density bump and vortex coexist in the same location. Moreover, a ridge of high specific angular momentum exists at the same angular position relative to the origin. This configuration can be explained as follows. The gravitational torque of the binary acts on the gas close to the inner edge of the gap to increase its specific angular momentum. The gas with enhanced angular momentum moves outward radially, exchanging the radial positions of the gas. This motion is observed as a vortex around the density bump.

This is similar to the situation where Rayleigh's criterion (e.g., Chandrasekhar 1961) is violated: a disk with a negative gradient of specific angular momentum is highly unstable. Since the initial state with a negative gradient of specific angular momentum is unstable, such a state may not appear in the numerical simulations. The violation of Rayleigh's criterion is sufficient to cause Rossby wave instability (RWI) in many cases, which causes vortices and pressure bumps (Li et al. 2000; Ono et al. 2016, 2018). RWI requires the vortensity, ${\left({\rm{\nabla }}\times {\boldsymbol{v}}\right)}_{z}/{\rm{\Sigma }}$, to have a local minimum for a barotropic gas (Lovelace et al. 1999; Ono et al. 2016). We confirm that the vortensity in our models has a local minimum at the location of the vortex.

The asymmetry does not result from keeping the angular momentum of the infalling gas constant. During several rotations of the binary, the gravitational torque produces an angular momentum distribution in the circumbinary disk. Consequently, the disks extend above the centrifugal radius for the original angular momentum, R_cent = j_inf²/(GM_tot). Note that the models with a rigidly rotating envelope exhibit asymmetry in the circumbinary disk. In these models, the infalling gas has a range of angular momenta.

For a binary system with a circumbinary disk, the stars rotate with higher angular velocity than that of the circumbinary disk. Therefore, the circumbinary disk receives angular momentum from the stars via gravitational torque. The spiral arms, gaps, and asymmetry in the circumbinary disk are natural outcomes of binary systems with circumbinary disks. We also note that a circumbinary disk has a higher angular momentum than the infalling envelope as a result of such gravitational torque (e.g., Figure 10).

In our models with very high specific angular momenta (j_inf ≳ 1.8), the centrifugal radii are large enough for the gravitational torque to act weakly on the circumbinary disks. When we consider the rigidly rotating envelope model, which may sometimes be realistic, the asymmetry appears as described in Section 4.9. Moreover, when we consider a binary orbit with high eccentricity, the binary may interact with the circumbinary disk effectively, even in high angular momentum cases (e.g., Dunhill et al. 2015; Price et al. 2018). This results in asymmetry in the circumbinary disk.

Recent high-resolution observations have revealed circumbinary disks with asymmetric structures. Takakuwa et al. (2017) observed an asymmetric circumbinary disk around the protobinary L1551 NE and compared their observations with numerical simulations based on the models described in this paper. Takakuwa et al. (2014) analyzed the density and velocity structures in the circumbinary disk and found evidence that they were influenced by gravitational torque. Such gravitational torque is a driving force behind the asymmetry of the circumbinary disk, as shown in Section 4.5.

The young binary UY Aur exhibits asymmetric brightness in the circumbinary disk, as seen in near-infrared scattered light (Hioki et al. 2007). Hioki et al. (2007) examined a model where the near side of the disk is brighter than the far side because of forward scattering by the dust. The model requires a thick circumbinary disk with a disk opening angle of Φ_open = 45°. Our 3D models reproduce thick circumbinary disks as shown in the lower panels of Figure 1 and the right panel of Figure 9, which are likely to be satisfied by the condition of the opening angle, depending on the disk surface for scattering light. The opening angle becomes small when light is scattered at high density. The large thickness of such circumbinary disks is attributed to the spiral arms, which disturb the disk and increase its vertical velocity.

The young binary HD 142527 has an asymmetric circumbinary disk showing a horseshoe dust structure (e.g., Fukagawa et al. 2013). Not only dust but also gas seems to have an asymmetric structure with a significant contrast level (Boehler et al. 2017). The gas contrast has been reported as 3.75 and is consistent with the contrast of the fiducial model, as shown in Section 4.4. However, Lacour et al. (2016) indicated that the companion is unlikely to be responsible for the gap of the circumbinary disk because of the small apocenter distance compared to the gap radius. The disk–companion interaction is likely to be too weak to induce asymmetry in the circumbinary disk via the mechanism described in this paper. Such a weak interaction is seen in the models with a high specific angular momentum (see lower right panel of Figure 16), where the asymmetry is unclear. Recently, Price et al. (2018) reproduced the asymmetric circumbinary disk by assuming that the companion with highly eccentric orbits interacts with the circumbinary disk.

Recent observations have resolved each of the circumstellar disks in young binary systems (e.g., Lim et al. 2016a, 2016b; Takakuwa et al. 2017). The masses of circumstellar disks are sensitive to the gas supply from the infalling envelopes, as are the masses of the circumbinary disks. When the infalling envelope consists of gas with a large angular momentum, the circumsecondary disk tends to have a larger mass than the circumprimary disk. In contrast, when the infalling envelope has a low angular momentum even in part, it contributes to the growth of the circumprimary disk because the gas falls onto the central region of the binary. The disk masses are also affected by the viscosity of the disk, which is not taken into consideration here.

All the circumstellar disks in the simulations have clear outer edges. Such circumstellar disks with sharp outer edges have been observed in the young binaries L1551 IRS 5 and L1551 NE by Very Large Array observations (Lim et al. 2016a, 2016b). The simulations indicate that the sizes are insensitive to the gas temperature and the angular momenta of the infalling envelopes. The sizes of the circumstellar disks reflect those of the Roche lobes. In other words, tidal truncation controls the disk sizes, as suggested by Lim et al. (2016a, 2016b).

Circumbinary planets are considered to be formed in circumbinary disks. The vortices in stellar disks can capture solid particles and initiate the formation of planetesimals (Barge & Sommeria 1995; Klahr & Henning 1997). Vortices and density and pressure bumps exist in our simulations. The vortex persists in the circumbinary disk, as seen in the simulations here, while it tends to migrate inward in the circumstellar disk around a single star (Richard et al. 2013). This is because the vortex shown here is continuously excited by the gravitational torque of the binary.

Turbulence hinders the formation of planetesimals. The turbulent feature in the circumbinary disk calms down after accretion onto the circumbinary disk is terminated (Figure 19). After that, the vortex in the circumbinary disk still exists. This suggests that the formation of planetesimals proceeds more efficiently after the main accretion phase.

Shi et al. (2012) also observed asymmetry in the circumbinary disks of an equal-mass binary using MHD simulations. They assumed a relatively cold gas with c = 0.05 and a weak magnetic field of β = 100. According to their simulations, the angular velocity of the asymmetry was Ω_p ≃ 3.2 × 10⁻³ Ω_⋆, which is significantly smaller than those of our models with moderate temperatures (c ≥ 0.05) but consistent with our cold models (c ≤ 0.02). Shi & Krolik (2015) doubled the sound speed (c = 0.1) and still observed asymmetry in the circumbinary disk. The angular velocity of the asymmetry had higher values of Ω_p = (0.17–0.21)Ω_⋆, which are roughly consistent with our moderate-temperature models. The difference probably comes from the fact that the gap in the circumbinary disk was located outside the computational domain in Shi & Krolik (2015) and the gas flow inside the gap was not reproduced.

Ochi et al. (2005) investigated accretion rates onto binary stars using 2D numerical simulations. They argued that the primary accretes more than the secondary even when more gas falls through the L2 point than through the L3 point. This means that the gas is transferred from the secondary to the primary through the bridge between the primary and secondary. In contrast, our simulations show that the secondary has a higher accretion rate than the primary in almost all the models examined here, and only the model with very low angular momentum ${j}_{\inf }\,=\,0.5$ shows a negative value of Γ (Table 1). Moreover, the gas tends to flow from the primary to the secondary through the bridge, as shown in Figure 4 (right) and Figure 5 (right). The results for the accretion rates in this paper are consistent with Bate & Bonnell (1997) and Young et al. (2015).

Satsuka et al. (2017) recently performed 3D simulations to investigate accretion onto the seeds of binary stars. In their models, an infalling envelope changes its density and angular momentum based on a collapsing cloud core model as time proceeds. Because they calculated the models for the early phase (t ≲ 20–50), some models exhibited spiral arms without forming circumbinary disks. Their results may reflect the initial conditions of short-term evolution.

The models examined here ignore the self-gravity of the gas. The self-gravity is important just after the fragmentation of collapsing cloud cores, where the gas envelope around the fragments is more massive than the fragments. The orbits of the fragments change as time proceeds because of the interactions between the fragments and the gas envelope and between the fragments themselves, as shown by Matsumoto & Hanawa (2003) and Matsumoto et al. (2015b), who performed numerical simulations that took self-gravity into account. In order to reproduce the evolution of binaries in the early phase, not only self-gravity but also its back-reaction on the orbits of the stars (or fragments) should be taken into consideration.

In the case of the protobinary L1551 NE, the mass of the circumbinary disk is 0.009–0.076M_☉, and the masses of the primary and secondary stars are 0.68M_☉ and 0.13M_☉, respectively (Takakuwa et al. 2017). The disk-to-star mass ratio is therefore M_CBD/(M₁ + M₂) = 0.01–0.09, indicating that the circumbinary disk is not likely to be gravitationally unstable (Kratter & Lodato 2016). However, the mass of the circumbinary disk could be comparable to that of the secondary star when it has the upper limit value. In this case, the orbital motion of the secondary star may be affected by the gravity of the circumbinary disk.

We assume the isothermal equation of state, ignoring temperature distribution in the gas. The temperature distribution of the gas should be determined by many factors: compression heating due to infalling gas and shock waves, radiation heating from protostars, cosmic-ray heating, and cooling by dust and molecular emission. Modeling the realistic temperature distribution increases a number of model parameters, and it also increases computational costs for simulations on long-term evolution. The dynamics of the circumbinary gas is significantly affected by the gas temperature as shown in Figure 11. The isothermal assumption is useful to investigate the dependence on the temperature.

We have examined models in which the infall gas is in the steady state. This results in an infalling envelope with a spatially constant specific angular momentum. Such infalling envelopes have been suggested by observations of the envelopes around protostars (e.g., Ohashi et al. 1997).

The observations have also suggested that infalling envelopes increase their specific angular momenta during their evolution (Yen et al. 2011, 2013, 2017). The models examined here do not include time-dependent infall into the system. When considering a binary with a separation of less than ∼100 au and masses of ∼1M_☉, the rotation period is less than 10³ yr, which is considerably shorter than the period of the accretion phase, i.e., 10^4–5 yr, thereby indicating that our quasi-steady-state models are applicable to such binary systems.

The magnetic fields redistribute the angular momenta in the circumbinary disk and circumstellar disks. If we take the magnetic fields into account, the accretion rates for both the stars will increase. The extent of the circumbinary disks can be influenced by the magnetic fields. Change in the angular momentum distribution can affect the orbital evolution of binaries (e.g., Zhao & Li 2013). The effects of the magnetic fields depend highly on magnetic diffusion, which is controlled by the degree of ionization of the gas. The models examined here correspond to the case where the degree of ionization is low. The effect of the MHD will be investigated in future work.