Gravitational Focusing and the Star Cluster Initial Mass Function

TL;DR

This paper explores how gravitational focusing and accretion processes can naturally produce the observed power-law initial mass function of star clusters, supported by numerical simulations.

Contribution

It demonstrates through N-body simulations that gravity-driven accretion alone can generate the cluster mass function without gas physics.

Findings

N-body simulations produce a power-law mass distribution.

Gravitational focusing can explain the cluster initial mass function.

Conditions for cluster formation timescales are identified.

Abstract

We discuss the possibility that gravitational focusing, is responsible for the power-law mass function of star clusters $N(\log M) \propto M^{-1}$. This power law can be produced asymptotically when the mass accretion rate of an object depends upon the mass of the accreting body as $\dot{M} \propto M^2$. While Bondi-Hoyle-Littleton accretion formally produces this dependence on mass in a uniform medium, realistic environments are much more complicated. However, numerical simulations in SPH allowing for sink formation yield such an asymptotic power-law mass function. We perform pure N-body simulations to isolate the effects of gravity from those of gas physics and to show that clusters naturally result with the power-law mass distribution. We also consider the physical conditions necessary to produce clusters on appropriate timescales. Our results help support the idea that gravitationally-dominated accretion is the most likely mechanism for producing the cluster mass function.

Tables

Table 1. Table 1 ChaNGa runs

run	N	geometry	seed	$\overline{{\ell }_i}$	$ϵ$	min ${\ell }_k$	${\alpha }_i$
40s1b	40000	sphere	1	0.047	0.001	0.016	0.01
40s2b	40000	sphere	2	0.047	0.001	0.016	0.01
40d1b	40000	disk	1	0.031	0.001	0.016	0.01
40d2b	40000	disk	2	0.031	0.001	0.016	0.01
400s1b	400000	sphere	1	0.021	7E-4	0.016	0.01
400s2b	400000	sphere	2	0.021	7E-4	0.016	0.01
400d1b	400000	disk	1	0.015	7E-4	0.016	0.01
400d2b	400000	disk	2	0.015	7E-4	0.016	0.01

Table 2. Table 2 Fits for the mass function slope of all ChaNGa runs

run	binning	slope $\pm$ error
40s1b	averaged	-0.89 $\pm$ 0.03
	min	-0.87 $\pm$ 0.03
	max	-0.90 $\pm$ 0.03
40s2b	averaged	-0.91 $\pm$ 0.08
	min	-0.80 $\pm$ 0.04
	max	-0.95 $\pm$ 0.04
40d1b	averaged	-0.99 $\pm$ 0.04
	min	-1.07 $\pm$ 0.05
	max	-1.01 $\pm$ 0.04
40d2b	averaged	-1.06 $\pm$ 0.03
	min	-1.06 $\pm$ 0.03
	max	-1.09 $\pm$ 0.05
400s1b	averaged	-0.99 $\pm$ 0.05
	min	-0.90 $\pm$ 0.03
	max	-1.00 $\pm$ 0.02
400s2b	averaged	-0.98 $\pm$ 0.06
	min	-0.85 $\pm$ 0.04
	max	-1.00 $\pm$ 0.03
400d1b	averaged	-1.18 $\pm$ 0.04
	min	-1.20 $\pm$ 0.03
	max	-1.16 $\pm$ 0.02
400d2b	averaged	-1.17 $\pm$ 0.04
	min	-1.10 $\pm$ 0.03
	max	-1.20 $\pm$ 0.03

Equations

$$\dot{M}=\frac{4\pi G^2{M}^2{\rho }_{\circ }}{(c_s^2+v_{\infty }^2{)}^{3/2}}\equiv \alpha M^2,$$

Full text

††slugcomment: Accepted to ApJ

1. Introduction

In recent years significant progress has been made on the form of the stellar cluster IMF (scIMF; Zhang & Fall, 1999; de Grijs et al., 2003; McCrady & Graham, 2007). As summarized by Fall & Chandar (2012), clusters in the Milky Way (Lada & Lada, 2003), Magellanic Clouds, M83, M51, and the Antennae exhibit a similar power-law mass function with $d\log N/d\log M\equiv\Gamma=-0.9\pm 0.15$. Others suggest a Schechter-type mass function is a better fit, with a truncation at high masses (Gieles et al., 2006; Bastian, 2008); however, all agree on the power-law behavior of $\Gamma\sim-1$ at lower masses.

