From large-scale to protostellar disk fragmentation into close binary stars

From large-scale to protostellar disk fragmentation into close binary stars

Leonardo Di G. Sigalotti    Fidel Cruz    Ruslan Gabbasov    Jaime Klapp    José Ramírez-Velasquez
January 1, 2018January 7, 2018July 1, 2019
January 1, 2018January 7, 2018July 1, 2019
January 1, 2018January 7, 2018July 1, 2019

Recent observations of young stellar systems with the Atacama Large Millimeter/submillimeter Array (ALMA) and the Karl G. Jansky Very Large Array (VLA) are helping to cement the idea that close companion stars form via fragmentation of a gravitationally unstable disk around a protostar early in the star formation process. As the disk grows in mass, it eventually becomes gravitationally unstable and fragments, forming one or more new protostars in orbit with the first at mean separations of 100 astronomical units (AU) or even less. Here we report direct numerical calculations down to scales as small as AU, using a consistent Smoothed Particle Hydrodynamics (SPH) code, that show the large-scale fragmentation of a cloud core into two protostars accompanied by small-scale fragmentation of their circumstellar disks. Our results demonstrate the two dominant mechanisms of star formation, where the disk forming around a protostar, which in turn results from the large-scale fragmentation of the cloud core, undergoes eccentric () fragmentation to produce a close binary. We generate two-dimensional emission maps and simulated ALMA 1.3 mm continuum images of the structure and fragmentation of the disks that can help explain the dynamical processes occurring within collapsing cloud cores.

circumstellar matter — methods: numerical — planets and satellites: formation — protoplanetary disks — stars: formation — stars: pre-main sequence
Corresponding author: Jaime\move@AU\move@AF\@affiliation

Área de Física de Procesos Irreversibles, Departamento de Ciencias Básicas,
Universidad Autónoma Metropolitana - Azcapotzalco (UAM-A), Av. San Pablo 180, 02200 Ciudad de México, Mexico


Instituto de Ciencias Básicas e Ingenierías, Universidad Autónoma del Estado de Hidalgo (UAEH),
Ciudad Universitaria, Carretera Pachuca-Tulancingo km. 4.5 S/N, Colonia Carboneras,
Mineral de la Reforma, C.P. 42184, Hidalgo, Mexico


Instituto Nacional de Investigaciones Nucleares (ININ), Departamento de Física
Carretera México-Toluca km. 36.5, La Marquesa, 52750 Ocoyoacac, Estado de México, México \move@AU\move@AF\@affiliationabacus-Centro de Matemática Aplicada y Cómputo de Alto Rendimiento,
Departamento de Matemáticas, Centro de Investigación y de Estudios Avanzados (Cinvestav-IPN),
Carretera México-Toluca km. 38.5, La Marquesa, 52740 Ocoyoacac, Estado de México, México


Instituto Venezolano de Investigaciones Científicas (IVIC), Centro de Física,
Apartado Postal 20632, Caracas 1020A, Venezuela \move@AU\move@AF\@affiliationESPOL Polytechnic University, Escuela Superior Politécnica del Litoral (ESPOL),
Physics Department, FCNM, Campus Gustavo Galindo km. 30.5, Vía Perimetral,
P.O. Box 09-01-5863, Guayaquil, Ecuador

1 Introduction

Stellar systems, consisting of two or more stars, are formed by two different mechanisms – large-scale fragmentation of cloud cores of gas and dust during their early isothermal collapse phase (Larson, 2001; Fisher, 2004), leading to widely separated companion protostars [ astronomical units (AU) or larger], and small-scale fragmentation of their circumstellar disks early in the collapse due to gravitational instabilities (Rodríguez et al., 2005; Connelley et al., 2008; Pech et al., 2010; Murillo & Lai, 2013; Tobin et al., 2013), leading to bound systems separated by tens of AU. This scenario of star formation seems to be consistent with multiplicity surveys of the stellar population in nearby star-forming regions (Connelley et al., 2008; Kraus et al., 2011; Tobin et al., 2016a). Circumstellar disks are then expected to play an important role at closer separations, implying that many young stars with widely separated companions should also have a close companion (Adams et al., 1989; Andalib et al., 1997). New clues fitting this idea have been recently revealed by Karl G. Jansky Very Large Array (VLA) and Atacama Large Millimeter/submillimeter Array (ALMA) observations of previously-unseen young protostars (Tobin et al., 2016a, b).

Numerical simulations of idealized protostellar disks with a central object have shown unstable behavior toward the development of long-wavelength, spiral shaped instabilities, with the eccentric mode (Bonnell & Bate, 1994; Woodward et al., 1994), when the disk-to-central-object mass ratio (Adams et al., 1989). However, these models lack physical reality since protostellar disks form and grow in mass within a collapsing environment as a result of the infall of material from the outer cloud envelope (Eisner, 2012; Tobin et al., 2012). On the other hand, if more protostars form within the same environment, the gravitational interaction among them may limit the disk size and affect the way angular momentum is transferred. However, the disk fragmentation model has succeeded in reproducing the constraints imposed by the observed statistical properties of low-mass objects, such as low-mass binary systems, brown dwarfs, and even planetary mass objects, which in turn have not yet been explained by means of other formation mechanisms (Stamatellos & Whitworth, 2009). Until now, no direct hydrodynamical collapse calculations have been reported in the literature that demonstrate the two dominant mechanisms for binary/multiple star formation mainly due to limited mass resolution.

