2.1. The Core and the Initial Setup of the Particles

We consider the gravitational collapse of a variant of the so-called “standard isothermal test case” which was first calculated by ^{Boss & Bodenheimer (1979)} and later calculated by ^{Burkert & Bodenheimer (1993)} and ^{Bate & Burkert (1997)}.

In this paper, the core radius is R₀ = 4.99 × 10¹⁶ cm ≡ 0.016 pc ≡ 3335 AU and its mass is M₀ = 5M_⊙. Thus, the average density and the corresponding free fall time of this core are ρ₀ = 1.91×10⁻¹⁷ gcm⁻³ and t_ff ≈ 4.8×10¹¹ s≡ 15244yr, respectively.

We set N_p particles on a rectangular mesh of side length 2 R₀ (the core diameter), to represent the initial core configuration. We then partition the simulation volume into small elements, each with a side length given by ∆x = ∆y = ∆z = R₀ /N_xyz where N_xyz is an integer given by 133, so that there are 133³ grid elements in the entire volume; at the center of each small grid we place a particle, which is next displaced from its initial location a distance of the order ∆x/4.0 in a random spatial direction. The boundary conditions applied to the simulation box, both in the initial setup and during the time evolution are named vacuum boundary conditions ( so that we have a constant volume with non-periodic boundary conditions).

It should be noted that we here consider only a uniform density core, so that its average density is given by ρ₀ or a number density of n₀ = 9.510⁶ particles/cm³. This is a very dense core, whose density is 10 times higher than that of the “standard isothermal test case.”

Thus, for all the simulations, the particles have the same mass, according to m_i = ρ₀ × ∆x∆y ∆z with i=1,...,N_p, where N_p=996 040. Then the mass is given by m_p=5.0 ×10⁻⁶ M_⊙.

2.2. Evolution Code

The gravitational collapse of our models has been followed by using the fully parallelized particlebased code Gadget2; see ^{Springel (2005)} and also ^{Springel et al. (2001)}. Gadget2 is based on the tree − PM method for computing the gravitational forces and on the standard smoothed particle hydrodynamics (SPH) method for solving the Euler equations of hydrodynamics. Gadget2 implements a Monaghan-Balsara form for the artificial viscosity; see ^{Monaghan & Gingold (1983)} and ^{Balsara (1995)}. The strength of the viscosity is regulated by the parameter αν=0.75 and βν=12×αv; see equations 11 and 14 in ^{Springel (2005)}. In our simulations we fixed the Courant factor to be 0.1.

2.3. Resolution and Thermodynamics Considerations

Following ^{Truelove et al. (1997)} and ^{Bate & Burkert (1997)}, in order to avoid artificial fragmentation, any code for collapse calculation must have a minimum resolution given in terms of the Jeans wavelength λ_J:

where c is the instantaneous speed of sound and ρ is the local density. To obtain a more useful form for a particle based code, the Jeans wavelength λ_J can be transformed into a Jeans mass, given by

The values of the density and speed of sound must be updated according to the following equation of state:

which was proposed by ^{Boss et al. (2000)}, where γ ≡ 5/3, and for the critical density we assume the value ρ_crit = 5.0 × 10⁻¹⁴ grcm⁻³.

The smallest mass particle that a SPH calculation must resolve in order to be reliable is expressed in terms of the particle mass, m_r, given by m_r ≈ M_J/(2N_neigh), where N_neigh is the number of neighboring particles included in the SPH kernel; see ^{Bate & Burkert (1997)}. Hence, a simulation satisfying all these requirements must satisfy mpmr<1

Thus, for the turbulent core under consideration we have M _j ≈ 1.1 x 10^-3 _⊙ and m _r ≈ 1.4 x 10^-5M_⊙, since we took N_neigh= 40. The ratio of masses is given by m_p/m_r=0.33 and then the desired resolution is achieved in our simulations.

2.4. The Turbulent Velocity of the Particles

To generate the turbulent velocity spectrum, we set up a second mesh, in which the partition is determined by the values of N_x, N_y and N_z, which will initially be given by 64 and later on they will be increased to 128, so that the total number of grid elements will change from 64³ to 128³. It should be noted that the velocity vector of a particle changes with the number of modes.

