ON THE INITIAL MASS FUNCTION OF LOW-METALLICITY STARS: THE IMPORTANCE OF DUST COOLING

With the calculations shown in this work are aimed to model the collapse of star-forming clouds at very low metallicity environments. For that, we make use of a version of Gadget 2 (Springel 2005; a smoothed particle hydrodynamics (SPH) code). In order to advance the calculations later in the collapse, when densities become very high, we make use of the sink particle technique (Bate et al. 1995; Jappsen et al. 2005). This approach enables us to carry the calculations even after the formation of the first contracting protostar.

These sink particles are created when a minimum of 100 SPH particles exceed the threshold number density of 5.0 × 10¹⁵ cm⁻³, which is reached when the gas becomes optically thick and further fragmentation is unlikely. For further discussion of the details of our sink particle treatment, we refer the reader to Clark et al. (2011a).

As discussed, e.g., by Omukai et al. (2005, 2010), the gas becomes optically thick at densities around 10¹⁵ cm⁻³. Beyond this value, further fragmentation is strongly suppressed and our choice of sink particle formation threshold ensures that single sink particles represent individual protostars rather then groups of objects (see also the discussion by Greif et al. 2012).

Furthermore, SPH with the smoothing we employ (i.e., the standard SPH cubic spline smoothing) prevents fragmentation below the resolution limit. As shown by many works (Bate & Burkert 1997; Whitworth 1998; Hubber et al. 2006), SPH does not suffer from artificial fragmentation, at the resolution limit if gravitational forces and pressure gradients are resolved equally well.

As for the dust grain properties and chemical network employed, we use the same recipe as in Dopcke et al. (2011).

The calculations performed in this work consisted in a set of four simulations, with metallicities Z/Z_☉ = 10⁻⁴, 10⁻⁵, 10⁻⁶, and the metal-free case. In each simulation we employed 40 million SPH particles, with a total mass of 1000 M_☉. The initial number density was taken to be 10⁵  cm⁻³ and an initial temperature of 300 K. We also included rotational and turbulent energy, which were taken to be 2% for rotation and 10% for turbulence, when compared to the gravitational potential energy.

The mass resolution is 2.5 × 10⁻³ M_☉, which corresponds to the mass of 100 SPH particles (see, e.g., Bate & Burkert 1997). The redshift chosen was z = 15, when the cosmic microwave background (CMB) temperature was 43.6 K. The dust properties were taken from Goldsmith (2001). In the calculations, the opacities vary linearly with Z, which means for instance that for the Z/Z_☉ = 10⁻⁴ calculations, the opacities were 10⁻⁴ times the original values.

Dust cooling is a consequence of inelastic gas–grain collisions, and thus the thermal energy loss from gas to dust vanishes when they have the same temperature. We therefore expect the cooling to cease when the dust reaches the gas temperature. In order to evaluate the effect of dust on the thermodynamic evolution of the gas and verify this assumption, we plot in Figure 1 the temperature and density for the various metallicities tested. With the simulations, we can compare the thermodynamic evolution of the dust and gas, at the point of time just before the first sink particle was formed (see Table 1). The dust temperature (shown in blue) can vary from the CMB temperature floor on the beginning of the collapse to the gas temperature (shown in red) at higher densities.

Zoom In Zoom Out Reset image size

Changes in metallicity influence the density at which dust cooling becomes efficient. For Z = 10⁻⁴ Z_☉, dust cooling begins to be efficient at n ≈ 10¹¹ cm⁻³, while for Z = 10⁻⁵ Z_☉, the density where dust cooling becomes efficient increases to n ≈ 10¹³ cm⁻³. For the Z = 10⁻⁶ Z_☉ case, dust cooling becomes important for n ≳ 3 × 10¹³ cm⁻³, preventing the gas temperature from exceeding 1500 K. For comparison, in the metal-free case the gas reaches temperatures of approximately 2000 K.

The efficiency of the cooling is also expressed in the temperature drop at high densities. The temperature of the gas decreases to a few hundred kelvin, a value much lower than the metal-free case. Such a temperature drop can significantly increase the instability of the gas to fragmentation. At much higher densities (n ≳ 10¹³ cm⁻³), the temperatures of dust and gas become coupled since compressional heating dominates over all of the cooling. After this, the thermal evolution of the gas is approximately adiabatic.

When we compare our results to the calculations of Omukai et al. (2010), we find good agreement with their one-dimensional hydrodynamical models, although we expect some small difference due to effects of the turbulence and rotation (see Dopcke et al. 2011) and also due to the use of different dust opacity models.