Here we present the results of high-resolution protostellar collapse simulations, using a consistent Smoothed Particle Hydrodynamics (SPH) code. The calculations capture the structure and fragmentation of protostellar disks after the large-scale fragmentation of a cloud core by spanning spatial scales ranging from several thousands of AU to about 0.1 AU. We started the simulations from a gas sphere of uniform density and temperature representative of an idealized molecular cloud core, with parameters corresponding to the so-called “standard isothermal test case” (Burkert & Bodenheimer, 1993), coupled to a barotropic pressure-density relation to mimic the transition from the isothermal to the nonisothermal collapse phase (Boss et al., 2000). We improve on mass resolution by allowing the calculations to work with unprecedented large numbers of particle neighbors, which translates into approximate second-order accuracy of the whole SPH discretization. Because of the lack of mass resolution, until now it was assumed that the collapse of the standard test case ends up with the formation of a wide binary with each core forming just a single object, since no fragmentation was seen in the simulations during the formation of the first protostars (Kitsionas & Whitworth, 2002; Arreaga-García et al., 2007; Riaz et al., 2014). However, fragmentation of a protostellar disk does not necessarily stop with the formation of a second protostar. Although our simulations do not show how the disks evolve much beyond the initial disk fragmentation, present-day disk simulations predict that as fragmentation proceeds in high-mass accretion disks, new protostars can form at increasingly larger radii, with the process resulting in a chaotic multiple stellar system (Krumholz et al., 2009; Peters et al., 2010; Bate, 2018). On the other hand, since gas accretion from the outer envelope stops only when its angular momentum equals that of the protostars, it is relatively easier for them to gain fresh gas through the accretion disk, driving the system toward equal mass protostars (Bate & Bonnell, 1997). The paper is organized as follows. In Section 2, we provide a brief description of the numerical methods along with the initial conditions employed for the collapse model calculations. The results of the simulations are described in Section 3 and the conclusions are given in Section 4.

2 Numerical Methods

2.1 The SPH Scheme

The calculations of this paper were performed using a modified version of the SPH-based GADGET-2 code (Springer, 2005), where the interpolation kernel was replaced by a Wendland C function to allow support of large numbers of neighbors (Wendland, 1995; Dehnen & Aly, 2012). Wendland functions have positive Fourier transforms and so they can work with arbitrarily large numbers of neighbors without suffering from a pairing instability, where particles come into close pairs and become less sensitive to small perturbations within the kernel support (Dehnen & Aly, 2012). Moreover, Wendland kernels do not allow particle motion on a sub-resolution scale, maintaining a very regular particle distribution even in highly dynamical tests (Rosswog, 2015). An improved scheme for the artificial viscosity relying on the method developed by Hu et al. (2014) was also implemented, which uses a novel shock indicator based on the total time derivative of the velocity divergence with a limiter that applies the same weight to the velocity divergence and vorticity. This allows the artificial viscosity to distinguish true shocks from purely convergent flows and discriminate between pre- and post-shocked regions. With this method viscous dissipation is effectively suppressed in subsonically convergent flows and in regions where the vorticity dominates over the velocity divergence, thus avoiding spurious angular momentum transport in the presence of vorticity.

To ensure formal convergence and first-order (i.e., ) particle consistency of the SPH equations, we adopt power-law dependences to set the smoothing length () and the number of neighbors () within the kernel support in terms of the total number of particles () (Zhu et al., 2015). In particular, particle consistency (or first-order accuracy), i.e., satisfaction of the normalization condition of the kernel function in discrete form can only be achieved when is sufficiently large for which the finite SPH sum approximation approaches the continuous limit. This is consistent with the results of an error analysis of the SPH representation of the continuity and momentum equations, which show that particle consistency is lost due to zeroth-order errors that would persist when working with a fixed () number of neighbors even though and (Read et al., 2010). If particle consistency is achieved, then particle consistency is automatically ensured because of the symmetry of the kernel function. Out of the family of possible curves describing the dependence of on , we choose the scaling relations and to set the initial values of and in terms of . These scalings were derived by requiring that , which accommodates large numbers of neighbors for given while still keeping reasonably large values of . According to the parameterization given by Zhu et al. (2015), an exponent of corresponds to , which is an intermediate value in the interval . This interval is valid for low-order kernels as the Wendland C function and is appropriate for quasi-random particle distributions for which the particle approximation error goes as . As is increased, these relations comply with the asymptotic limits , , and with for full consistency of the SPH equations (Rasio, 2000). The above scalings have important implications on the minimum resolvable mass , where is the particle mass. Since and , then , implying improved mass resolution when the number of neighbors is increased. This feature makes a key difference with previous SPH collapse calculations since it allows resolution of small-scale structures in the flow through many orders of magnitude increase in density and pressure. Convergence and consistency testing on three-dimensional flow problems have demonstrated the second-order accuracy (i.e., the -particle consistency) of the code (Gabbasov et al., 2017).


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefEquatorial density evolution in a 266-AU region around the cloud center for model U2 after large-scale fragmentation into a binary protostar: (a) 44.14 kyr, (b) 44.39 kyr, (c) 44.43 kyr, and (d) 44.68 kyr. The logarithm of the density is plotted. For (a), (b), (c), and (d) the maximum density is: , , , and ( in g cm), respectively. The insets show a magnification of the binary cores and their whirling disks. The disks develop a two-armed spiral pattern associated with a linear growth of a gravitational instability. In (b) and (c), the gravitational instability in disk B enters a nonlinear growth phase and so one of its spiral arms is seen to fragment into a second protostar.

