The Stellar IMF from Isothermal MHD Turbulence

The origin of the stellar initial mass function (IMF) is still not fully understood. Many processes, such as magnetic support, radiative and mechanical feedbacks from young stars, density enhancements or pressure support from supersonic turbulence, dynamical interactions between accreting stars, or competitive accretion, may affect the mass distribution of stars. Their relative importance varies with environment and is still disputed. Numerical simulations of star formation have started to reveal the mechanisms controlling the star formation rate (SFR), thanks to systematic parameter studies based on large sets of simulations with a modest dynamic range of scales (e.g., Padoan & Nordlund 2011a; Federrath & Klessen 2012; Padoan et al. 2012), or, more recently, with large-scale simulations where many star-forming regions are formed ab initio and the time evolution and scatter of the SFR can also be investigated (Padoan et al. 2017). While the SFR is already converged in these simulations, a much larger dynamic range is needed to achieve a numerically converged IMF covering the whole spectrum of stellar masses, from brown dwarfs to massive stars. The computational cost of such experiments is beyond the reach of most studies of star formation, so stellar IMFs from numerical simulations are scarce and generally of low statistical significance. Typically limited to ∼100 stars, these simulations barely constrain the IMF turnover and do not yield enough massive stars to probe the Salpeter range. Unable to pursue the necessary range in space and timescales, numerical studies often opt for increasing the physical complexity, for example, including radiative and mechanical feedbacks, in the hope of discerning the effect of new processes on the IMF, despite the small sample size and dubious numerical convergence.

The first numerical IMF with a large enough number of sink particles to probe both the IMF turnover, including brown dwarfs, and the Salpeter range was presented in Padoan et al. (2014b), where we simulated randomly driven, supersonic MHD turbulence in a 4 pc region with a total mass of 3000 M_⊙ and a maximum resolution of 50 au. After a long initialization phase without self-gravity, the simulation was evolved with self-gravity, generating 1288 sink particles over a period of 3.2 Myr (2.7 free-fall times). The mass distribution of the sink particles was consistent with a Chabrier IMF (Chabrier 2005) at low masses and a power law with Salpeter's slope (Salpeter 1955) above 1–2 M_⊙. In previous numerical studies, the size of the simulated region, L, the duration of star formation (between the first and the last sink particles), t_SF, and the total number of sink particles, N_*, were respectively ∼10, ∼100, and ∼10 times smaller than in Padoan et al. (2014b). The radiation hydrodynamic (HD) smoothed particle hydrodynamics (SPH) simulations by Bate (2012) and Bate (2014) described a 500 M_⊙ region with L = 0.4 pc and yielded at most N_* = 183 over a time t_SF ∼ 0.09 Myr. The HD grid-based simulations by Krumholz et al. (2012), representing a 1000 M_⊙ region, had L = 0.46 pc, N_* = 158, and an even shorter star formation time, t_SF ∼ 0.02 Myr. The more recent MHD simulations by Myers et al. (2014), modeling a region of the same size and mass, yielded at most N_* = 92 over a time t_SF ∼ 0.05 Myr.

The earlier barotropic SPH simulations by Bate (2009), with 500 M_⊙ in a region with L = 0.4 pc, achieved a large number of sink particles, N_* = 1254, and a slightly longer star formation time, t_SF ∼ 0.15 Myr. However, that IMF peaked at 0.02 M_⊙, an order of magnitude below the peak of the Chabrier IMF, so most of the sink particles had very low masses, which are created only in the absence of a magnetic field or radiation feedback (Bate 2012), and would not exist in nature. An even earlier SPH simulation by Bonnell et al. (2003), modeling 1000 M_⊙ in a region with L = 1 pc, reached N_* ≈ 400 and t_SF ∼ 0.35 Myr. This IMF peaked around 0.3 M_⊙ and produced a realistic Salpeter range. However, lacking both magnetic field and radiation feedback, the realistic peak was almost certainly the result of the very limited number of SPH particles, 5 × 10⁵ compared with 3.5 × 10⁷ in Bate (2009). With such a low number of SPH particles, the IMF was incomplete just below its peak and, more importantly, the turbulent velocity field could not be resolved well enough to describe the smallest scales of the turbulent fragmentation process responsible for the IMF turnover. Similar considerations apply to the more recent SPH simulations by Ballesteros-Paredes et al. (2015), modeling a region of the same size, same mass, and similar duration of star formation as in Bonnell et al. (2003).

SPH and grid-based simulations without a magnetic field have been used to claim that the radiative feedback from accreting protostars (perhaps even the mechanical feedback from stellar outflows) is necessary to reproduce the observed IMF (Bate 2009, 2012, 2014; Krumholz et al. 2012). MHD simulations focusing on the effect of the magnetic field strength on star formation have included the radiative feedback as well, without exploring the isothermal case (Myers et al. 2014) and thus could not establish the relative importance of the magnetic field and radiation feedback in controlling the IMF. According to the turbulent fragmentation models (Padoan et al. 1997; Padoan & Nordlund 2002, 2011b; Hennebelle & Chabrier 2008, 2009; Hopkins 2012), the IMF originates primarily as the consequence of supersonic turbulence and can be reproduced under an isothermal gas approximation. Thus, it is important to carefully test the case of supersonic MHD turbulence of an isothermal gas, using as realistic as possible initial conditions and pursuing a large dynamic range in space and timescales in order to probe the whole IMF with a large statistical sample.

As mentioned above, in Padoan et al. (2014b) we already succeeded in reproducing the observed IMF with a simulation of isothermal MHD turbulence that represented a significant step forward in terms of the size of the simulated region, the length of the star formation time, and the number of sink particles. This model has since been refined with slightly better hydrodynamics and more optimal parameters for the sink particle model (Frimann et al. 2016b; Jensen & Haugbølle 2017). In this work, we carry out a more systematic study of the isothermal MHD case, using a large number of simulations to test the dependence of the results on the numerical parameters of the sink particle model (Sections 2 and 3 and Appendix C), to verify the numerical convergence of the IMF (see Section 4), and to test the predicted variability of the IMF due to variations of physical parameters (virial parameter or mean density) or resulting from the early time evolution of the IMF (see Section 5). We also use observed properties of molecular clouds (MCs) to constrain the environmental dependence of the IMF peak (Section 5.2) and stress that the IMF time evolution predicted by our turbulent fragmentation model (Padoan & Nordlund 2002, hereafter PN02), a consequence of the relatively long timescale of massive star formation (Padoan & Nordlund 2011b; Padoan et al. 2014b), is not necessarily predicted by other turbulent fragmentation models (Hennebelle & Chabrier 2008, 2009; Hopkins 2012), where massive stars result from the collapse of massive cores, in line with the scenario of massive star formation by McKee & Tan (2002, 2003). However, we do not see an accelerated accretion rate as the stars gain mass, so the results of our simulations are also at odds with the predictions of the competitive accretion scenario (Bonnell et al. 2001; Bonnell & Bate 2006).

This work focuses on star formation under physical conditions typical of Galactic MCs, i.e., turbulent regions of cold molecular gas following the observed velocity–size relation (with a large scatter). Different conditions in more extreme environments may require a careful consideration of detailed processes that are neglected here, such as radiation and mechanical feedbacks.

The simulations are carried out with a locally developed version of the public adaptive mesh refinement (AMR) code ramses (Teyssier 2002) that includes random turbulence driving and a robust algorithm for sink particles. Compared to the public version, it has been heavily modified to scale better on supercomputers with large number of cores per node, using an OpenMP/MPI hybrid parallelization, improved MPI load balancing, and with special attention to minimizing and bundling MPI communication. We have also improved the stability of the HLLD in high Mach number flows (see Appendix A), optimized the conjugate gradient method used for solving self-gravity, and improved the consistency of gravitational forces and the stability of sink particle orbits when using subcycling in time (see Appendix B). For the models considered in this paper we use periodic boundary conditions and an isothermal equation of state.

To initialize the turbulent state, we first run without self-gravity for ∼20 dynamical times, starting with uniform density and magnetic fields and a random acceleration with power only at wavenumbers 1 ≤ k ≤ 2 (k = 1 corresponds to the computational box size). The driving force keeps the rms sonic Mach number, ${{ \mathcal M }}_{{\rm{s}}}\equiv {\sigma }_{{\rm{v}}}/{c}_{{\rm{s}}}$, at an approximate constant value, where σ_v is the three-dimensional rms velocity and c_s is the isothermal speed of sound.

The driving is implemented as in Padoan & Nordlund (2004, 2011a) and Padoan et al. (2012). We use purely solenoidal driving, because of the large separation between the typical driving scale of supernovae (∼70 pc) and the size of the simulated region (4 pc). In Pan et al. (2016) we showed that, in supernova-driven turbulence, the compressive ratio (compressive over solenoid power) at MC scales is consistent with that obtained in the inertial range of randomly driven turbulence with purely solenoidal driving.