During the collapse, the thermodynamic evolution of the gas takes different paths depending on the metallicity, as shown in the density–temperature diagram (Figure 1). In order to explain this behavior, we take a closer look at the cooling and heating processes involved. In Figure 2 we show the main cooling and heating rates divided into four panels for the different metallicities. These rates were calculated by averaging values of individual SPH particles in one density bin, where the total density range was divided into 500 bins in log space.

Zoom In Zoom Out Reset image size

At densities below n ≈ 10¹⁰ cm⁻³, dust cooling is unimportant in all of the runs. At these densities, the dominant coolant is H₂ line emission, while the heating is dominated by compressional (pdV) heating at n ≲ 10⁸ cm⁻³, and by three-body H₂ formation heating at higher densities.

At higher densities, dust cooling starts to play a more important role. In the Z = 10⁻⁴ Z_☉ simulation, dust cooling exceeds pdV heating at n ≈ 10¹⁰ cm⁻³, although it does not exceed the H₂ formation heating rate until n ≈ 10¹¹ cm⁻³. Once this occurs, and dust cooling dominates, the gas temperature drops sharply. In the Z = 10⁻⁵ Z_☉ simulation, on the other hand, dust becomes the dominant coolant only at n ≈ 10¹³ cm⁻³, and so the temperature decrease happens later and is smaller. Finally, in the Z = 10⁻⁶ Z_☉ case, dust cooling becomes competitive with pdV heating only at the very end of the simulation, and so the effect on temperature evolution is less evident.

The other thermal processes play a minor role during the collapse. For example, H₂ dissociation cooling only becomes important in the runs with Z = 10⁻⁶ Z_☉ and 0, and only for n > 10¹³ cm⁻³. At very high densities (n > 10¹⁴ cm⁻³), H₂ collision-induced emission (CIE) cooling also begins to be important. For more details on H₂ heating and cooling processes in this very high density regime, we refer to Clark et al. (2011a).

The angular momentum is transported to smaller scales during the collapse. This process leads to the formation of a rotation supported disk-like structure. However, the rotation is not sufficient to hold the collapse of the disk, and finally it fragments forming multiple objects.

In Figure 3, we show the density structure of the gas. These snapshots were taken immediately before the first protostar was formed. Observe that at large scales the gas cloud properties are the same for all metallicities. Differences in the thermodynamic evolution appear only at n ≳ 10¹¹ cm⁻³ (see Figure 1). As a consequence, we observe variations in the cloud structure only in the high-density regions.

Zoom In Zoom Out Reset image size

One way to study the effect of dust cooling on the fragmentation behavior and the final stellar IMF is to look at the changes in the number of Bonnor–Ebert (M_BE) masses contained in this central dense region. Using the definition from Bromm et al. (2009),

$$\begin{equation} M_{{\rm BE}} = 500\ M_{\odot } \left(\frac{T}{200\ \rm K}\right)^{3/2} \left(\frac{n}{10^4 \ \rm cm^{-3}}\right)^{-1/2}, \end{equation} \tag{ 1 }$$

for an atomic gas with temperature T and number density n, we have computed the number of Bonnor–Ebert masses contained within a series of concentric radial spheres centered on the densest point in each of our four simulations. The results are shown in Figure 4.

Zoom In Zoom Out Reset image size

At the beginning of the simulation, the cloud had ∼3 M_BE. During the collapse, the gas cools and reaches ∼6 M_BE in all cases. Cooling and heating are different depending on the metallicity, and this difference is seen for distances smaller than ∼400 AU. The Z = 10⁻⁴ Z_☉ case, for instance, has twice the number of M_BE for distances smaller than ∼10 AU, when compared to the other cases. This will have direct consequences for the fragment mass function as we will see in the next section.

When the conditions for the formation sink particle are met (see Section 2.1), such particles begin to form in the density peaks (Figure 5). Then, a disk is built up in these regions, where fragmentation also occurs (Tohline 1980). During further collapse, this dense region creates spiral structures. For Z = 10⁻⁵ Z_☉ and 10⁻⁴ Z_☉, density waves build up spiral structures, which become locally gravitationally unstable and go into collapse. The formation of binary systems by triple encounters (Binney & Tremaine 2008) transfers kinetic energy to some sink particles, causing them to be ejected from the high-density region. For Z = 10⁻⁵ Z_☉, when the star formation efficiency (SFE) is 0.5%, fragmentation has already occurred in a secondary dense center, at a distance of ∼20 AU from the first dense region.

For Z = 10⁻⁶ Z_☉ and 0, the formation of spiral structures is not observed. In these two runs, star formation occurs mainly in the central clump.

The simulations were stopped at a point in time when 4.7 M_☉ of gas has been accreted into the sink particles. This is sufficient to identify the two fundamental modes of fragmentation discussed in Section 3.5.