One explanation for the scIMF has been that it simply reflects the mass distribution of parent giant molecular clouds (e.g., Elmegreen & Efremov, 1997), perhaps with some accounting for feedback by stellar energy input (Fall et al., 2010). However, the observed mass functions of inner Milky Way giant molecular clouds (GMCs) yield $\Gamma\sim-0.5$ (Williams & McKee, 1997); the Antennae also show similarly flat mass functions (Wilson et al., 2003; Wei et al., 2012), although outer Milky Way clouds may exhibit $\Gamma\lesssim-1.1$ and LMC clouds $\Gamma\sim-0.7$ (Rosolowsky, 2005). Moreover, while GMC masses obviously provide an upper limit to cluster masses, it isn’t obvious why the global cloud mass function should be reflected in cluster masses that are one to two orders of magnitude smaller. Even if the fraction of dense gas which can form clusters is similar among molecular clouds (and the evidence for this is mixed; see Bastian, 2008; Lada et al., 2010; Burkert & Hartmann, 2013; Battisti & Heyer, 2014), the manner in which the dense gas component might be broken up into individual protocluster clouds is not necessarily the same (see Hennebelle, 2012, for a model of modified turbulent acccretion).

Long ago, Zinnecker (1982) showed that if bodies accrete mass at a rate proportional to the square of their mass, $\dot{M}\propto M^{2}$, a population with a narrow initial mass range will develop a power-law mass distribution which asymptotically approaches $\Gamma=-1$. Zinnecker pointed out that “Bondi-Hoyle-Lyttleton” (BHL) accretion, with a mass accretion rate characterized by

[TABLE]

(see discussion in Edgar, 2004) has the requisite dependence on the central gravitating mass $M$.

In the standard picture of BHL accretion, a central gravitating mass, $M$, travels through an infinitely large ambient medium $\rho_{\circ}$ which has constant sound speed $c_{s}$ at constant relative velocity $v_{\infty}$. The mass’s gravity focuses material into a wake bound to the object from which it accretes material at the rate in Eq. 1. But the molecular clouds within which stars and clusters form exhibit density and velocity fields that are far from uniform; in addition, the cloud itself is self-gravitating. Nevertheless, using isothermal SPH simulations with decaying turbulence, Ballesteros-Paredes et al. (2015) (BP15) demonstrated the development of power-law sink mass functions with $\Gamma=-1$. These results suggest that gravitational focusing - which is at the heart of the BHL accretion process - is able to operate despite the complex, time-variable environment.

In this contribution, we examine the conditions for which gravitational focusing can produce star clusters with the required mass function. We show that observations of cluster-forming clouds are consistent with the requirements provided the observed supersonic velocity distributions are understood as being gravitationally-generated rather than turbulence driven by an external agent. We suggest that the signature of the dominance of gravitational over thermal physics is the power-law behavior of the mass function thereby linking the cluster and stellar IMF.

2. Numerical Simulations

The power-law sink mass functions found by BP15 were produced in SPH calculations with an initial supersonic velocity field but allowing the turbulence to decay. To minimize the role played by thermal physics, an isothermal equation of state was adopted. If, as suggested by the possible connection with BHL accretion gravity dominates, a similar result should be found with gravitational accretion without including any gas physics. We therefore explore results using a pure N-body code.

2.1. Numerical Setup

To test the scenario of purely gravitational accretion, we use the ChaNGa (Charm N-body GrAvity solver) code (Jetley et al., 2008; Menon et al., 2015). Since ChaNGa is a cosmological code, using the “dark matter” particle implementation ensures that all interactions are purely gravitational, with no gas physics included. With a standard $\Lambda$-CDM cosmology, we choose a set of units in which the code length unit corresponds to $1$ $\mathrm{kpc}$ of proper distance, the mass unit is $2.22\times 10^{5}$ $\mathrm{M_{\odot}}$ and the time unit is $1$ $\mathrm{Gyr}$. For generality, the simulation can be rescaled to representative values. Particles are originally distributed so as to have a uniform density within the chosen geometry. We generate a three dimensional initial turbulent velocity field with the Kolmogorov power law spectrum $P(k)\propto\mathbf{k}^{-11/3}$ with maximum wavenumber $k_{max}=64$, so that each run is seeded with turbulent velocity fluctuations. The primary function of the initial turbulence is to ‘stir’ up random density fluctuations.