2.2 Initial Conditions

The initial conditions for protostellar collapse were chosen to be those of the standard isothermal test case (Burkert & Bodenheimer, 1993), coupled to a barotropic pressure-density relation to mimic the transition from the isothermal to the nonisothermal collapse phase (Boss et al., 2000). The models were assumed to remain isothermal at an initial cloud temperature of 10 K up to a critical density g cm, which separates the isothermal from the nonisothermal collapse. Two highly resolved calculations were considered (models U1 and U2), which differed only in the total number of particles. Model U1 started with million particles and neighbors, while model U2 had million particles and neighbors. The particles were initially distributed in a glass configuration (Couchman et al., 1995) within a sphere of mass , uniform density g cm, radius pc, and solid-body rotation s. The uniform-density background was perturbed azimuthally with a 10%, variation. With this choice of the initial parameters, the ratios of the thermal and rotational energies to the absolute value of the gravitational energy are, respectively, and . The models are close to virial equilibrium (), and apart from the assumption of uniform density, they could represent cloud cores that begin to infall after ambipolar diffusion has lessened the magnetic field support. The initial free-fall time is kyr. Although the standard isothermal model is an idealization of a real cloud core, it provides a simple model from which to learn how nonaxisymmetric perturbations grow from a structureless medium. In addition, when the transition from isothermal to adiabatic collapse is prolonged, convergence is more difficult to attain and resolution requirements become significantly more demanding. These initial conditions are therefore suitable for testing the resolution performance of the code from several thousands of AU down to scales of 1 AU or less.

3 Results

Models U1 and U2 collapsed in a similar form. Figure 1 shows the density evolution of model U2 at selected times after large-scale fragmentation into a binary protostar. A bar connecting two overdense fragments separated by AU is shown in Fig. 1a. Up to this point the results of the simulations are very similar to those reported elsewhere in the literature (Kitsionas & Whitworth, 2002; Arreaga-García et al., 2007; Riaz et al., 2014), except that now the nascent binary protostar is bridged by a centrally condensed, thicker bar. At this time, the binary cores are within the nonisothermal phase of contraction with maximum densities of orders of magnitude higher than the initial density. Spinning of the fragments causes the bar to fan out close to them and develop well-pronounced spiral arms (not shown in Fig. 1), which extend outward for about 350 AU. As the cores accrete low angular momentum mass from the bar and the outer spiral arms, their orbital separation decreases to around 191 AU and the bar dissipates (Figs. 1b-d). As a result of accretion, well-defined protostellar disks of size about 50 AU in diameter have formed around the primaries (disks A and B). After about 250 years (Fig. 1b), the disks had increased in size and developed a strong two-arm spiral pattern. The larger insets in each frame show a magnification of the binary cores and their surrounding disks. Such spiral arms, other than providing the main source of angular momentum transfer, are a clear signature of self-gravitating disks. At this time, the gas in one of the spiral arms of disk B has become locally unstable and started to condense at a distance of 15 AU from the primary, forming a second protostar (Figs. 1b and c). By 44.68 kyr (Fig. 1d), the bar has almost completely dissipated, leaving two linearly aligned condensations of very low mass which move in opposite directions attracted by the larger and more massive protostellar cores. At the time of Fig. 1d, the fragment formed from disk B has already completed a full revolution around the primary.


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefTime variation of Toomre’s parameter along the radial structure of disks A (left) and B (right). The color bar in the middle of both plots indicates the values of . Disk A shows values of and hence it does not fragment. In contrast, is achieved in one of the spiral arms of disk B as highlighted by the undulating yellow region between 10 and 20 AU from the primary. This agrees with a recent revision of the Toomre’s criterion which indicates that in the spiral arms is a necessary condition for fragmentation of the disk (Takahashi et al., 2016). As a result, the disk undergoes eccentric () fragmentation into a second protostar, forming a close binary.

3.1 The Toomre Parameter

The standard reference for estimating the importance of gravitational instability (GI) in young protostellar disks is the parameter , defined by (Toomre, 1964)


