© ESO 2018

1 Introduction

The formation of stars inside a cluster depends on local as well as on global conditions. Since the majority of stars forms inside clusters (e.g., Lada & Lada 2003), clarifying these effects will lead to an advanced understanding of the star formation processes. It is commonly accepted that the high-mass end of the initial mass function (IMF) follows some power law dN∕d logM ∝ M^Γ, where Salpeter (1955) derived that Γ ≃−1.35 with the field star mass spectrum. The shape of the IMF around and below solar masses is quite different. While Kroupa (2001) suggested a power-law behavior with lower values of Γ at lower masses (see also Hillenbrand 2004), Chabrier (2003) proposed a lognormal distribution that peaks at about 0.2–0.3 M_⊙. Bastian et al. (2010) compiled IMFs observed in clusters and reported strong variations of Γ between –0.4 and –2 (see also, e.g., Moraux et al. 2007), and this value is, furthermore, sensitive to the mass range considered.

Many numerical studies have been dedicated to the understanding of the IMF. An important difficulty arises from the need to resolve sufficiently small spatial scales while at the same time considering relatively large spatial scales that are necessary for adequate statistics. Earlier cluster formation simulations (e.g., Bate et al. 2003; Bate & Bonnell 2005; Bate 2005; Clark et al. 2008; Offner et al. 2008, 2009), although valuable in the pioneering illustration of some IMF physics, are limited in either statistics or resolution and practically do not produce well-defined IMFs. More recently, thanks to the increase in computing power, larger simulations that provided more reliable statistics have been performed. Bonnell et al. (2003) simulated an isothermal 1000 M_⊙ cloud of 1 pc diameter, with Mach ~10, and produced an IMF that was broadly consistent with Γ ~−1. Bate (2009a, 2012) simulated 500 M_⊙ clouds of 0.404 radius with Mach number , with either a polytropic equation of state (eos) or radiation hydrodynamics. The mass spectra they presented have Γ slightly higher than –1 in the mass range 0.1–3 M_⊙ in both cases. Maschberger et al. (2010, 2014) typically found Γ ≳−1 in sub-clusters of a 10⁴ M_⊙ simulation inside a cylinder of 10 pc length and 3 pc diameter (Bonnell et al. 2008, 2011). These simulations used a piecewise polytropic eos, with the gas being isothermal at density 5.5 × 10⁻¹⁹−5.5 × 10⁻¹⁵ g cm^-3. Krumholz et al. (2011) simulated 1000 M_⊙ clouds at 2.4 × 10⁵ cm^-3 and studied the effect of radiation hydrodynamics compared to the isothermal condition. Their mass spectra also showed Γ > −1 in the mass range 0.1–2 M_⊙, either with or without radiation. Girichidis et al. (2011) started with various density profiles for 100 M_⊙ molecular clouds at a density 4.6 × 10⁵ cm^-3, with , and generally found Γ ~ −1. When a density profile ρ ∝ r⁻² is employed, an initial Mach number leads to the formation of a single star. Ballesteros-Paredes et al. (2015) simulated a piece of isothermal molecular clouds of 1000 M_⊙ in a 1 pc box for several Mach numbers and obtained an IMF compatible with Γ ~−1, whose shapeis established quite early. When we examine the mass spectra of these studies in detail, statistical fluctuations remain, and it is not obvious howthe value of Γ is determined. It is common to see slopes that are shallower at the intermediate-mass range than the canonical –1.35 inferred by Salpeter (1955).

This work aims to examine the stellar mass spectrum shape as a result of the initial conditions by covering a wide range of initial parameters. We perform a series of numerical simulations of cluster formation inside molecular clouds of 1000 M_⊙ with varying initial conditions, namely density and turbulence level. To overcome the above-mentioned computational obstacles, a compromise is made between the need for high resolution and the need for sufficient statistics to obtain meaningful mass spectra. The first Larson core scale (a few AU) is resolved to ensure a good description of the gravitational fragmentation. Sink particles are employed as a subgrid modeling of stars. The choice of the cloud mass (1000 M_⊙) allows us to produce a statistically significant cluster. With this simple setup, we show that the stellar distribution is environment dependent and is universal only under certain circumstances.

The molecular gas was initially confined in a sphere of radius 0.042−1 pc, with the corresponding number density between 10³ and 10⁷ cm^-3. While noting that this configuration is rather dense compared to general molecular clouds, the very wide density range allows us to illustrate the influence of the initial density. To isolate the effect of initial density, we intentionally left out magnetic field, cooling, radiative transfer, and all stellar feedback effects, while representing the thermodynamics of the gas with a simple smoothed two-slope polytropic eos. These numerical experiments, therefore, are not intended to represent fully realistic molecular clouds.

In the second section, we present the numerical setup and the initial conditions. The results are qualitatively described in the third section. The fourth section is devoted to the stellar mass spectra inferred from the various simulations, while in the fifth section, the model of Hennebelle & Chabrier (2008) is adapted to the case of a collapsing cloud, and we compare theory to simulation results. The sixth section discusses the results, and the seventh concludes the paper.

2 Simulations

2.1 Initial conditions

The simulation box was initialized with a spherical molecular cloud of mass M = 1000 M_⊙ and a density profile , where r is the distance to the cloud center, and ρ₀ and r₀ are the density and size of the central plateau, respectively. The density contrast between the cloud center and edge was a factor of ten, and the radius consequently was 3r₀. The simulation box was twice the size of the cloud, and the surrounding space was patched with diffuse medium of density ρ₀ ∕100. The turbulence was initialized from a Kolmogorov spectrum with random phases. The seed was the same for all simulations, while the amplitude was scaled to match the assigned turbulent energy level and box size. The temperature, T, is given by a smoothed two-slope polytropic eos, such that the gas was at 10 K at low density and followed the dependence T ∝ ρ^2∕3 at a number density higher than 10¹⁰ cm^-3.

Two effects of the cloud initial conditions were studied: the density, and the turbulence level. We performed a series of simulations by varying the compactness of the cloud, defined by the ratio between the free-fall time and the sound-crossing time t_ff ∕t_sc = 0.15,  0.1,  0.05,  and 0.03. The four setups increased in density, while the virial parameter remained the same. This was done by fixing the ratio between the free-fall time and the turbulence-crossing time t_ff ∕t_tc = 1.1, which is close to virial equilibrium. On the other hand, we also set t_ff∕t_tc = 1.5,  0.5,  0.3,  and 0.1, while keeping t_ff∕t_sc = 0.05. The level of turbulence is important since turbulence provides dynamical support against self-gravity and creates local overdensities through shocks. The simulation parameters are listed in Table 1. The number density n = ρ∕(μm_p), where μ = 2.33 is the mean molecular weight and m_p is the atomic hydrogen mass, is presented in the table instead of the volumetric density ρ.