We include both spherical and disk geometries, run with different random seeds. The initial conditions are sub-virial for both cases, with the initial virial parameter, $\alpha_{i}=2|K|/|U|=5\sigma^{2}R/GM$, given in Table 1, leading to an overall collapse. All runs have particle distributions of total mass $M=1$, total radius $R=1$, and we use $G=1$ in code units; disk runs start off with a thickness, $h=0.1$. We employ the default periodic boundary conditions, with a box size set to $10R$, to avoid interactions with the boundary. Force smoothing to avoid tight binaries is implemented in the form of a softening parameter $\epsilon$, based on Dehnen (2001). Table 1 lists each run and its relevant parameters.

2.2. Cluster Finding with FOF

To evaluate the formation of star clusters over time in the simulation, we implement a cluster finder based on the Friends of Friends (FOF) algorithm. FOF has been widely used in the cosmology community as a halo-finder instrumental in extracting halo mass functions from simulations (Knebe et al., 2011). The algorithm has only one free parameter, the linking length $b$, and operates on a set of particle positions in three dimensions. To locate clusters, the algorithm iteratively determines all the particles in a group that are within the linking length $b$ away from at least one other particle in the group.

The FOF algorithm offers a few advantages as a cluster finder: it is easier to capture irregular structures, an advantage when looking for stellar clusters, unlike some density threshold methods that construct halos/clusters out of spheres. With only one free parameter, there is relatively little fine tuning needed to produce mass functions, as long as there is some basis for the selection of the linking length.

To generate complete mass functions that probe substructure across different size scales, we use a hierarchical FOF. This approach utilizes a range of linking lengths to generate several sets of mass functions, with smaller values of the linking length finding smaller sub-structures. The mass functions generated at each scale can be summed, removing duplicate structures where applicable to create one hierarchical mass function. We empirically determine the base linking length for our setup by using the value that optimizes both the amount of particles in groups and the number of groups created. For most of the setups, this value runs close to $0.2\bar{\ell}_{i}$ , where $\bar{\ell}_{i}$ is the initial average interparticle distance, which guarantees groups that are at least 125 times denser than the initial state and an overdensity of at least 3 above the expected interparticle density if the final state were homogeneous. With hierarchical FOF, we probe down to $0.08\bar{\ell}_{i}$ at maximum .

The final step for our cluster finder is to filter out structures that are clusters of clusters, which we don’t want to include in our analysis. In the case of FOF, clusters that are adjacent to one another can be grouped together into one single cluster, even if they are not representative of a single object. To filter out these cases, we adopt a criterion based on the behavior of the average cluster density with mass. In a mass accretion scenario that produces large scale structures by the accretion of smaller scale structures onto a larger gravitating bound structure, the mean densities of more massive structures tend to increase. This will not be true for clusters of clusters, which are grouped solely on the basis of their adjacency. The average density of these objects flattens out and decreases, consistent with groups that are growing more in radius than they are in mass.

In Figure 1, we see that there exists a rough truncation mass above which the clusters are likely to be clusters of clusters due to a turnover in average cluster density. By only including clusters less massive than the truncation mass in our analysis, we remove clusters of clusters, but keep all sub-structures that comprised them due to the hierarchical nature of the data sets.

2.3. Results

2.3.1 Cluster Growth Through Gravitational Focusing

In Figure 2, we present a series of three epochs of cluster formation in run 40s1b. Run 40s1b is an initially spherical distribution of 40000 particles. In Figure 3, we show the time evolution of its mass function, where the right panels of the Figure correspond to the three epochs in Figure 2. (For the sake of generality, we will use the number of particles in a cluster as a proxy for the mass of the cluster for all figures, i.e. $M=N[\bar{m}]$). Figure 2 frame a) shows an initial phase where most groups are of a similarly low mass. At this point, over half of the particles in the simulation can be assigned to a group. This distribution is a result of the initial turbulent mixing in the simulation. By the epoch in Figure 2 frame b), some intermediate mass groups have been created. We see this in the corresponding frame in Figure 3 d), where the tail of the distribution has started to grow, creating a shallower slope. In the last epoch shown in Figure 2 c) and the last panel of Figure 3, the intermediate mass clusters have rapidly accreted new members, such that many of them are now some of the most massive clusters and the mass function has grown to approach the asymptotic $\Gamma=-1$ slope. While at this stage, all groups found by FOF are bound by virtue of the virial parameter $\alpha<1$, there is a subset of particles that were not assigned to groups and remain unbound. We include the final power law fits of the mass function for all runs in Table 2, as well as sample values for fits with small and large numbers of bins. The “averaged” slope values are a weighted average of fits over a range of bin sizes and different density refinement algorithms.