where is the sound speed, is the epicyclic frequency (which is equal to the angular velocity in a Keplerian disk), is the gravitational constant, and is the disk surface density. As , self-gravity in the disk becomes increasingly important and when , the disk becomes unstable to the growth of long-wavelength, spiral density waves. A recent revision of the Toomre’s criterion indicates that spiral arm formation occurs when , whereas disk fragmentation due to nonlinear instability of the spiral arms sets in when (Takahashi et al., 2016).

The behavior of self-gravitating disks is complex and depends on a range of disk properties including size, surface density, temperature, and thermal physics. The Toomre parameter at any position within the disk can be calculated using the observationally more convenient expression given by (Kratter & Lodato, 2016)


where is the mass of the central protostar, is the disk mass, is the angular velocity, is the disk scale height, and -4. For , the disk remains stable. As , self-gravity dominates and the disk becomes marginally unstable to the growth of spiral density waves. When , the disk may ultimately fragment due to nonlinear growth of the GI. Figure 2 shows the time evolution of along the radial structure of disks A and B. The dynamics of both disks is dominated by the formation of spiral arms when , whereas fragmentation of disk B occurs at kyr when in one of the spiral arms (see Figs. 1b and c). This result is in very good agreement with a recent revision of the condition for self-gravitational fragmentation of protoplanetary disks (Takahashi et al., 2016). Evidently, the condition defines the time and location within the disk where a GI sets in.

In order to evaluate the parameter plotted in Fig. 2, we first calculate the mass of the central core, , by summing up the mass of all particles lying within a sphere of radius 2 AU around the center of mass of the core. Similarly, the disk mass, , is calculated by summing up the mass of all particles contained within a cylindrical pillbox of radial extent 2 AU 50 AU and width AU centered with the core. The upper frame in Fig. 3 shows the time evolution of and for disks A and B of model U2. Similar plots were obtained for the lower resolution case U1. For comparison, the solid red line depicts the mass of the gravitationally unstable region (i.e., the second fragment) in disk B as a function of time for which according to Fig. 2. This mass grows steeply from kyr to kyr and then approximately linearly (at a much slower rate), reaching a mass of by kyr. As the second protostar grows to a comparable mass to the primary (i.e., against for the primary), the disk shrinks and eventually dissipates (Fig. 1d). During this stage, the secondary orbits around the primary on a very short timescale compared to the stellar lifetime, completing two orbits in less than 1 kyr by the time the simulations were terminated (see Section 3.4 and Fig. 8 below for more details). The bottom frame of Fig. 3 shows the time variation of the ratio for both disks. We note that disk B has , as expected for the growth of the mode (Adams et al., 1989). At kyr, for disk B, whereas at this point disk A is not yet massive enough with values of that are barely above 0.3.

The scale height is calculated as a function of radius by assuming hydrostatic equilibrium of the disk in the vertical -direction. A rectangle of size AU coinciding with the plane of the cylindrical pillbox is chosen on which all particles within the box are projected. The rectangular domain is first partitioned into 50 bins of radial width 1 AU and height 40 AU. The bins are then subdivided into 40 square cells of area 1 AU, where cellwise values of the density are calculated by averaging the contribution from all particles lying within each cell. In this way, density profiles are constructed in the -direction (i.e., along successive bins) for each radius and fitted to a Gaussian distribution of the form (Dong et al., 2016)


to determine and , where is the disk midplane () density. Part of this procedure is illustrated in Fig. 4, where the projected positions of particles on the plane of the pillbox are shown in the left panels for disks A and B. For clarity purposes only three bins bounded by the vertical red lines are depicted at selected radial positions. The panels on the right show the resulting Gaussian-like density distributions along the -coordinate corresponding to the selected bins as indicated by the colored arrows. The maximum values of the distributions for each successive bin define the midplane density, , and the widths of the distributions define the scale height, . After formation of the protostellar disks, the values of , , , and are determined after each timestep and employed in Eq. (2) to calculate the parameter displayed in Fig. 2 as a function of disk radius and time.

3.2 Two-dimensional Emission Maps


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHref(a) Time evolution of the masses of the primary protostars (thin solid blue and dashed green lines) and disks (thick solid blue and dashed green lines) produced in model U2 for systems A and B. (b) Disk mass over central protostar mass ratio, , for disks A and B. The thick solid red line in (a) depicts the mass involved in the gravitationally unstable region within disk B as a function of time. The horizontal dotted line in (b) marks the analytically derived limit above which eccentric () disk fragmentation is expected to occur (Adams et al., 1989). Disk B grows to , while disk A exhibits values of that are barely above the limit for eccentric fragmentation.

Two-dimensional (2D) emission maps of the protostellar disks of Fig. 1 can be generated as follows. As a first step, we calculate the temperature for all particles by combining the barotropic pressure-density relation with the ideal gas law to give


