Does turbulence determine the initial mass function?

Aside from the initial velocity fields, our setup is identical to that in BBB03: We set up a series of turbulent, spherical clouds, with $50\text{ M}_{\odot}$ of gas of uniform density with diameter 0.375 pc. The corresponding initial free-fall time is $t_{\mathrm{ff}}=1.90\times 10^{6}$ years. The minimum Jeans mass at the opacity limit is $M_{\mathrm{min}}\approx 0.0011\text{ M}_{\odot}$. We use 3.5 million SPH particles, consistent with BBB03, who showed that about 75 particles are required per $M_{\mathrm{min}}$ (see also Bate & Burkert 1997). Particles were distributed in a uniform random distribution. We adopt code units with a length unit of $0.1{\rm pc}$, mass unit of $1{\rm M}_{\odot}$ and time units such that $G=1$.

We adopt a barotropic equation of state $P=K\rho^{\gamma}$. Following BBB03, we prescribe $\gamma=1$ (i.e. isothermal) for densities lower than the opacity limit for fragmentation ($\rho=10^{-13}\text{ g cm}^{-3}$), $\gamma=7/5$ for $10^{-13}\text{ g cm}^{-3}<\rho<10^{-10}\text{ g cm}^{-3}$ and $\gamma=1.1$ for $\rho>10^{-10}\text{ g cm}^{-3}$. We define the constant $K$ to be such that the sound speed is $c_{s}=1.84\times 10^{4}$ cm s^-1 during the isothermal phase (i.e. 10 K assuming a mean molecular weight $\mu=2.46$) and in the $\gamma=7/5$ regime such that the pressure remains continuous when $\gamma$ changes. As discussed by Bate (2009a), using a barotropic equation of state over-produces low mass stars and brown dwarfs compared to observations, since the cold gas surrounding the protostars fragments too readily (c.f. Fig. 6). Several groups (Bate, 2009b, 2012; Offner et al., 2009; Krumholz et al., 2010; Commerçon et al., 2010) showed that this can be solved by modelling radiation in the flux-limited diffusion approximation. However, simulations with radiation are expensive, precluding the kind of statistical study we perform here, the radiation algorithm is not yet implemented in phantom, and a barotropic equation of state is sufficient to answer the question of whether the PDF influences the IMF. We also ignore magnetic fields which change the star formation rate and perhaps also the IMF (Ostriker et al., 1999; Heitsch et al., 2001; Vázquez-Semadeni et al., 2005; Tilley & Pudritz, 2007; Price & Bate, 2008, 2009; Myers et al., 2014).

We impulsively drive turbulence in each cloud, as in BBB03, by imposing an initial supersonic turbulent velocity field. The amplitude of the velocity fluctuations follow a power spectrum $P(k)\propto k^{-4}$ where $k$ is the wavenumber, in order to be consistent with Larson’s scaling relation. We generate each field via a Fourier transform on a $64^{3}$ grid, which is then interpolated onto the SPH particles. The coefficient of each Fourier mode is drawn from a Rayleigh distribution with each mode also given a uniform random phase between $[-\pi,\pi]$. This is equivalent to sampling from a cylindrical bivariate Gaussian (Dubinski, Narayan & Phillips, 1995).

To obtain a purely solenoidal velocity field, we take the curl of a vector field to produce a divergence-free velocity field. Similarly for a purely compressive velocity field, we take the gradient of a scalar field to produce a curl-free field. We compute the gradients in Fourier space. Velocities are normalised so that the initial kinetic energy is equal to the gravitational potential energy, giving an initial RMS Mach number of $\mathcal{M}=6.4$. We performed simulations using 7 realisations of the initial velocity field for each case (solenoidal or compressive), realised by changing the seed in the random number generator for the phases and amplitudes.

Following BBB03, we introduce sink particles (Bate et al., 1995) when the central density of pressure-supported fragments reaches $\rho_{s}=10^{-11}$ g cm^-3, two orders of magnitude higher than the opacity limit. Once $\rho_{s}$ is exceeded and sink formation conditions are satisfied, we replace gas particles within 5 au with a sink particle. Gas particles within 5 au are accreted if they pass checks for angular momentum and boundness, with their mass and momentum added to the sink. Gravity between sinks is softened within 4 au; gas particles are accreted without checks within this radius.

Figure 1 shows the evolution of column density from $t=0$ to $t=0.3t_{\rm ff}$ (left to right) in two representative calculations, using solenoidal driving (top, as in BBB03) and compressive driving (bottom). Shocks form quickly in both cases, due to the impulsive supersonic velocity field, but are stronger in the compressive case, driving the formation of large scale filaments after only $0.3t_{\mathrm{ff}}$. For the solenoidal case, $\nabla\cdot\mathbf{v}=0$ initially by definition, so there are no regions that initially promote collapse.

Figure 2 shows the subsequent small-scale fragmentation in the compressive cloud, with the first protostar formed after just $0.2t_{\rm ff}$. The process in all other clouds appears visually very similar. Gas flows into dense cores along filaments (e.g. Gómez & Vázquez-Semadeni, 2014; Smith et al., 2016; Federrath, 2016; Klassen et al., 2016), feeding young protostars via accretion discs. The process is chaotic and dynamical, with close encounters between stars resulting in the destruction of accretion discs, and the ejection of smaller mass objects. Bound systems form and get destroyed by interactions on a very short time-scale. The stars live in a competitive environment, where those that grow in mass quickly stay in the dense regions and accrete further material, whilst ejecting lower mass objects.