Figure 6 shows the mass distribution of sink particles at that time. We typically find sink masses below 1 M_☉, with somewhat smaller values in the 10⁻⁴ Z_☉ case compared to the other cases. No sharp transition in fragmentation behavior was found, but rather a smooth and complex interaction between kinematic and thermodynamic properties of the cloud.

Table 1 lists the main sink particle properties. It shows that the time taken to form the first sink particle is slightly shorter for higher metallicities. This shorter time is a consequence of the more efficient cooling by dust, which decreases the thermal energy that was delaying the gravitational collapse. In Table 1 we also observe that the star formation rate is lower for Z = 10⁻⁴ Z_☉. This is because star formation started at an earlier stage of the collapse, when the mean density of the cloud was lower and there was less dense gas available to form stars.

To better understand whether the resulting protostellar cluster was affected by varying the metallicity, we plot the final sink mass distribution in Figure 6. It shows that for the simulations with Z ⩽ 10⁻⁵ Z_☉, the resulting sink particle mass function is relatively flat. There are roughly equal numbers of low-mass and high-mass stars, implying that most of the mass is to be found in the high-mass objects. This mass function is consistent with those found in other recent studies of fragmentation in metal-free gas (Greif et al. 2011; Smith et al. 2011). If the sink particle mass function provides a reliable guide to the form of the final stellar IMF, it suggests that at these metallicities, the IMF will be dominated by high-mass stars.

The resolution limit in the simulations is 2.5 × 10⁻³ M_☉. All sink particles mass function presented in Figure 6 have their masses well above this value by construction. These simulations were intended to study the influence of the metallicity on the mode of star formation. Such a goal could be achieved when the sink particles have accreted a total of approximately 5 M_☉. However, these systems should not be taken as star formation is completed, since the sink particles would continue accreting. Other mechanisms should also influence the mass function, and we list them in Section 4.

One way to explain the final mass distribution of the fragments is to look at the timescales for mass accretion and fragmentation. The degree of gravitational instability inside a volume can be represented by the number of Bonnor–Ebert masses contained in this volume. We can therefore estimate the fragmentation timescale by computing the time taken for the central dense region to accrete one Bonnor–Ebert mass. In other words, we have $t_{\rm frag} \equiv {M}_{{\rm BE}}/\dot{M}$, where $\dot{M}$ is the accretion rate. This value is shown as a function of the enclosed gas mass in Figure 7, where the values are calculated for particles in spherical shells, and the center is taken to be the densest SPH particle. $\dot{M}$ is obtained by summing up the mass of the particles (m_p) inside a shell times their radial velocity (v_r), $\dot{M} \equiv \sum _{{\rm shell}} m_{p}v_{r}$. This way of estimating the expected mass accretion in the sink particles is motivated by the results in Clark et al. (2011b, see their Figure 3(b)). They demonstrate that the mass accretion rate based on the radial mass infall profile is a good estimate of that actually measured by the sink particles.

For comparison, we also plot the accretion timescale, here defined as the time taken by the gas to accrete the mass enclosed by that radius, $t_{\rm acc} \equiv {M}_{\rm enc}/\dot{M}$. When the fraction t_frag/t_acc > 1, one expects that the gas enclosed by this shell is going to be accreted faster than it can fragment, favoring high-mass objects. Conversely, for t_frag/t_acc < 1, the gas will fragment faster than it can be accreted by the existing fragments, and the final mass distribution is expected to have more low-mass objects. Note that as defined here, the timescale on which new fragments form t_frag/t_acc is the inverse of the quantity M_gas/M_BE plotted in Figure 4.

In Figure 7, the simulation with Z = 10⁻⁴ Z_☉ has the lowest values for t_frag/t_acc, over a wide range of M_enc. This indicates that more low-mass fragments are expected to form in this case, leading to a steeper fragment mass function.

Now we can compare the predicted values before sink formation started with the final accretion and fragmentation timescales. These values are designed to represent the characteristic timescales on which the mass histogram changes: the fragmentation timescale (τ_frag) is the time in which the number of fragments change by a significant amount, while the accretion timescale (τ_acc) represents the time in which the existing fragments grow in mass. We therefore define $\tau _{\rm frag} \equiv n/\dot{n}$, and $\tau _{\rm acc} \equiv M/\dot{M}$, where n is the number of sink particles and M is the total mass incorporated into sink particles.

Both timescales should increase over time, since the number of sink particles and the total mass also increase. This reflects the fact that it takes longer to change the shape of the mass histogram in frequency when the number of elements is higher, and in mass when the total mass is high. The key point is that these timescales are related to the shape of the histogram. If the timescale for accretion is shorter than the timescale for fragmentation, the histogram will tend to be dominated by high-mass objects. Conversely, if the fragmentation timescale is shorter than the accretion timescale, the low-mass part will be populated before the objects can grow and occupy higher mass bins. Therefore, a comparison between τ_acc and τ_frag helps us to understand the shape of the mass histogram.