where K, g cm, and . At temperatures below 100 K, a value of is appropriate because the rotational and vibrational degrees of freedom of molecular hydrogen are frozen out, and so only translational degrees of freedom need be considered (Boss et al., 2000). However, precise knowledge of the dependence of temperature on density at the transition from isothermal to nonisothermal collapse will require solving the radiative transfer problem coupled to a fully self-consistent energy equation. This transition was studied for the collapse of an initially centrally condensed core obeying a spherically-symmetric Gaussian density profile, using nonisothermal thermodynamics and solving the mean intensity equation in the Eddington approximation with detailed equations of state (Boss et al., 2000). The collapse was seen to remain strictly isothermal up to g cm. At higher densities the collapse was near isothermal, with the temperature rising very slowly from 10 K to K when g cm. At densities higher than this, the temperature was seen to rise sharply with the density and the collapse became nonisothermal. Therefore, a value of the critical density g cm would be more representative of the near isothermal phase. However, a value of g cm, while prolonging the isothermal phase of collapse to higher densities, makes convergence to be more difficult to attain and resolution requirements to be significantly more demanding. Larger and more massive disks than those calculated here can therefore be obtained by anticipating the near isothermal phase using lower transition densities in barotropic collapse calculations. Figure 5 shows the projection on the disk midplane of the disk particle densities (normalized to the critical density g cm) and temperatures [calculated using Eq. (4)] as functions of radius for disks A and B at kyr. The prominent peaks in density and temperature at represent the central protostars, while those at AU correspond to the secondary fragment produced in disk B. The upper portions of the density peaks are well above the line , implying that the protostars are in a nonisothermal phase of collapse. The small peaks present at about 10 AU from the central protostar in Figs. 5a and c clearly indicates that by this time disk A is in the process of fragmenting into a new secondary as can be better seen from the emission maps of Fig. 6a.

The images displayed in Fig. 6 correspond to disk surface density and midplane temperature maps at selected times during the evolution of disks A and B. These images are generated from the numerical data by first drawing a cubic box of sides 56 AU that encloses the entire disk, where the geometrical center () and the plane of the box are made to coincide with the center-of-mass of the primary core and the aproximate equatorial plane of the disk, respectively. The box is divided into a number of parallel planes above and below the plane, with separations of 0.35 AU between each other. Each successive plane is in turn partitioned into square pixels, with each pixel corresponding effectively to an area of AU. This results in a box composed of regular cubic cells of sides AU. A value of the density is assigned to each cell center by simply averaging the density of all those particles within a cell. This procedure permits drawing density profiles from cellwise densities along the -direction over the entire volume of the box, which are then fitted to a Gaussian distribution to determine the disk midplane density, , and scale height, (see Section 3.1). The disk surface density (in g cm) is then calculated from the estimates of using the expression , where is the gas volume density (Dong et al., 2016). The disk midplane temperature is finally calculated using Eq. (4) with replaced by as determined from the maximum of the Gaussian -density profiles.


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHref(a) Projected particle positions on the plane of a cylindrical pillbox of radius 50 AU and height 40 AU enclosing disk A. The vertical red lines binds three bins of radial width 1 AU at selected radii from the disk center. (b), (c), and (d) Gaussian-like density distributions corresponding to the bins in (a) as indicated by the colored arrows. (e) The same of (a) for disk B. (f), (g), and (h) The same of (b), (c), and (d) for disk B. The maximum and width of the distributions define the disk midplane density and scale height used in Eq. (3), respectively.


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefParticle density distribution as a function of radius for (a) disk A and (b) disk B. The horizontal solid line marks the transition from isothermal to nonisothermal collapse. Also shown is the particle temperature distribution as a function of radius for (c) disk A and (d) disk B. The time shown in all plots is kyr. The density and temperature peaks centered at represent the primary protostars, while those centered at AU correspond to the second fragment produced in disk B.

The radiative flux emerging from the disk surface is finally calculated using the following relation (Dong et al., 2016):


where is the Stefan-Boltzmann constant, is the disk midplane temperature, and is the gas surface density from the midplane of the disk to its surface. Equation (5) relates the emergent radiative flux from the disk surface to the disk midplane temperature. Here and are, respectively, the Planck and Rosseland optical depths from the midplane to the disk surface, with and being the Planck and Rosseland mean opacities, respectively. The Planck and Rosseland mean opacities are calculated from the opacity tables reported by Semenov et al. (2003). We may see from Fig. 6 that disk A, which was less prone to fragmentation, shows clear signs of being forming a second protostar at a distance of around 10 AU from the primary by kyr. This feature complies with theoretical expectations of the role of disks in the formation of close stellar systems (Kratter & Lodato, 2016).


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHref(a) Gas surface density (top row) and midplane temperature (bottom row) of disk A viewed face-on at four different times. (b) The same of (a) for disk B. By kyr, when the calculation was terminated, disk A shows signs of fragmentation into a second protostar. The logarithm of the surface density (in g cm) and the temperature (in K) are plotted. Each row of images uses the color bars on their right sides. In all images the central protostar is drawn with a black (surface density) and dark blue (temperature) circle to increase the surface density and temperature contrasts.