In the final frame (f), the mass function appears to be a little flatter than -1 at lower masses, but this is inevitable, the reasons for which are twofold: the initial conditions and the limits on the cluster finder. First, the $\Gamma=-1$ limit is an asymptotic one, and as the starting cloud is finite, eventually smaller mass clumps will simply cease to form due to a lack of replenishing material. Second, while a hierarchical mass function ought to probe the substructure of clumps, we are limited by the smallest linking length we set for the cluster finder: an arbitrary limit due to computational constraints.

In Figure 4, we find that the flatter feature does not exist in the disk geometry; the mass functions produce a relatively constant slope across the entire mass range. Disks are prone to an edge effect, where the outside of the disk tends to collect a ring of particles as it contracts (Burkert & Hartmann, 2004). This thin layer of excess pile-up creates a region with an overabundance of lower mass groups which, in this case, compensates for the underpopulation of lower mass clusters we see in the spherical geometry. Despite these differences, slopes between the different geometries vary by $<0.2$.

2.3.2 Mass Accretion History

For additional insight into the development of the mass function, we look at the mass accretion history of representative groups. Selecting groups from the power law tail of the distribution, we characterize the fashion in which the group attains mass by tracing the former groups its members belonged to over time. We can demonstrate that for a typical group at early times there exists a dominant sub-group which will attract smaller groups and accrete them (Figure 5). Combined with the demonstration that the $\Gamma=-1$ power law is not present at early times, it is unlikely that the mass distribution we see is seeded in the initial conditions, but instead is a function of gravitationally-driven accretion from a random distribution.

3. Gravitational focusing vs. BHL accretion

As mentioned in the Introduction, Zinnecker (1982) showed that if the accretion of an initial population of similar masses scales directly as $M^{2}$ (as is the limiting case of BHL accretion), $\Gamma\rightarrow-1$ asymptotically. As shown in Fig. 6, $\dot{M}$ vs $M$ starts off as roughly $M^{2}$, but later shows bending at low and high masses. The same behavior was seen in BP15, in which the SPH simulation nevertheless produced a $\Gamma=-1$ power law in the sink mass function. BP15 attributed this to “starving” the small masses by competition with higher mass sinks, and slowing down the high masses due to depletion of their environments. But if the strict BHL accretion formula, or a pure $\dot{M}\propto M^{2}$ is not strictly applicable, why does $\Gamma\rightarrow-1$?

In the case of the SPH runs with sink formation, BP15 suggested that the answer was that the gas densities and velocities were uncorrelated with the final sink mass; thus on average the $M^{2}$ dependence wins out. They found that specific local groups of sinks, embedded in similar environments, exhibited accretion rates proportional to $M^{2}$ but with different values of $\alpha$ (equation 1) for differing groups. BP15 then argued that adding up the individual groups with similar power-law distributions results in a common power law.

The difficulty with attempting to find an analytic or schematic understanding of the development of the mass function is that the environment is not static, but instead is strongly perturbed by the formation of local centers of gravitational attraction. A naive application of the BHL formula (equation 1) might suggest that accretion over large enough scales to make massive clusters is difficult, given the observed increase in cloud velocity dispersions with increasing size (Larson, 1981). However, as Heyer et al. (2009) showed, the Larson velocity-size relation for molecular clouds is consistent with near-virial motions, as would be expected if gravity were the driving force (see also Ballesteros-Paredes et al., 2011). In this case it is difficult to identify the appropriate value of $v$ to use in equation 1, and so formally the BHL accretion picture doesn’t apply. However, the power-law mass function that results and its development with time ( Figures 3, 4) is quite similar to simulations of pure BHL accretion (see, e.g., Figure 10 in Hsu et al., 2010)). Perhaps a schematic way of understanding the results is that the initial random velocity fluctuations average out so that overdensities can accrete matter at a rate $\propto M^{2}$. Whatever the interpretation, the common results of both the SPH simulations in Ballesteros-Paredes et al. (2015) and the N-body calculations in this paper provide strong evidence for gravitational accretion producing power law distributions with $\Gamma\sim-1$.