Let us denote the side length of this second mesh by L₀, so that it is proportional to the core radius R₀

where C_R is a constant, the value of which will also determine the collapse model under consideration, see Table 1.

Thus, the size of each grid element is given by δx = L₀/N_x, δy = L₀/N_y and δz = L₀/N_z. In Fourier space the partition is given by δK_x = 1.0/(N_x × δ_x) , δK_y = 1.0/(N_y × δ_y) and δK_z = 1.0/(N_z × δ_z).

Each Fourier mode has the components K_x = i_KxδK_x, K_y = i_Ky δK_y and K_z = i_Kz δK_z, where the indices i_Kx, i_Ky and i_Kz take values in [−N_x/2,N_x/2], [−N_y/2,N_y/2] and [−N_z/2,N_z/2], respectively. The magnitude of the wave number is K=Kx2+Ky2+Kz2 . Then K_max is proportional to 3Nx2L0 . A Fourier mode can equally be described by a wavelength λ = 2π/K. Then we see that

We will consider two types of turbulence, so that we will be able to compare how the nature of the initial turbulent spectrum affects the core collapse. The initial power of the velocity field will be given by:

where the spectral index is a constant.

Let us first consider the Fourier transform F- of the velocity field, so that v→r→=Fv→K→ or explicitly:

It can then be shown that the velocity v→K→ can be written as:

where the scalar and vector potential functions are given by:

so that their Fourier transforms are Φr→=FΦK→ and A→r→=FA→K→, respectively. According to the Helmholtz decomposition theorem, the velocity field in physical space will then be determined in general by:

2.4.1. Divergence-free Turbulence

Let us consider the case Φ( K→) ≡ 0, which implies that K→∙v→K→=0. The velocity field can then be written in terms of a vector potential A→ alone. According to equation 10, ∇∙v→r→ = 0; this is the divergence free turbulence spectrum, whose Fourier transform can be approximated by the following discrete summation:

In order to have the power described in equation 6, we choose the vector potential A→K→ given by:

where is C→ a vector whose components are denoted by (CKx,CKy,CKz), and take values obeying a Rayleigh distribution. Hence, the magnitude of C→ is calculated by means of the formula C=σ×-2.0×log1.0-u, where u is a random number in the interval (0,1) and σ = 1.0 is a fixed parameter. There must be one phase Φ_K for each component of the vector C→ so that each of (ΦKx,ΦKy,ΦKz) takes random values in the interval [0,2π].

Thus, following ^{Dobbs et al. (2005)}, the components of the particle velocity are given by the above equation 13.

2.4.2. Curl-free Turbulence

On the other hand, let us consider the case in which A→K→ = 0, so that K→× v→K→= 0. According to equation 10, this velocity field satisfies ∇× v→= 0 and for this reason we call it the curl-free turbulence spectrum.

In order to have the power spectrum shown in equation 6, we choose the scalar potential function given by:

where there is now only one wave phase random function Φ_K that takes random values in the interval [0,2π]. Therefore, the velocity field will be determined by:

where the spectral indexwill be fixed in our simulations to be n = −1 and thus we will have v² ≈ K−³; see ^{Arreaga & Klapp (2014)}.

In the two types of turbulence spectra, the SPH particles have initially a Gaussian distribution of velocity. Subsequently, by using all the simulation particles, we obtain that the average Mach number is, remarkably, almost the same for all our simulations, M=1.5. The velocity dispersion is also very similar for the two types of spectrum: σ_v= 0.21 km/s. Following the definitions of skewness (or third moment) and kurtosis (or fourth moment) of a given distribution, see ^{Press et al. (1992)}, we observe differences in the values of the skewness and kurtosis of the velocity distribution spectra: 0.4 and 0.003 for the divergence-free type and 0.38 and 0.03 for the curl-free type, respectively.

2.5. Initial Energies

It is well known that the global dynamical evolution of the core is determined by the ratio of the thermal energy to the gravitational energy, denoted by α and the ratio of the kinetic energy to the gravitational energy, denoted by β; see, for instance, ^{Miyama et al. (1984)}, ^{Hachisu & Heriguchi (1984)} and ^{Hachisu & Heriguchi (1985)}, who obtained a criterion of the type αβ < 0.2 to predict the collapse and fragmentation of a rotating isothermal core.