3.3 ALMA Images


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefALMA 1.3 mm continuum images of (a) disk A at 44.97 kyr, (b) disk B at 44.47 kyr, both at a viewing angle of , and (c) disk A at 44.97 kyr, (d) disk B at 44.47 kyr, both at a viewing angle of . The side of the disks nearest the observer is oriented toward the bottom of each frame. The source distance is assumed to be at 100 pc from the observer. The angular resolution of each image is shown with the blue circle drawn in the lower left corner, corresponding to a synthetic beam size of (9 AU at this distance). The simulation images were produced with an integration time of 1 hour. The color bar and numbers on the right indicate the flux intensity (in mJy beam).

Disks with fragment very quickly on orbital timescales (i.e., within kyr). This process reduces the disk mass and stabilizes the system against further activity (Stamatellos et al., 2011). Therefore, the growth of the GI is a short-lived phenomenon that is difficult to observe. In this sense it is useful to simulate ALMA images to help interpret the observations. Figure 7 shows simulated ALMA images of disks A and B at viewing angles of and and at two different times during the evolution. These correspond to full resolution images obtained by transforming the 2D emission maps of Fig. 6 with the aid of the simobserve and simanalyze tools under Common Astronomy Software Application (CASA) for a wavelength of 1.3 mm (230 GHz; ALMA band 6, configuration C43-8 consisting of a full array of 50 antennas of size 12 m). A reference distance of 100 pc from the source was assumed. These images bear a remarkable resemblance to the ALMA images of the Class 0 triple protostar L1448 IRS3 (consisting of a tight binary and a more widely separated tertiary) as observed when the tertiary, which has the largest peak intensity in the system, is removed (Tobin et al., 2016b). However, compared to this protostellar system, the length scale in our results for disk B is smaller by a factor of .

3.4 Orbital evolution of second protostar in disk B

The orbital evolution of the second protostar around the primary in disk B is shown in Fig. 8. The secondary emerges at a distance of AU from the primary protostar at about kyr, reaching a maximum separation of AU by 44.43 kyr. During this stage the fragment accretes mass from the disk at a high rate as shown by the steep increase of , losing orbital angular momentum during this stage and migrating to a minimum separation of AU by kyr when it has already completed half of an orbital period. Thereafter, as the disk dissipates the accretion rate declines and the secondary moves away in its orbit until it reaches a maximum separation of AU at kyr by the time it has completed a full revolution about the primary. The secondary continues in its elliptic orbit and completes a second period by kyr when the calculation is terminated. At this time, the mass of the secondary has become comparable to the primary. During this second orbital period the minimum and maximum distances from the primary are and AU, respectively. This gives an orbital period of years, which corresponds to a much shorter timescale compared to the Class 0 lifetime of kyr (André et al., 2000). The eccentricity was found to decrease from for the first period to for the second period, suggesting that the secondary is approaching a more circular orbit as it comes close to the primary. However, a definite conclusion on the orbital stability of this double-star system (i.e., against mergers) cannot be reached because it will demand pursuing the calculation farther in time up to stellar densities.


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHref(a) Time evolution of the masses of the primary (solid blue line) and secondary (solid magenta line) protostars in disk B. (b) Orbital evolution of the secondary core around the primary protostar. The curve displays the instantaneous position of the secondary after completion of two orbital periods.

4 Discussion and Conclusions

Fragmentation of the kind we have seen in our simulations should be a normal process in the formation of stellar systems. However, at the stage of the evolution shown in Fig. 1d, slightly more than 10% of the gas has been converted into protostars, and therefore the system is very much in its infancy. Beyond the epoch covered by our simulations, heating from the protostars may become strong enough to start dissociation of the H available in the disk and so it will lose its ability to further fragment as most of the coolant is removed. As accretion from the cloud envelope continues, some protostars may become hot enough to initiate ionization of the atomic hydrogen, terminating the accretion. Other than the host star irradiation, disks are also heated by two additional sources: accretion energy and external radiation. However, for deeply embedded cores as may be the case of Class 0 objects, the envelope may insulate the disk from external sources, providing a somewhat lower temperature thermal bath (Kratter & Lodato, 2016).