4. Discussion

The simulations in this paper and in Ballesteros-Paredes et al. (2015) refer to collections of particles and sink particles, respectively. Protostar clusters involve a combination of gas and stars. To apply our results to real clusters, we must assume that the protocluster gas $+$ star mass function maps directly into the final scIMF. This seems reasonable because the total mass in stars is generally thought to be $\gtrsim 0.3-0.5$ of the original protocluster cloud mass for the cluster to remain gravitationally bound (unless gas loss is extremely slow; Hills, 1980; Mathieu, 1983; Geyer & Burkert, 2001; Lada & Lada, 2003). Star formation regions with lower star formation efficiencies will likely not be able to become clusters at all (Kroupa et al., 2000). In addition, observations of nearby star clusters, like the ONC, show that clusters themselves form in regions of dense gas at high efficiencies, without significant contributions from diffuse gas components of molecular clouds (Hillenbrand & Hartmann, 1998). Moreover, the similarity of cluster mass functions both younger and older than 10 Myr in the interacting galaxy system of the Antennae and the LMC, suggests that feedback does not alter the shape of the cluster mass function (Fall et al., 2010).

Our identification of the cluster mass function as a result of gravitationally-driven accretion implies that clusters form subvirially. Generally speaking, the formation of a massive star cluster must occur in $\lesssim 3-10$ Myr, to avoid disruption by the energy input from the massive stars. This requires that the free-fall time be $\sim 10$ Myr or less, which requires molecular hydrogen densities $\gtrsim 100{\mathrm{c}m^{-3}}$, a plausible value for molecular clouds. The other requirement is that expansion velocities or dispersions be at most roughly virial (assuming some dissipation of supersonic motions). While observed velocity dispersions of molecular clouds increase with increasing scales (Larson’s first law; Larson, 1981), recent observations suggest that these motions are roughly virial (Heyer et al., 2009) (see also Larson (1981)). If the supersonic velocities are largely gravitationally-driven, avoiding the problem of too rapid dissipation (Ballesteros-Paredes et al., 2011), this second criterion is automatically satisfied.

Testing the hypothesis of gravitationally-driven accretion and/or subvirial initial conditions directly with the currently available observations is difficult (Proszkow et al., 2009; Kuznetsova et al., 2015), but the advent of Gaia may make it possible to search for collapse in statistically-significant samples of stars, especially combined with radial velocity measurements (Tobin et al., 2009; Kounkel et al., 2016; Da Rio et al., 2016). The gravitationally-driven accretion seen in simulations tends to form infalling filamentary streams, and this may be observable in the gas using appropriate molecular tracers.

Morphological tests are possible and even easier to apply. As Burkert & Hartmann (2004) showed, gravitational focusing tends to produce concentrations of mass in finite clouds near regions of smaller radii of curvature at the cloud edge; the simplest example of this is formation of dense gas near the ends of filamentary clouds. While this picture is consistent with the spatial structure of the Orion A cloud (Hartmann & Burkert, 2007), it might be possible to put this on a firmer statistical basis, using large-scale maps from the Spitzer and Herschel Space Telescopes and other facilities (Churchwell et al., 2009; André et al., 2010; Mairs et al., 2016).

Our picture of gravitationally-focused accretion naturally explains the tendency of stars to form in groups and clusters. It is also consistent with the idea that the supersonic motions in dense molecular clouds tend to be driven by gravity (Ballesteros-Paredes et al., 2011). Finally, this picture suggests a close connection with the upper-mass slope of the stellar IMF (Ballesteros-Paredes et al., 2015), with departures from the limiting $\Gamma=-1$ slope to that of the Salpeter value arguably the result of feedback from high-mass stars halting accretion.

Detailed numerical simulations of star and cluster formation in galaxy models with complicating effects such as stellar feedback and magnetic fields could probe the limits of this simple picture and help constrain the expected masses of star clusters as a function of environment for comparison with observations.

AB acknowledges important conversations with Clare Dobbs, LH with Mike Fall, Javier Ballesteros-Paredes, and Mark Heyer. This work was supported in part by NASA grants NNX16AB46G and NNH15ZDA001N-XRP,by the University of Michigan, and in part through computational resources and services provided by Advanced Research Computing at the University of Michigan, Ann Arbor.