Since the smooth density profile and the seeded turbulent field are not fully self-consistent, we ran a relaxation phase before the actual simulation. For each run, the system was evolved without self-gravity during 0.3 t_tc at the lowest level of refinement (2⁸) to prepare amore realistic density field with local fluctuations. The relaxation effects are discussed in Appendix A.

2.2 Numerical setup

The simulations were run with the adaptive mesh refinement (AMR) magnetohydrodynamics (MHD) code RAMSES (Teyssier 2002; Fromang et al. 2006). The system was evolved following ideal MHD equations. Although these simulations were not magnetized, the MHD solver was employed to avoid numerical uncertainties when the magnetic field might be included in follow-up studies. All simulations were run on a base grid of 2⁸ , corresponding to different physical resolutions for different density cases, except for the densest case D, which was run on 2⁷ . After we switched on gravity, the refinement required that the local Jeans length always be resolved by ten cells. Canonical runs were refined to level 14, and we varied the resolution to check for numerical convergence over up to 4 AMR levels (see Table 1). The lowest resolution in all runs was 38 AU, and the best resolution reached 2 AU.

2.3 Sink particles

Sink particles were used in our simulations as a subgrid model of stars (Krumholz et al. 2004; Federrath et al. 2010a). The algorithm we used wasdeveloped by Bleuler & Teyssier (2014), and it proceeds as follows: first, density peaks are identified. Over a given density threshold (n_sink = 10¹⁰ cm^-3), a sink particle is placed when the local gas properties satisfy the criteria of maximum refinement level, local virial boundedness, and flow convergence. The sink particle then interacts with the gas through gravity and continues to accrete mass from its surrounding. A sink accretes 75 percent of the mass that exceeds n_sink from surrounding cells with the accretion radius, r_acc, that is, four times the smallest cell size. The angular momentum of the accreted gas is added to the sink. The effect of angular momentumin regulating mass accretion or fragmentation may modify the resulting sink particle mass spectrum in the presence of a magnetic field, and it is currently ignored in the this simplified study as the magnetic field is not considered, and even if existed, it is thought to have no strong influence on the core at such small scales. The pre-stellar core (a few tens of AU) is resolved in the simulations, and this guarantees that the sink particles represent individual stars.

2.4 Missing physics

We performed an ensemble of numerical experiments to quantify the influence of initial conditions and spatial resolution, and the physics used here is greatly simplified. Many other processes not included here are nonetheless important, and may even be dominant in some circumstances. For example, the role of the accretion luminosity that emanates from the protostars has been stressed by Krumholz et al. (2007), Bate (2009b), and Commerçon et al. (2011). Accretion luminosity provides a substantial heating source to the gas, which leads to greater thermal support and may reduce or even suppress the fragmentation, that is, the formation of several objects. The magnetic field also likely plays an important role (Hennebelle et al. 2011; Peters et al. 2011; Myers et al. 2013). In particular, due to magnetic braking and magnetic support, it tends to reduce the fragmentation within clusters. This discussion may become even more crucial in the context of massive disks, which, as we show, form in some of the simulations presented below. A magnetic field may drastically modify the picture there (Hennebelle & Teyssier 2008; Machida et al. 2008), even when non-ideal MHD processes are included (Masson et al. 2016; Hennebelle et al. 2016; Wurster et al. 2017). The same is true for the accretion luminosity, which can also efficiently heat the disk and stabilize it against fragmentation (Offner et al. 2009; Commerçon et al. 2010), although Stamatellos et al. (2012) argued that intermittent accretion may limit the effect of radiative heating.

3 Overall evolution

3.1 Qualitative description

3.1.1 Large scales

Figure 1 shows the column density along the z-axis for three snapshots of simulations A1++ (top panels)and C1+ (bottom panels). The first snapshot corresponds to an early time where less than 20 M_⊙ of gas has been accreted onto sink particles, while for the later times, 200–300 M_⊙ of gas has been accreted.

At time 0.155 Myr, a network of dense filaments is already present in run A1++. The sink particles, represented by red dots with sizes proportional to the particle masses, are located within these filaments and sometimes even at the intersection between filaments. This trend becomes clearer at 0.203 Myr, while the correlation appears less tight at 0.424 Myr. This is because the sink particles have accreted a substantial fraction of the initial cloud (about one third of the cloud mass) and have also undergone N-body interactions.

At time 0.003 Myr, run C1+ displays significantly more filaments than run A1++ because the thermal support is lower in run C1+. The sink particles are also located inside filaments, and they tend to be more broadly distributed as a result of the widespread filament formation. This is even more obvious at time 0.004 Myr.

3.1.2 Small scales

Figure 2 shows three zooms for simulations A1++ (top panels) and C1+ (bottom panels) at various times.

The zooms corresponding to simulation A++ reveal two massive fragmenting disks at time 0.155 Myr (top left panel, located at x ≃ 1.195, x ≃ 1.21 and y ≃ 1.4 pc) seen face-on and one disk at time 0.424 Myr (top right panel, located at the center) seen edge-on. Their typical size is a few hundreds of AU, and they likely exhibit a rotating pattern to which several sinks are associated. In the two cases, the disks are relatively isolated without other objects in their immediate vicinity (see also Fig. 1). On the other hand, the zoom at time 0.291 Myr (top middle panel) shows the presence of a many sink particles, with no sign of rotating and fragmenting structures. Small disks can be seen around some sink particles. The trend observed here is typical of the situation in run A. Massive fragmenting disks are relatively common in regions that are not too crowded, but rare, if not absent, in regions where numerous sink particles are found. This suggests that in run A1++, disk fragmentation is a significant mode that leads to the formation of new objects (however, see Sect. 2.4).

The zooms corresponding to run C1+ (bottom panels) also show disks (for example, three disks are located at the centers of each of the three panels). These disks, however, show no signs of fragmentation, even though the left and middle panels are relatively isolated. This trend appears to be generic for run C. Disk fragmentation does not seem to be an important mode. The reason is that disks are much smaller, with a typical size on the order of a few tens of AU. The absence of large disks is most probably due to the high number density of sink particles, which produces dynamical interactions that lead to efficient transport of angular momentum. Disk truncation by stellar encounter has indeed been recognized as an efficient process (Clarke & Pringle 1993; Breslau et al. 2014; Jílková et al. 2016).

Fig. 1 Column density maps of runs A1++ at 0.155, 0.203, 0.424 Myr (top row from left to right) and C1+ at 0.003, 0.004, 0.008 Myr (bottom row from left to right). Run A1++ is more globally collapsing, and the sinks, initially forming inside the filamentary network, are more concentrated in a central cluster. Run C1+ is initially much denser and the evolution of the filamentary structure is more pronounced, with more widespread sink formation.