Radiative hydrodynamics simulations of fragmenting disks were seen to reproduce the constraints set by the observed statistical properties of brown dwarfs and low-mass hydrogen-burning stars (Stamatellos & Whitworth, 2009). More specifically, disk fragmentation can explain the shape of the mass distribution of low-mass stars, the lack of brown dwarfs as close companions to solar-mass stars, the presence of disks around brown dwarfs, the statistics of low-mass binary stars, and the formation of free-floating planetary objects. However, the main argument against the formation of brown dwarfs and low-mass stars by this mechanism is whether massive enough disks are actually realizable in nature. In view of the angular momentum content of their parental cores and the gas accretion from the cloud envelope, the formation of such disks appears to be inevitable. On the other hand, massive () and extended ( AU) disks that do not fragment are expected to dissipate on a viscous timescale (of order Myr) (Stamatellos et al., 2011). Simulations that ignore the effects of magnetic fields find that such disks form frequently in turbulent and/or rotating molecular cloud cores (Attwood et al., 2009). Simulations of protostellar disk formation including turbulence and magnetic fields also show that protostellar disks usually do not extend over more than about 100 AU in size (Myers et al., 2013; Joos et al., 2013; Seifried et al., 2015). However, their rarity suggests that either they form and quickly fragment or their formation is suppressed by the effects of magnetic braking (Mellon & Li, 2009). Alternatively, if typical disk sizes are not larger than AU, their non-detection could be a resolution issue. Using ideal magnetohydrodynamics (MHD), Hennebelle & Teyssier (2008) finds that disk formation is suppressed if the magnetic field is strong enough. On the other hand, if the magnetic field is misaligned, protostellar disks may form but fragmentation is suppressed (Hennebelle & Ciardi, 2009; Commerçon et al., 2010). However, resistive MHD calculations reveal that disk formation is possible because Ohmic dissipation dominates over other processes at relatively high densities ( cm) (Nakano et al., 2002). Further resistive MHD calculations by Machida et al. (2008) and Machida et al. (2010) showed that disk formation and fragmentation after the formation of the first core are possible, while Krasnopolsky et al. (2010) found that the formation of disks with sizes of AU will require resistivities of the order of cm s.

In general, the masses and radii of protostellar disks formed at early stages of the star formation process are much debated because of uncertainties in these properties (Stamatellos et al., 2011). The uncertainties arise because the disks around newly formed protostars (i.e., Class 0) are in almost all cases hidden within the infalling cloud envelope. Therefore, radiative transfer models are needed to distinguish the submillimeter and millimeter emission of the disk from the emission of the envelope. Although our simulations were not evolved much beyond the initial disk fragmentation, comparison with observations indicates that many detected disks around Class I and Class II objects may be remnants of Class 0 disks that have fragmented at a very early stage. Disks around young protostars ( kyr) with masses and sizes consistent with the masses (0.01–0.03) and sizes (20–50 AU) predicted by our simulations have been observed (Rodríguez et al., 2005; Eisner, 2012).

Our results demonstrate that protostellar disks can experience gravitational instability and fragmentation at very young ages, leading to the formation of hierarchical multiple systems as was recently detected by ALMA observations (Tobin et al., 2016b). The advent of mathematically consistent resistive MHD, radiative hydrodynamics models, together with the discovery of larger samples of Class 0 objects with unstable disks using facilities such as ALMA, will help to assess the contribution of disk fragmentation to the production of close binary/multiple stars.

