Astron. Astrophys. 319, 547-560 (1997)

3. Results

3.1. Comparison between the FD and the TREESPH evolution

The sequence of equatorial density contours of Fig. 1 shows the time evolution undergone by model FDh. For all of the model calculations the time is measured in units of the initial free-fall time for collapse (, where is the initial reference density). By comparison, Fig. 2 displays a similar sequence of contour plots for the same evolution, as calculated with code TREESPH using particles. An essentially identical sequence was obtained using particles. Since no assumptions about symmetry were made in the TREESPH calculations, the plots of Fig. 2 refer to the -plane containing the particle of maximum density. The contours are then obtained by projecting onto that plane the positions of all particles within the range , where is defined to be a small fraction f of the initial cloud radius. We choose for the calculation and for the lower resolution run .

Fig. 1a-d. Equatorial density contours for the standard isothermal collapse test, as obtained for model FDh. The contour lines correspond to the logarithm of the density. The high-density value is labeled. Counterclockwise rotation is assumed. A sequence in time is presented: a (upper left) , b (upper right) , c (lower left) , and d (lower right) . The entire cloud is shown in a -c. A region of radius cm is shown in d.

Fig. 2a-d. Equatorial density contours for the standard isothermal collapse test, as calculated with code TREESPH using 13997 particles. The contour lines correspond to the logarithm of the density. The high-density value is labeled. Counterclockwise rotation is assumed. A sequence in time is presented: a (upper left) , b (upper right) , c (lower left) , and d (lower right) . The entire cloud is shown in a -c. A region of radius cm is shown in d. An essentially identical sequence was obtained using 7123 particles.

In spite of differences in the numerical methods and set up of the initial conditions, both calculations agree that up to the time they could be compared, the result of the evolution is the formation of a protostellar binary system. The intermediate forms of the collapse are also seen to compare qualitatively well. Due to the rapid rotation assumed for the standard isothermal test case, the cloud collapses to form a rotationally flattened disk about the equatorial plane. Meanwhile, direct fragmentation into a binary system occurs as a result of the gravitational growth of the initial density perturbation. Fig. 3 shows the evolution of the maximum amplitude of the first four even modes in a Fourier expansion of the equatorial density for model FDh. The mode amplifies slowly during the first half of the free-fall time and thereafter it grows more rapidly reaching amplitudes of by . Substantial growth of the higher order modes is seen to occur only after about with maximum amplitudes by the end of the calculation. As expected, the odd modes always kept down at a level of noise with maximum amplitudes .

Fig. 3. Growth of the maximum amplitude of modes in the equatorial plane for the calculation of model FDh. Odd modes kept down at a level of noise with maximum amplitudes .

Contrary to previous first-order calculations, no intermediate bar configuration forms in the present models. This feature was first noted by MB who compared results for the standard isothermal test case from two different second-order FD codes: one based on Cartesian coordinates and the other on spherical coordinates. Their results with the spherical coordinate code, however, indicate that in spite of the second-order accurate methods, an intermediate bar may eventually form if the allowed equatorial resolution is not adequate to solve the rotationally flattened disk forming from collapse. In particular, using a -grid with fairly equidistant points for , the evolution was seen to pass through an intermediate bar configuration as in previous first-order calculations while, for the same , the intermediate bar disappeared when the calculation was repeated using a stretched (non-uniform) -grid about the equatorial plane, thereby allowing finer resolution there at the expense of much coarser resolution near the rotational axis. The non-uniform -grid (with ) used in MB gives an angular spacing, , near the equatorial plane which is a factor of 6 smaller than the corresponding one allowed in the calculation of model FDh (with fairly equidistant points for ). In spite of this difference in equatorial resolution, the contours of Fig. 1 closely resemble those given in MB at comparable maximum densities. This result clearly indicates that the disappearance of the intermediate central bar is a consequence of the smaller diffusive truncation errors associated with code FD, and that enhanced equatorial resolution plays a minor role in determining that feature. This last statement is further confirmed by the results obtained for model FDl, with fairly equidistant points for as in the spherical coordinate model in MB that produced an intermediate bar. The evolution of this model is illustrated in the sequence of contour plots of Fig. 4. Although this calculation was made with lower initial radial and azimuthal resolution (i.e., and ) compared to model FDh and the calculations in MB, we see that the intermediate forms of collapse are very similar to those shown in Fig. 1 and that again no bar configuration is seen to form. The growth of the maximum equatorial amplitude of the first four even modes closely follows that given in Fig. 3, with the mode reaching amplitudes of by the end of the calculation (). Furthermore, comparison of the results from the TREESPH calculations shows that varying the number of particles has also little effect on the intermediate forms of the evolution. In both of these model calculations, the protostellar binary formed by direct fragmentation without the occurrence of an intermediate central bar (see Fig. 2).