Fig. 2 Zoomed column density maps of runs A1++ (top row) and C1+ (bottom row). Red dots represent the sink particles, and their sizes correspond to the sink mass. Disk formation around sink particles is observed in both cases, while the disks in run A1++ are more extended, showing likely signs of rotation and fragmentation, and the disks in run C1+ are significantly smaller, with no sign of fragmentation.

3.2 Density PDF

Figure 3 shows the probability distribution function (PDF) of the density field at four time steps and for six simulations, corresponding to four initial densities and three initial turbulence levels. The PDFs are very similar to one another. They present a peak that is roughly n₀∕400, with n₀ being the central density, and then a power law at high densities with a slope of about –1.5. In Table 1, we list the position of the density PDF peak, n_{pdf peak}, because it may be more representative than the initial peak density, n₀, in particular, because we let the cloud relax before running the simulation with self-gravity. For most cases, the PDF is almost invariant in time. A slight decrease in turbulence level does not make much difference (C1t03), while with significantly lowered turbulence (C1t01), the global collapse shifts the PDF toward high densities in time. This PDF shape is indeed expected and has been inferred in previous calculations (Kritsuk et al. 2011; Hennebelle & Falgarone 2012; Girichidis et al. 2014). The n^−1.5 power law stems from the n ∝ r⁻² density profile that develops during the collapse (e.g., Shu 1977). At very high density, that is to say, n > 10¹⁰ cm^-3, the PDF becomes shallower. This is a consequence of the eos, which becomes adiabatic. At low density, that is, below the mean density, we see a behavior that broadly looks lognormal with a plateau near the peak and then a stiff decrease.

Fig. 3 Density PDF (normalized in random units) for runs A1++, B1++, C1+, D1, C1t03, and C1t01. At low densities, the shape of the PDF is broadly lognormal, while at high density, it is ∝ n−1.5 . The PDF is almost invariant in time given that the cloud is sufficiently supported, as it is in most cases. In run C1t01, the initial turbulent support is very low and thus the cloud global collapse is strong, giving a density PDF that shifts toward high densities in time.

3.3 Sink particle analysis

3.3.1 Density profiles around sink particles

Figure 4 shows the mean density profile around sink particles in simulation A1++ (top panels) and C1+ (bottom panels). Theprofile is averaged among all radial directions and among sinks, starting from the distance of twice the smallest cell size. When there exists a nearby sink, the density profile is considered only up to half of the distance in between. The solid lines represent the mean value while the shaded areas show the standard deviation. Three ranges of mass are considered, below 0.01 M_⊙, between 0.01 and 0.1 M_⊙, and between 0.1 and 1 M_⊙. For A1++, results are presented for sinks younger than 1000 yr (left), 5000 yr (middle), and all sinks (right), while for C1+, 300 yr (left) and 1000 yr (middle) were used for age filtering.

Several interesting trends can be inferred. First of all, at younger ages (i.e., left and middle panels), the density profiles are broadly ∝ r⁻² and proportional, by a factor of a few, to the density profile of the singular isothermal sphere (SIS), , where the isothermal sound speed c_s = 200 m s⁻¹ at 10 K (e.g., Shu 1977, n_SIS = ρ_SIS∕(μm_p) with the black solid line). Second, the density at early times (left panels) is significantly above n_SIS , particularly between 20–50 AU, where it can be as high as ~20 n_SIS. This value will be used in the companion paper (hereafter paper II) to analyze the gas surrounding the first Larson core and estimate the expected peak of the stellar distribution. At later times (middle and right panels), the density drops to lower values that are comparable to or even lower than n_SIS. There is large dispersion both due due the age dispersion of sinks and the asymmetric environment in which older sinks are found. This is particularly the case in the right panels that include sinks of all ages, since the sinks decorrelate from the local overdensity in which they form at later evolutionary stage. The sink particle environment changes rapidly, that is to say, in a few kyr.

3.3.2 Time-dependence of accretion

Figure 5 shows sink mass against accretion time for runs A1++, B1++, C1+, and D1. For the four runs, the final sink mass, M_f , and the accretion timescale, t_acc,60, taken as the time that a sink needs to accrete 60 percent of its mass, although with a large dispersion, are correlated, and t_acc,60 increases with M_f on average.

In the case of run A1++, M_f and t_acc,60 are nearly proportional. Interestingly, there is an excess of points at M_f ≃ 10⁻² M_⊙ and t_acc,60 ≃ 3 × 10³ yr. The abundance of low-mass sinks is a possible signature of the disk fragmentation discussed above. For run B1++, there is no clear sign of bi-modality. The relationbetween M_f and t_acc,60 is still broadly linear, but a deviation seems to occur at high mass. Finally, runs C1+ and D1 present similar behaviors, while the latter exhibits a slightly larger dispersion. There is a flat dependence of M_f on t_acc,60 for M < 0.1 M_⊙ and a much stiffer one for larger M_f with . This behavior is quantitatively interpreted in Sect. 5.5.1.

4 Mass spectra of sink particles

Figures 6 and 7 show the mass functions of sink particles at total accreted mass of 20,   50,  100,  150,  200,  and 300 M_⊙, where applicable, since some runs are less evolved. To assess the numerical convergence of the runs, we performed runs at several spatial resolutions for most cases, and they are presented in the same row. This allows us to compare the various runs at more or lessthe same physical resolution, therefore asserting that the differences are not mere consequences of different resolutions. The stellar distribution from the simulations were fitted with lognormal distributions, and their characteristic mass and width of distribution are listed in Appendix B. The purpose was to measure the position of the peak in a systematic way. We also drew power-law distributions, which constitute good fits above ~ 0.1 M_⊙.

As a general remark, the shape of the mass spectrum is determined rather early for all runs, that is to say, even when only 20 − 50 M_⊙ have been accreted. This is compatible with the analysis in Sect. 3.3, where it has been found that the accretion time is short (<10⁵ yr for model A and <500 yr for model C in Fig. 5) and that once the sinks have decorrelated from their environment (sinks in the low-density environment in Fig. 4), they no longer accrete significantly.

Fig. 4 Density profiles around sink particles for runs A1++ at 0.26 Myr (top panels) and C1+ at 0.008 Myr (bottom panels). From left to right, sinks aged under 1000 yr, 5000 yr, or all sinks are used to produce the average profiles for run A1++, while 300 yr and 1000 yr are used as thresholds for run C1+. Sink particles are separated into three mass bins <0.01 M⊙ (blue), 0.01–0.1 M⊙ (green), and 0.1–1 M⊙ (red), average age of selected sinks is shown in the legend. Shaded regions represent the standard deviation of the density profiles for all sinks and all radial directions. The SIS density profile, , is plotted in black for reference. The density profile at younger ages is always broadly ∝ r−2 , with a factor~10 higher than the SIS, while sinks of older ages have probably accreted all mass from the core in which it is embedded or have been ejected, thus decorrelating from the r−2 profile. Thisis reflected by the large dispersion that is particularly evident with more massive sinks.