We acknowledge support from the abacus-Centro de Matemática Aplicada y Cómputo de Alto Rendimiento of Cinvestav-IPN under CONACYT grant EDOMEX-2011-C01-165873. The calculations of this paper were performed using the abacus supercomputer facilities.


  • Adams et al. (1989) Adams, F. C., Ruden, S. P., & Shu, F. H.  1989, ApJ, 347, 959
  • Andalib et al. (1997) Andalib, S. W., Tohline, J.W., & Christodoulou, D. M.  1997, ApJS, 108, 471
  • André et al. (2000) André, P., Ward-Thompson, D., & Barsony, M.  2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell (Tucson: University of Arizona Press), 59
  • Arreaga-García et al. (2007) Arreaga-García, G., Klapp, J., Sigalotti, L. Di G., & Gabbasov, R.  2007, ApJ, 666, 290
  • Attwood et al. (2009) Attwood, R. E., Goodwin, S. P., Stamatellos, D., & Whitworth, A. P.  2009, A&A, 495, 201
  • Bate (2018) Bate, M. R.  2018, MNRAS, 475, 5618
  • Bate & Bonnell (1997) Bate, M. R. & Bonnell, I. A.  1997, MNRAS, 285, 33
  • Bonnell & Bate (1994) Bonnell, I. A. & Bate, M. R.  1994, MNRAS, 271, 999
  • Boss et al. (2000) Boss, A. P., Fisher, R. T., Klein, R. I., & McKee, C. F.  2000, ApJ, 528, 325
  • Burkert & Bodenheimer (1993) Burkert, A. & Bodenheimer, P.  1993, MNRAS, 264, 798
  • Commerçon et al. (2010) Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R.  2010, A&A, 510, L3
  • Connelley et al. (2008) Connelley, M. S., Reipurth, B., & Tokunaga, A. T.  2008, AJ, 135, 2526
  • Couchman et al. (1995) Couchman, H. M. P, Thomas, P. A., & Pearce, F. R.  1995, ApJ, 452, 797
  • Dehnen & Aly (2012) Dehnen, W. & Aly, H.  2012, MNRAS, 425, 1068
  • Dong et al. (2016) Dong, R., Vorobyov, E., Pavlyuchenkov,Y., Chiang, E., & Liu, H. B.  2016, ApJ, 823, 141
  • Eisner (2012) Eisner, J. A.  2012, ApJ, 755, 23
  • Fisher (2004) Fisher, R. T. A.  2004, ApJ, 600, 769
  • Gabbasov et al. (2017) Gabbasov. R., Sigalotti, L. Di G., Cruz, F., Klapp, J., & Ramírez-Velasquez, J. M.  2017, ApJ, 835, 287
  • Hennebelle & Ciardi (2009) Hennebelle, P. & Ciardi, A.  2009, A&A, 506, L29
  • Hennebelle & Teyssier (2008) Hennebelle, P. & Teyssier, R.  2008, A&A, 477, 25
  • Hu et al. (2014) Hu, C.-Y., Naab, T., Walch, S., Moster, B. P., & Oser, L.  2014, MNRAS, 443, 1173
  • Joos et al. (2013) Joos, M., Hennebelle, P., Ciardi, A., & Fromang, S.  2013, A&A, 554, A17
  • Kitsionas & Whitworth (2002) Kitsionas, S. & Whitworth, A. P.  2002, MNRAS, 330, 129
  • Krasnopolsky et al. (2010) Krasnopolsky, R., Li, Z.-Y., & Shang, H.  2010, ApJ, 716, 1541
  • Kratter & Lodato (2016) Kratter, K. & Lodato, G.  2016, ARA&A, 54, 271
  • Kraus et al. (2011) Kraus, A. L., Ireland, M. J., Martinache, F., & Hillenbrand, L. A.  2011, ApJ, 731, 8
  • Krumholz et al. (2009) Krumholz, M. R., Klein, R. I., McKee, C. F., Offner, S. S. R., & Cunningham, A. J.  2009, Science, 323, 754
  • Larson (2001) Larson, R. B.  2001, in The Formation of Binary Stars, ed. H. Zinnecker & R. D. Mathieu (IAU Symposium), Vol. 200, 93
  • Machida et al. (2008) Machida, M. N., Tomisaka, K., Matsumoto, T., & Inutsuka, S.-I.  2008, ApJ, 677, 327
  • Machida et al. (2010) Machida, M. N., Inutsuka, S.-I., & Matsumoto, T.  2010, ApJ, 724, 1006
  • Mellon & Li (2009) Mellon, R. R. & Li, Z.-Y.  2009, ApJ, 698, 922
  • Murillo & Lai (2013) Murillo, N. M. & Lai, S.-P.  2013, ApJL, 764, L15
  • Myers et al. (2013) Myers, A. T., McKee, C. F., Cunningham, A. J., Klein, R. I., & Krumholz, M. R.  2013, ApJ, 766, 97
  • Nakano et al. (2002) Nakano, T., Nishi, R., & Umebayashi, T.  2002, ApJ, 573, 199
  • Pech et al. (2010) Pech G., et al.  2010, ApJ, 712, 1403
  • Peters et al. (2010) Peters, T., Klessen, R. S., Mac Low, M.-M., & Banerjee, R.  2010, ApJ, 725, 134
  • Rasio (2000) Rasio, F. A.  2000, Progr. Theoret. Phys. Suppl., 138, 609
  • Read et al. (2010) Read, J. I., Hayfield, T., & Agertz, O.  2010, MNRAS, 405, 1513
  • Riaz et al. (2014) Riaz, R., Farooqui, S. Z., & Vanaverbeke, S.  2014, MNRAS, 444, 1189
  • Rodríguez et al. (2005) Rodríguez, L. F., Loinard, L., D’Alessio, P., Wilner, D. J., & Ho, P. T. P.  2005, ApJL, 621, L133
  • Rosswog (2015) Rosswog, S.  2015, Living Rev. Comput. Astrophys., 1, 1
  • Semenov et al. (2003) Semenov, D., Henning, Th., Helling, Ch., Ilgner, M., & Sedlmayr, E.  2003, A&A, 410, 611
  • Seifried et al. (2015) Seifried, D., Banerjee, R., Pudritz, R. E., & Klessen, R. S.  2015, MNRAS, 446, 2776
  • Springer (2005) Springel, V.  2005, MNRAS, 364, 1105
  • Stamatellos & Whitworth (2009) Stamatellos, D. & Whitworth, A. P.  2009, MNRAS, 392, 413
  • Stamatellos et al. (2011) Stamatellos, D., Maury, A., Whitworth, A., & André, P.  2011, MNRAS, 413, 1787
  • Takahashi et al. (2016) Takahashi, S. Z., Tsukamoto, Y., & Inutsuka, S.  2016, MNRAS, 458, 3597
  • Tobin et al. (2012) Tobin, J. J., et al.  2012, Nature, 492, 83
  • Tobin et al. (2013) Tobin, J. J., et al.  2013, ApJ, 779, 93
  • Tobin et al. (2016a) Tobin, J. J., et al.  2016a, ApJ, 818, 73
  • Tobin et al. (2016b) Tobin, J. J., et al. 2016b, Nature, 538, 483
  • Toomre (1964) Toomre, A.  1964, ApJ, 139, 1217
  • Wendland (1995) Wendland, H.  1995, Adv. Comput. Math., 4, 395
  • Woodward et al. (1994) Woodward, J. W., Tohline J. E., & Hachisu, I.  1994, ApJ, 420, 247
  • Zhu et al. (2015) Zhu, Q., Hernquist, L., & Li, Y.  2015, ApJ, 800, 6
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description