Fig. 4a-d. Equatorial density contours for the standard isothermal collapse test, as obtained for model FDl. The contour lines correspond to the logarithm of the density. The high-density value is labeled. Counterclockwise rotation is assumed. A sequence in time is presented: a (upper left) , b (upper right) , c (lower left) , and d (lower right) . The entire cloud is shown in a -c. A region of radius cm is shown in d.

The variation of the maximum density with time is shown in Fig. 5. The results indicate that the maximum density is much less sensitive to resolution than was previously found by ML. We see that the TREESPH sequence with 13997 particles closely follows that with 7123 particles, the former showing slightly higher maximum densities only after the first free-fall time. The effects of resolution are made more evident by intercomparing the FD sequences, where at comparable phase the differences in maximum density are larger. Although the maximum density in model FDh is always above that of the TREESPH sequences, the trends of the time variation are very similar. The major differences between the FD and the TREESPH sequences are seen to occur after about when the evolution of the fragments is dominated by their internal forces. These differences, however, cannot be explained solely in terms of differences in the final resolution achieved by the two methods. Indeed we find that the TREESPH calculations produce forming-fragment regions with a relatively higher rotational energy ( by ) compared with the FD models. Thus the central regions of the TREESPH fragments may become rotationally unstable and bounce, giving rise to a temporary decrease of the maximum density. The excess of rotation at the centre of the FD fragments was less severe but still enough to produce a slow collapse. This explains the much milder increase of the maximum density for the FD sequences after . During this phase, the gravitational energy at the centre of the fragments has sufficiently increased to induce a more rapid collapse, as evidenced by the further sudden rise of the maximum density. A similar behaviour is also observed to occur in the evolution of the TREESPH fragments but on a relatively longer time-scale. This is because the fragments must dispose of their excess of rotational energy before trying to collapse upon themselves. Transfer of angular momentum to the spiral arms is the likely mechanism through which the TREESPH fragments may eventually reduce the importance of the rotational energy and collapse to higher densities.

Fig. 5. The variation of maximum density with time, as calculated with code FD (curve 1: model FDh; curve 2: model FDl) and with code TREESPH using 13997 (curve 3) and 7123 (curve 4) particles.

Another interesting question regarding the standard isothermal collapse test case is whether or not the evolution will result in a run-away collapse of the fragments. The results with code TREESPH show that the fragments may evolve farther in time without undergoing rapid collapse upon themselves, in good agreement with ML calculations. This aspect was in contrast with the results from early first-order FD calculations, where as suggested by ML, run-away collapse occurred because the description of the fragments could have been degraded by the excessive numerical diffusion of angular momentum from them. However, the results with code FD indicate that in spite of the decreased numerical diffusion associated with the second-order methods, the description of the fragments can also critically depend on spatial resolution. This aspect is made evident in Fig. 5, where in the case of model FDh the maximum density suddenly stops increasing at and thereafter starts oscillating as in the TREESPH calculations, while at comparable values, the maximum density in the sequence of model FDl continues increasing, indicating a run-away collapse of the fragments as in all previous first-order calculations. Unfortunately, the calculation of model FDh was halted shortly after (at ) because the time steps became extremely small. On the contrary, the calculation of model FDl was continued up to when . By this time, however, the accuracy of the evolution has already been degraded due to the lack of resolution as the fragments rapidly concentrated in a few grid cells. The TREESPH calculations were terminated at to allow comparison with ML results for the long term evolution.