4.1 Dependence on density

The first row of Fig. 6 shows results for runs A1. In run A1++, several features appear. First of all, the mass spectrum is essentially flat between 0.1 and 3 M_⊙. For masses >3 M_⊙, there is a stiff decrease. At masses ≃ 0.01 M_⊙, we see an excess of objects, as we noted in Sect. 3.3, that is likely a consequence of disk fragmentation and is possibly a spurious effect in the simulation that may disappear with more complete physics. At lower resolutions, runs A1+ and A1 show mass spectra that are compatible with the spectrum of run A1++, although they are noisier.

The second row of Fig. 6 shows results for runs B1. In run B1++, we can distinguish a flat mass spectra for 0.05 M_⊙ < M < 0.3 M_⊙ and then a power-law type distribution at higher masses with an exponent ≃ −0.5 to –1. At masses M > 2 M_⊙, there might even be another regime, described by a stiffer power law (this remains to be confirmed), as shown by the yellow and black histograms. The shapes of the mass spectra of run B1+ are broadly compatible with those of B1++, but the flat part is absent, and instead a peak appears at 0.1 M_⊙. The mass spectra of run B1 are significantly flatter than those of B1+ and B1++, illustrating the need for resolution and careful numerical convergence tests.

The third row of Fig. 6 shows results for runs C1 (see also the right panel of the bottom row). In run C1+, a peak appears at 0.1 M_⊙, while at high mass, a power law with an exponent equal to about –0.75 is inferred. Run C1 exhibits very similar mass spectra, with the peak and the power law almost identical. Even the total number of objects is comparable (when comparing at 200 accreted M_⊙). Thus it seems reasonable to claim that numerical convergence has been reached for runs C1+ and C1. Runs C1– and C1– – also show a peak at 0.1 M_⊙, while the power-law behavior is noisier. The slope may be slightly flatter for C1– –.

The bottom left panel of Fig. 6 shows results for run D1. The mass spectra are very similar to those of run C1 with a peak at 0.1 M_⊙ and a power law with exponent about –0.75 to –1. There are possibly more low-mass objects (<0.1 M_⊙) than in run C1+.

4.2 Dependence on turbulence

The runs with various levels of turbulence are presented in Fig. 7. Run C1t15 has twice more kinetic energy than run C1, while runs C1t03 and C1t01 have 10 and 100 times less kinetic energy, respectively. Increasing the turbulence level to slightly super-virial (C1t15) does not make any visible difference the mass spectra, which present the same peak and power-law behavior as run C1.

For run C1t03, the peak of the mass spectra seems to be shifted to slightly higher masses, say, tentatively ≃ 0.2 M_⊙ and there are slightly fewer objects. The exponent of the power law is unchanged. When the turbulence is lowered even more (C1t01), one dominating heavy star is formed at the center. This means that less mass is available to form other objects, and their number is significantly reduced. The mass spectra are now very flat over almost two orders of magnitude in mass.

5 Analytical modeling and physical interpretation

5.1 General formalism

The theory developed by Hennebelle & Chabrier (2008, 2009, 2013) consists of identifying the mass that at scale R is gravitationally unstable, that is, the mass contained within the clumps in which gravity dominates other supports. This is achieved by writing that at all scales the mass at a density that exceeds a scale-dependent threshold (obtained through the virial theorem) is equal to the mass contained in the structures of mass lower than or equal to the (turbulent) Jeans mass.

As discussed in Hennebelle & Chabrier (2013), a difficult question is to what extent the distribution should take into account the crossing timescale of the density perturbations. Small perturbations can indeed be rejuvenated several times while the large perturbations are still collapsing. In the case of cold dark matter halos, there is no such contribution because the fluctuations are imprinted at the beginning of the Universe and are then amplified by gravity. In the case of a turbulent molecular cloud, where a turbulent cascade proceeds, it seems that the small scales are continuously reprocessed.

In the present case, the clouds we considered collapse in about one free-fall time. Moreover, it is unclear whether a turbulent cascade truly develops because the crossing time is typically a few free-fall times. In this context, it therefore seems reasonable not to weight the mass spectra by local free-fall time. As we show below, this is entirely consistent with the numerical results.

As in Hennebelle & Chabrier (2008), we write that (1)

where V _c is the cloud volume, is the mean density, is the density PDF, and the local fluctuation is . This equation is the same as the one inferred in Hennebelle & Chabrier (2008). The first expression for M_tot(R) represents the mass contained within structures of mass . It is equal to the mass of the gas, which, smoothed at scale R, has a logarithmic density higher than a critical threshold . The second expression signifies that the number of structures of mass is , where is the number density of structures ofmass between M^′ and M^′ + dM^′ and P(R, M^′) is the probability to find a gravitationally unstable cloud of mass M^′ embedded inside a cloud of gas, which at scale R has a logarithmic density higher than . It has been demonstrated by Hennebelle & Chabrier (2008) that P(R, M^′) = 1 is a reasonable assumption.

Taking the derivative of Eq. (1) with respect to R, we obtain (2)

We note that we have ignored here any dependence of on R, which would introduce a second term that does not play a significant role in most circumstances.

When the expression of the critical density threshold, , is specified, the mass spectrum of the self-gravitating pieces of fluid can be inferred using the geometrical expression (3)

where typically C _m = 4π/3.

5.2 Density probability function

As suggested by Eq. (2), the density PDF plays a major role, and it is necessary to use the right PDF to obtain reliable results. In particular, most analytical models have been using a lognormal PDF, which is not a valid approximation for collapsing clouds.

5.2.1 Lognormal distribution

Numerical simulations of turbulent isothermal clouds have revealed (Vuez-Semadeni 1994; Padoan et al. 1997; Kritsuk et al. 2007; Federrath et al. 2008) that a reasonable approximation of is the lognormal distribution (4)

where and

In this expression, is the Mach number and b a non-dimensional coefficient that depends on the forcing, which typically varies between 0.25 when the forcing is purely solenoidal and almost 1 when the forcing is applied on compressible modes (Federrath et al. 2010b). For simplicity and because it is usually not so important, we neglect here the scale dependence of σ₀ (see, e.g., Hennebelle & Chabrier 2013, for a discussion).

5.2.2 Power-law distribution

When the collapse proceeds, the density PDF that develops in the inner part of the clump is simply (Kritsuk et al. 2011; Hennebelle & Falgarone 2012; Girichidis et al. 2014) (5)

where and ρ₀ are constants of normalization.

5.3 Instability criterion

The last step is to specify the density threshold by writing the virial theorem. As shown in Hennebelle & Chabrier (2008, 2013), the condition for a density fluctuation at scale, R, to be unstable is (6)