We computed the density PDFs by binning the particles into 2000 bins equally spaced between $-10<\log_{10}(\rho)<10$ in code units. We then computed the standard deviation, $\sigma_{\ln\rho}$ by fitting a log-normal distribution to the PDF (using scipy.optimize.curve_fit in Python). Note that the PDF computed in this way is mass, rather than volume-weighted. Both volume and mass-weighted PDFs are expected to be log-normal when the equation of state is approximately isothermal (Padoan et al., 1997b; Passot & Vázquez-Semadeni, 1998; Scalo et al., 1998; Nordlund & Padoan, 1999; Ostriker et al., 2001).

Figure 3 shows the time evolution of the (mass weighted) RMS Mach number, $\mathcal{M}$, (top panel) and standard deviation, $\sigma_{\ln\rho}$, (centre panel) for our entire set of calculations, with the solid lines showing the mean from the 7 different simulations for each type of velocity field and the shaded region shows the 1$\sigma$ standard deviation. The bottom panel shows the evolution in the $\mathcal{M}$-$\sigma_{\ln\rho}$ plane. In both the solenoidal and compressive clouds $\mathcal{M}$ decays with time due to the dissipation of energy by shocks, reaching a minimum before rising again once bound structures have formed.

Comparison of PDFs in decaying turbulence simulations is complicated by the time evolution of the velocity field. In our calculations the initial density field is uniform and the PDF thus develops in response to the initial turbulent velocity field. Since the clouds evolve on different time-scales, it is not particularly meaningful to compare their PDFs at the same time. Rather — for the purposes of our study — Equation 1 suggests that they should be compared at the same RMS Mach number $\mathcal{M}$ so that the only difference is from the different mixing parameters $b$.

The lower panel of Figure 3 shows that the initial collapse of the cloud roughly corresponds to $\sigma_{\ln\rho}\lesssim 2$. Once $\sigma_{\ln\rho}$ reaches this value $\mathcal{M}$ rises again once fragmentation begins. Also, the PDF is no longer log-normal. We thus use the time interval where $\sigma_{\ln\rho}<2$ to compare the density PDFs prior to the onset of star formation. The standard deviation of the PDFs is different not only at the same time early in the evolution of the cloud, but also at the same RMS Mach number.

Figure 4 shows the resultant PDFs computed at the time when all calculations have the same RMS Mach number of $\mathcal{M}=5.5$, which is when $\sigma$ differs most between the simulations. The difference in the PDF produced by compressive vs. solenoidal driving is similar to that shown by (e.g.) Federrath et al. (2008) and Federrath et al. (2010), except that we show the mass-weighted version. Compressive driving produces a broadening of the PDF caused by the collision of stronger shocks which in turn create larger variations in the density field. This demonstrates that our different choices of impulsive driving indeed drive significant differences in the density PDF prior to star formation.

Figure 5 shows the total stellar mass as a function of time, measured by the mass in sink particles. The onset of star formation occurs at $t\approx 0.2t_{\rm ff}$ in the compressive case, compared to $t\approx 1.1t_{\rm ff}$ in the solenoidal case. Once star formation starts in each calculation, the rate at which material is converted to stars is higher by a factor of $\sim 2$ in the compressive clouds compared to the solenoidal cloud. The overall efficiency of star formation is similar in both types of calculation over the time we have continued the simulations, with $\approx 15$% of the gas mass converted to stars. However, the efficiency is higher on an absolute scale since this occurs over a shorter timescale in the compressive case. Also, the end of the simulations does not mark the end of the star formation process since the mass in stars continues to increase.

Figure 6 shows the IMFs from our simulations, combining all 7 realisations with solenoidal (left) and compressive driving (right), with the cumulative IMFs shown in Figure 7. The IMFs of stars that have finished accreting (235 of 388 and 298 of 533 sinks for solenoidal and compressive, respectively) are shown in red, while the IMF of all stars is shown in blue. The lowest mass possible in our calculations is $\approx 0.005\text{ M}_{\odot}$ from the opacity limit for fragmentation, which sets the low-mass cutoff. The IMFs appear similar to those shown in BBB03 but with better statistics because of our multiple realisations. Our IMFs are also similar to those found by Bate (2009a) from one calculation of a 500 $\text{M}_{\odot}$ cloud. In particular, we observe the statistically significant excess in low mass stars and brown dwarfs compared to the Kroupa (2001) and Chabrier (2005) IMFs (dashed and solid lines, respectively) that occurs when a barotropic equation of state is employed (e.g. Bate, 2009a, b).

There is no obvious difference between the IMFs produced by the different types of driving. Statistics confirm this — a Kolmogorov-Smirnov test gives a p-value of 0.71 between the two distributions when considering all sink particles, and a p-value of 0.98 when considering only sinks that have finished accreting. This means we cannot reject the hypothesis that the samples come from the same underlying distribution. Thus, while the type of driving changes the density PDF, the resultant IMFs are indistinguishable.