The density contours of Fig. 1d (at ) and 2d (at ) show the protostellar binary resulting from the high-resolution calculations with code FD and code TREESPH, respectively. The maximum density in Fig. 1d is and that in Fig. 2d is . The general shape of the density contours in Fig. 1d is very similar to that given in MB. Because of the higher azimuthal resolution used here ( against in MB and for model FDl), the binary fragments evolved on a slightly shorter time-scale. The fragments obtained with code TREESPH appear to have an elongated shape and resemble the ML fragments more closely than those found with the FD codes. The presence of well-defined spiral arms just outside the main fragments is also evident in Fig. 2d. These features are also present in Fig. 1d but in a less prominent fashion. We further note that despite these differences, the time-scale and the maximum densities for the TREESPH evolution closely agree with those given in MB using their Cartesian FD code.

The final shape of the fragments may be influenced by the accretion of matter from the outer cloud regions. The accretion process is clarified in Figs. 6 and 7, which show the velocity field in the equatorial plane at the times of Figs. 1d and 2d, respectively. We see that the accretion is highly anisotropic with strong pressure gradients forming preferentially in the azimuthal direction. These gradients involve a sharp transition from supersonic to subsonic velocities as shown more clearly in Fig. 6. The matter which approaches the fragments azimuthally, i.e., in the sense of the general cloud rotation (downstream matter), loses angular momentum due to the gravitational attraction exerted by the fragment and then swings inwards. Most peripheral mass with high angular momentum is seen to swing further outwards and then to approach the fragment on curved trajectories (upstream matter). While these features are essentially the same in Figs. 6 and 7, there is still a basic difference between them - accretion of mass onto the TREESPH fragments occurs preferentially through the spiral arms. This would explain their more elongated shape, which is also evident in Fig. 7 as shown by the regions where the velocity arrows are heavily superimposed. Furthermore, the fragments produced here look quite flattened about the equatorial plane and surrounded by a planar shock through which they accrete low angular momentum mass in the z -direction.

Fig. 6. Velocity field in the equatorial plane at , for model FDh. The rotation axis is in the centre of the plot. The numbers on the box sides are given cm. The maximum velocity is 1.55 km s-1.

Fig. 7. Velocity field in the equatorial plane at , as calculated with code TREESPH using 13997 particles. The numbers on the box sides are given cm. The maximum velocity is 1.26 km s-1. Only particles within a fraction of the initial cloud radius above and below the equator are shown.

3.2. Integral properties of the fragments

In numerical calculations of protostellar collapse and fragmentation we are always concerned with the properties of the resulting fragments. Here the fragment region is defined as the cell (or particle) of maximum density, plus all surrounding cells (or particles) with densities higher than a factor of the maximum value. For simplicity, in the sequel we shall refer to these cells (FD case) and particles (TREESPH case) as the elements making up the fragment. The total mass of the fragment is then determined by summing up the masses of the individual elements. The internal velocity field of the fragment, say , is here defined with respect to a spherical coordinate system whose origin coincides with the instantaneous location of the element of maximum density (centre of the fragment). In this frame, the fragment's rotation axis is assumed to be perpendicular to the equatorial plane of the cloud. These velocities should not be confused with those describing the internal motion of the cloud with respect to its centre of mass (see Figs. 6 and 7). In the FD calculations, was always seen to coincide with the origin of the grid coordinate system . In the TREESPH calculations, equatorial symmetry is well maintained and is found to lie approximately in the -plane containing the particle of maximum density, with during the evolution. Therefore we can always define a spherical coordinate system with origin at and with respect to which the actual particle velocities are assumed to be defined. Thus the location of an element relative to the fragment's centre can be determined from its location with respect to the cloud's centre by means of the transformations