with a_J being a dimensionless geometrical factor of order unity. Taking for example the standard definition of the Jeans mass, the mass enclosed in a sphere of diameter equal to the Jeans length, we obtain a_J = π^5∕2∕6. With Eq. (3), this implies (7)

where is the critical massat scale R. In Eq. (6), η is the exponent that enters the velocity dispersion, and its value is typically expected to be 0.4–0.5, and V ₀ is the velocity dispersion at 1 pc (and not at the cloud scale).

After normalization, Eq. (7) becomes (8)

where . The Jeans mass , Jeans length , and Mach number at the cloud scale are given by (9)

5.4 Asymptotic analysis

5.4.1 Expected behavior during collapse: two support regimes

In the case of a the power-law density PDF, we obtain the mass spectrum as a function of mass with Eqs. (2), (5), and (8). To gain further insight, we performed an asymptotic analysis by writing (10)

which leads to (11)

If the thermal support is dominant, then η = 0, while η ≃ 0.5 if the turbulent dispersion is dominant. From Eq. (2), we obtain (12)

Therefore, two asymptotic behaviors are expected. If the thermal support is dominant, we expect , while if the turbulent support is dominant, then (assuming η = 0.5). This leads to dN∕d logM ∝ M⁰ and dN∕d logM ∝ M^−3∕4, respectively.

This qualitatively matches the numerical results very well since simulations A1, which have the lowest initial density and therefore the highest thermal support, present a flat mass spectrum, that is, d N∕d logM ∝ M⁰, while runs C1 and D1, which have the densest initial conditions and therefore the lowest thermal support, present mass spectra entirely compatible with dN∕d logM ∝ M^−0.75 for M > 0.1 M_⊙.

If the weighting by the free-fall time were taken into account, that is to say, if we were to account for the possibility that several generations of low-mass stars might form while massive stars are forming, the mass spectrum would need to be multiplied by . Then we would infer dN∕d logM ∝ M⁻¹, regardless of the value of η and of the relative importance between thermal and turbulent support. From our results, this does not seem to be the case.

5.4.2 Lognormal case

For comparison, it is worth recalling the asymptotic behaviors that have been inferred for the case of lognormal PDF. Hennebelle & Chabrier (2008) inferred (in the region where the PDF is relatively flat, that is to say, where δ² inside the exponential of Eq. (4) does not dominate) that (13)

which leads to dN∕d logM ∝ M⁻² if the thermal support dominates the turbulent dispersion and dN∕d logM ∝ M^−1.25 otherwise. Taking into account the weighting by the free-fall time would modify these conclusions. Hennebelle & Chabrier (2013) inferred that (14)

leading to to dN∕d logM ∝ M⁻³ if the thermal support dominates and dN∕d logM ∝ M^−1.5 if the turbulent support dominates.

It might be surprising at first that for a lognormal PDF, thermal support leads to spectra that are stiffer than that with turbulent support, while the reverse is true for the power-law PDF ∝ ρ^−1.5. The reason lies in the dependence of on ρ (see Eqs. (12)–(14)) and the density-mass relation inferred from the virial theorem (Eq. (11)). If the PDF decreases too stiffly with density, the dense material is rare, and since for a specific clump mass, thermal support leads to clumps more diffuse than turbulent dispersion, thermal support tends to make larger clumps. In this respect, the critical exponent for the PDF is –1 if free-fall time is not taken into account and –1.5 otherwise.

5.5 Comparison with the simulation results

To compare the simulation results with the analytical predictions, we need to specify a mean density and a Mach number. For the first, we use a radius,R_c, 30 percent larger than the radius indicated in Table 1 since this latter value corresponds to the radius before the relaxation phase. A velocity dispersion also has to be specified in order to derive the Mach number. A difficulty arises because a substantial part of the motions is due to the collapse itself, and it is not straightforward to distinguish infall and turbulence. For the cloud velocity dispersion, we take . For a virialized cloud, we may have chosen , while infall should be excluded for our present purpose of estimating the turbulent dispersion. This leads to the normalized velocity dispersion at 1 pc (15)

5.5.1 Comparison of the accretion times

We first compare the accretion time-scale displayed in Fig. 5 and the prediction of the formalism. For this purpose we use Eqs. (3) and (8). Since the free-fall time is given by , we have a link between the accretion time, assumed to be the free-fall time, and the mass. Before a quantitative comparison, it is enlightening to consider the asymptotic behavior. Writing again M ∝ R^1+2η, we obtain (16)

Two regimes appear: if thermal support dominates turbulent dispersion, we obtain M ∝ τ_ff, while is obtained with η = 0.5 in the opposite case (). This qualitatively matches the behaviors observed in Fig. 5 very well since for run A1++, M ∝ t_acc,60, while for run C1+ and D1, a much stiffer relation is obtained.

Using the complete Eq. (8), while applying 1.3 R_c and V ₀ from Eq. (15), we plot M against 0.6 τ_ff for the fourruns in Fig. 5 (solid black lines). It is justified with a ρ ∝ r⁻² density profile that the fraction of accreted mass is proportional to the collapse time. Sinks of mass close or below 0.1 M_⊙ are strongly affected by the adiabatic equation of state, as described in paper II, and the physics that governs their formation and timescale is not the one described by the present analytical model. Therefore, the present comparison is meaningful only for masses ≳ 0.1 M_⊙. With this in mind, the model provides a very reasonable fit for the four runs at masses ≳ 0.1 M_⊙.

The stiff dependence of the mass on the accretion time for run C1+ and D1 is well described by the analytical model. This is particularly interesting because this behavior is directly linked to the turbulent dispersion (which in a loosely sense can be called turbulent support). This implies that the origin of the massive stars is indeed a consequence of the coherence induced by the support, which can be either thermal or turbulent, and is not due to purely competitive accretion.

5.5.2 Comparison of the mass spectra

The analytical prediction is obtained by combining Eqs. (2), (5) and (8) with V ₀ given by Eq. (15). For the sake of comparison, we also computed the predicted mass spectra employing a lognormal density PDF. Figure 8 shows the results for four runs. Mass spectra at the time when the total sink mass equals to 150 M_⊙ were used for comparison. Since the shape of the mass spectrum varies little with time, we normalized the number of objects from the simulations to match the number predicted by the analytical model. The analytical model with the power-law density PDF stated by Eq. (5) and the mass spectra inferred from simulations, solid and dashed blue lines, respectively, generally agree pretty well for M > 0.1 M_⊙, where the analytical models are applicable (see paper II), except for models B1, C1, and D1 at high mass, where the lognormal PDF seems to provide better fits. The changing slope of the mass spectrum around ≃ 1 M_⊙ may be a consequence of the PDF, which shows a transition between a lognormal behavior at low density and a power law at high density (see Fig. 3).