We calculate the average τ_acc and τ_frag at each point in time. The value for the mean τ_frag and τ_acc is calculated by the time-weighted individual times as in Equations (2) and (3):

$$\begin{eqnarray} \tau _{\rm frag} \equiv & \frac{1}{t} \sum \frac{1}{\dot{n_i}} \Delta t = \frac{1}{t} \sum \frac{t_i}{n_i} \Delta t, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \tau _{\rm acc} \equiv & \frac{1}{t} \sum \frac{\left\langle M \right\rangle }{\dot{M_i}}\Delta t = \frac{1}{t} \sum \frac{t_i \left\langle M \right\rangle }{M_i}\Delta t , \end{eqnarray} \tag{ 3 }$$

where i refers to each snapshot of the simulation, and it varies from two to the total number of snapshots. n_i, M_i, t_i, and t are the number of sink particles, the total mass in sinks, the time at the snapshot i, and the total time of sink particle formation in each simulation, respectively. 〈M〉 is the average mass of a sink particle in the simulations including all metallicities. This average mass (〈M〉) is used to make τ_acc in units of time, and comparable to the accretion time. Other values were tried, such as the total average mass for each metallicity. However, the accretion time was not affected considerably. Thus, the comparison between accretion and fragmentation leads to qualitatively equal results, independent of the value used to get the correct dimension for τ_acc.

Important to note that τ_{frag, merging} is calculated in the same way as τ_frag, with the difference that the number of sink particles (n_i) is always lower or equal to the non-merging case. That is why τ_{frag, merging} has higher values when compared to τ_frag.

Figure 8 shows the average timescales for fragmentation and accretion for different metallicities. We also plot the values considering the sink particles that could have merged due to collisions (see Section 4). These results explain the difference in the sink particle mass distribution in Figure 6. For Z ⩽ 10⁻⁵ Z_☉, the fragmentation time is always larger than the accretion time, indicating that the sink particles will accrete faster than they can be generated, resulting in a flatter mass distribution. On the other hand, when the fragmentation time is longer than the accretion time (for Z = 10⁻⁴ Z_☉), the gas fragments, rather than moving to the center and being accreted. As a consequence, the low-mass end of the protostellar mass function grows faster than the high-mass end, and the slope of the mass function steepens. This behavior agrees well with the predictions from before fragmentation started, as shown in Figure 7.

Note that the values in Figure 7 were calculated before the formation of the first sink particle, while the values in Figure 8 were calculated using the sink particle properties. A comparison between them is useful to evaluate whether the gas cloud properties from before star formation started could be used to predict its star formation behavior. The trend in both figures is that the fragmentation timescale is normally shorter than the accretion timescales for Z = 10⁻⁴ Z_☉. From these results we conclude that the mass distribution in Figure 6 can be explained by the timescales in Figure 8, in particular the fact that the Z ⩽ 10⁻⁵ Z_☉ simulations have more high-mass objects. The last finding is that the transition from Pop III to Pop II star formation mode is not abrupt, in the sense that there is no metallicity below which the gas cannot fragment. The transition is rather in the stellar IMF, and gas clouds with Z ≲ 10⁻⁵ Z_☉ form a more flat IMF, while gas clouds with Z ≳ 10⁻⁴ Z_☉ produce a cluster with more low-mass objects (see also Clark et al. 2008).

Another property of the star-forming cloud that we observed to vary in our calculations is the spatial mass distribution. The dependence of the enclosed gas and sink mass on the distance from the center of mass is shown in Figure 9. The Z = 0 case has almost all of the sink particle mass concentrated within r < 8 AU. The gas density for this case is also higher in this region, when compared to the other metallicities, showing that the gas and sink particle mass densities follow each other. The mass in sink particles exceeds the gas mass for small radii, being the most important component in the gravitational potential. For r > 150 AU, the gas becomes the most massive component, for all Z. Girichidis et al. (2012) also reported this behavior, but in their case the sink particles already started to dominate the potential below r ≈ 10³ AU.

This higher concentration of gas and sinks at the center occurs because for the Z = 0 case, the gas had higher temperatures in the central region. For high temperatures, the criterion for gravitational instability requires higher densities, which are achieved only very close to the center. As a consequence, the sink particle formation criteria are met just for short distances from the center.

Consequently, the dominance of sink particles mass in the gravitational potential over the gas mass, for radii smaller than 150 AU, shows the importance of treating gas and stars together in this sort of problem. It also suggests that N-body effects, such as ejections and close encounters, should play an important role in the formation of these dense star clusters, even in the very earliest stages of their evolution (see Smith et al. 2011; Greif et al. 2012).