The refinement strategy is based on the overdensity. We refine the root grid when the density reaches a certain level, typically 10 times the average density, and then refine using a Truelove refinement criterion (Truelove et al. 1997), where the grid is refined every time the density increases by a factor of four. This refinement strategy is retained when self-gravity is turned on, though the number of AMR levels is increased to better resolve the gravitational collapse. Additionally, we expand the most refined grid so that it is at least eight cells across, to counteract the creation of many small island grids around density peaks. In contrast to Kritsuk et al. (2006) and Schmidt et al. (2009), we do not refine based on velocity or pressure gradients. However, we verified that the tails of the density probability density functions (pdf's) are sufficiently converged in the turbulence simulations used as initial conditions (see Figure 1).

Zoom In Zoom Out Reset image size

Gravitational collapse is central to simulations of star formation. As cloud cores in the turbulent flow become unstable and start to collapse, the density rises by many orders of magnitude. With current computational capabilities it is very challenging to capture the collapse to stellar densities even with very deep AMR hierarchies, and the correspondingly small time step would lead to prohibitively expensive simulations (though see Nordlund et al. 2014; Kuffmeier et al. 2017, for examples where stellar densities are almost reached). In this work, we are not interested in the detailed stellar physics, and to correctly capture the physics at larger scales of tens of astronomical units, we use sink particles to represent collapsed-gas regions. These should be formed robustly where the gas has unequivocally started to collapse gravitationally, but the conversion from gas in a grid cell to a sink particle representation cannot happen at an arbitrarily high density. An isothermal gas will fragment while collapsing, and if the Jeans length is not sufficiently resolved, the gas can undergo numerical fragmentation (Truelove et al. 1997, 1998). The normalized Jeans length at a given density and numerical resolution, the Jeans number, is

$$\begin{eqnarray}&&{L}_{{\rm{J}}}=\displaystyle \frac{{\lambda }_{{\rm{J}}}}{{\rm{\Delta }}x}=\sqrt{\displaystyle \frac{{c}_{{\rm{s}}}^{2}\pi }{G\rho }}\displaystyle \frac{1}{{\rm{\Delta }}x},\end{eqnarray} \tag{ 1 }$$

and it has been argued that one should have L_J ≥ 4 (Truelove et al. 1997) everywhere.

To detect gravitational collapse in the gas and to convert the gas into a sink particle, we use a number of criteria:

• 1.  
  Sink particles can only be created in cells where the gas density is above a threshold value ρ_s. We generally make sure that this only happens at the highest AMR level and that the Jeans length is resolved. In most of the presented runs the Jeans length is resolved with at least two cells at ρ_s (L_J,s = 2), though tests have shown that even if the Jeans length is slightly under resolved we do not generate artificially a higher amount of sinks (in accordance with Gong & Ostriker 2013). This is likely because even if the gas has started to fragment a bit, the different fragments will afterward be accreted by the sink particle. A typical value for ρ_s in the runs presented here is 10⁵–10⁶ times the average density.
• 2.  
  The gravitational potential is required to have a local minimum at the cell where the sink particle will be created. This is evaluated by smoothing with a 2 × 2 × 2 average and comparing with the potential in the 26 neighboring smoothed cells.
• 3.  
  The velocity field has to be converging in the cell, ∇ · v < 0.
• 4.  
  No other previously created sink particle can be present within an exclusion radius, r_ex, of the cell where the new particle is created.
• 5.  
  In runs with an active energy equation, we also disallow creation of sinks from very warm gas. This can otherwise spur the creation of sink particles in, e.g., the high-density expanding shell of a supernova.

This sink particle recipe has already been used in a number of papers (Padoan et al. 2012, 2014b, 2016a, 2016b, 2017; Vasileiadis et al. 2013; Nordlund et al. 2014; Frimann et al. 2016a, 2016b; Kuffmeier et al. 2016, 2017; Pan et al. 2016; Jensen & Haugbølle 2017). In the past, we had used a simpler recipe based only on density criteria in lower-dynamic-range simulations aimed at deriving the SFR (Padoan & Nordlund 2011a). Due to the much larger resolution in following works and the added goal of predicting robustly the stellar IMF, we have switched to using all the above criteria. The individual conditions are similar to what has been used by Federrath et al. (2010), Bate et al. (1995), Krumholz et al. (2004), Gong & Ostriker (2013), and Bleuler & Teyssier (2014). A major difference between our models and many of the cloud-scale models in the literature is the very deep AMR grids we can afford. This makes it possible to use a high-density threshold, ρ_s, while still resolving the Jeans length and ensures that the gas is already undergoing a physical collapse, at the point where the sink particle is created.

A sink particle is born without any mass but will immediately acquire mass by accretion. Sink particles can accrete gas from cells that are within the accretion radius, r_acc, if the gas density is higher than a threshold value, ρ_acc. In the runs presented below we are using two different accretion recipes. A simple recipe is used for the test runs, while a more complicated recipe, capable of tracing the accretion down to very low densities, is used for the convergence runs. The stellar masses obtained with the two recipes are practically identical (see Figure 14 in Appendix C), and the only difference is that the second recipe allows us to get a more precise picture of the accretion history at accretion rates down to 10⁻⁹ M_⊙ yr⁻¹.

In the simple recipe, the amount of gas accreted is such that the density in the cell is left to be just below ρ_acc. Typically, we set ρ_acc to be ρ_s/2. The sink particles only accrete mass, momentum, and, if present, passive scalars. Thermal energy is removed from the gas in proportion to the amount of mass removed, while magnetic field cannot be deposited on the sink particles. To avoid a pileup of high-density gas just outside the accretion radius, r_acc, it is necessary to have ρ_acc ≤ ρ_s. The accretion distance is typically set to one to two Jeans' lengths at r_acc around the sink particle. This simple accretion recipe of removing gas to bring it below a critical density is similar to what we have used in the past (Padoan & Nordlund 2011a; Padoan et al. 2012, 2014b; Vasileiadis et al. 2013; Nordlund et al. 2014) and what is used by Federrath et al. (2010). At the very high densities and small r_acc we are considering, due to the deep AMR hierarchy, the Bondi–Hoyle accretion radius is larger than r_acc for most of the lifetime of a sink particle, and correct accretion rates are therefore naturally enforced, and we do not need the more complicated accretion recipes employed by other groups (e.g., Krumholz et al. 2004; Gong & Ostriker 2013).

In the more complicated recipe, we set ρ_acc to a much lower value, allowing for accretion from gas at lower densities, but we also add the condition that the relative speed between the sink particle and the cell is lower than $\sqrt{2}$ the Kepler velocity, to avoid artificial accretion of unbound gas. We also taper off the efficiency of the accretion depending on how bound the gas in a cell is to the sink particle. In a simple picture where a sink particle has an effective geometric cross section of πR², the accretion rate should be ${\dot{m}}_{\mathrm{gas}}=\pi {R}^{2}\rho v$, where v is the relative speed between the gas and the sink particle. This formula is our starting point, but we modify it to take into account that unless R is the grid spacing in our code this would correspond to accretion at a distance. Let d be the distance between the center of the cell and the nearest sink particle and v_K = (Gm_sink/d)^1/2 be the Kepler velocity. The mass accretion rate in a single time step Δt from a cell with density ρ is then

$$\begin{eqnarray}{\dot{m}}_{\mathrm{gas}}=\displaystyle \frac{{\rm{\Delta }}\rho {\rm{\Delta }}V}{{\rm{\Delta }}t}=\left\{\begin{array}{ll}{\alpha }_{\mathrm{rate}}\,\rho {v}_{K}{\rm{\Delta }}{x}^{2}{f}_{v} & \mathrm{if}\,\rho \leqslant {\rho }_{\mathrm{th}}\\ 0.5\rho {\rm{\Delta }}V/{\rm{\Delta }}t & \mathrm{if}\,\rho \gt {\rho }_{\mathrm{th}}\end{array}\right.,\end{eqnarray} \tag{ 2 }$$

where Δρ is the change in the gas density in the cell volume ΔV, and ρ_th is a threshold density that we normally set to 2ρ_s as a safety valve to avoid a large amount of gas piling up faster than it is accreted inside the accretion radius. α_rate = 0.2 controls how efficient accretion proceeds compared to rotation. In the unrealistic case of spherical accretion with zero rotation α_rate = 1. f_v is a tapering function that limits the rate depending on how far the cell is from the sink and how large the relative speed is,

$$\begin{eqnarray}{f}_{v}=\left[1-{\left(\displaystyle \frac{d}{{r}_{\mathrm{acc}}}\right)}^{2}\right]\times \left\{\begin{array}{ll}0 & \mathrm{if}\,v\geqslant \sqrt{2}{v}_{K}\\ 2-{\left(\tfrac{v}{{v}_{K}}\right)}^{2} & \mathrm{if}\,{v}_{K}\lt v\lt \sqrt{2}{v}_{K}\\ 1 & \mathrm{if}\,v\leqslant {v}_{K}.\end{array}\right.\end{eqnarray} \tag{ 3 }$$

An optimal accretion algorithm should work as a transparent boundary condition such that no discontinuity is seen across the accretion radius. We have tested the formula both in the low-resolution limit, where a sink particle is streaming through a low-density gas, and in the opposite limit with very high resolution, where the sink is accreting from an accretion disk. In both cases, the above numerical parameters, α_rate = 0.2 and ρ_th = 2 ρ_s, result in a smooth transition between gas inside and outside the accretion radius, r_acc.

A large fraction of the mass in protostellar systems is lost through winds and jets (Matzner & McKee 2000; Alves et al. 2007; Könyves et al. 2010), launched from small scales not resolved in the current simulations. To account for this mass loss, we apply an efficiency factor, _acc = 0.5, to the accreted mass and momentum of the sink particles in the simulation, so that only ${\epsilon }_{\mathrm{acc}}{\dot{m}}_{\mathrm{gas}}{dt}$ of the mass is accreted to a given sink particle in a time interval dt (the nonaccreted gas fraction is simply removed from the simulation). Without a more detailed knowledge of the typical mass loss in outflows, we prefer to stay with a single value of _acc, rather than expand the parameter space with yet another dimension. In an older version of the code, used for some of the smaller runs in Appendix C, we accreted the full amount of mass (_acc = 1), and the accretion rates and sink particle masses were readjusted after the simulation had finished running. Since this overestimates the gravitational pull of the individual stars, one cannot simply rescale after the fact with the same factor. Two otherwise-identical test runs have shown that to rescale from a run with _acc = 1 to one with _acc = 0.5 the appropriate factor is 2.4 (see Figure 14 in Appendix C). As discussed by Kuffmeier et al. (2017) in the context of deep zoom-in simulations and by Federrath et al. (2014) in the context of jet and outflow models as a function of resolution, _acc depends on the maximum resolution in the model. In our case with a maximum resolution of 50 au for the high model, we do not resolve disks, except for the most massive stars, and we are therefore not accounting for outflows. Thus, it is appropriate to use _acc = 0.5 in all runs.

This paper is based on a large suite of simulations that describe the evolution of a generic MC piece of size L_box = 4 pc, using isothermal supersonic turbulence, self-gravity, magnetic field, and a subgrid sink particle model for the gravitational collapse of the gas. The simulations are characterized by three nondimensional numbers: the sonic Mach number, ${{ \mathcal M }}_{{\rm{s}}}={\sigma }_{{\rm{v}}}/{c}_{{\rm{s}}}$, the Alfvénic Mach number, ${{ \mathcal M }}_{{\rm{a}}}={\sigma }_{{\rm{v}}}/{v}_{{\rm{a}}}$, and the virial number, ${\alpha }_{\mathrm{vir}}=5{\sigma }_{{\rm{v}}}^{2}R/(3{GM})$, which measure the relative strength of kinetic energy to thermal, magnetic, and gravitational energies, respectively. In the expression for the Alfvénic Mach number, v_a is the Alfvén speed corresponding to the mean magnetic field and the mean density in the computational domain. R is a characteristic size we use to define the virial number of a simulation, which we take to be R = L_box/2, and M is the total mass in the simulation, M = M_box.

The basic parameters of our simulations cannot be chosen arbitrarily. They are constrained by observations and the physical conditions in the interstellar medium (ISM), and by what is feasible with current computational capabilities. At large scales of hundreds of parsecs, the ISM has a complex thermal structure. On parsec scales, high-density cold molecular gas is shielded from UV radiation and cooled by dust and atomic lines balanced by cosmic-ray and shock heating, forming an effectively isothermal medium from where stars are born. Thus, to be realistic, our model should have a size smaller than ∼10 pc. On the other hand, because we aim at studying the emergence of the IMF including the high-mass Salpeter range, the model should contain at least several thousand solar masses of gas. Numerically, the resolution needed to properly characterize the density fluctuations induced by supersonic turbulence limits the sonic Mach number corresponding to the rms velocity and therefore the size of the box we can consider (because of the velocity–size relation followed by MCs (Larson 1981), albeit with a significant scatter). Another constraint is the desire to resolve small enough scales that we marginally resolve large disks around the sink particles. In summary, the approximation of an isothermal gas sets an ultimate upper limit for a realistic box size, while numerical constraints set a more stringent limit. The requirement of sufficient mass sets a lower limit.

Taking the above constraints into consideration, we have chosen ${{ \mathcal M }}_{{\rm{s}}}=10$, ${{ \mathcal M }}_{{\rm{a}}}=5$, and α_vir = 0.83 for the reference simulations. To convert to physical units, we assume an isothermal sound speed of 0.18 km s⁻¹, corresponding to a temperature of T ≈ 10 K, and a mean molecular weight μ = 2.37, appropriate for cold MCs. Given the nondimensional parameters, the rest depends on the physical box size. Choosing L_box = 4 pc (in line with a characteristic velocity–size relation in MCs) corresponds to setting a total mass M_box = 3000 M_⊙, a mean density of 795 cm⁻³, and a mean magnetic field strength of 7.2 μG. The resulting mean column density is consistent with MC observations (e.g., Heyer et al. 2001; Roman-Duval et al. 2010), and the magnetic field strength with OH Zeeman splitting measurements (e.g., Crutcher 2012), once density and magnetic fluctuations arising in the super-Alfvénic turbulence are properly modeled (Lunttila et al. 2008, 2009). This gives a free-fall time, ${t}_{\mathrm{ff}}=\sqrt{3\pi /32G{\rho }_{0}}$, of 1.18 Myr and a dynamical or crossing time of t_dyn = 1.08 Myr. The virial parameter is known to control the SFR (Krumholz & McKee 2005; Padoan et al. 2012, 2014a, 2017) and is expected to influence the peak of the IMF (see Section 5.2). Thus, we complement our base model, high, with three additional models, light, heavy, and massive, with M_box = 1500 M_⊙, M_box = 6000 M_⊙, and M_box = 12,000 M_⊙ respectively.

The main properties of the five convergence models with a total mass of 3000 M_⊙ and the three models light, heavy, and massive are listed in Table 1. All models use six levels of refinement on top of the root grid. These are the reference models discussed in the main sections of the paper. An additional 36 models are discussed in Appendix C, where we explore the model dependence on the numerical parameters of our sink particle implementation. Computationally, the project has required more than 50 million CPUh at four different HPC centers and has resulted in ∼100 TB of data.

Table 1.  Numerical Parameters of the Large-scale Runs Exploring Convergence and the Dependence on α_vir

	Run Parameters										Creation of Sinks				Accretion to Sinks
Run	Root	NAMR	Δx	LJ	ρref	Mbox	αvir	tend	SFE	Nsink	LJ,s	ρs	ρs	rex	ρacc	ρacc	ρth	racc	acc
	Grid		au		$\langle \rho \rangle$	M⊙		Myr				cm−3	$\langle \rho \rangle$	Δx	cm−3	$\langle \rho \rangle$	ρs	Δx
16	163	6	800	2.0	2	3000	0.83	1.6	13%	108	2	6.6 × 106	8.3 × 103	8	4227	5.3	2	4	0.5
32	323	6	400	2.5	5	3000	0.83	1.8	13%	169	2	2.6 × 107	3.3 × 104	8	4227	5.3	2	4	0.5
low	643	6	200	3.6	10	3000	0.83	2.4	13%	279	2	1.1 × 108	1.3 × 105	8	4227	5.3	2	4	0.5
med	1283	6	100	7.2	10	3000	0.83	2.5	13%	363	2	4.2 × 108	5.3 × 105	8	4227	5.3	2	4	0.5
high	2563	6	50	14.4	10	3000	0.83	2.5	13%	410	2	1.7 × 109	2.1 × 106	8	4227	5.3	2	4	0.5
light	2563	6	50	14.4	20	1500	1.67	4.0	5%	86	2	1.7 × 109	4.3 × 106	8	4227	5.3	2	4	0.5
heavy	2563	6	50	14.4	5	6000	0.42	0.7	5%	614	2	1.7 × 109	1.1 × 106	8	4227	5.3	2	4	0.5
massive	2563	6	50	14.4	2.5	12,000	0.21	0.3	3%	1223	2	1.7 × 109	5.3 × 106	8	4227	5.3	2	4	0.5

Download table as:  ASCIITypeset image

The results from running the numerical tests discussed in Appendix C have guided us toward a definitive set of parameters for the sink particle model. We require a density threshold for sink particle creation corresponding to resolving the Jeans length with at least two grid cells, L_J,s = 2. Very close to a sink particle the flow will be artificially disturbed and we require an exclusion radius of r_ex = 8 cells to avoid the creation of spurious sink particles, the occurrence of which can be tested with a very powerful method, the nearest-neighbor histogram (see Appendix C.5). In addition, it is important that the density at sink formation ρ_s is significantly above the highest densities reached by the turbulence alone, which sets a minimum requirement of ≈10⁵ for the gas overdensity at the highest AMR level, which in this case requires six levels of refinement and a gas density of ≈10⁸ cm⁻³ (in our reference model high, we actually have ρ_s = 1.7 × 10⁹ cm⁻³ and a minimum cell size of Δx = 50 au). Finally, the model needs to have a high enough resolution that small-scale high-density cores, which are the progenitors of low-mass stars, are resolved and can collapse; otherwise, a numerical IMF turnover would be created, dictated by the resolution. This requires both that the turbulence is sufficiently well resolved to sample well the high-density tail of the gas density pdf and that the Jeans length is sufficiently resolved to allow the gravitational collapse and suppress numerical fragmentation (Truelove et al. 1997).

To test the numerical convergence of the IMF, we have carried out five runs, 16, 32, low, med, and high, with a varying root grid size of 16³, 32³, 64³, 128³, and 256³, respectively, and six levels of refinement (reaching a minimum cell size Δx = 50 au in the reference simulation high). The corresponding minimum numerical values of the Jeans length (in units of Δx of each simulation) are L_J = 2, 2.5, 3.6, 7.2, and 14.4, while the parameters in the sink particle model are identical in all five runs (see Table 1). Notice that instead of using the same overdensity threshold ${\rho }_{\mathrm{ref}}=10\langle \rho \rangle$ for refinement for all the runs, giving a factor of two increase in L_J between runs when the spatial resolution is doubled, we have adopted a slightly more generous refinement strategy (lower density thresholds) for the runs 16 and 32 to avoid resolving the Jeans length with less than two cells (see column with ρ_ref in Table 1). However, the threshold density for sink creation, ρ_s, does increase by a factor of four between runs with increasing resolution by a factor of two, in order to keep the numerical value of the Jeans length at the density of sink creation, L_J,s, independent of resolution (as all other sink particle model parameters), L_J,s = 2. In these simulations, we have also imposed a small maximum time step size of 80 days at the highest level of refinement. The small time step size facilitates an accurate integration of sink particle orbits over the long integration time of most of these runs.

To illustrate the general evolution of the simulations, in Figure 2 we show the column density in the high run at different times. At the time when gravity is turned on, the density field is well mixed and well described by a lognormal pdf (see Figure 1). At later times, the densest gas decouples from the turbulent flow and starts to collapse. Some effect of self-gravity can also be seen at larger scales, such as the formation of a large (∼1 pc), elongated dense region that has assembled most of the high-density gas and most of the sink particles by the end of the simulation. This region is itself embedded in a looser, but still discernible, overdense structure that stretches the full width of the computational volume (see the bottom right panel of Figure 2).

Zoom In Zoom Out Reset image size

The star formation efficiency (SFE) is the fraction of mass turned into stars at a given time t: SFE(t) = M_sink(t)/M_box, where M_sink(t) is the total mass in sink particles at time t. The SFR is the time derivative of SFE and can be expressed in nondimensional units if multiplied by a characteristic timescale. We follow the convention introduced by Krumholz & McKee (2005) to adopt the free-fall time of the mean density as the characteristic time, so the nondimensional SFR is defined as SFR_ff(t) = dSFE(t)/d (t/t_ff). Figure 3 shows the time evolution of SFE and SFR_ff in the convergence runs. The time dependence is well fitted by a power law, SFR_ff ≈ 0.06 (t/t_ff)^0.5, a slow increase of SFR_ff from 0.04 at 0.5 Myr (0.4 t_ff) to 0.08 at 2.4 Myr (2 t_ff). This increase in SFR_ff is probably related to the formation of a dominant dense cluster toward the end of the run (see Figure 2). The increasing density in the cluster-forming region decreases the local virial parameter, increasing the SFR (Padoan et al. 2012, 2017).

Zoom In Zoom Out Reset image size

We do not expect this trend of increasing SFR_ff to continue for a much longer time under more realistic conditions. To keep the evolution realistic for a longer timescale, past SFE = 0.13, external forcing from scales beyond the 4 pc box size would be needed to properly account for the interaction of larger-scale turbulence with the cluster-forming region, as observed in larger-scale simulations with supernova driving (Padoan et al. 2017). With supernova driving, the typical disruption time of a MC is ≈2 t_dyn (Padoan et al. 2016b), equivalent to ∼2.5 Myr in the simulations of this work, meaning that external feedback from larger scales should play an important role in the evolution of a 4 pc region after approximately 2 Myr. Protostellar feedback from the high-mass stars in the dense cluster could also become important as a local feedback agent at that stage, but it is neglected here.

Compared to the high run, SFR_ff is converged in the med run and nearly converged in the low run, while at lower resolutions the runs are clearly not converged, and possibly also affected by numerical fragmentation, given the low minimum Jeans number of L_J = 2. There is a period between 0.7t_ff and 1.2t_ff where the low run has a higher SFR_ff, but otherwise the three runs low, med, and high show converged SFR_ff. This result is in agreement with previous studies of the SFR (e.g., Padoan & Nordlund 2011a; Padoan et al. 2012; Federrath & Klessen 2012), where it was found that SFR_ff converges at a relatively low numerical resolution, which allowed those works to explore a broad parameter space with many relatively low resolution simulations. On the other hand, the numerical convergence of the stellar IMF is much more demanding and, in our opinion, has never been convincingly achieved. It has not even been tried, so far, in the case of isothermal MHD turbulence, which is the main goal of this work.

The convergence test of our numerical IMFs is shown in Figures 4 and 5. It has been carried out for a single snapshot of each simulation, near the end of each run, to take advantage of the highest value reached by SFE and thus the larger statistical sample. All five runs are compared at SFE = 0.13, corresponding to a time of 1.61, 1.91, 2.46, 2.62, and 2.64 Myr after the formation of the first star, in order of increasing resolution (SFE_ff is higher in the low-resolution runs, so the same value of SFE is reached in a shorter time). The IMF histograms computed at those times are plotted in Figure 4. To avoid the confusion generated by noisy overlapping histograms, the IMFs have been vertically shifted, except in the case of the reference run high, where N(m) is indeed the number of sink particles in each logarithmic mass interval of the histogram.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The figure shows a clear shift of m_peak toward lower values as the resolution increases, although the shift is quite small between the two highest-resolution runs. To quantify the result of this comparison, we have fitted all the IMFs with a lognormal function for sink masses m ≲ 2 M_⊙ (dotted lines), and with a single power law for masses m > 2 M_⊙ (solid straight lines). Both models clearly provide very good representations, in their respective mass ranges, of the shape of the IMFs. In the case of the reference run high, the IMF is complete over more than three orders of magnitude in mass, with two orders of magnitude covered by the lognormal fit, and over one order of magnitude by the power-law fit. The smallest mass bin where the IMF is assumed to be complete is taken to be the one just above a sharp drop in N(m), at approximately 0.02 M_⊙ in the case of the run high (and larger values by approximately a factor of two for each consecutive step of decreasing resolution). The cutoff is much more apparent when the histograms are computed with narrower mass bins, which is what we have done in order to determine these approximate completeness limits.

The dashed line in Figure 4 shows the Chabrier IMF for single stars (Chabrier 2005) below 2 M_⊙, continued by the Salpeter IMF (Salpeter 1955) at larger masses. It fits almost exactly our highest-resolution IMF, except that the Salpeter range is slightly higher in the case of the simulation, partly because the turnover region has a slightly larger width than in the Chabrier IMF. The values of the best-fit parameters are plotted in Figure 5, showing that all IMFs above 16 have a width, σ_m, consistent with that of the Chabrier IMF, except for the slight increase in the case of the run high, and a power-law slope, Γ, consistent with the Salpeter value. In the case of the lowest-resolution run, 16, σ_m cannot be measured and Γ is slightly smaller than Salpeter's value. The IMF peak, m_peak, shows a clear dependence on resolution (top panel of Figure 6), but also a trend toward numerical convergence (which may or may not have been already reached at the resolution of the run high).

Zoom In Zoom Out Reset image size

An alternative way to test for numerical convergence of the IMFs with increasing resolution is to use the cumulative IMFs. Although the value of m_peak is best defined by a lognormal fit to the IMF, as we did above, the procedure has some dependence on the choice of bin size and location, while the cumulative IMF is immune to such choices. In Figure 7, we show the cumulative IMFs of the convergence-test simulations at SFE = 0.13, expressed as the number of sink particles above the mass m as a function of m. The figure shows a clear tendency toward numerical convergence. The two highest-resolution runs, med and high, have essentially the same cumulative IMF down to a mass of order 0.1 M_⊙, lower than the IMF peak. The rate of convergence is illustrated in Figure 6, where we plot the number of stars above a given mass, m_min, as a function of the root grid size of the simulation. Even in the case of the total number of stars, m_min = 0 M_⊙, the cumulative IMFs are clearly converging with increasing resolution, consistent with the convergence of m_peak in the top panel of Figure 5. For stars with m > 0.1 M_⊙, the convergence is achieved at the resolution of the run med, and all runs have essentially the same number of intermediate- and high-mass stars (m_min = 1.0 M_⊙).

Zoom In Zoom Out Reset image size

Based on the above convergence tests, we can conclude that we have found, for the first time, clear evidence of the convergence of the IMF turnover (with a nearly converged value of m_peak in agreement with the observations) in the case of isothermal MHD turbulence, as predicted by the turbulent fragmentation models of the IMF (PN02; Padoan et al. 1997; Hennebelle & Chabrier 2008, 2009; Padoan & Nordlund 2011b; Hopkins 2012).

The universality of the stellar IMF is hotly debated. While most works emphasize the apparent invariance of the IMF (e.g., Chabrier 2005; Bastian et al. 2010; Massey 2011; Offner et al. 2014; Weisz et al. 2015), some stress compelling evidence of IMF variability (e.g., Kroupa 2001; Dib et al. 2010, 2017; Marks et al. 2012; Kroupa et al. 2013; Scholz et al. 2013; Dib 2014). Observational determinations of the stellar IMF have been approximated with different models, such as a multicomponent power-law function (Kroupa 2001, 2002), a tapered power-law function (Parravano et al. 2011), or a lognormal function below 1 M_⊙ continued by a power law at larger masses (Chabrier 2005). The slope of the power law at large masses is usually assumed to be Γ = 1.35, as first estimated by Salpeter (1955), or slightly shallower (e.g., Γ = 1.3 in Kroupa's and Chabrier's IMFs). However, in the most complete and homogeneous study to date, based on 85 resolved clusters in M31, the IMF is actually found to be somewhat steeper, with ${\rm{\Gamma }}={1.45}_{-0.06}^{+0.03}$ (Weisz et al. 2015). Besides this well-defined mean value, the value of Γ exhibits variations from cluster to cluster that are found to be within the error bars, with only a few outliers (e.g., Weisz et al. 2015), or interpreted as indications of intrinsic IMF variations (e.g., Dib et al. 2017). The IMF peak is found to be at approximately 0.2 M_⊙, but statistically significant variations from cluster to cluster are suggested by Dib (2014). These variations, if confirmed, may reflect both the environment and the age of the observed stellar populations.

Physical models for the origin of the IMF predict a dependence on the average physical parameters of the star-forming environment, which may in principle result in a larger IMF variability than observed. One may explore extra processes that would keep the predicted IMF invariant, such as the radiative feedback from the accretion luminosity of protostars (Krumholz et al. 2016), which appears to be crucial in simulations neglecting the magnetic field (e.g., Bate 2012; Krumholz et al. 2012). However, it is also possible that the variability predicted by the theory does not violate the observational constraints. In this section, we consider the turbulent fragmentation model by PN02 and use our simulations to test its prediction for the dependence of the IMF turnover on physical parameters. We then apply the model to the physical parameters of MCs derived from large MC surveys, showing that the predicted IMF variations are within the observational constraints. The PN02 model also implies an early time evolution of the IMF, which we show to be qualitatively confirmed by the simulations.

5.1. The IMF Turnover from Turbulent Fragmentation

The origin of the characteristic stellar mass, essentially the turnover and peak of the IMF, is arguably the most fundamental question in star formation. Padoan et al. (1997) proposed that the IMF turnover is the direct result of the pdf of gas density in a turbulent MC and modeled the turnover as a probability distribution of Jeans masses in the isothermal gas with a lognormal density pdf. Although that early model did not account for the power-law tail of the IMF at large masses, its original explanation for the origin of the IMF turnover has been essentially retained in following turbulent fragmentation models that can predict the full IMF (PN02; Hennebelle & Chabrier 2008, 2009; Hopkins 2012).

A numerical derivation of the turnover mass from the PN02 model yields ${m}_{\mathrm{peak}}\sim {M}_{\mathrm{BE},0}{{{ \mathcal M }}_{{\rm{s}}}}^{-1.1}$, where ${M}_{\mathrm{BE},0}$ is the Bonnor–Ebert (BE) mass with the external density equal to the average density of the star-forming region and ${{ \mathcal M }}_{{\rm{s}}}$ is the rms Mach number of the turbulent flow (Equation (7) in Padoan et al. 2007). This result can be easily derived as the characteristic BE mass in the turbulent flow, meaning the BE mass with external density equal to the characteristic post-shock density, or, equivalently, the BE mass with external pressure equal to the characteristic dynamic pressure of the turbulent flow. The standard BE mass confined by a thermal pressure ${P}_{\mathrm{th},0}$ is given by

$$\begin{eqnarray}&&{M}_{\mathrm{BE}}=1.182\displaystyle \frac{{\sigma }_{\mathrm{th}}^{4}}{{G}^{3/2}{P}_{\mathrm{th},0}^{1/2}}.\end{eqnarray} \tag{ 4 }$$

Including the dynamic pressure of the turbulence, the external pressure is given by ${P}_{0}={P}_{\mathrm{th},0}+{P}_{\mathrm{dyn},0}={P}_{\mathrm{th},0}(1+{{{ \mathcal M }}_{{\rm{s}}}}^{2})$. Substituting into the previous expression of M_BE, we get a modified turbulent BE mass:

$$\begin{eqnarray}&&{M}_{\mathrm{BE},{\rm{t}}}\approx \displaystyle \frac{1.182\,{\sigma }_{\mathrm{th}}^{4}}{{G}^{3/2}{P}_{0}^{1/2}}=\displaystyle \frac{{M}_{\mathrm{BE},0}}{{(1+{{{ \mathcal M }}_{{\rm{s}}}}^{2})}^{1/2}}\approx \displaystyle \frac{{M}_{\mathrm{BE},0}}{{{ \mathcal M }}_{{\rm{s}}}},\end{eqnarray} \tag{ 5 }$$

which is a good approximation to the turnover mass in the turbulent fragmentation models mentioned above, providing an intuitive explanation of the origin of the IMF peak. To test the validity of this prediction, we express the IMF peak as

$$\begin{eqnarray}&&{m}_{\mathrm{peak}}\equiv {\epsilon }_{\mathrm{BE}}\,{M}_{\mathrm{BE},{\rm{t}}},\end{eqnarray} \tag{ 6 }$$

where _BE is a local efficiency parameter analogous to _acc in the sink particle accretion model, and use the simulations to verify whether it provides a good fit to the numerical IMFs.

For this purpose, we use the four simulations light, high, heavy, and massive with a root grid of 256³ cells and six AMR levels, with four different values of the virial parameter (see Table 1). The virial parameter is varied by leaving the rms velocity constant and increasing or decreasing the mean density (total mass) in the computational volume by a factor of two or four relative to the reference run high (see Section 3 and Table 1). The overdensity threshold at which the root grid is refined is changed from ${\rho }_{\mathrm{ref}}=10\langle \rho \rangle$ in run high to ${\rho }_{\mathrm{ref}}=20\langle \rho \rangle ,5\langle \rho \rangle ,\,\mathrm{and}\,2.5\langle \rho \rangle$ in light, heavy, and massive, respectively, to keep the minimum Jeans number constant, at 14.4. The IMFs from these four simulations are shown in Figure 8, where the histograms are shifted vertically by a factor of four between consecutive runs, except for the top histogram, to minimize the confusion of overlapping plots. The IMFs are all sampled at SFE = 0.024, corresponding to a time of 2.07, 0.83, 0.46, and 0.23 Myr after the formation of the first star, for the runs light, high, heavy, and massive, respectively. We have chosen a rather low SFE for this comparison because the run light has a very low SFR_ff, such that to reach a much higher SFE the simulation should be integrated for much longer than 2 Myr. As commented above, on a scale of 4 pc the influence of larger-scale feedbacks should become quite significant after approximately 2 Myr, making this idealized setup driven by a random force somewhat questionable at later times. Despite the short timescale of the higher α_vir runs at SFE = 0.024, we have found that the value of m_peak (and the ratios of its values from different runs) is already reasonably stable to allow this comparison.

Zoom In Zoom Out Reset image size

The dotted lines in Figure 8 are lognormal fits of the IMFs (the power-law fit at larger masses is not possible in this case because the high-mass tail is not developed yet at this early time). The lowest-mass bin for the fit is based on the approximate IMF completeness limit judged as in the numerical convergence test, while the highest-mass bin is approximately 10 × m_peak, assuming that the beginning of the power-law tail is also shifted to higher masses as the mean density decreases. The IMF peak clearly shifts toward smaller values as the mean density increases, as predicted by the isothermal turbulent fragmentation model of the IMF. The best-fit values of the lognormal peaks are shown in Figure 9, plotted as a function of the α_vir value of each run (a proxy for the inverse of the mean density at a fixed rms velocity and size). The prediction of Equation (6) is shown by the open circles, after normalizing the relation by the measured value of m_peak in the run high. The normalization corresponds to the choice _BE = 0.64, quite close to the related local efficiency parameter set in the sink particle accretion model, _acc = 0.5.

Zoom In Zoom Out Reset image size

Figure 9 shows that the measured variation of m_peak with the mean density is approximately consistent with the prediction of Equation (6). Although the slight discrepancy in the case of the run massive may seem significant, it is not significant if one takes into account the uncertainty in the measured value for the run high. Furthermore, because we have established that the value of m_peak in the run high may not be fully converged (see Figure 5), it is also possible that the value of m_peak in the run massive is even less converged, as the total mass in this run is larger and the peak smaller than in the run high. The increasingly higher lack of numerical convergence with increasing mean density could then explain the observed deviation from the prediction of Equation (6).

We set the system rms velocity assuming a temperature of 10 K and the system size (or total mass) based on a standard Larson velocity–size relation (see Section 3). If we chose not to follow the observed Larson relation, both the rms velocity and the size (or total mass) of the system could be rescaled, as long as the nondimensional parameters of the simulation, ${{ \mathcal M }}_{{\rm{s}}}$ and α_vir, were not changed. Thus, one may suspect that the predicted IMF peak is consistent with the numerical IMFs only for specific values of gas temperature or system size, but it can be easily shown that this agreement is immune to the rescaling of the simulation. The virial parameter can be expressed as

$$\begin{eqnarray}&&{\alpha }_{\mathrm{vir}}\propto {{{ \mathcal M }}_{{\rm{s}}}}^{2}{{TM}}_{\mathrm{tot}}^{-2/3}{\rho }_{0}^{-1/3}.\end{eqnarray} \tag{ 7 }$$

Because both ${{ \mathcal M }}_{{\rm{s}}}$ and α_vir are fixed in the simulation, the mass can only be scaled according to ${M}_{\mathrm{tot}}\propto {T}^{3/2}{\rho }_{0}^{-1/2}\propto {M}_{\mathrm{BE},0}$. This shows that imposing a value for both ${{ \mathcal M }}_{{\rm{s}}}$ and α_vir in the simulation implies a fixed value of the ratio ${M}_{\mathrm{tot}}/{M}_{\mathrm{BE},0}$, and thus a fixed value of ${M}_{\mathrm{tot}}/{M}_{\mathrm{BE},{\rm{t}}}$. Thus, rescaling the temperature or size (or total mass) of the system does not affect our comparison of the predicted IMF peak with the IMF peak from the simulations.

To fully test the prediction of the turbulent fragmentation model with respect to the IMF turnover (and the width of the IMF as well), we should also consider the dependence of m_peak on the sonic and Alfvénic rms Mach number. Because all the simulations of this work have the same Mach number, this important test will be addressed in a separate study.

5.2. Variability of the IMF Turnover with Environment

The theoretical and numerical prediction that the IMF peak scales with ${M}_{\mathrm{BE},{\rm{t}}}$ implies an environmental dependence of the IMF. In previous works, we have already stressed that if the virial parameter does not vary significantly in star-forming regions, and assuming standard velocity–size and mass-size relations, the predicted IMF peak should have only mild variations (Padoan et al. 2007). Here, we try to quantify the expected scatter of m_peak based on the scatter in the observed properties of star-forming regions. We can express m_peak as a function of the nondimensional parameters of the simulation and the total mass:

$$\begin{eqnarray}&&{m}_{\mathrm{peak}}\approx 1.124\,{M}_{\mathrm{tot}}{{{ \mathcal M }}_{{\rm{s}}}}^{-4}{\alpha }_{\mathrm{vir}}^{3/2},\end{eqnarray} \tag{ 8 }$$

which shows that for constant α_vir and for standard Larson relations, M_tot ∝ L² and σ_v ∝ L^1/2, m_peak is constant. However, observed MCs have a range of values of α_vir and yield Larson relations with a significant scatter and with exponents in general different from those standard values. Thus, our IMF model should predict non-negligible IMF peak variations from cloud to cloud.

In order to quantify the observational scatter in m_peak predicted by the model, we consider two of the largest Galactic MC samples available: the MC catalog by Heyer et al. (2001), extracted from a decomposition of the ¹²CO FCRAO Outer Galaxy Survey (Heyer et al. 1998), and the MC catalog by Roman-Duval et al. (2010), extracted from the UB–FCRAO Galactic Ring Survey (Jackson et al. 2006). To limit the distance and mass uncertainties, Heyer et al. (2001) consider only MCs with circular velocity <−20 km s⁻¹, which yield a sample of 3901 clouds. Roman-Duval et al. (2010) provide an estimate of the error in the mass determination of each of the 750 MCs in their catalog. We select a subsample of their clouds with a mass error <20%, in order to minimize the scatter in cloud properties due to observational errors instead of intrinsic cloud differences. Finally, we retain only MCs with mass >10³ M_⊙ (smaller clouds would not yield a well-sampled IMF), resulting in 720 MCs from the Outer Galaxy Survey and 174 MCs from the Galactic Ring Survey.

Figure 10 shows the estimated value of m_peak for the clouds in the two observational samples. In the case of the Outer Galaxy, m_peak shows a tendency to decrease with increasing cloud mass, while m_peak is essentially independent of cloud mass in the case of the Galactic Ring. On the average, the expected IMF peak is more than twice larger for clouds in the Outer Galaxy than for those in the Galactic Ring, because of the larger values of α_vir in the Outer Galaxy clouds. For the most massive clouds (few ×10⁵ M_⊙), where α_vir is relatively low also in the case of the Outer Galaxy Survey, the two samples give approximately the same value, m_peak ≈ 0.25, consistent with the peak of the Chabrier IMF. In order to estimate a characteristic value of the peak, we consider the clouds with α_vir < 3.0, because of the strong suppression of star formation at larger values of the virial parameter (e.g., Padoan & Nordlund 2011a; Padoan et al. 2012, 2017), and with mass M_cl > 10⁴ M_⊙, because most of the mass is in the most massive clouds, based on the cloud mass distribution. With these subsets of clouds from the two surveys, the mean and standard deviations are m_peak = 0.6 ± 0.25 M_⊙ and m_peak = 0.26 ± 0.09 M_⊙ for the outer and inner Galaxy, respectively, with over 90% of these star-forming clouds yielding values in the range 0.1 M_⊙ < m_peak < 1.0 M_⊙.

Zoom In Zoom Out Reset image size

This scatter in the peak of the stellar IMF predicted for different MCs is the consequence of the scatter in the velocity–size and mass-size relations, or, equivalently, the scatter in the relation between virial parameter and mass (see Figures 31, 33, 34, and 35 in Padoan et al. 2016b and Figures 5–7 in Padoan et al. 2016a). We have recently shown that supernova-driven turbulence generates MCs with properties consistent with the observations (Padoan et al. 2016a, 2016b; Pan et al. 2016). Because of this successful comparison between MCs selected from our simulation and the observations, we can use the simulation to infer that most of the scatter in the observational Larson relations may originate from true physical variations from cloud to cloud, rather than be dominated by statistical uncertainties in the observational measurements. Thus, we conclude that the predicted variations of the IMF peak from cloud to cloud, illustrated by Figure 10, are realistic. This result is consistent with the recent finding that the IMFs of young nearby stellar clusters show significant variations from region to region. Using a Bayesian analysis of the IMFs of eight young Galactic clusters, Dib (2014) has demonstrated that the posterior probability distribution functions of the IMF parameters of different clusters do not generally overlap within the 1σ uncertainty level. In the case of the Chabrier plus power-law fit, he derives IMF peak values in the range 0.29–0.69 M_⊙; in the case of the fit with Parravano's tapered power law, the range is even larger, 0.14–0.80 M_⊙.⁴ These observed IMF peak values are consistent with the ones predicted by our model applied to the star formation conditions of typical Galactic MCs. Thus, their scatter is consistent with that expected as a consequence of cloud-to-cloud variations in ${{ \mathcal M }}_{{\rm{s}}}$ and α_vir at fixed cloud mass (essentially the scatter in the Larson relations).

5.3. Variability of the IMF from Time Evolution

As explained in Padoan & Nordlund (2011b), the PN02 turbulent fragmentation model implies a time evolution of the IMF, because more massive stars are the result of converging motions from larger scales in the turbulent flow, requiring larger time to assemble the stellar mass (the turnover time of turbulent eddies increases with their size) than lower-mass stars. Therefore, at very early times, massive stars are still not fully formed, as they require a timescale comparable to the turnover time of the largest turbulent scales in the flow, of order of 1 Myr in typical MCs. This is much longer than the formation time of 100 kyr in the model of massive star formation of McKee & Tan (2002, 2003).

It should be stressed that the mechanism of massive star formation (and thus the origin of the Salpeter slope of the IMF tail) in the turbulent fragmentation models of Hennebelle & Chabrier (2008) and Hopkins (2012) is quite different than in PN02 and, unlike PN02, may lead to the McKee and Tan scenario of massive star formation. In these models, massive stars originate from massive cores that manage to exceed their Jeans mass. The reason why a large mass is needed to exceed the Jean mass is that the turbulence is included as a source of pressure support defining the Jeans mass, despite the fact that such a generalization of the Jeans mass is actually valid only in the case in which the turbulent outer scale is much smaller than the core size and the turbulent velocity is much smaller than the speed of sound (Chandrasekhar 1951), both conditions being largely violated in the context of these models. The collapse of such a massive core cannot occur until it is fully formed, meaning until it has exceeded this turbulent Jeans mass. Once that happens, the core collapses and forms a massive star essentially in a free-fall time, similarly to the scenario of the McKee and Tan model. However, massive prestellar cores as predicted by these turbulent fragmentation models may have too low gas density (too large sizes), on average, compared with observed cores, or even with the initial conditions of the McKee and Tan model, because they only need to be mild density fluctuations in the turbulent flow, rather than post-shock regions.

Turbulent pressure support against self-gravity plays no role in PN02, where the turbulence is only viewed as a source of density enhancement through shocks. Prestellar cores are assumed to emerge in the post-shock gas, where the turbulence has been largely dissipated. The inertial converging flows feeding such post-shock cores can accumulate enough mass to form a massive star, over a characteristic turnover time on the scale of such flows, much longer than the free-fall time in the post-shock gas. At the post-shock density, such mass would be many times larger than the Jeans mass (excluding support from turbulent motions that is not important in the post-shock gas), so the core cannot be supported against collapse for the whole time necessary to gather all the available mass. As soon as the critical mass for collapse in the post-shock gas has been reached, a protostar of intermediate mass is formed by the collapse of the core, and the rest of the mass has to be accreted through a circumstellar disk fed by the same converging flows that had assembled the prestellar core. In other words, the stellar mass predicted by the PN02 model should be seen as the total mass available to form a star, while the actual mass of a prestellar core (prior to its collapse into a protostar) could be significantly smaller, at least in the case of massive stars (see Figure 1 in Padoan & Nordlund 2011b). This results in a difference between the prestellar core mass function (MF) and the stellar IMF, with the prestellar core MF having a steeper high-mass tail than the Salpeter IMF (Padoan & Nordlund 2011b). In the case of low-mass stars, the stellar mass is not much larger than the characteristic BE mass in the post-shock gas, so most of the core mass is assembled before the core collapses.

Earlier turbulence simulations without self-gravity and sink particles have already demonstrated the post-shock origin of prestellar cores as assumed in PN02 (Padoan et al. 2001, 2007), at odds with the scenario of Hennebelle & Chabrier (2008) and Hopkins (2012). Using a clump-find algorithm (instead of sink particles), Padoan et al. (2007) identified dense post-shock cores containing many Jeans masses (and not supported against self-gravity by their turbulent pressure) in very large simulations of supersonic MHD turbulence without self-gravity. They also found that the core mass distribution was consistent with a power law with the Salpeter slope, proving that such cores could contain the mass reservoir responsible for the formation of massive stars. Evidently, if self-gravity had been present in the simulations, those massive cores would have collapsed much before gathering their total mass, and the rest of their mass would have been accumulated over many free-fall times, as indeed shown by more recent simulations with self-gravity and sink particles, such as in the work by Padoan et al. (2014b). In Padoan et al. (2014b), using a simulation with almost identical physical and numerical parameters as the model high in this work, we obtained nearly 1300 sink particles over a time of 3.2 Myr, with a mass function closely following a Chabrier IMF at small masses and a Salpeter IMF at masses larger than 1–2 M_⊙. We used that simulation to argue that the large-scale infall from the turbulent inertial flows feeding the protostars through an accretion disk could explain the observed luminosity distribution of protostars. We also showed that, on average, the time to gather 95% of the final stellar mass, t₉₅, increases with increasing final stellar mass, M_f, according to t₉₅ = 0.45 Myr × (M_f/1 M_⊙)^0.56, so it takes on average over 1 Myr to form a 10 M_⊙ star (see Figure 13 in Padoan et al. 2014b). However, we did not see an accelerated accretion rate as the stars gain mass, so our results are at odds with the predictions of the competitive accretion scenario (Bonnell et al. 2001; Bonnell & Bate 2006).

The dependence of the formation time on the final stellar mass is confirmed by the simulations of this work. Figure 11 shows the dependence of t₉₅ on M_f, for our reference simulation high with 256³ root grid. Only stars that have practically stopped accreting by the end of the simulation are included in the plot, to ensure that the value of t₉₅ is not artificially truncated by the finite integration time. This is enforced by selecting only the sink particles whose accretion rate averaged over the final 100 kyr of the simulation is less than 10% of their accretion rate averaged from their birth time to the time they reach 95% of their final mass. This selection retains 72.5% of the sink particles. We have verified that a much more stringent selection, where the accretion rate in the last 100 kyr has dropped to less than 0.005% of its lifetime average, retains only 28% of the sink particles but yields the exact same power-law fit given below, albeit with larger uncertainties of the slope and intercept. Even the power-law fit obtained by retaining all the sink particles yields the same parameters as in Equation (9) below, within the 1σ uncertainties, as long as the two stars more massive than 20 M_⊙ (and still actively accreting) are excluded from the fit. The apparent insensitivity of the relation between t₉₅ and M_f is due to the long integration time of the simulation relative to the t₉₅ values of even the most massive stars.

Zoom In Zoom Out Reset image size

Figure 11 shows that the values of t₉₅ increase with increasing M_f, with a lower envelope approximately consistent with a constant infall rate of approximately 0.7 × 10⁻⁵ M_⊙ yr⁻¹. A power-law fit of the average values of t₉₅ in logarithmic bins of M_f gives

$$\begin{eqnarray}&&{t}_{95}=0.51\pm 0.02\,\mathrm{Myr}\times {({M}_{{\rm{f}}}/1{M}_{\odot })}^{0.58\pm 0.02},\end{eqnarray} \tag{ 9 }$$

consistent with our previous result in Padoan et al. (2014b). Although the fit would be barely affected by including the two lowest-mass bins, those bins were excluded because the relation seems to flatten at the lowest masses. Furthermore, we do not expect t₉₅, as we measure it, to scale with M_f for masses below the IMF peak, as such stars (mainly brown dwarfs) do not result from characteristic turbulent compressions at a certain (small) scale, but from very rare compression events biased toward very large density (necessary for the BE mass to reach brown-dwarf values; Padoan & Nordlund 2004).

The relatively long formation timescale should result in a time dependence of the IMF. Using the same simulation high, we quantify this time evolution through a lognormal fit of the IMF between 0.03 and 1.5 M_⊙ and a power-law fit of the IMF above 1.5 M_⊙. Figure 12 shows the IMF and the two fits at seven select times. Every time includes stars of ages between 0 and Δt, so it is equivalent to observe a star-forming region of age Δt. At very early times, only a few stars have been formed, so the IMF would be poorly defined. To increase the statistical sample, for every value of Δt we consider all the adjacent time intervals of size Δt included between the beginning and the end of the star formation process (2.5 Myr in total) and measure the mass of all the sink particles formed in each of those intervals (hence with age between 0 and Δt). In this way, the stellar sample is always of the order of the total number of sink particles at the end of the simulation, for every value of Δt. Under the reasonable assumption that the average physical conditions in the simulated star-forming region are stationary, this procedure should provide the correct IMF also for very small values of Δt. For clarity, the IMF histograms in Figure 12 have been shifted vertically by an arbitrary value, except for the case of 0.03 Myr, which shows the actual number of stars in each mass bin.

Zoom In Zoom Out Reset image size

Figure 12 shows a clear time evolution of the IMF, with the lognormal part at low masses becoming gradually broader over the first megayear of evolution, and a gradual buildup of the power-law tail. The time dependences of the best-fit parameters are plotted in Figure 13, where the horizontal dashed lines mark the value of mean and standard deviation of the observational IMF from Chabrier (2005) (middle and top panels) and the Salpeter IMF slope (bottom panel). Both the peak and the width, σ_m, of the IMF grow with time and settle after approximately 1 Myr.

Zoom In Zoom Out Reset image size

While σ_m settles at a value indistinguishable from the observational one, m_peak remains a bit larger than Chabrier's value. The time-averaged value of the IMF peak after the first 1 Myr is $\langle {m}_{\mathrm{peak}}\rangle =0.28\pm 0.02$. The exponent of the power-law tail has a rather large uncertainty, due to the limited number of intermediate- and high-mass stars in the simulation. However, it is clearly decreasing over time as the high-mass tail of the IMF develops. The gradual buildup of the IMF tail takes nearly 2 Myr, after which its power-law slope, Γ, is essentially constant and consistent with Salpeter's value.

This time dependence of the IMF is very hard to constrain observationally, because a very young star-forming region, say, less than 0.5 Myr, would typically have very few stars and thus a poorly defined IMF (we cannot generate many Δt intervals as we have done with the simulation). Furthermore, the variations of m_peak due to variations in the physical parameters of the star-forming regions according to Equation (8) may be difficult to disentangle from the effect of time evolution, particularly if the physical parameters of a given region also vary with time.

We have presented results from a large set of simulations of randomly driven supersonic MHD turbulence, meant to represent a characteristic MC region of 4 pc and 3000 M_⊙. To match the assumptions of the turbulent fragmentation model of the IMF, we have adopted an isothermal equation of state, neglecting both radiative and mechanical feedbacks from protostars. Our main results are as follows:

• 1.  
  Thanks to a large number of test runs, we have found a definitive set of optimal parameters of our sink particle model: number of cells per Jeans' length at the sink particle creation density L_J,s = 2, exclusion radius for creation r_ex = 8 Δx, ratio of accretion threshold density to sink creation density ρ_th/ρ_s = 2, and accretion radius r_acc = 4 Δx.
• 2.  
  We have found clear trends toward a numerical convergence of the IMF and have shown that isothermal MHD turbulence can reproduce both the IMF turnover, consistent with the Chabrier IMF, and the power-law tail at high masses, consistent with the Salpeter slope.
• 3.  
  The dependence of the IMF peak on the mean gas density predicted by the turbulent fragmentation model has been confirmed using three simulations with different values of the mean density.
• 4.  
  We have estimated the expected variations of the IMF peak due to variations of virial parameter and Mach number in MCs, based on observed cloud properties derived from large Galactic MC surveys. The predicted variations are in line with the observational evidence of IMF variations in young stellar clusters.
• 5.  
  The timescale of star formation increases with stellar mass, with 10 M_⊙ stars requiring on average over 1 Myr to reach 95% of their final mass. This results in a time dependence of the IMF. It takes of order 1 Myr for the width and peak of the IMF to grow to their final values, and nearly 2 Myr for the power-law high-mass tail to get established and for the slope to settle to the Salpeter value.
• 6.  
  The time evolution of the IMF is expected from the turbulent fragmentation model of PN02, but not necessarily from other turbulent fragmentation models. It is also in contradiction to the idea that massive stars form from the rapid collapse of massive cores, as well as to the idea that the growth of massive stars is controlled primarily by their own mass through competitive accretion.

We conclude that isothermal MHD turbulence yields a stellar IMF consistent with the observations. Other physical processes, such as radiation and outflows from protostars, competitive accretion, or dynamical interactions, may not play a dominant role. Because we have focused on typical star-forming conditions in Galactic MCs, it is still possible that these other processes become dominant in more extreme environments, such as near the Galactic center or in starburst galaxies. The important implications for the formation of massive stars mentioned in this work will be addressed in detail in a separate work.

T.H. is supported by a Sapere Aude Starting Grant from the Danish Council for Independent Research. P.P. is supported by the Spanish MINECO under project AYA2014-57134-P. Å.N. is supported by a large framework grant from the Danish Council for Independent Research. Research at the Centre for Star and Planet Formation is funded by the Danish National Research Foundation (DNRF97). Computing resources for this work were partly provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. We acknowledge PRACE for awarding us access to the computing resource JUQUEEN based in Germany at Forschungszentrum Jülich and access to the computing resource CURIE based in France at CEA, where part of the simulations where carried out. The astrophysics HPC facility at the University of Copenhagen, which is supported by a research grant (VKR023406) from VILLUM FONDEN, was used for carrying out part of the simulations, the analysis, and long-term storage of the results.

Performing simulations of magnetized and highly supersonic flows is a challenge for Godunov solvers, which is only amplified when coupled with self-gravity, sink particles, and fine-coarse transitions of an adaptive mesh. The average Mach number in the presented simulations is 10, but at exceptional points in the flow the combined advection and fast-mode speed can be more than 500 times the sound speed. To stabilize the solver in these extreme conditions, we have made several changes to the public MHD solver included in ramses. First of all, we have contributed a small bug fix to how the magnetic fluxes are updated at fine–coarse interfaces that occasionally would result in nonphysical solutions. Second, we have made a modification of how states in the HLLD solver are computed when in principle a "0/0" situation arises (Miyoshi & Kusano 2005) in the computation of the transverse velocities and magnetic fields in the Riemann fan. The equations for the intermediate states are

$$\begin{eqnarray}&&{e}^{* }=\rho (S-v)(S-{v}^{* })-{B}_{\parallel }^{2}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{v}_{\perp }^{* }={v}_{\perp }-{B}_{\parallel }{B}_{\perp }({v}^{* }-v)/{e}^{* }\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{B}_{\perp }^{* }={B}_{\perp }e/{e}^{* },\end{eqnarray} \tag{ 12 }$$

where ρ, v, B, e, and S are the density, velocity, magnetic field, internal energy, and combined advection and fast-mode speed, respectively. Here ∥ and ⊥ indicate the parallel and perpendicular components, respectively, with respect to the normal of the interface over which the Riemann problem is solved. This is indeterminate when e* = 0. Instead, the limiting expressions (Miyoshi & Kusano 2005)

$$\begin{eqnarray}&&{e}^{* }=\rho (S-v)(S-{v}^{* })-{B}_{\parallel }^{2}=0\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{v}_{\perp }^{* }={v}_{\perp }\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{B}_{\perp }^{* }={B}_{\perp }\end{eqnarray} \tag{ 15 }$$

should be used. In practice, numerical round-off introduce large errors when e* is close to 0, compared to the quantities involved. We have therefore modified the condition, such that we use the limiting expressions already when ${e}^{* }\lt {10}^{-4}{B}_{\parallel }^{2}$. This has removed cases where the solver crashed immediately from one time step to the next. Both fixes have been propagated to the public version of ramses.

In a finite-volume code, unphysical ringing and striping in the solution can arise near strong shocks. The principle mechanism to counteract the problem is using a slope-limiting procedure that lowers the order of interpolation from centered to interface values in the vicinity of the shocks introducing diffusion. If the slope-limiting procedure is only based on variations in variables across the interface, in a dimensionally split MHD code plane-parallel shock fronts can trigger the so-called carbuncle instability (Sutherland et al. 2003). To counteract this, we have introduced a new three-dimensional slope-limiting procedure, which has the added benefit that it introduces diffusion isotropically close to shocks irrespective of the direction of propagation maintaining, e.g., spherical symmetry in supernova-driven shock fronts. Letting f(i, j, k) be the variable that is to be slope limited, we then compute all 26 slopes

$$\begin{eqnarray}&&{\rm{\Delta }}{f}_{\alpha \beta \gamma }(i,j,k)=f(i+\alpha ,j+\beta ,k+\gamma )-f(i,j,k),\end{eqnarray} \tag{ 16 }$$

where α, β, γ = −1, 0, +1, and also the face-aligned slopes

$$\begin{eqnarray}&&{{df}}_{x}(i,j,k)=\displaystyle \frac{1}{2}(f(i+1,j,k)-f(i-1,j,k))\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{{df}}_{y}(i,j,k)=\displaystyle \frac{1}{2}(f(i,j+1,k)-f(i,j-1,k))\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{{df}}_{z}(i,j,k)=\displaystyle \frac{1}{2}(f(i,j,k+1)-f(i,j,k-1)).\end{eqnarray} \tag{ 19 }$$

This is used to construct an isotropic slope limiter

$$\begin{eqnarray}&&\mathrm{lim}\,=\,\min \left(1,\delta \displaystyle \frac{\min (| \min ({\rm{\Delta }}f)| ,| \max ({\rm{\Delta }}f)| )}{\sqrt{{{df}}_{x}^{2}+{{df}}_{y}^{2}+{{df}}_{z}^{2}}}\right),\end{eqnarray} \tag{ 20 }$$

where 1 ≤ δ ≤ 2 is an input parameter that smoothly interpolates the slope limiter between behaving as a min-mod and a monotonized-central limiter. The limiter is then applied to the final slopes equally in all directions: ${{df}}_{x,y,z}(i,j,k)\,\to \mathrm{lim}\,\times \,{{df}}_{x,y,z}(i,j,k)$. We have typically used a value of δ = 1.5 for the simulations presented in this paper.

With the above modifications the HLLD solver in ramses has turned out to be exceptionally stable. Even though the solution is stable, sometimes the time step will decrease dramatically owing to cells with strong magnetic fields and low densities typically in the vicinity of sink particles. To finish the runs in a reasonable time, we calculate the combined advection and fast-mode speed, c_B, at all cell centers. If ${c}_{B}\gt {c}_{B}^{\max },$ we use a diffusive Lax-Friedrichs (LLF) solver on all connecting edges and interfaces of the cell to calculate the Riemann problem. This is effective in diffusing the solution sufficiently for the time step to stabilize at reasonable values. Typically we have used ${c}_{B}^{\max }=500$, which activates LLF in approximately 1 out of 10⁴ cells.

In ramses the Poisson equation is solved on the adaptive mesh using a combination of multigrid and conjugate gradient methods (Guillet & Teyssier 2011). The code can use not only spatial but also temporal subcycling, with a factor of two resolution difference between each AMR level. The latter is required when using 10 or more AMR levels for efficient evolution of the model. The code is structured such that the gravitational potential is calculated on the coarse levels first, with a recursive call midway through the time step to the time evolution routine on the finer level. The time evolution of the cells is done after this recursive call, resulting in the finest AMR levels to be evolved forward first in a W-like evolution cycle for the different AMR levels. Therefore, the boundary conditions on individual patches in the AMR hierarchy, computed by interpolating the potential at coarser levels to boundary cells constructed on the fly, can be out of temporal sync at the second step in a subcycle. If two sink particles are placed inside individual patches with coarser grid cells in between, this will lead to a delayed transmission of the gravitational force and a complete decay of, e.g., binary systems in less than 10 orbits. In a periodic box, tests with a few particles show that the delayed transmission of gravity through lower levels subjects moving particles to self-interaction and will slow them down on relatively short timescales. Notice that self-interaction is in practice not a problem in our simulation since a single particle will not dominate the potential. While these effects are most obvious for the particles (both sink and dark matter), similar effects happen for the gas dynamics. Fortunately, the gravitational potential is the second integral of the density, and as such in general smooth, both spatially and in time. In the case of time subcycling we have implemented extrapolation forward in time in the boundary cells between different levels. Letting Δt_c be the past time step at the coarse level and Δt_f be the current time step at the fine level, the gravitational potential at the boundary then becomes ${\phi }_{t+{\rm{\Delta }}t/2}=(1+{\rm{\Delta }}{t}_{f}/{\rm{\Delta }}{t}_{c}){\phi }_{t}-({\rm{\Delta }}{t}_{f}/{\rm{\Delta }}{t}_{c}){\phi }_{t-{\rm{\Delta }}t}$. This extrapolation has effectively removed the decaying orbits and self-interaction induced by time subcycling and is now also integrated in the public version of the ramses code.

Federrath et al. (2010) use a combination of a grid-based potential for the gas and a direct computation of the gravitational forces of the sink particles, which has a cost that is proportional to (N_cell + N_sink) × N_sink. They also subcycle the sink particles in time to resolve the orbits of individual sink orbits. While this method is arguably very precise in evolving in particular multiple stars in close orbits, it rapidly becomes prohibitively expensive when considering large-scale simulations, such as those presented in this paper, with up to 1000 sink particles. Gong & Ostriker (2013) present an alternative implementation where the combined gravitational potential of the gas and sink particles is computed on the grid, and the force on individual particles is found through interpolation of the forces on the nearest grid cell centers. While this method is scalable, it is not precise at the subcell scale, and long-term stable orbits require binaries to be separated by at least 10 grid cells at all times, and even then the orbits are prone to precession and some level of decay.

To couple gas and sink particles gravitationally, we have developed two methods for evaluating the force on individual particles. We use either CIC (eight neighbors) or TSC (27 neighbors) interpolation from the sink particle position to deposit the mass on the grid. The combined density of the gas and the sink particles is used as input to the Poisson equation, and the resulting gravitational potential ${\phi }_{\mathrm{sink}+\mathrm{gas}}$ is used for finding the gravitational acceleration of the gas.

In the first method to calculate the gravitational acceleration on the sink particles we use a matching CIC or TSC interpolation of the gradients of ϕ_sink+gas to the sink particle position. This corresponds closely to the method described in Gong & Ostriker (2013). In the second method, which is the default for the runs presented here, only the gas is updated by calculating the gravitational acceleration from ${\phi }_{\mathrm{sink}+\mathrm{gas}}$. To find the gravitational acceleration of the sink particles, we solve the Poisson equation again, but this time for the gas only. The gravitational acceleration of the gas on the sink particles can then be computed from ϕ_gas, while the gravitational force of the sink particles instead is evaluated by direct summation of the forces from individual particles. We use a cubic spline softening of the gravitational force (Federrath et al. 2010) with a typical softening length of 0.3Δx when directly computing the forces between different sink particles.

Both of the methods are scalable to thousands of sink particles, and the second method gives as precise orbits in multiple-star systems as the method of Federrath et al. (2010) at the added cost of requiring an additional solution to the Poisson equation.

The smaller-scale 64³ root grid test models in Appendix C, RUN1–RUN27, are using the first method, while all other runs have been evolved with the second method. We have performed test runs at ${{ \mathcal M }}_{{\rm{s}}}=10$ with both the first and the second method. Even though individual binary systems are much better described by the second method, the overall impact on the IMF in the test runs is negligible and significantly smaller than the Poisson scatter in the IMF distribution.

To move the sink particles through the volume, we use a simple leapfrog Kick–Drift–Kick algorithm, identical to the method used to integrate the orbits of dark matter particles in ramses. The global Courant condition is modified to take into account both the velocity v_x,y,z and acceleration a_x,y,z of the sink particles. Let Δx^s be the minimum of the distances between sink particles and the softening length

$$\begin{eqnarray}&&{\rm{\Delta }}{x}^{s}=\min [{ \mathcal S }{\rm{\Delta }}x,\mathrm{sink}\mbox{--}\mathrm{sink}\,\mathrm{distances}],\end{eqnarray} \tag{ 21 }$$

where ${\rm{\Delta }}x={2}^{-{\ell }}$ is the cell size at refinement level ℓ and ${ \mathcal S }=0.3$ is the softening number. Then the limit on the time step imposed by each sink particle i is

$$\begin{eqnarray}&&{{dt}}_{{\rm{s}}}=\min \left[\displaystyle \frac{{C}_{{\rm{\Delta }}t}{\rm{\Delta }}{x}^{s}}{2\max [| {v}_{x}^{i}| ,| {v}_{y}^{i}| ,| {v}_{z}^{i}| ]},{\left(\displaystyle \frac{{C}_{{\rm{\Delta }}t}{\rm{\Delta }}{x}^{s}}{| {a}^{i}| }\right)}^{1/2}\right],\end{eqnarray} \tag{ 22 }$$

where ${C}_{{\rm{\Delta }}t}$ is the Courant factor, which is set to 0.5 in the runs.

The sink particle algorithm contains a number of density and distance thresholds for creation and accretion of gas to the sink particles, and it is important to understand the impact of different choices, and in particular the impact on the SFR, and the IMF. To test the algorithm, we have created a large grid of models that enables us to study the dependence of the IMF on all parameters. The test runs are prepared using the same physical parameters and driving as the convergence runs and are listed in Table 2. Our grid of models mainly affects three aspects of the simulation.

Table 2.  Numerical Parameters of the Test Runs

	Run Parameters								Creation of Sinks				Accretion to Sinks
Run	Root	NAMR	Δx	LJ	Mass	αvir	tend	SFE	LJ,s	ρs	ρs	rex	ρacc	ρacc	Nbody	racc	acc
	Grid		(au)		(M⊙)		(Myr)			(cm−3)	($\langle \rho \rangle$)	(Δx)	(cm−3)	($\langle \rho \rangle$)		(Δx)
RUN1	643	6	200	8.0	3000	0.83	1.5	25%	8	6.6 × 106	8.3 × 103	16	3.3 × 106	4.1 × 103	no	2	1.0
RUN2	643	6	200	8.0	3000	0.83	1.6	25%	8	6.6 × 106	8.3 × 103	32	3.3 × 106	4.1 × 103	no	2	1.0
RUN3	643	6	200	8.0	3000	0.83	1.6	25%	8	6.6 × 106	8.3 × 103	64	3.3 × 106	4.1 × 103	no	2	1.0
RUN4	643	6	200	8.0	3000	0.83	1.5	25%	4	2.6 × 107	3.3 × 104	32	1.3 × 107	1.6 × 104	no	2	1.0
RUN5	643	6	200	8.0	3000	0.83	1.5	25%	4	2.6 × 107	3.3 × 104	32	1.3 × 107	1.6 × 104	no	2	1.0
RUN6	643	6	200	8.0	3000	0.83	1.6	25%	4	2.6 × 107	3.3 × 104	64	1.3 × 107	1.6 × 104	no	2	1.0
RUN7	643	6	200	8.0	3000	0.83	1.5	25%	2	1.0 × 108	1.3 × 105	16	5.2 × 107	6.6 × 104	no	2	1.0
RUN8	643	6	200	8.0	3000	0.83	1.5	25%	2	1.0 × 108	1.3 × 105	32	5.2 × 107	6.6 × 104	no	2	1.0
RUN9	643	6	200	8.0	3000	0.83	1.5	25%	2	1.0 × 108	1.3 × 105	64	5.2 × 107	6.6 × 104	no	2	1.0
RUN10	643	7	100	8.0	3000	0.83	1.5	25%	8	2.6 × 107	3.3 × 104	16	1.3 × 107	1.6 × 104	no	2	1.0
RUN11	643	7	100	8.0	3000	0.83	1.6	25%	8	2.6 × 107	3.3 × 104	32	1.3 × 107	1.6 × 104	no	2	1.0
RUN12	643	7	100	8.0	3000	0.83	1.5	25%	8	2.6 × 107	3.3 × 104	64	1.3 × 107	1.6 × 104	no	2	1.0
RUN13	643	7	100	8.0	3000	0.83	1.6	25%	4	1.0 × 108	1.3 × 105	16	5.2 × 107	6.6 × 104	no	2	1.0
RUN14	643	7	100	8.0	3000	0.83	1.6	25%	4	1.0 × 108	1.3 × 105	32	5.2 × 107	6.6 × 104	no	2	1.0
RUN15	643	7	100	8.0	3000	0.83	1.5	25%	4	1.0 × 108	1.3 × 105	64	5.2 × 107	6.6 × 104	no	2	1.0
RUN16	643	7	100	8.0	3000	0.83	1.5	25%	2	4.2 × 108	5.2 × 105	16	2.1 × 108	2.6 × 105	no	2	1.0
RUN17	643	7	100	8.0	3000	0.83	1.5	25%	2	4.2 × 108	5.2 × 105	32	2.1 × 108	2.6 × 105	no	2	1.0
RUN18	643	7	100	8.0	3000	0.83	1.5	25%	2	4.2 × 108	5.2 × 105	64	2.1 × 108	2.6 × 105	no	2	1.0
RUN19	643	8	50	8.0	3000	0.83	1.6	25%	8	1.0 × 108	1.3 × 105	16	5.2 × 107	6.6 × 104	no	2	1.0
RUN20	643	8	50	8.0	3000	0.83	1.6	25%	8	1.0 × 108	1.3 × 105	32	5.2 × 107	6.6 × 104	no	2	1.0
RUN21	643	8	50	8.0	3000	0.83	1.6	25%	8	1.0 × 108	1.3 × 105	64	5.2 × 107	6.6 × 104	no	2	1.0
RUN22	643	8	50	8.0	3000	0.83	1.6	25%	4	4.2 × 108	5.2 × 105	16	2.1 × 108	2.6 × 105	no	2	1.0
RUN23	643	8	50	8.0	3000	0.83	1.6	25%	4	4.2 × 108	5.2 × 105	32	2.1 × 108	2.6 × 105	no	2	1.0
RUN24	643	8	50	8.0	3000	0.83	1.6	25%	4	4.2 × 108	5.2 × 105	64	2.1 × 108	2.6 × 105	no	2	1.0
RUN25	643	8	50	8.0	3000	0.83	1.6	25%	2	1.7 × 109	2.1 × 106	16	8.3 × 108	1.0 × 106	no	2	1.0
RUN26	643	8	50	8.0	3000	0.83	1.5	25%	2	1.7 × 109	2.1 × 106	32	8.3 × 108	1.0 × 106	no	2	1.0
RUN27	643	8	50	8.0	3000	0.83	1.6	25%	2	1.7 × 109	2.1 × 106	64	8.3 × 108	1.0 × 106	no	2	1.0
RUN28	5123	3	200	8.0	3000	0.83	2.0	25%	1	4.2 × 108	5.2 × 105	8	2.1 × 108	2.6 × 105	yes	2	1.0
RUN29	5123	3	200	8.0	3000	0.83	2.0	25%	2	1.0 × 108	1.3 × 105	4	5.2 × 107	6.6 × 104	yes	2	1.0
RUN30	5123	3	200	8.0	3000	0.83	1.9	25%	2	1.0 × 108	1.3 × 105	8	5.2 × 107	6.6 × 104	yes	2	1.0
RUN31	5123	3	200	8.0	3000	0.83	2.0	25%	2	1.0 × 108	1.3 × 105	16	5.2 × 107	6.6 × 104	yes	2	1.0
RUN32	5123	3	200	8.0	3000	0.83	2.0	25%	4	2.6 × 107	3.3 × 104	8	1.3 × 107	1.6 × 104	yes	2	1.0
RUN33	643	6	200	8.0	3000	0.83	1.5	25%	2	1.0 × 108	1.3 × 105	8	5.2 × 107	6.6 × 104	yes	2	1.0
RUN34	10243	2	200	8.0	3000	0.83	2.1	25%	2	1.0 × 108	1.3 × 105	4	5.2 × 107	6.6 × 104	yes	2	1.0
test	1283	6	100	7.2	3000	0.83	2.5	13%	2	4.2 × 108	5.3 × 105	8	2. 1 × 108	2.6 × 105	yes	4	1.0
acc	1283	6	100	7.2	3000	0.83	2.5	13%	2	4.2 × 108	5.3 × 105	8	4227	5.3	yes	4	1.0

Note. the Test and acc runs are identical to the med run, except that the acc run is using an accretion efficiency _acc = 1 and the test run, in addition to that, is only using a density threshold for accretion like the rest of the test runs. We have set in boldface the main parameters that are explored in each set of runs.

Download table as:  ASCIITypeset image

C.1. Creation

C.2. Accretion

C.3. Numerical Resolution

C.4. Star Formation Efficiency and Average Stellar Mass

The models have run with self-gravity and sink particles for up to 1.5 free-fall times or 1.8 Myr, similar to the lifetimes of nearby star-forming regions. At the end of the runs, the SFE reaches 25%. Notice that since these models have an accretion efficiency _acc = 1, their SFEs have to be rescaled with a factor of 2.4, corresponding to a rescaled SFE of 10%, to compare them with the convergence runs as shown in Figure 14 and discussed in Section 2. Figure 15 shows that the overall SFE is extremely robust between runs, which is a consequence of gravitational collapse: if gas reaches a certain threshold density, it will almost unavoidably collapse. Depending on the details of the sink recipe, this will result in either accretion to an existing sink particle or creation of a new sink particle, but the total accreted mass remains nearly constant.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The average mass of the sink particles produced in each run is shown in Figure 16. We have grouped the runs in terms of the Jeans number at the threshold density for sink particle creation, L_J,s. There are both common trends and a significant dependence on the numerical details of the average mass: there is a clear dependence on the exclusion radius, the number of AMR levels, and the density threshold at which sink particles are created. Given that the SFE is very similar for all runs, a smaller average mass corresponds to a higher number of sink particles and vice versa. A smaller exclusion radius, a higher number of AMR levels, and a lower threshold for sink particle creation all result in a lower average mass and therefore a higher number of sink particles. The difference between different models diminishes with increasing threshold for sink creation, showing convergence at the highest resolution. The parameter with the smallest effect on the average mass is the exclusion radius. Whether adding a higher number of AMR levels, or changing the sink threshold, both have a larger impact.

Zoom In Zoom Out Reset image size

As we will see below, the reasons for the significant dependence of the average stellar mass on the sink particle creation parameters and N_AMR may be numerical: with a small root grid of 64³ cells, the threshold for sink creation is too low, resulting in the spurious creation of extra particles. Thus, besides being a clear illustration of the effect of numerical parameters, Figure 16 is also a clear indication that the resolution of these test runs is too low to study the stellar IMF, though probably sufficient to study the SFR (as indicated by Figure 15). Although we cannot afford a similar systematic set of test simulations with a much larger root grid size, we expect that the dependence on numerical parameters shown in Figure 16 would be significantly weaker in simulations with a root grid size of 256³ cells or larger. The higher-resolution runs can resolve an increasing number of close binaries, which should not be confused with the artificial creation of spurious sink particles. Finally, another important effect that is explored below is how well the turbulent fragmentation and the smallest collapsing cores are resolved.

A common feature of all models is that the average sink particle mass becomes nearly stationary after roughly 0.8 free-fall times. This can be understood in terms of the time it takes for the most massive stars in a cluster to accrete. Only after ∼0.8 free-fall times is enough of the high-mass Salpeter slope developed to keep the average mass essentially constant, as demonstrated by Figures 12 and 13 in Section 5.3, which are based on our highest-resolution simulations.

C.5. Nearest-neighbor Separation at Formation

To investigate the effects of different parameter choices, a higher-order statistics than the SFE or the average stellar mass is needed. In this section we will introduce the nearest-neighbor separation at formation, and we use it below, together with the IMFs of the different runs, to better understand the differences between the test runs. The nearest-neighbor separation at formation, ${ \mathcal D }$, or "neighbor statistics" for brevity, is the histogram of distances from a set of sink particles, at the time of their formation, to the nearest other particle in the simulation volume. It is both very simple to extract from a numerical simulation and a powerful statistic to judge the fidelity of a sink particle distribution, if we have prior knowledge about how it should behave.

In this subsection we will argue that ${ \mathcal D }$ should increase as a function of distance for reasonable assumptions about the stellar distributions. This can then be used to detect when models are polluted by sink particles generated as a consequence of nonphysical fluctuations close to other sink particles.

Assume that there is no difference between the statistical distribution of newly formed stars and more evolved stars. Then ${ \mathcal D }$ can be defined in terms of the nearest-neighbor function of a point-particle set. For a point-particle set with N particles the probability of finding another point at a distance r of a given point is

$$\begin{eqnarray}&&{D}^{N}(r)=-\displaystyle \frac{{dP}(V(r))}{{dr}},\end{eqnarray} \tag{ 23 }$$

where P(V(r)) is the conditional probability (averaged over all points) that there is another point inside the volume V centered on a given point, and the derivative gives the probability that the nearest neighbor is in the interval [r, r + dr]. Since in a star-forming region stars are formed one at a time, ${ \mathcal D }$ can then be found by forming the cumulative sum

$$\begin{eqnarray}&&{{ \mathcal D }}^{N}(r)=\displaystyle \frac{1}{N}\displaystyle \sum _{n=1}^{N}{D}^{n}(r).\end{eqnarray} \tag{ 24 }$$

This is a highly nontrivial distribution related to the n-point correlation functions. As an example, if the neighbor diagram has a strong peak close to the exclusion radius for creation of new particles, r_ex, it could indicate that

• 1.  
  the overall spatial resolution is so low that clumps that should fragment into many stars are not well resolved;
• 2.  
  the minimum density for sink formation, ρ_s, is low enough that disturbances in the fluid dynamics around other sink particles are sufficient to induce the formation of new sink particles; and
• 3.  
  the Jeans length is not sufficiently resolved, and therefore numerical fragmentation results in a large number of artificial "multiple systems."

To elucidate the properties of this function, we may consider the simplest theoretical example of a stationary Poisson process with a volume point density of n, where we have

$$\begin{eqnarray}&&{D}^{n}(r)=-\displaystyle \frac{d(1-{e}^{-{nV}(r)})}{{dr}}=4\pi {r}^{2}{{ne}}^{-{nV}(r)}.\end{eqnarray} \tag{ 25 }$$

If we assume a starting density of one star per unit volume appropriate for our periodic box, we can then write the neighbor statistics for N stars as

$$\begin{eqnarray}&&{{ \mathcal D }}^{N}(r)=\displaystyle \frac{4\pi {r}^{2}}{N}\displaystyle \sum _{n=1}^{N}n\exp [-{nV}(r)].\end{eqnarray} \tag{ 26 }$$

To test for numerical artifacts, the solutions at small distances are most relevant, where NV(r) < 1. There we have ${ \mathcal D }(\mathrm{ln}r)\propto {r}^{3}$ with a binning in terms of $\mathrm{ln}r$, as done in the histograms from the numerical simulations. At larger scales the function has a relatively well-defined peak close to the point where all terms in the sum in Equation (26) still contribute to the total histogram, i.e., where NV(r) ∼ 1, corresponding to the average separation between particles of the final point set.

Define the two-point correlation function, ξ(r), as the excess probability, compared to a homogeneous point set, of finding another particle at distance r from a given point

$$\begin{eqnarray}&&{dP}=n(1+\xi (r)){dV}.\end{eqnarray} \tag{ 27 }$$

Assuming for a general distribution that all higher-order correlation functions are negligible, we can then repeat the analysis above to find that

$$\begin{eqnarray}&&{{ \mathcal D }}^{N}(r)=\displaystyle \frac{4\pi {r}^{2}[1+\xi (r)]}{N}\displaystyle \sum _{n=1}^{N}n\exp [-{{nV}}_{\xi }(r)],\end{eqnarray} \tag{ 28 }$$

where ${V}_{\xi }(r)\equiv {\int }_{V}[1+\xi (r)]{dV}$. For large N,

$$\begin{eqnarray}{{ \mathcal D }}^{N}(r) & \approx & \displaystyle \frac{4\pi {r}^{2}[1+\xi (r)]}{{V}_{\xi }{(r)}^{2}}\{[1+{V}_{\xi }(r)]\exp [-{V}_{\xi }(r)]\\ & & -[1+{{NV}}_{\xi }(r)]\exp [-{{NV}}_{\xi }(r)]\}.\end{eqnarray} \tag{ 29 }$$

For a clustered distribution such as those observed in young stellar clusters, and assuming a power law ξ(r) ∝ r^β, we find that at small separations, where ξ(r) ≫ 1,

$$\begin{eqnarray}&&{ \mathcal D }(r)\propto {r}^{3+\beta }d\mathrm{ln}r.\end{eqnarray} \tag{ 30 }$$

Observationally, the two-point correlation functions in the nearby Taurus, ρ Oph, and Orion regions are found to be a power law with exponents in the range −2 < β < −1 (Gomez et al. 1993; Larson 1995; Gomez & Lada 1998). This corresponds to a slowly increasing trend for ${ \mathcal D }$ as a function of distance, but with some uncertainty. In the numerical models we are applying the statistics at the birth of the star, while observationally we see a time-evolved distribution. Spectroscopic multiple systems formed through, e.g., disk fragmentation may eject stellar companions to the field, and tight multiple stellar systems may have formed dynamically from systems with larger separations. In a clustered system, ${ \mathcal D }$ may also have contributions from higher-order correlation functions. Thus, the observational evidence seems to favor a constant or weakly increasing nearest-neighbor separation at formation as a function of the logarithmic distance binning, though there is a considerable uncertainty for the innermost scales, below ∼1000 au. A decreasing ${ \mathcal D }$ would require a highly peaked two-point correlation function with a slope less than ≈−3, which seems to be at odds with the observations.

Because ${ \mathcal D }$ should increase as a function of distance (up to the mean particle separation), we do not expect it to peak at the smallest possible distance in the simulation, that is, at the exclusion radius, r_ex (see below). Thus, a sharp peak in ${ \mathcal D }$ at r_ex is to be interpreted as clear evidence of spurious sink particles emerging from numerical artifacts (e.g., unphysical density fluctuations near the accretion radius).

C.6. Impact of Sink Formation Parameters

In this subsection we discuss how the different numerical parameters in the sink formation model impact the IMF and the neighbor statistics, ${ \mathcal D }$, discussed above.

Exclusion radius (r_ex). The purpose of the exclusion volume is to avoid the formation of spurious sink particles near the accretion radius of a previously created sink particle. However, secondary density peaks may have a physical origin and lead to the formation of bona fide stellar companions, so the exclusion radius should be kept as small as possible in order to resolve as many physical companions as possible. We use the neighbor statistics, ${ \mathcal D }$, as our tool of choice to detect spurious sink particles. As r_ex is a purely numerical radius, the physics of star formation should not yield any special feature in ${ \mathcal D }$ at a distance equal to r_ex. A sharp peak at r_ex, for example, is certainly not expected, because ${ \mathcal D }$ is predicted to slowly decrease with decreasing distance, as concluded in the previous subsection. On the other hand, the numerical models for sink creation and accretion may easily cause a peak in ${ \mathcal D }$ at r_ex (for example, if the gas reservoir within r_acc is depleted too rapidly, the density may temporarily peak at a distance around r_ex, if r_ex ∼ r_acc). Thus, we assume that a peak in the neighbor statistics at a distance equal to r_ex is evidence for spurious sink particles and hence too small r_ex; conversely, we assume that if no obvious peak is present at a distance equal to r_ex, the value of r_ex is large enough and the contamination from spurious sink particles is small.

Figure 17 shows ${ \mathcal D }$ for several test runs with 64³ root grid and six to eight levels of refinement (top to bottom rows of panels). Each panel shows the histogram for three different values of r_ex, with the smallest value, r_ex = 16Δx, indicated by the vertical dashed line. The panels show the expected behavior: small values of r_ex may result in a peak of ${ \mathcal D }$ at a distance equal to r_ex, while the peak decreases and even disappears as the value of r_ex is increased. The effect of r_ex clearly depends on the parameter L_J,s and hence on the density threshold where sink particles are created, ρ_s, which is discussed below. One can see that with L_J,s = 4, shown in the middle column of Figure 17, only the largest value of r_ex is marginally successful at removing the unphysical peak (one remains at the intermediate resolution with seven refinement levels, shown in the central panel), while with L_J,s = 2 the peak is completely removed, also for the smallest value of r_ex.

Zoom In Zoom Out Reset image size

The panels in Figure 18 show the corresponding IMFs, where it is clear that the spurious sink particles corresponding to the peaks in the neighbor statistics are almost exclusively low-mass ones. At masses larger than their peak mass, the IMFs are essentially insensitive to the choice of r_ex. Furthermore, as for the neighbor histogram, the IMF is practically independent of r_ex when the lowest value of L_J,s is adopted, as shown in the right column of Figure 18.

Zoom In Zoom Out Reset image size

As discussed below, L_J,s = 2 is the optimal value of the numerical Jeans length at the density of sink creation, so we can conclude that, in these test runs, values of r_ex as small as r_ex = 16Δx do not result in significant numbers of spurious sink particles. In principle, we could adopt an exclusion radius as small as r_ex = 2 r_acc (a smaller value, comparable to r_acc, is guaranteed to generate spurious sinks). To test whether this is acceptable, and also to test whether the evidence collected from the low-resolution test runs holds true at higher resolution, we also consider three runs with a 512³ root grid and three levels of refinement and three different values of r_ex (RUN29, RUN30, and RUN31). The dependence of the IMF and ${ \mathcal D }$ of these runs on r_ex is shown in Figure 19. The neighbor statistics does not show any peak at r_ex, not even at the smallest value, r_ex = 2r_acc = 4Δx, and the IMFs are barely affected by the choice of r_ex. Thus, we conclude that, as long as L_J,s = 2 (see below), we can adopt the smallest possible value of r_ex, that is, r_ex = 2r_acc. In our convergence-test runs, where r_acc = 4Δx, we have adopted r_ex = 8Δx.

Zoom In Zoom Out Reset image size

Density at sink formation (L_J,s). The density at sink formation, ρ_s, has to be as high as possible, so that only bona fide collapsing regions are chosen, but not so high that the Jeans length is not resolved at the available spatial resolution. Because the turbulence alone can generate very large densities that are not necessarily collapsing regions, ρ_s must be at least several times larger than the highest post-shock density in the absence of self-gravity at the spatial resolution of the simulation. In our case, a turbulent flow with a sonic Mach number of 10, a characteristic high density without self-gravity, is roughly ∼10⁵ times the average density (depending on spatial resolution). However, a very large value of ρ_s is not a sufficient condition to avoid spurious sink particles. In the converging flow in the vicinity of an already-formed sink particle, the gas may reach densities higher than ρ_s owing to the gravitational potential around the sink, possibly triggering the formation of spurious sinks. To address this concern, we parameterize the dependence of the IMF on the density at sink formation through the value of the numerical Jeans length at the density of sink formation, L_J,s. As ρ_s is increased (with constant Δx), the value of L_J,s decreases.

The ideal value of L_J,s is a compromise between the requirement of spatially resolving the Jeans length at the end of the AMR hierarchy (high L_J,s) and the need to adopt a large value of ρ_s to avoid the creation of spurious sink particles in secondary density peaks that should not collapse into separate stars (low L_J,s). Too low values of L_J,s would cause overmerging of density fluctuations, meaning that a single star would emerge out of a fragmented region that should have yielded a multiple system (in our experience, this outcome is more likely than the risk of artificial fragmentation due to a violation of Truelove's condition). Too high values of L_J,s may result in spurious sink particles.

The dependence of the neighbor statistics and the IMF on L_J,s is shown in Figures 17 and 18. Each column of panels in those figures has a constant value of L_J,s, with L_J,s = 8 on the left, L_J,s = 4 in the middle, and L_J,s = 2 on the right. Figure 17 shows that an artificial peak of spurious sinks at a distance equal to r_ex is always present for the largest values of L_J,s, at any value of r_ex. As L_J,s is decreased, the peak decreases, particularly for larger values of r_ex. It is only with L_J,s = 2 that the peak is completely gone, independent of r_ex. As mentioned above in relation to the dependence on r_ex, the peak in the neighbor diagram is reflected by an increase in the number of low-mass stars, as shown by the IMFs of Figure 18.

To make sure that these results hold true at higher resolution, we consider again the simulations with a 512³ root grid, specifically RUN28, RUN30, and RUN32, where L_J,s = 1, 2, and 4, respectively. Figure 20 shows that the neighbor histogram yields a peak at r_ex when L_J,s = 4. The largest value where the peak is suppressed is L_J,s = 2, as in the lower-resolution runs. Thus, in order to resolve the Jeans length while also avoiding spurious sink particles, the best choice is L_J,s = 2, which is the value we have adopted in the main simulations of this work. A similar conclusion was reached by Gong & Ostriker (2013).

Zoom In Zoom Out Reset image size

C.7. Impact of Numerical Resolution

In this subsection we discuss how the numerical resolution impacts the IMF and the nearest-neighbor statistics.

Root grid size. The impact of changing the root grid size while keeping all other parameters fixed can be seen in Figure 21. The root grid resolution impacts how well the turbulence is resolved throughout the volume, but if the number of cells on the refined levels is determined by the same Jeans number, the resulting neighbor statistics and IMF are almost identical. We only see a minor difference with a slight increase in the number of low-mass stars when using a higher-resolution root grid. This makes sense, since a better resolution of the turbulence facilitates the creation of rare peaks with both high density and low mass.

Zoom In Zoom Out Reset image size

Number of refinement levels. The number of refinement levels controls the minimum cell size, and correspondingly the maximum density that can be reached, if everything else is kept fixed. The most realistic cases for our test runs are shown in the right column of Figure 18, corresponding to RUN7–RUN9 (N_AMR = 6), RUN16–RUN18 (N_AMR = 7), and RUN25–RUN27 (N_AMR = 8), where in all cases relatively high densities are reached with a Jeans number at sink formation of L_J,s = 2. The increase in the number of levels of refinement results in a better resolution of low-mass stars. This is because higher densities are reached before sink formation, as required for small cores to start collapsing. However, the turbulence is not well sampled in these runs because of a too-small root grid resolution.

As found with the convergence runs in Table 1, by simultaneously increasing the root grid size and the maximal resolution, we can properly populate the low-mass end of the IMF. A very high number of AMR levels is in itself not enough to do it.

To summarize this section, we list a guiding set of rules that we have established for setting the numerical parameters of the sink particle model and the numerical resolution:

• 1.  
  The numerical Jeans length at sink formation should be as high as possible, but in practice it is limited to a rather low ideal value, L_J,s = 2, to avoid the creation of spurious sink particles.
• 2.  
  The exclusion radius should be kept as small as possible. As long as L_J,s = 2, the exclusion radius can be set equal to twice the accretion radius, r_ex = 2r_acc (r_ex = 8Δx in our main simulations where r_acc = 4Δx).
• 3.  
  A sufficiently high resolution root grid and a large enough minimum Jeans number L_J are needed to resolve both the turbulence and the gravitational collapse. A high number of levels of refinement makes it possible to resolve small, high-density cores that are the natal sites of low-mass stars. Therefore, to populate the IMF toward the low-mass end, it is necessary to increase all three quantities, N_root, N_AMR, and L_J, simultaneously.

These principles have been followed to design the convergence runs listed in Table 1, and they have made it possible to demonstrate the trend toward a convergence of the numerical model.