To conclude, simulations and analytical models both predict the power-law change from d N∕d logM ∝ M⁰ in the regime dominated by thermal energy (model A1) to dN∕d logM ∝ M^−3∕4 in the turbulence-dominated regime (model C1 and D1), without any adjustment of free parameters.

6 Discussions

6.1 Comparison with previous works

Through the series of simulations A1 through D1 with increasing initial density, we have demonstrated that there are indeed two regimes of stellar distribution formation in massive collapsing clouds. First, the stellar distribution is flat in the low-density regime, where thermal energy is the dominant support against self-gravity. Second, d N∕d logM ∝ M^−3∕4 in the high-density regime, where turbulence is dominant. Alternatively, varying the turbulence strength in runs C1 at fixed density also shows this effect: the stellar distribution becomes flat as the initial turbulent energy is lowered. Table 2 lists a summary of simulation parameters from the literature and the regimes investigated in these works.

Recently, two studies by Bertelli Motta et al. (2016) and Liptai et al. (2017) concluded that in their numerical experiments, the mass spectrum is not influenced by the initial turbulence, which is expected to have an effect on the IMF, according to the reservoir-based theories proposed by Hennebelle & Chabrier (2008), Schmidt et al. (2010) and Hopkins (2012). Our results are partly compatible with their conclusion because we see from runs C1, C1t15, and C1t03 that varying the turbulence by a factor 5 does not lead to very significant change in mass spectrum. The reason is that self-gravity triggers the development of a power-law tail that is not considered by Hennebelle & Chabrier (2008) and Hopkins (2012). When the density PDF relevant for gravitational collapse is considered, the comparison between the reservoir-based theory and the simulation results becomes very good for M > 0.1 M_⊙ (near and below this value, the physics to be considered is different, as described in paper II), as shown in Fig. 8. However, the mass spectra produced in our numerical simulations tend to be too shallow compared to the Salpeter exponent. Therefore, while turbulence could possibly be less important in setting the low-mass part of the IMF than originally proposed by Padoan et al. (1997), Hennebelle & Chabrier (2008), and Hopkins (2012), it is likely playing a determining role for the massive stars.

A numerical experiment comparable to ours has been conducted by Ballesteros-Paredes et al. (2015), who simulated a 1000 M_⊙ cloud at density 1.7 × 10⁴ cm^-3 with Mach number varying between 4 and 16. The sink mass spectrum from their study is rather flat at the beginning, while a power law of Γ ~−1 develops very quickly and there is no obvious difference between different initial conditions. Their conditions are broadly compatible with our runs B1 (or probably would stand between our runs B1 and C1). Given that the slope that they inferred is not very accurately determined, Γ ≃−1 is close enough to the value − 3∕4 we that obtain, especially because for a high mass, we obtain some slightly steeper mass spectra (see Fig. 8).

Another similar set of parameters has been used by Bate (2009a), which again falls between runs B1 and C1. From Fig. 3 of Bate (2009a), the mass spectrum between 0.1 and 1 M_⊙ seems to be shallower than the Salpeter exponent of − 1.35 and compatible with Γ = −0.75. For higher masses, the statistics are less significant, but broadly compatible with Γ ≃−1.

Furthermore, Girichidis et al. (2011) also find mass spectra that are compatible with Γ larger than or near –1 (particularly run TH-m-1 in their Fig. 8). Since they consider 100 M_⊙, it is not easy to compare their results with our runs, however.

Finally, one important question is why the regime Γ ≃ 0 has not been reported in the works listed in Table 2. The conditions for Γ ≃ 0 are not typical of present-day molecular clouds. However, this regime may have been observed in the context of primordial stars. The collapse of minihaloes, thought to be the progenitors of Pop III stellar systems, has been studied and flat mass spectra have been inferred with Γ ≃−0.2 (Stacy & Bromm 2013; Stacy et al. 2016). Since the temperature of the gas in these minihaloes is on the order of 1000 K, the thermal support is much higher.

6.2 Astrophysical interpretation: competitive accretion and mass reservoir

A long-standing debate is concerned with the way stars acquire their mass, and based on collapse calculations similar to those performed in this work, the role played by accretion has been emphasized, leading to the model named competitive accretion (e.g., Bonnell et al. 2001, 2004). The underlying idea is that massive stars tend to attract gas through their gravity and might accrete as d M∕dt ∝ M².

Another type of models has been emphasizing the role of the mass reservoir (Padoan et al. 1997; Hennebelle & Chabrier 2008; Hopkins 2012). The fundamental difference with the competitive accretion scenario, which is stochastic in nature, is that the mass of the reservoir largely determines that of the star. Therefore, the mass of the stars is determined before their formation. The analytical model developed in this paper, identical to the model proposed by Hennebelle & Chabrier (2008) except for the density PDF, belongs to this category. The mass reservoir that predetermines the mass of the star is a direct consequenceof the thermal support and turbulent dispersion. The comparisons between our simulations and the analytical model (Figs. 5 and 8) show very good agreement (we recall that no free parameter had to be adjusted). This reinforces the role of the mass reservoir and suggests that it may be dominant.

This does not mean that the conversion from the gas reservoir to the star is entirely deterministic and that there is no competition at all. The timescales displayed in Fig. 5 present significant dispersions as a consequenceof all fluctuations in the system, which certainly include some kind of competitive accretion. The mean valuesand the general trends, however, can be well reproduced by a deterministic model that emphasizes the role of apreexisting coherent reservoir of gas. Another important aspect is the choice of which particular star will accrete most of the mass of the reservoir in which it is embedded. It is likely the case that inside a turbulence-dominated reservoir, several thermally dominated reservoirs can be identified, leading to several low-mass stars that may eventually accrete the mass of the larger turbulence-dominated reservoir and become massive stars.

Fig. 7 Mass spectrum evolution for models C1t15, C1t03, and C1t01 corresponding to three levels of initial turbulence. Except for C1t01, which initially has a very low turbulent energy, the mass spectra are not significantly different.

Fig. 8 Mass spectra as predicted from analytical model for four sets of parameters corresponding to runs A1++, B1++, C1+, and D1. Simulation results at 150 M⊙ accreted (blue dashed, same as the magenta histogram in Fig. 6 with absolute value shifted to compare to models), model with power-law (blue solid) and lognormal (red dashed) density PDF are plotted, as well as the power-law relation (–3/4 in cyan and –1 in black) for reference. The vertical dot-dashed line indicates a mass lower limit above whichthe model is applicable (see paper II).

It is worth stressing that massive stars form in the same way as low-mass stars. They simply result from larger reservoirs. It is true, however, that in many circumstances, the reservoir of the massive stars will be predominantly determined by the turbulent dispersion (the velocity dispersion appearing in the virial criterion expressed by Eq. (6)) instead of by the thermal support. Since turbulence presents many fluctuations, this reservoir is less clearly determined than the reservoir dominated by the thermal support.