The gravitational and kinetic energies can be approximated in terms of the physical parameters of the core considered in this paper:

where G is Newton’s universal gravitation constant, M₀ is the mass, R₀ is the radius, and 〈V〉 the average velocity. The thermal energy E_therm (kinetic plus potential interaction terms of the molecules) can be approximated by:

where N is the total number of molecules in the gas, k_B is the Boltzmann constant, and T is the temperature of the core.

These energies can be calculated in terms of the SPH particles as follows:

where p _i is the pressure and Φ_i are the values of the pressure and gravitational potential at the location of particle i, with velocity given by v_i and mass m_i; the summations include all the simulation particles.

Now, the values of the speed of sound c₀ and the level of turbulence, which is adjusted by multiplying the velocity vector by an appropriate constant, are chosen so that the velocity field fulfills the following energy requirements:

where the energies entering in these ratios have been calculated using the relations of equation 18.

The virial theorem states that if a cloud is in thermodynamical equilibrium, then the dimensionless energy ratios satisfy the following relation:

which will be used in a plot of the next section. It should be noted that we do not include a turbulence term, γ_turb, as all the turbulent energy has already entered in the kinetic ratio, β.

In order to compare our simulations with other papers, we consider Equation 1 of ^{Walch et al. (2009)}, in which the turbulent energy of a core is approximated by the sum of two terms: the first one contains the FWHM velocity width reported by ^{André (2007)} for the diazenylium molecule N₂H⁺, and the second one depends on the thermal energy k_B T.

As we mentioned at the end of § 2.4, our simulations have an average Mach number of 1.5, so by using our value of c₀=35 820 cm/s, we get that the average velocity 〈v〉 = 0.53 km/s. However, an estimate of the velocity dispersion used in observations, denoted here by σve, can be obtained from the empirical relation Ekin=M0σve2∕2, and using equation 16 and the second relation of equation 19, we obtain σve=β6GM05R0=0.57 km/s, which is of the same order as 〈v〉, as expected. By using a mass of 29 uma for the diazenylium molecule and the value of σve in equation 1 of ^{Walch et al. (2009)}, we obtain for our simulations a ratio of the turbulent energy to the gravitational energy of γ_turb = 0.11. In the simulations reported by ^{Walch et al. (2009)}, the ratio γ_turb = 0.010 was used.

2.6. The Models

The models considered in this paper are summarized in Table 1, whose entries are as follows. Column 1 shows the model number; Column 2 shows the value of the size constant C_R as defined in equation 4; Column 3 shows the partition used in the Fourier mesh, so that the total number of grid elements of the cubic mesh is the cube of that number; we mentioned in § 2.5 that the average Mach number of the simulations are very similar; however, less than 1 percent of the simulation particles can attain very high velocities: Column 4 shows the peak Mach number of the simulations at the initial time (the snapshot zero); Column 5 gives the peak velocity found in the collapsed region where the resulting protostar is formed; Column 6 shows the number of the figure in which the model configuration is shown; Column 7 shows the type of turbulence used to generate the initial velocity spectrum where the label D-F means divergence-free and C-F means curl-free; Column 8 shows the resulting configuration, where the label S means a single protostar and B means a binary system of protostars.

In order to further characterize our models, we calculate the total angular momentum of each model by using all the simulation particles, so that L→=∑i=1Npr→×p→, with the linear momentum p~ given by p→=mpv→, where m_p is the mass of the particle and v→ its velocity; see § 2.1. Figure 1 shows the specific angular momentum L/M₀ (the total angular momentum | L→ | of the initial distribution of particles divided by the core mass) against the model number. By looking at Figure 1, one can see that (i), when the Fourier mesh size increases in terms of the core diameter, so does the | L→ |; and (ii), when the number of Fourier modes increases from 64³ to 128³, | L→ | decreases. The three horizontal lines in Figure 1 indicate the observed values for real cores with radii varying within 1-5 ×10¹⁷ cm; see ^{Goodman et al. (1993)}; ^{Bodenheimer (1995)}. In fact, the observed values of L/M₀ are higher than those of our models. The two vertical lines in Figure 1 indicate the models where we will observe binary formation; see § 3.