where , and are the coordinates of the element of maximum density. Differentiating with time the above relations gives the velocity of the fragment element in terms of its components and the velocity of the element of maximum density as measured in the main cloud frame. Let now consider the equatorial plane of the cloud, as shown in Fig. 6, and assume that the binary fragment produced is such that its centre (point of maximum density) has . In this case, we recover the identity , and if we further assume that the elements of the fragment are only undergoing rotational motion about its axis of spin , then and . We see from Fig. 6 that this is not the case and that deriving as described above is necessary to improve on the accuracy of the rotational properties. The gravitational, thermal, rotational, and kinetic non-rotational energies as well as the spin angular momentum of the fragment are then evaluated by summing the volume-element integrated values of , , , , and , respectively, over all elements composing the fragments. This determines the values of , , and which denote, respectively, the ratios of the thermal, rotational, and kinetic non-rotational to the gravitational energy for the clumps. Finally, the spin specific angular momentum of the fragment is found by dividing its spin angular momentum by its total mass .

We have calculated the properties of the binary fragments for two different choices of the factor . In Table 1, we use and compare our results (S32 and S80: code FD; K7123 and K13997: code TREESPH) with those from previous first-order donor-cell calculations by Bodenheimer & Boss (1981) and de Felice & Sigalotti 1992 (B32: Boss' code with ; T32: Tohline's code with ; FS32: Sigalotti's code with ), second-order FD calculations by MB (M64: Myhill's Cartesian code with 64 points; B64: Boss' spherical code with ), and most recent SPH calculations by Bate et al. (1995) using a version of Willy Benz's code with 8024 particles (BBPi: calculation made using an isothermal equation of state; BBPp: calculation made using a polytropic equation of state with for g cm^-3). In Table 2, we re-calculate the integral properties by setting according to the actual shape of the fragments as given in Fig. 8. The last two columns in Table 1 and rows in Table 2 list the binary separation in terms of the initial cloud radius and the time at which the fragment properties were evaluated. The specific angular momentum of the initial cloud is compared with that of the resulting fragments to determine the reduction factor of the angular momentum. Since the initial conditions produced symmetric binaries, Tables 1 and 2 show the properties of only one fragment.

Table 1. Fragment properties for : comparison with previous calculations

Table 2. Integral properties of the actual fragments

Fig. 8a and b. Equatorial density contours showing the binary fragment region: a (left) as in Fig. 1d at (code FD), b (right) as in Fig. 2d at (code TREESPH). A region of radius cm is shown in a and b. The high- and low-density values are labeled. The fragmentary appearance of the fragments in b is only a transient feature.

3.3. Discussion of the results

We note from Table 1 that the second-order calculations resulted in larger binary separations compared with the first-order ones. This is expected because in all of the second-order calculations the cloud is seen to fragment directly into a binary system. On the contrary, the first-order calculations were largely affected by numerical diffusion which acted primarily to inhibit direct fragmentation. In this case, the initial density perturbation was seen to decay into a central bar configuration. Subsequent collapse and fragmentation of the bar then resulted in a binary system of separation smaller than that found by direct fragmentation. Numerical diffusion of angular momentum from the growing fragments could also explain the lower final values of along with the increased efficiency of angular momentum redistribution compared with the second-order calculations.

Although Table 1 provides a direct comparison of the fragment properties with previous calculations of the standard isothermal collapse test case, we find that our results can be compared more closely using Table 2, where the properties account for the actual shape of the fragments (see Fig. 8). Considering the differences between the numerical methods, we see that there is a general good agreement between the FD and the TREESPH results. The fragments which form have low values of () and relatively high values of (). Significant non-rotational motion is also present in their internal structure as suggested by the moderate values of (). Note that the fragments produced in model FDl have much lower values of () and higher (), which is consistent with these clumps being undergoing run-away collapse. All four of our model calculations agree that about 20 percent of the total cloud mass is in the form of fragments, and that their spin angular momentum is not reduced by factors larger than with respect to the initial cloud value.

The fragments are seen to form at such low values of that they could experience rapid collapse and eventually sub-fragment in the further evolution. A rapid collapse of the fragments would be accompanied by an equally rapid rise of the maximum density. With the exception of model FDl, our results indicate that the fragments undergo such a rapid collapse until a certain value of the maximum density is reached. Thereafter the collapse slows down and the maximum density may start oscillating. Since we were unable to follow the calculation of model FDh beyond , this trend is barely shown in Fig. 5. Previous SPH calculations by Gingold & Monaghan (1981) and ML, and more recently by Bate et al. (1995), have shown that the long term evolution could be controlled by the gravitational interaction between the two fragments. In ML a clear tendency for the fragments to get closer together was observed to occur up to when they apparently started to coalesce. A tendency for the fragments to fall towards each other is also seen to occur in our calculations as shown in Fig. 9. In this figure, we plot the time variation of the binary separation (i.e., one-half of the distance between the points of maximum density divided by the initial cloud radius) up to the time the calculations were compared. Although the binary separation is always nearly larger for the FD models, it is seen to decrease at a rate similar to that found in the TREESPH calculations.

Fig. 9. Variation of the binary separation distance with time, as calculated with code FD (circles: model FDl; plus signs: model FDh) and code TREESPH (asterisks: model with 7123 particles; crosses: model with 13997 particles). The binary separation is measured as one-half of the distance between the points of maximum density divided by the initial cloud radius.

The long term evolution was also studied by Bate et al. (1995) who followed the calculations well beyond the point where ML estimated the fragments to start coalescing. Following the initial fragmentation, the growing fragments were also seen to get progressively closer together. However, by when the two fragments reach a minimum separation distance, the binary system have had time to acquire enough orbital angular momentum to prevent merging and force the fragments to pass each other at periastron. Thereafter the clumps continue to orbit and the disc surrounding one of them fragments into one more object, yielding a triple protostellar system by . The sequence of contour plots in Fig. 10 shows the long term evolution of the binary system up to 1.80 , as calculated with code TREESPH using 13997 initial particles. A similar sequence was also obtained using 7123 particles. As the two clumps approach each other, coalescence at the log contour level is seen to occur by (Fig. 10b). At this time, the binary separation is a factor of 2 smaller than the value given in Table 2 (Fig. 8b). This trend continues up to as shown in Fig. 10c and d, where coalescence at the four outermost contour levels is evident. Although these results are qualitatively in agreement with the general trends described in ML, we note that by , the fragments are still distinct and seem to approach a point of minimum separation without merging, as described in Bate et al. (1995). The question of whether the fragments will definitively merge or survive eventually sub-fragmenting, will require following the evolution farther in time than we were able to do here. In any case it seems likely that the protostellar binary is merely a transient feature.

Fig. 10a-d. Equatorial density contours showing the long term evolution of the binary system, as calculated with code TREESPH using 13997 particles. The contour lines correspond to the logarithm of the density. The high- and low-density values are labeled. Counterclockwise rotation is assumed. A sequence in time is presented: a (upper left) , b (upper right) , c (lower left) , and d (lower right) . A region of radius cm is shown. The binary separation distance in terms of is 0.22, 0.16, 0.12, and 0.09 in a, b, c, and d, respectively.

In ML the fragments were estimated to be in approximate equilibrium after in the long term evolution. We find from Table 2 that at the times when the properties were evaluated, for the fragments formed in model FDh, and and 0.47 for those formed in the TREESPH calculations with 7123 and 13997 particles, respectively. While virial equilibrium requires that , it may well be the case that by these times the outer regions of the fragments are already close to equilibrium and that only a small amount of matter at their centres will continue to collapse. This view is clarified in Fig. 11, where the variation of the maximum density with time is shown up to for both of our TREESPH calculations. We see that the centre of each fragment undergoes further collapse reaching a maximum density peak by . Thereafter the fragments bounce at their centres and the maximum density begins to oscillate about an equilibrium configuration. During this phase, the evolution of the fragments is clearly influenced by the strong gravitational interaction between them. After , the oscillations of the maximum density produce peaks that are slightly larger in the calculation with 13997 particles. However, later on (by ), the oscillations in both sequences start to converge more closely, indicating that the results with respect to run-away collapse of the fragments are reaching the correct continuum limit.

Fig. 11. The variation of maximum density with time as calculated with code TREESPH. Curve 1 refers to the sequence with 13997 particles and curve 2 to that with 7123 particles. The long term evolution is shown up to .

© European Southern Observatory (ESO) 1997

Online publication: July 3, 1998
helpdesk@link.springer.de