Fig. 5 Final sink mass against accretion time (tacc,60 is used) for runs A1++, B1++, C1+, and D1. The solid black lines represent the free-fall time of the mass reservoir derived from the analytical model (see Sect. 5.5.1). The behavior of the distributionand the analytical model are very similar except for short accretion times (corresponding to mass ≲0.1 M⊙), where the physics to be considered is different. This strongly suggests that the mass that built the stars essentially comes from a preexisting deterministic reservoir.

Fig. 6 Sinkmass spectra for models A1, B1, C1, and D1 at various time steps, corresponding to a comparable amount of accreted mass. Various spatial resolutions are explored for the purpose of verifying numerical convergence and direct comparison between models. For models B1, C1, and D1, a peak at ≃ 0.1 M⊙ is obtained independently of the numerical resolution. While model A1 presents a flat mass spectrum, models C1 and D1 present at high mass a power-law behavior that is compatible with M−3∕4 (see also Fig. 8).

6.3 Different regimes leading to the IMF

The results presented in this paper and its companion suggest that the IMF is a consequence of several different physical regimes. Here we try to clarify the situation by introducing some terminology and classification. Unfortunately, a difficulty arises because two independent types of processes play a role and complicate the situation.

First of all, we have to distinguish between the regimes for which the thermal support is dominant and the regime for which the turbulent dispersion dominates. Mathematically speaking, this is the respective amplitude of the two terms appearing in the right-hand side of Eq. (7). Second, there are also three density regimes that need to be distinguished. They broadly correspond to three ranges of mass, the most massive stars, intermediate-mass stars, and low-mass stars, or density regimes I, II, and III. Regime I corresponds to the lognormal part of the density distribution, regime II the power-law density distribution, and regime III to the density at which the gas is adiabatic (formation of the first Larson core). Therefore, density regimes I and II must be subdivided into thermally and turbulence-dominated regimes, or regimes Ia, IIa, and Ib, IIb. Below, the term density regime refers to regimes I, II, or III, while the term regimerefers to regimes Ia, Ib, IIa, IIb, or III (i.e., a combination of the density PDF and support against gravity).

Density regimes I and II have in common that they present a preexisting deterministic mass reservoir that is due either to thermal support or to turbulence dispersion and gravity. Since the density PDF of density regime I is set by turbulence, while the density PDF of density regime II is due to collapse and gravity, they are turbulently and gravitationally driven, respectively. Regime III is different and is described in paper II. It relies on the physics of the so-called first Larson core, which is related to the transition from isothermal to adiabatic phase and to the density and temperature at which the second collapse occurs. In general, there is no clear relation to a pre-existing coherent mass reservoir in density regime III. We stress that simulations A1 essentially present regimes IIa and III. Simulations C1 and D1 are essentially IIb and III with a possible transition to regime Ib at high masses. Simulation B1 is likely to represent regimes IIb, IIa, and III.

The existence of these different density regimes may be a strong clue to understanding the apparent universality of the IMF (Bastian et al. 2010). Regime III is largely universal, except in terms of metallicity, which may affect the properties of the first Larson core (Masunaga & Inutsuka 1999). Regime Ib, in which the lognormal density PDF is set by turbulence and turbulent dispersion dominates thermal support, is likely the relevant regime (at least in the present-day Universe, and it could be different for Pop III, see Stacy & Bromm 2013) in a wide range of conditions and is also expected to be universal because gravity and turbulence are universal processes. The transition between density regimes I and II and the existence of regimes Ia and IIa is not universal, however. Depending on the local density and Mach numbers, the masses for which the reservoir is built from alognormal density PDF and those for which it is built from a power-law density PDF may vary. However, since the difference between the mass spectra expected in regimes Ib and IIb is small, they may not be easy to distinguish observationally, especially because the statistics on a particular cluster usually present some noise.

Finally, we stress that the existence of these various regimes is compatible with observations. For example, in Fig. 11 of Alves de Oliveira et al. (2013), the mass spectrum is compatible with a cut at a mass lower than 0.1 M_⊙ (regime III), a power law with Γ ≃−3∕4 between 0.1 an 1 M_⊙ (regime IIb), and then a power law with a Salpeter exponent at higher mass (regime Ib).

7 Conclusions

We investigated the fragmentation of a collapsing 1000 M_⊙ clump and the resulting stellar mass distribution. In this paper we focused on the power-law part at the high-mass end of the distribution,which is related to the gas isothermal phase, while in paper II, we investigated the origin of the peak of the distribution,which is is a consequence of the adiabatic eos at high density. We performed a series of numerical simulations in which we varied the initial density by four orders of magnitude and the turbulence energy by two orders of magnitude. We verified numerical convergence by systematically varying the resolution. To assess the physical processes at play in the simulations, we developed an analytical model that we compared to the simulation results. The model was adapted from Hennebelle & Chabrier (2008) and considered a power-law density PDF that was appropriated for a collapsing clump instead of a lognormal one, valid for a non-collapsing turbulent medium.

We found two different regimes. When thermal support is dominant, it is inferred that d N∕d logM ∝ M⁰, while when the local turbulent dispersion dominates, we obtain dN∕d logM ∝ M^−3∕4. As explained in paper II, this is valid for masses higher than about ten times the mass of the first Larson core, M_L , around which the thermodynamics and the tidal forces due to the first Larson core are essential. For masses lower than a few M_L in our simulations, the mass spectrum starts to decrease with the mass, except for the most diffuse clouds, where disk fragmentation leads to the formation of objects down to a mass of M_L , that is, a few 10⁻² M_⊙. We stress, however, that the physics included in the present simulations is too simplistic regarding disk formation and fragmentation.

While the mass spectra of the densest clouds qualitatively resemble the observed IMF, for masses higher than 0.1 M_⊙ they exhibit an exponent that is shallower than the Salpeter exponent of − 1.35. Nonetheless,we observe a possible transition toward a slightly steeper value that is broadly compatible with the Salpeter exponent for masses above a few solar masses. This change in behavior is associated with the change in density PDF, which switches from a power-law distribution to a lognormal. Therefore, our results suggest that while gravitationally induced fragmentation may be important for low masses, turbulently induced fragmentation is likely responsible for setting up the IMF above a few solar masses.

Acknowledgements

The authors thank the anonymous referee for a report that improved the manuscript. We thank Gilles Chabrier for a critical reading of the manuscript that improved the clarity of the paper significantly. This work was granted access to HPC resources of CINES under the allocation x2014047023 made by GENCI (Grand EquipementNational de Calcul Intensif). This research has received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013 Grant Agreement No. 306483).

Appendix A Influenceof the initial fluctuations