Fig. 1 The initial specific angular momentum L/M0 against the model number. The horizontal lines indicate the observed values of L/M0 for real cores; see ^{Goodman et al. (1993)}; ^{Bodenheimer (1995)}. 

3.1. Models 1-4

By looking at the left and middle panels of Figures 4, 5, 6, and 7, which correspond to models 1, 2, 3, and 4, we note that as the wavelength of the initial perturbation increases, so do the sizes of the initial over-density clumps, and the number of clumps that form across the core decreases.

Fig. 4 Iso-density plots for model 1. The color figure can be viewed online. 

Fig. 5 Iso-density plots for model 2. The color figure can be viewed online. 

Fig. 6 Iso-density plots for model 3. The color figure can be viewed online. 

Fig. 7 Iso-density plots for model 4. The color figure can be viewed online. 

Due to the loss of homogeneity in the initial distribution of the over-dense clumps, the protostar in formation accretes mass from the surroundings in an anisotropic manner; as a consequence, the protostar moves slightly from the center, in the south and southwest directions, as can be seen in the right panels of Figures 6 and 7, which correspond to models 3 and 4.

Figure 2 shows that for models 1-4, the larger the initial perturbation size, the earlier the local peak in the density curve occurs. Besides, we note that the peak velocity in the collapsing region decreases as the size of the initial perturbation increases; see Column 5 of Table 1.

3.2. Models 5-8

Here we re-calculate models 1-4, keeping all their parameters unchanged except for the number of grid elements of the Fourier mesh; see Column 3 of Table 1. We first note that the number of small overdensity clumps formed at the initial snapshots increases, but they are still homogeneously distributed across the core volume; see the left panels of Figures 8 and 9 and compare them with those of models 2-3.

Fig. 8 Iso-density plots for model 5. The color figure can be viewed online. 

Fig. 9 Iso-density plots for model 6. The color figure can be viewed online. 

The loss of homogeneity in the initial distribution of over-density clumps is noticeable in model 7, as can be seen in Figure 10. A significant reduction in the number of the initial over-density clumps can also be seen in the initial snapshot of model 8; see the left panel of Figure 10. However, by comparing this panel with that of model 4, shown in Figure 7, we see that there are slightly more over-density clumps than those formed in model 4.

Fig. 10 Iso-density plots for model 7. The color figure can be viewed online. 

3.3. Models 9-11

In these models, the wavelength of the initial perturbation far exceeds the core radius; see Column 2 of Table 1. It can be seen in the left panel of Figure 12 that only two over-density clumps are clearly visible in the equatorial slice of model 9. These perturbations are coupled to deform the central core simultaneously in two orthogonal directions, see the panel in the middle of Figure 12. By looking at the right panel of Figure 12, we see the appearance of small spiral arms around the protostar, indicating that the protostar has gained a net angular momentum as a consequence of the small number of coupled perturbations that acted upon it initially. For model 10, we see that only one over-density is visible at the initial snapshot shown in the left panel of Figure 13. For this reason, the still forming protostar can be seen to be highly deformed, mainly along a single direction; see the left and middle panels of Figure 13. Thus, the right panel of Figure 13 shows that the still forming protostar drags a long tail of gas along this single direction. The same behavior is mainly observed for model 11; see Figure 14.

Fig. 11 Iso-density plots for model 8. The color figure can be viewed online. 

Fig. 12 Iso-density plots for model 9. The color figure can be viewed online. 

Fig. 13 Iso-density plots for model 10. The color figure can be viewed online. 

Fig. 14 Iso-density plots for model 11. The color figure can be viewed online. 

3.4. Models 12-14

For these models we use more Fourier modes to re-calculate the models of § 3.3. Let us examine model 12 shown in Figure 15 and compare it with model 9 shown in Figure 12. We first note that there is only one perturbation mode visible in the left panel of Figure 15; despite this, there are a few dominant directions along which the initial over-density grows; see the middle panel of Figure 15. The velocity of the collapsing peak is significantly greater in model 12 than that observed in model 9; see Column 5 of Table 1. Because of this, it is very likely that the protostar formed in model 12 has gained a net angular momentum larger than that of the protostar of model 9, as is suggested by looking at the right panel of Figure 15.

Fig. 15 Iso-density plots for model 12. The color figure can be viewed online. 

The size of the initial over-density clump induced by the perturbation mode of model 13 looks smaller than that induced in model 10; compare Figure 16 with Figure 13.

Fig. 16 Iso-density plots for model 13. The color figure can be viewed online 

In model 14 we observe the occurrence of fragmentation of the central core region, so that a binary protostar system is formed.

3.5. Models 15-18

We consider now the set of models generated with the curl-free turbulence spectrum. In the left panel of Figure 18, we see that the number of over-density clumps formed initially in model 15 is smaller than the number observed in model 1, shown in the left panel of Figure 4. The clumps seen in Figure 18 are also in fact larger that those seen in Figure 4. Despite this, we still see that the curl-free turbulence spectrum creates a homogeneous distribution of over-density clumps. At time 0.67t_ff, when the strong collapse is yet to begin, the spherical symmetry of the collapsing core is already visible in the middle panel of Figure 18. For this reason, we see that the resulting protostar accretes gas from the remaining core almost isotropically; see the right panel of Figure 18.

Fig. 17 Iso-density plots for model 14. The color figure can be viewed online. 

Fig. 18 Iso-density plots for model 15. The color figure can be viewed online. 

When the size of the Fourier mesh begins to increase in terms of the core radius, which is shown in models 16-18, the initial over-density clumps are formed thicker and their number is smaller than those seen in model 15; see the left panels of Figures 19, 20 and 21. As a consequence, the homogeneity and isotropy of the over-density clump distribution is slightly lessened, and then some clumps grow forming large filaments, as can be observed in the right panels of Figures 20 and 21. All these models finish with a single protostar formed in the collapsed central region.

Fig. 19 Iso-density plots for model 16. The color figure can be viewed online. 

Fig. 20 Iso-density plots for model 17. The color figure can be viewed online. 

Fig. 21 Iso-density plots for model 18. The color figure can be viewed online. 

By comparing the left panels of Figures 19, 20 and 21, corresponding to models 16-18, with those of models 2-4, shown in Figures 5, 6 and 7, we conclude that when the Fourier mesh begins to increase, the models generated with a divergence-free turbulence lose the homogeneity and isotropy of the initial overdensity clump distribution earlier than those models generated with a curl-free turbulence.

3.6. Models 19-22

When the number of Fourier modes increases, then the number of over-density clumps induced initially also increases, so that an homogeneous distribution of small length over-density clumps is formed; see the left panel of Figure 22 and compare it with the left panel of Figure 18. When the size of the Fourier mesh increases, which is the case of models 20-22, then we observe the formation of a smaller number of over-density clumps; see the left panels Figures 23 and 24 and compare them with those of Figures 19 and 20. However, for model 22, the length of the over-density clumps begins also to increase; see the left panel of Figure 25 and compare it with that of Figure 21 of model 18.

Fig. 22 Iso-density plots for model 19. The color figure can be viewed online. 

Fig. 23 Iso-density plots for model 20. The color figure can be viewed online. 

Fig. 24 Iso-density plots for model 21. The color figure can be viewed online. 

Fig. 25 Iso-density plots for model 22. The color figure can be viewed online. 

We thus conclude that with a larger number of modes, the loss of the homogeneity and isotropy of the initial over-density clump distribution is more delayed; see the left panels of Figures 23, 24, and 25.

We do not see any significant difference between the initial panels of models 19-20 shown in Figures 22-23 and those of models 5-6 shown in Figures 8-9, so that the two types of turbulence spectra are almost indistinguishable. However, by comparing the initial over-density distribution of clumps for models 21 and 22, shown in Figures 24 and 25, with those of models 7 and 8, shown in Figures 10 and 11, we see that the over-density clumps are in general larger for the divergence-free turbulence type than those generated with the curl-free turbulence type.