One conventional method to initialize a molecular cloud in simulations is to give a density profile, usually smooth, and introduce randomly seeded turbulence, in which case the cloud is artificially constructed such that the velocity and density fields are not self-consistent. We generate relaxed initial conditions by evolving the cloud with the hydrodynamic equations without considering self-gravity during a fraction of the system turbulence-crossing time. This relaxation generates local density fluctuations that create a self-consistent correlation between density and velocity fields. Since it is important to verify that the results do not depend on the initial fluctuations and on this procedure to set up the initial conditions, we have explored the dependence on the initial conditions using virialized clouds of 0.75 and 0.084 pc radius (see parameters of run series A and C in Table A.1).

Fig. A.1 Mass spectra of runs A0, A1d, A2, and A2d, with the same color coding as in Fig. 6. The cloud with a radius of 0.75 pc is evolved without gravity on levels 8 or 9 during 0, 0.3, and 0.6 turbulence-crossing times. With longer relaxation, more massive stars develop, and the IMF is slightly more top-heavy.

Fig. A.2 Mass spectra of runs C0, C1, and C2 with the same color coding as in Fig. 6. The cloud with a radius of 0.084 pc evolved directly without relaxation or evolved without gravity on levels 8 and 9 during 0.3 turbulence-crossing time. There is no significant dependence of the IMF on the initial density fluctuations.

We performed simulations after three different relaxation times, 0, 0.3 or 0.6 t_tc . To verify that the initial fluctuations are not determinant, we also varied the initial resolutions by using a configuration with either 256³ (level 8) or 512³ (level 9). The relaxation was made on the uniform base grid, and the AMR refinement was activated only for the second phase and was allowed down to level 14. The resulting mass spectra for runs A are shown in Fig. A.1 and runs C in Fig. A.2. There are no noticeable trends for either runs A or Cwith varied resolution of relaxation, except for runs A1d and A2d, which are relaxed during 0.6 t_tc and have probably lost too much kinetic energy and thus yield fewer low-mass stars.

Appendix B Gaussian fit to the IMF

The IMFs in Figs. 6 and 7 are fitted to a Gaussian distribution. The purpose of these fits is to provide a systematic way to estimate the peak position and width of the distribution. We note that while reasonable, these fits are not always very good. In particular, the fitted distribution tends to underestimate the number of objects around the peak. In Table B.1, the peak mass and the width of the Gaussian fit are listed.

References

All Tables

All Figures

Fig. 1 Column density maps of runs A1++ at 0.155, 0.203, 0.424 Myr (top row from left to right) and C1+ at 0.003, 0.004, 0.008 Myr (bottom row from left to right). Run A1++ is more globally collapsing, and the sinks, initially forming inside the filamentary network, are more concentrated in a central cluster. Run C1+ is initially much denser and the evolution of the filamentary structure is more pronounced, with more widespread sink formation.

In the text

Fig. 2 Zoomed column density maps of runs A1++ (top row) and C1+ (bottom row). Red dots represent the sink particles, and their sizes correspond to the sink mass. Disk formation around sink particles is observed in both cases, while the disks in run A1++ are more extended, showing likely signs of rotation and fragmentation, and the disks in run C1+ are significantly smaller, with no sign of fragmentation.

In the text

Fig. 3 Density PDF (normalized in random units) for runs A1++, B1++, C1+, D1, C1t03, and C1t01. At low densities, the shape of the PDF is broadly lognormal, while at high density, it is ∝ n−1.5 . The PDF is almost invariant in time given that the cloud is sufficiently supported, as it is in most cases. In run C1t01, the initial turbulent support is very low and thus the cloud global collapse is strong, giving a density PDF that shifts toward high densities in time.

In the text

Fig. 4 Density profiles around sink particles for runs A1++ at 0.26 Myr (top panels) and C1+ at 0.008 Myr (bottom panels). From left to right, sinks aged under 1000 yr, 5000 yr, or all sinks are used to produce the average profiles for run A1++, while 300 yr and 1000 yr are used as thresholds for run C1+. Sink particles are separated into three mass bins <0.01 M⊙ (blue), 0.01–0.1 M⊙ (green), and 0.1–1 M⊙ (red), average age of selected sinks is shown in the legend. Shaded regions represent the standard deviation of the density profiles for all sinks and all radial directions. The SIS density profile, , is plotted in black for reference. The density profile at younger ages is always broadly ∝ r−2 , with a factor~10 higher than the SIS, while sinks of older ages have probably accreted all mass from the core in which it is embedded or have been ejected, thus decorrelating from the r−2 profile. Thisis reflected by the large dispersion that is particularly evident with more massive sinks.

In the text

Fig. 5 Final sink mass against accretion time (tacc,60 is used) for runs A1++, B1++, C1+, and D1. The solid black lines represent the free-fall time of the mass reservoir derived from the analytical model (see Sect. 5.5.1). The behavior of the distributionand the analytical model are very similar except for short accretion times (corresponding to mass ≲0.1 M⊙), where the physics to be considered is different. This strongly suggests that the mass that built the stars essentially comes from a preexisting deterministic reservoir.

In the text

Fig. 6 Sinkmass spectra for models A1, B1, C1, and D1 at various time steps, corresponding to a comparable amount of accreted mass. Various spatial resolutions are explored for the purpose of verifying numerical convergence and direct comparison between models. For models B1, C1, and D1, a peak at ≃ 0.1 M⊙ is obtained independently of the numerical resolution. While model A1 presents a flat mass spectrum, models C1 and D1 present at high mass a power-law behavior that is compatible with M−3∕4 (see also Fig. 8).

In the text

	Fig. 7 Mass spectrum evolution for models C1t15, C1t03, and C1t01 corresponding to three levels of initial turbulence. Except for C1t01, which initially has a very low turbulent energy, the mass spectra are not significantly different.
In the text

Fig. 8 Mass spectra as predicted from analytical model for four sets of parameters corresponding to runs A1++, B1++, C1+, and D1. Simulation results at 150 M⊙ accreted (blue dashed, same as the magenta histogram in Fig. 6 with absolute value shifted to compare to models), model with power-law (blue solid) and lognormal (red dashed) density PDF are plotted, as well as the power-law relation (–3/4 in cyan and –1 in black) for reference. The vertical dot-dashed line indicates a mass lower limit above whichthe model is applicable (see paper II).

In the text

	Fig. A.1 Mass spectra of runs A0, A1d, A2, and A2d, with the same color coding as in Fig. 6. The cloud with a radius of 0.75 pc is evolved without gravity on levels 8 or 9 during 0, 0.3, and 0.6 turbulence-crossing times. With longer relaxation, more massive stars develop, and the IMF is slightly more top-heavy.
In the text

	Fig. A.2 Mass spectra of runs C0, C1, and C2 with the same color coding as in Fig. 6. The cloud with a radius of 0.084 pc evolved directly without relaxation or evolved without gravity on levels 8 and 9 during 0.3 turbulence-crossing time. There is no significant dependence of the IMF on the initial density fluctuations.
In the text