3.7. Models 23-25

When the size of the perturbation mode increases further, as for models 23-25, then it is evident that there are only two main directions along which the initial over-density clump distributions grow: see the left panels of Figure 26, 27 and 28. For this reason, the deformation in the core shape is very significant; see the right panels of Figure 26, 27 and 28. Despite this, the final outcome of these models is still a single protostar, which is displaced away from the central region, as a consequence of the highly anisotropic accretion of gas.

Fig. 26 Iso-density plots for model 23. The color figure can be viewed online. 

Fig. 27 Iso-density plots for model 24. The color figure can be viewed online. 

Fig. 28 Iso-density plots for model 25. The color figure can be viewed online. 

By comparing the initial clump distributions of models 23-25 with those of models 9-11, shown in Figures 12, 13 and 14, we note that the over-density grows in orthogonal directions for these two sets of models. This is expected according to equations 11 and 15, as the velocities are orthogonal vectors.

3.8. Models 26-28

When the number of Fourier modes increases, as in models 26-27, which are shown in Figures 29, 30 and 31, in a set of models with a large number of modes, as in models 23-25 shown in Figures 26, 27 and 28, then we observe that the core is more deformed, so that even core disruption can be achieved, as was the case with model 27, in which a binary protostars system is formed.

Fig. 29 Iso-density and plots for model 26. The color figure can be viewed online. 

Fig. 30 Iso-density plots for model 27. The color figure can be viewed online. 

Fig. 31 Iso-density plots for model 28. The color figure can be viewed online. 

The same behavior was observed in models 12- 14, shown in Figures 15, 16 and 17, in which we see again the orthogonal directions of the growth of the over-density clumps when compared with those observed in models 26-28.

3.9. Integral properties

We present here some integral properties of the resulting protostars, such as the mass and the values of the energy ratios α_f and β_f. These properties are calculated by using a subset of the simulation particles, which is determined by means of the following procedure. We first locate the highest density particle in the collapsed central region of the last available snapshot for each model; this particle will be the center of the protostar. We then find all the particles which (i) have density above some minimum value, given by log₁₀ (ρ_min/ρ₀) = 1.0 for all the turbulent models; (ii) are also located within a given maximum radius r_max from the protostar’s center.

All the calculated integral properties are reported in Table 2, whose entries are as follows. The first column shows the number of the model; the second column shows the parameter r_max in terms of the initial core radius R₀; the third column shows the mass of the protostar given in terms of M_⊙; the fourth and fifth columns give the values of α_f and β_f, respectively.

There are two lines in Table 2 for models 14 and 27, as their outcomes are binary protostar systems, so that each line indicates the properties of each component of the binary protostar. For these resulting binary systems, we simply define the binary separation as the distance between the centers associated with each protostar. We find the the separation is 163 AU for each binary system.

We find that the average mass of the protostars for models 1-13, excluding model 14, whose outcome is a binary system, is 0.94 M_⊙, with a standard deviation of 0.47 M_⊙. Analogously, the protostar average mass for models 15-26 and 28, excluding the mass of the binary of model 27, is 0.67 M_⊙, with a standard deviation of 0.34 M_⊙. It must be emphasized that these mass relations do not change much if we include the mass of the protostars in the binaries.

In Figure 32 we show the mass of the protostars in terms of the simulation model number. It seems that there is nothing that allows us to clearly distinguish the turbulent spectrum used initially; neither is there any trend to assess the effect that the increase in the number of Fourier modes can have on the resulting protostar mass; see also Table 2. However, there seems to be a tendency to low protostar masses for higher model numbers in both turbulence spectra. If this is true, then it would imply that the larger the perturbation mode, the lower the protostar mass.

Fig. 32 Proto-stellar mass in terms of the simulation model number. 

In Figure 33 we show that most protostars are near or on their way towards virialization, as they are close to the virial line; see § 2.5.It should be noticed that this tendency of collapsing cores to evolve toward the virial line was first pointed out by ^{Boss (1980)} and ^{Boss (1981)}. We also see that higher values of α_f and β_f are obtained for the curl-free turbulence type.

Fig. 33 Dimensionless energy ratios of the resulting protostars.