Note to: A&A, 611, A88 (2018), https://doi.org/10.1051/0004-6361/201731522

© ESO 2018

1 Introduction

Star formation is believed to have major consequences on the structure of our Universe. In spite of sustained efforts, the origin of the initial mass function (IMF) in particular is still a matter of debate (e.g., Offner et al. 2014). While the gravo-turbulent fragmentation of cold dense molecular gas is a widely accepted scenario of star formation, it is in apparent contradiction with the observed IMF exhibiting no significant variation among different environments (e.g., Kroupa 2001; Chabrier 2003; Bastian et al. 2010; Offner et al. 2014), although Dib et al. (2017) found more variations than are commonly assumed.

In particular, the physical origin of the peak of the IMF, around 0.3 M_⊙, is still uncertain. Several ideas have been proposed, which fall into three categories that are not mutually exclusive. First, it has been proposed that the thermodynamics of the gas could provide a particular density to which a Jeans mass is associated. For example, Larson (1985) computed the gas temperature in dense core conditions and concluded that owing to molecular cooling, the effective adiabatic exponent, γ_ad, is expected to be on the order of 0.7 for densities below ≃10⁵ cm⁻³, while at higher densities, γ_ad ≃ 1 as a result of gas-dust coupling. Since the fragmentation sensitively depends on γ_ad, the question arises as to whether the change of the effective equation of state at ≃10⁵ cm⁻³ could lead to a particular Jeans mass, which for a temperature of 10 K is ~0.3 M_⊙. A similar, though not identical, line of explanations has been put forward by Elmegreen et al. (2008): the thermodynamics of the gas in various environments could result in a Jeans mass that weakly depends on the gas density. The second category of explanations has been proposed by Bate (2009b), Krumholz (2011), and Guszejnov et al. (2016). The idea is that because of the heating from the accretion luminosity, the gas temperature increases in the vicinity of the protostar, and therefore the local increase of the Jeans mass leads to a possible self-regulation. Again, these studies infer that the Jeans mass has a weaker dependence on the gas density than in an isothermal gas. Finally, the third category of explanations invokes a compensation between the density and the velocity dispersion, which vary in opposite directions with clump masses through Larson relations (Hennebelle 2012; Lee & Hennebelle 2016b). Considering clumps of increasingly higher mass, the density decreases with the mass, which in turn increases the Jeans mass, while the increasing Mach number tends to create denser gas through shocks, yielding a slowly varying characteristic mass. None of these explanations is well established yet. For example, while it has been suggested by various simulations (e.g., Jappsen et al. 2005; Bate 2009b) that thermodynamics influences the mass spectrum of the objects that form, there is no well-established set of simulations demonstrating the invariance of the peak of the IMF for a broad range of physical conditions.

In the companion paper (Lee & Hennebelle 2018, hereafter Paper I), we conduct a systematic exploration of the initial conditions by varying the initial density and turbulence of the molecular cloud. The physics is greatly simplified at this stage because our goal is to identify the physical processes at play. At relatively high densities (≳ 10⁵   cm⁻³), the shape of thestellar distribution, including the peak position and the high-mass end slope, becomes no longer dependent on the global density of the parent molecular cloud. This means that some physical mechanisms are operating at local scales and do not depend on the large ones, such as the collapsing clump itself. This result is surprising at first because the mean Jeans mass changes by more than one order of magnitude in our series of simulations. Moreover, the thermodynamics treated in these simulations assumes that the gas remains isothermal up to a density of a few 10¹⁰   cm⁻³, where the gas becomes adiabatic. While in the context of a collapsing clump, the isothermal regime does not provide any characteristic mass, the adiabatic regime does, because the Jeans mass increases with increasing density. The lowest Jeans mass, estimated at the density where the gas becomes adiabatic, is a few 10⁻³   M_⊙ (Rees 1976; Whitworth et al. 2007) and therefore is almost two orders of magnitude lower than the peak position of the stellar distribution around 0.1 M_⊙ observed in the simulations. Consequently, the origin of a characteristic mass in our simulations is not obvious and needs to be elucidated. For this purpose, we perform in this paper a series of simulations with initial conditions identical to those in Paper I, while altering the gas equation of state (eos). With the variation of the stellar distribution peak position (around 0.1 M_⊙) with respect to the varying eos, we try to shed light on the mechanisms that determine the IMF peak mass.

In the second section, we present the numerical setup we employed, including the initial conditions and the various eos used throughout the paper. In the third section, we present the results of the numerical simulations. In the fourth section, we infer the mass of the first Larson cores for various eos, and we study its correlation with the mass spectrum peak in the simulations. In the fifth section, we develop an analytical model that accounts for the factor, on the order of ten, between the mass of the first Larson core and the peak of the stellar distribution obtained from numerical simulations. In the sixth section, we provide a discussion of our results. The seventh section concludes the paper.

2 Numerical setup and initial conditions

2.1 Numerical setup and runs

We used the adaptive mesh refinement (AMR) magnetohydrodynamics code RAMSES (Teyssier 2002; Fromang et al. 2006) to evolve the hydrodynamical equations. Simulations were initialized with a Bonner-Ebert-like spherical molecular cloud with a density profile , where r₀ is the size of central plateau and r the distance to the cloud center, as described in Paper I (see more details therein).

Run C1in Paper I is used as the canonical run here and is labeled C10h, as we present below. The cloud has 1000 M_⊙ and 0.084 pc radius, with the simulation box twice the size of the cloud. The central density n₀ = ρ₀∕(μm_p) ≃ 6 × 10⁷ cm⁻³, where μ = 2.33 is the mean molecular weight and m_p is the atomic hydrogen mass. We note that this density is very high and may not be realistic to describe most Milky Way clouds. However, our goal here is to investigate the physical processes rather than performing a detailed comparison with observations. A turbulent velocity field following a Kolmogorov spectrum with random phases was seeded, and a relaxation phase without self-gravity was run during ~ 30% of the turbulence-crossing time on the base grid (2⁸) to prepare for coherent density and velocity fluctuations. Initially, the turbulence was a mix of compressive and solenoidal modes with an energy ratio 1:2. The AMR scheme requires that the local Jeans length is always resolved by ten cells. The canonical run had a maximum refinement of 14 levels, corresponding to 4 AU resolution, while higher and lower resolution runs were also performed. The simulation parameters are listed in Table 1.

2.2 Equation of state

In a dense core, the temperature is about 10 K, except at high density, where the gas is optically opaque to its own radiation and is nearly adiabatic (e.g., Masunaga et al. 1998; Vaytet et al. 2012; Vaytet & Haugbølle 2017). The temperature T, which depends on the gas density ρ, is described with the eos (1)

where the isothermal temperature T₀ at low density is set to 10 K. We studied the effect of the eos by varying the critical turnover density ρ_ad, which designates the change between isothermal behavior and adiabatic heating, and the polytropic index γ, which describes the heating rate. The polytropic index of molecular hydrogen is γ₁ =5/3 and γ₂ =7/5 before and after the excitation of rotation modes of hydrogen molecules at around 100 K. At ~ 1000 K, the H₂ molecules start to dissociate and γ drops to γ₃ ~ 1.1. In this study, we aim to examine the effect of thermodynamics on the final mass a star can reach. We therefore also considered a more complete description of a full eos, (2)

where n = 3 is a smoothing parameter. This barotropic eos is employed for convenience. It signifies effectively T ∝ ρ^γ−1 with γ = 1 for ρ < ρ_ad, γ = γ₁ for ρ_ad < ρ < ρ_ad,2, γ = γ₂ for ρ_ad,2 < ρ < ρ_ad,3, and γ = γ₃ for ρ > ρ_ad,3. In the following studies, we fixed ρ_ad,2 = 30 ρ_ad and ρ_ad,3 = 1000 ρ_ad,2.

As shown in Table 1, n_ad = ρ_ad∕(μm_p) is varied between 10⁹ cm⁻³ and 10¹⁰ cm⁻³, while using a hard eos γ = 5∕3 or a soft one γ = 4∕3. We also study the effects of a full eos that describes the transition of γ from 1 to 5/3, 7/5, and then 1.1. The two latter transitions correspond roughly to T = 100 and 1000 K. We also set n_ad = 10¹³ cm⁻³ in run I, that is practically isothermal since such high density is hardly reached at this resolution.

2.3 Sink particle algorithm

We used the sink particle algorithm from Bleuler & Teyssier (2014) (see also Krumholz et al. 2004; Federrath et al. 2010), and we briefly describe it here. Sink particles are formed at the highest refinement level. The algorithm first identifies local density concentrations, which are referred to as clumps, above a given density threshold (we typically used 1/10 of the density threshold for sink particles n_sink). Connected clumps are merged if the density contrast between the peak and the saddle point is less then two. The final clumps are then passed on for a sink formation check.

The clump peak undergoes several criteria before a sink formation site is flagged. First, a density threshold is imposed such that the peak must have a density higher than n_sink, which we varied between 10¹⁰ and 3 × 10¹¹ to investigate its influence. Second, additional criteria are used to check whether the clump is virially bound and has a converging flow. A sink particle is placed at the density peak when all criteria are satisfied. Simple density thresholding was used in runs C1hd–o to study the numerical effects of the algorithm. The sink then interacts gravitationally with the gas component as well as other sinks and accretes from the surrounding with a threshold scheme. When the cells surrounding the sink particle exceeds n_sink, 75 percent (this numerical factor seems not to have a crucial effect on our results, see Appendix A) of the excessive gas mass in the cells is accreted onto the sink particle.

2.4 Missing physics

The physics was deliberately simplified because we tried to conduct a systematic set of simulations to clarify the influence of the initial conditions, eos, and resolution. Other processes not included here nonetheless play significant, possibly dominant roles. This is particularly the case for the accretion luminosity that emanates from the protostars (Krumholz et al. 2007b; Bate 2009b; Commerçon et al. 2011) and heats the gas, which in turn increases the thermal support. The magnetic field is another important process (Hennebelle et al. 2011; Peters et al. 2011; Myers et al. 2013) that likely affects the fragmentation of clusters through magnetic braking and magnetic support; we do not consider this here.

3 Results: mass spectra in numerical simulations

3.1 Isothermal runs

The sink distributions displayed in Fig. 1 present a well-defined peak and a power law at high masses that is roughly ∝ M⁻¹. The isothermal runs show that the stellar distribution peak moves from 0.05 to 0.01 M_⊙ when we change the resolution from 17 to 4 AU. As the collapse of isothermal gas is self-similar, there is no characteristic scale or mass in such simulations at high density. Roughly speaking, when we double the resolution, the highest density that can be described increases by a factor of about 2³ = 8. We therefore may expect that the peak position shifts by when two more refinement levels are introduced. This is compatible with the trend observed in Fig. 1.

The mass of sink particles therefore depends on the numerical resolution. This is likely what occurs in simulations that recovered the stellar distribution in a collapsing massive and isothermal cloud (e.g., Bonnell et al. 2006; Girichidis et al. 2011; Ballesteros-Paredes et al. 2015). Typically, the dependence of the peak position is expected once the collapsing cloud is sufficiently cold. It is worth emphasizing that this is not a statement about the absence of convergence of isothermal calculations in general. As shown in Girichidis et al. (2011) and in Paper I, initial conditions matter greatly. For example, Gong & Ostriker (2015) did obtain numerical convergence in their colliding flow calculations.

Fig. 1 Stellar distribution of isothermal runs with varied resolutions of 17, 8, and 4 AU in colored histograms, with time and total accreted mass shown in the legend. Lognormal fits are presented with dashed thin curves. Power laws of M−3∕4 and M−1 are shown for comparison with black and gray lines, respectively. The characteristic mass of the stellar distribution shifts to lower values with increasing resolution.

3.2 Runs with γ = 5∕3

With γ = 5∕3, the fiducial run (C10h) in Fig. 2 very quickly shows a well-defined stellar distribution that peaks at ≃ 0.1 M_⊙. The runs with four resolutions, 2, 4, 8, and 17 AU, show that the peak position is robust and does not change significantly with the resolution. This is a noticeable difference to the isothermal case. There is nevertheless some significant evolution of the mass spectra when the spatial resolution improves. First of all, the number of objects near the peak increases when the resolution improves from 17 to 4 AU (from log₁₀(dN∕dlogM) ≃ 1.6 to ≃ 2.1 at the time when 200 M_⊙ have been accreted), while this number does not seem to change significantly between the two highest resolutions (two bottom panels). As in the isothermal simulations, the high-mass part is also described by a power law ∝ M⁻¹, although this is much clearer in the highest resolution runs.

Figure 3 reveals that n_ad has an effect on the peak position. With n_ad = 10⁹, sinks are completely prevented from forming at the resolution of 4 AU (C1h, no mass spectra available). Decreasing the resolution allows sink particles to form. The run at 8 AU resolution (C1h–) produces top-heavy spectra with a loosely defined peak at ≃ 1 M_⊙, while a coarser resolution at 17 AU (C1h– –) shifts the peak to ≃0.3 M_⊙ and more sinks are formed. This effect results from the excessive heating at high densities, where the approximated expression of one single γ value stops being physical, and the high-temperature gas is thus always thermally supported against collapse at the finest resolution scale and sink particles hardly form. This behavior is interpreted more quantitatively in Sect. 4. The run with n_ad = 3 × 10⁹ at 8 AU resolution (C3h–) produces mass spectra that peak at lower mass than that of run C1h– at same resolution.

These results suggest that the eos has a strong effect on the characteristic mass of the stellar distribution. Moreover, we find that the stellar distribution peak is significantly higher than the Jeans mass at which the gas becomes opaque to its radiation (typically on the order of 10⁻³ M_⊙) or even the mass of the first Larson core (a few 10⁻² M_⊙). The link between the two is addressed in Sect. 4 and Sect. 5.

3.3 Runs with γ = 4∕3

As revealed by Fig. 4 for runs using γ = 4∕3, the peak shifts to slightly lower masses and the stellar distribution is slightly wider than with γ = 5∕3. The run with n_ad = 10⁹ cm⁻³ (C1s–) yields a slightly higher peak mass than that with n_ad = 10¹⁰ cm⁻³ (C10s–)since the heating starts at lower density, and thus the higher temperature at the same density provides more support against gravity. With γ = 4∕3, the gas heats up less quickly as density increases than at γ = 5∕3, and the high-density region suffers less from overheating (more sinks form in C1s– than C1h–).

3.4 Runs with a full eos

The first panel of Fig. 5 shows a very similar behavior to the fiducial run (C10h) and suggests that the γ = 5∕3 part of the eos matters the most. To test the influence of the sink particle scheme, we have performed a run for which n_sink = 3 × 10¹¹ cm⁻³ (C10fd, second panel). Applying higher threshold for sinks does not affect the peak mass, while the low-mass end becomes less populated. This shows that the exact moment where the sink particles are introduced does not control the peak position, but has some influence on the distribution of low masses. This is compatible with the eos being the most important process that determines the peak position.

Finally, we also considered n_ad = 10⁹ cm⁻³ using the full eos. The shift of the peak is as obvious with a full eos (C10f and C1f) as that with γ = 5∕3. We note that in run C1f we overcome the overheating problem at high resolution that we encountered in run C1h. Gas is allowed to collapse andform sink particles without difficulty. This confirms that the determining part of the eos corresponds to γ = 5∕3.

Fig. 2 Stellar distribution of runs with a smoothed two-slope polytropic eos with γ = 5∕3 and nad = 1010 cm−3, with varying resolutions of 17, 8, 4, and 2 AU. The characteristic mass of the stellar distribution is ≃ 0.1 M⊙, regardless of numerical resolution.

4 Mass of the first Larson core and peak of the sink distribution

As the previous section suggests that the thermodynamics is essential in determining the peak of the stellar distribution, we calculate the mass of the first Larson core and study its correlation with the peak mass obtained in the simulations.

4.1 First Larson core

As gravitational collapse proceeds, the density increases, and at some point, the dust becomes opaque to its own radiation. At this point, the gas becomes essentially adiabatic, and a so-called first Larson core forms (Larson 1969; Masunaga et al. 1998). While the gravitational energy can overcome the thermal support through contraction of an isothermal gas, gas with polytropic index γ > 4∕3 can resist the collapse. The first Larson core grows in mass by accretion, leading to the increase of the central density and temperature. When the temperature reaches about 1500 K, molecular hydrogen starts dissociating and therefore efficiently cools the gas, and the gravitational collapse resumes, which eventually leads to the formation of the second Larson core, namely the protostar. To obtain a protostar, it is probably necessary to accumulate sufficient mass to trigger the second collapse, and this process is essential in setting the peak of the stellar distribution. First Larson cores with lower masses cannot form protostars (or will do so through cooling, but on a longer timescale). While this process is not explicitly resolved in our simulations, the introduction of the sink particles to some extent mimics this sequence since they are introduced only when enough mass has been accumulated for the gas to be gravitationally unstable.

By integrating the hydrostatic equations of the gas, we deduce the density profile that a core develops before undergoing the secondary collapse. From the density profile, the mass can be inferred, provided that some radius can be estimated or specified. The equations of hydrostatic equilibrium are

where M, ρ, P, r, and G are the mass contained within a given radius, density, pressure, radius, and the gravitational constant. The derivative of Eq. (4) with respect to r yields (5)

where the second equality is obtained by substituting with Eqs. (3) and (4). The pressure is related to the density by the eos that we used in the simulations, (6)

where k_B is the Boltzmann constant, and the temperature T(ρ) is given by Eqs. (1) or (2).

This results in a second-order differential equation of ρ that we can solvenumerically using standard Runge-Kutta methods: (7)

Before the second collapse occurs, there is no singularity at the center, and thus we use the boundary conditions (8)

where ρ₀, the central density, remains a free parameter to be specified.

Figure 6 shows some examples of the density profile (upper panel) and integrated core mass (lower panel). The cases γ = 5∕3 (blue), γ = 4∕3 (cyan), and full eos (red) are plotted with n_ad = 10¹⁰ cm⁻³ (solid) and 10⁹ cm⁻³ (dashed), where the legends correspond to the labels of the simulations. The black curves display the isothermal case. We also show the profile calculated with the eos used by Bate et al. (2003) with γ = 7∕5 and n_ad = 2.45 × 10¹⁰ cm⁻³. The central density used as boundary condition for the integration was chosen to be n₀ = 100 n_ad, which is a reasonable value for illustrative purpose, while the actual value of n₀ that is reached in simulations depends on the resolution and sink formation algorithms.

In general, the result of the hydrostatic equilibrium without singularity is a central density plateau with a decreasing envelope. Depending on the central density, the solution may decrease infinitely or reach negative values at large radius, while it should be connected to a confining ambient pressure in physical conditions. For sake of convenience, we defined the boundary of this object where the density reaches the ambient density, namely 10⁷   cm⁻³ in this study. With a sharply decreasing density profile, this selection of the density cutoff does not have a strong effect on the derived mass. As long as the outer density profile has a radial dependence that is much steeper than r⁻², the object is clearly distinguishable from the isothermal envelope and its mass is well defined.

Figure 7 displays the size of the core truncated at ambient density 10⁷   cm⁻³ (top), the local power-law exponent of the density profile (middle) α = −rdn∕dr∕n, such that locally n ∝ r^−α, with the gray linetracing α = 2, and the integrated mass (bottom), as a function of the central density, n₀, with several n_ad and γ values corresponding to the simulations. Above a certain value of n₀ (when the adiabatic regime begins), α starts to be significantly larger than 2 and a core is well defined, providing a clear definition of the first Larson core mass, M_L. Exceeding this mass, the gas becomes gravitationally bound, and this corresponds to the formation of sink particles in simulations. With low values of n₀, on the otherhand, the integrated mass is sensitive to the ambient density cutoff and is almost constant regardless of n₀. This corresponds to the isothermal regime of the eos, and the density profile approaches that of a singular isothermal sphere (SIS; Shu 1977). The mass depends on the density cutoff, and no obvious core could be identified. The size is therefore physically meaningful only at high n₀, while at low n₀, it is merely a truncation of the SIS. In practice, we used α ≫ 2 to infer the radius of the Larson core. As can be seen, the size of the first Larson core is typically a few tens of AU, well resolved in our simulations. When a non-isothermal eos is used, the density is almost flat at the center and significantly drops at some radius. The higher the value of n₀, the smaller the size of the core, such that the integrated mass does not increase strongly with increasing n₀, and even reaches a characteristic value, except for the cases with γ = 5∕3, where the mass increases with n₀ at rate of roughly .

An interesting question is whether rotation could modify these conclusions. In particular, since rotation adds another support, it might lead to more massive first Larson cores. To investigate this possibility, we used the ESTER code (Espinosa Lara & Rieutord 2013; Rieutord & Espinosa 2013), which computes bidimentional equilibria of rotating polytropes in solid-body rotation. We found that the maximum mass of a rotating polytrope increases only slightly (tens of percents). The reason is that the rotation support is significant only in the outer part, while most of the mass lies in the central region.

Fig. 3 Stellar distributions of runs with a smoothed two-slope polytropic eos with γ = 5∕3, lowered nad at 3 × 109 (C3h–) and 109 cm−3 (C1h–, C1h– –). The characteristic mass of the stellar distribution increases with decreasing nad. At nad = 109 cm−3, the stellar distribution becomes top-heavy with improved resolution from 17 to 8 AU.

Fig. 4 Stellar distributions of runs with a smoothed two-slope polytropic eos with γ = 4∕3, nad = 1010 cm−3 (C10s–) and 109 cm−3 (C1s–). The characteristic mass of the stellar distribution increases by a factor ~ 3 when nad is lowered by a factor ten.

Fig. 5 Stellar distributions of runs with a full eos with nad = 1010 cm−3 (C10f, C10fd) and 109 cm−3 (C1f). The transitions to γ = 5∕3, 7∕5, and 1.1occur at nad, 30 nad, and 30 000 nad. Using a higher sink density threshold (C10fd) results in fewer low-mass stars but does not change the peak of the stellar distribution. Decreasing nad by a factor ten increases the peak mass by ~10.

Fig. 6 First Larson cores for various equations of state. The hydrostatic equilibrium equations are integrated with n0 = 100 nad as boundary condition. We present cases with γ = 5∕3 (blue), γ = 4∕3 (cyan), full eos (red), and isothermal (black). The legends correspond to the labels of the simulations, with solid lines for nad = 1010 cm−3 and dashed lines for 109 cm−3. We also show the case γ = 7∕5 with nad = 2.45 × 1010 cm−3 by Bate et al. (2003, green). The values of n0 are chosenfor illustrative purpose, and the actual values reached in simulations depend on the resolution (highest density resolved, see Fig. 7). The circle indicates the radius at which the density reaches 107   cm −3, a reference ambient density that truncates the core. As long as the density profile decreases steeply enough, the choice of this truncation density has no strong effect on the resulting mass.

Fig. 7 Properties of the first Larson core plotted against the central density n0. The hydrostatic equilibrium Eq. (7) is solved with various nad and γ values, and truncated at the ambient density n = 107 cm−3. Three quantities are shown: the radius at which the density decreases to n = 107   cm−3, the local power law exponent α = −rdn∕dr∕n, and the integrated mass. Top left: γ = 5∕3; top right: γ = 4∕3; bottom left: full eos. Bottom right: eos taken from Bate et al. (2003, solid), Bate (2005, dashed), and Bonnell et al. (2011, dotted). Each time the isothermal case is plotted in black for comparison. A core is well defined only when α ≫ 2, and thus giving a characteristic mass of the first Larson core. When a core is well defined, its mass is almost constant regardless of n0 when γ = 4∕3 or full eos is used. In contrast, the core mass increases with increasing n0 when γ = 5∕3, which explains the peak mass shift between runs C1h– and C1h– –.

4.2 Interpreting the numerical simulation results

When we assume, as we show below, that the peak of the stellar distribution is simply proportional to the mass of the first Larson core, M_L, these calculations can be used to interpret the results obtained from the numerical simulations.

For γ = 5∕3 and n_ad = 10¹⁰ cm⁻³, Fig. 2 shows that the peak of the stellar distribution does not change significantly with the resolution. Since increasing the resolution allows us to probe higher densities, we can interpret this behavior as a consequence of the weak dependence of M_L on the core central density, n₀, which, as Fig. 7 shows (blue solid lines) remains below 2 × 10⁻² M_⊙ for n₀ < 10¹², while noting that a self-gravitating core is well defined only when n₀ ≳ 2 × 10¹⁰ (α ≫ 2). With even higher resolution, we can expect that the peak of the stellar mass spectrum shifts to higher values, as seen with the C1h runs at varied resolutions.

For n_ad = 10⁹ cm⁻³ and γ = 5∕3 (Fig. 3), the self-gravitating first Larson core is always well defined, with a mass that increases with n₀. This is reflected by the shift in peak mass between runs C1h– – and C1h–. The eos with n_ad = 10⁹ cm⁻³ gives a mass of the first Larson core that is higher than that with n_ad = 10¹⁰ cm⁻³ at the same value of n₀, which is also in agreement with simulation results.

With γ = 4∕3, the effect of varied n_ad on the peakof the stellar distribution was found to be less significant (Fig. 4). This is again in good agreement with Fig. 7 (cyan lines), where at the same value of n₀, the mass varies less with n_ad than in the γ = 5∕3 cases.

Finally, the runs with full eos (Fig. 5) present a similar behavior to the runs with γ = 5∕3 (Fig. 3) regarding the relation between peak mass and n_ad. This is also expected as the mass dependence of M_L is relatively similar to that of γ = 5∕3 for n₀ below 10¹²   cm⁻³ (Fig. 7, red lines).

4.3 Correlation between the first Larson core and the stellar distribution peak

A link between the mass of the first Larson core, M_L, and the mass at which the stellar distribution peaks has been suggested in previous discussions. We now compare this in a more systematic way. For this purpose, we have compiled our simulations with other results from the literature, and we present in Fig. 8 the correlation between the peak mass of the stellar distribution and the mass of the first Larson core deduced from the eos we employed in the corresponding runs. The mass of the first Larson core is deduced from Fig. 7 to lie around the point where α starts to be larger than two. Since the definition is not entirely straightforward and also because M_L depends on n₀, we estimated error bars by considering that n₀ may vary by a factor 10. To estimate the peak position, we took the mass corresponding to the maximum values in the mass spectra. The error bars are given by the binning size.

We also included results from literature. Bate et al. (2003) simulated a 50 M_⊙ cloud with γ = 7∕5 above 10⁻¹³ g cm⁻³ and obtained a median sink mass of 0.7 M_⊙. Increasing the density by decreasing the size (Bate & Bonnell 2005) or increasing the mass and velocity dispersion while keeping the same density (Bate 2009a) gave similar median masses of 0.023 and 0.02 M_⊙ that are significantly lower than in the previous case. Alternatively, changing the critical density to 1.1 × 10⁻¹⁴ g cm⁻³ (Bate 2005) gave a median mass of 0.054 M_⊙. The median mass was used from these studies since there is not enough statistics and a peak was not always well defined. The eos used in the simulation series by Bate et al. is a piece-wise temperature. To calculate the density profile of a mass concentration, the temperature as well as its derivative are required to be continuous, therefor we used (9)

where n = 100 ≫ 1, to replace their eos when we inferred the mass of the first Larson core.

Figure 8 shows that the first Larson core mass is in good correlation with the stellar distribution peak, and we approximately obtain M_peak ≃ 10M_L. We offer an explanation for this additional factor of about ten in Sect. 5.

5 Stabilizing effects of tidal forces around the first Larson core: a model

To link the mass of the first Larson core to the mass of the forming star, we present a simple model to discuss the effect of a point mass in a density field. The questions we address are how the typical mass of a self-gravitating object embedded in a dense collapsing cloud is limited in general to a small fraction of that of this cloud, and how this fraction depends on the mass of the accreting object. Here we propose that the accretion onto the already formed self-gravitating object competes with the tendency to form newobjects from gravo-turbulent fragmentation. We therefore wish to estimate the typical radius at which the probability to form another self-gravitating fragment equals one. The mass enclosed within this radius is expected to be accreted by the central object, while the mass outside is expected to be distributed between forming fragments.

We discuss the conditions under which a density perturbation located near the point mass can become self-gravitating and form another core. The idea is that the first Larson core will “screen” itself by shearing out fluctuations through tidal forces. The modeling essentially proceeds in four steps:

• In Sect. 5.1 we describe a point mass that is surrounded by a density profile, with a perturbation at distance r_p and of size δr.

• In Sect. 5.2 we develop the mass condition. The perturbation needs to have a mass of at least M_L such that its collapse can result in the formation of a stellar object.

• In Sect. 5.3 we infer the energy condition. The gravitational field generated by the point mass, the density field, and the perturbation itself are used to calculate the gravitational energy of the perturbation. The idea is to compare the self-gravity of the perturbation to the destructive tidal field created by the point mass and the envelope. The gravitational energy of the perturbation is required to be negative such that it can contract to form another star.

• In Sect. 5.4 we estimate the mass that is expected do eventually be accreted onto the central object. The mass and energy constraints provide a condition on the density perturbations as a function of distance. Assuming a lognormal distribution for the density fluctuations (in addition to the local mean density, which itself is ∝ r⁻²), we can infer at which radius the probability of finding a self-gravitating perturbation of mass at least equal to M_L is one.

Fig. 8 Peak mass of the sink mass spectrum plotted against the mass of the first Larson core calculated with the eos we used in the simulations. The peak mass is read from the stellar distribution, with the bin size as error. The error bars from literature results are larger as they have lower statistics and higher fluctuations in the spectrum. The mass of the first Larson core is read from Fig. 7, with error bars indicating the mass range with the probable n0 range. Whenwe examine this by eye, there is a correlation Mpeak ~ 5MLarson. We also performed a simple linear regression, without considering the errors, to determine the scaling factor, and this led to .

5.1 Step 1: density profile around sink particles and perturbation

By extracting gas profiles around young sinks particles that have not yet decoupled from their natal gas, we have shown in Paper I that the cores that form stars have indeed some characteristic density profiles ρ ∝ r⁻², as expected from theory (Shu 1977). We therefore take the density field ρ_e (r_e) as a function of distance to the point mass (10)

where A is the amplitude of the SIS density profile (Shu 1977). In Paper I, the density profiles around young sink particles were inferred from simulations, and a typical value of A ~ 10 was found at early time when the sinks are still actively accreting (see the left panels of Fig. 4 of Paper I). For the sake of simplicity, we consider a uniform density spherical perturbation ρ_p at distance r_p: (11)

where r_e and r_p are the vectors pointing from the central object to a point in the envelope and to the center of perturbation, and η is a constant for the local density contrast of the perturbation.

To determine whether the perturbation is prone to self-gravitating collapse, we use the virial theorem to probe its stability (12)

The virial energy E_vir is a combination of the gravitational energy E_g and the thermal energy E_ther, with M_p being the mass of the perturbation. The gravitational acceleration g = g_L + g_e + g_p has three components, coming from the central first Larson core that is described with a point mass, the power-law density envelope, and the perturbation itself. We note that to correctly calculate the virial energy, only the tidal tensor part of the gravity should be considered and the mean value exerted on the center of the object should be removed. In the following, we therefore always use the relative gravitational acceleration with respect to the center of the perturbation. The collapse criteria are therefore

5.2 Step 2: the mass condition

The mass contained inside the perturbation radius δr_p around the center of perturbation is (15)

We designate θ the angle between δr and −r_p. With the trigonometric relation (16)

we can integrate (17)

where the normalized quantity u_p = δr_p∕r_p. The normalized mass is expressed as m_p(u_p, η) and is given in Appendix B.

As explained above, the mass of the perturbation must typically be equal to M_L since no star can form below this value. By requiring M_p(r_p, δr_p, η) = M_L, we obtain a relation between the size u_p and amplitude η of a perturbation at distance r_p.

5.3 Step 3: the energy condition

Now turning to the calculation of the virial parameter, we first show the importance of tidal forces in the vicinity of the sink particles, and then, to gain physical intuition, we perform a simplified 1D estimate before carrying out a full 3D calculation.

5.3.1 Tidal forces

Tidal forces have been advocated in various contexts to either limit or quench star formation (Bonnell & Rice 2008; Ballesteros-Paredes et al. 2009), but also, conversely, to trigger it (Jog 2013; Renaud et al. 2014). Ntormousi & Hennebelle (2015) discussed the influence of the density distribution and showed that in a collapsing cloud, the tidal forces tend to be initially compressive, which promotes fragmentation, and then the forces become stabilizing. The key quantity to study the influence of tidal forces is the gravitational stress tensor, ∂_i g_j, which is characterized by the three eigenvalues, λ_i. From the Poisson equation, we know that ∑ λ_i = −4πGρ. Thus there is at least one negative eigenvalue, say λ₁, and in most cases, two of them are negative. It is well known that certain gravitational fields present a value of λ₃ that is positive and that these fields render the formation of self-gravitating perturbations less likely (see Appendix C for the tidal forces in the simulation). Below, the tidal effect is taken into account by including the corresponding term in the virial theorem.

5.3.2 1D calculation of gravitational energy

The direction along r_p is the relevant dimension as the gravity of the central mass, and the envelope generates diverging tidal forces and prevents theoverdense perturbation from self-gravitating collapse. Based on the reasonable illustrative results and for simplicity, we first perform a virial integration along this direction on the axis of the perturbation, ignoring thermal energy, and try to obtain some physical insight from the simplified calculations. We define the position δr with respect to the center of perturbation at r_p, with the positive direction pointing away from the central mass. The gravitational acceleration, relative to the center of the perturbation, along the radial direction is (18)

By removing the acceleration in the moving coordinate of the perturbation, the tidal force generated by the gravity is adequately considered. The first two terms on the right-hand side represent the gravitational acceleration generated by the central star, the third and fourth terms are the acceleration by the envelope, and the last term is the gravity of the perturbation itself. The density is (19)

The virial integration therefore is (20)

where is the non-dimensionalized energy (a detailed derivation is given in Appendix D). This is a line integration, and the quantity has the dimension of energy per unit surface. We define a normalizing radius (21)

such that (22)

where is the mass of the envelope within r_p. The energy can be expressed as (23)

with contributions from the three components.

The solutions are found by simultaneously solving for all (24)

The first condition is the condition stated by Eq. (13), while the normalized expression stated by Eq. (17) is used.

We plot the critical density contrast η_c and the corresponding size of the perturbation, , against in Fig. 9, showing that with increasing distance from the central star, a lower density contrast (smaller η) and a lower level of mass concentration (larger δr_p) are needed for the perturbation to be self-gravitating. As η increases at fixed mass, u_p and decrease. At any given distance to a first Larson core, we can therefore derive a critical local density contrast η_c above which the perturbation is self-gravitating and is prone to collapse. Within this distance, all density fluctuations lower than η_c are prevented from collapsing by the tidal field of the central mass and its envelope. The value of η_c is typically higher than a factor of a few, while the radius of the perturbation must be significantly smaller than the distance from the central object. These numbers indicate that the conditions for a perturbation to become unstable are not easy to satisfy, and most perturbations that in the absence of tidal forces would be unstable are therefore rendered stable. This process favors further accretion onto the central object other than distributing the mass between new fragments.

Fig. 9 Critical value of ηc (blue dashed) and corresponding (black solid) as functions of the normalized distance . The gravitational potential energy is calculated in 1D.

5.3.3 Full 3D calculation

The SIS density profile generates a positive tidal force in the radial direction, while the field is compressing in the other two directions. A full 3D calculation can therefore be different from the 1D case, and it is necessary to clarify this effect. For conciseness, the full integration of Eq. (12) in 3D space is presented in Appendix E, and here we simply discuss the results: (25)

where is the normalized gravitational energy. The results are qualitatively similar to the 1D calculations. The values of η_c and corresponding are plotted in Fig. 10 as a function of the normalized distance (magenta). At a short distance, they are very similar to the 1D approximation and η_c must be larger than several, but at a large distance, η_c drops to ~1, implying that the perturbation is easily unstable.

However, we also need to consider the thermal support against self-gravity. We can infer the condition from Eqs. (12) and (17). Including the thermal support does not change the conclusions at small r_p much, while at larger distance, where the envelope density is low, thermal support becomes dominant and tidal forces are unimportant. In this context, η_c becomes an increasing function of r_p, thus imposing a lower limit on η_c (black lines in Fig. 10).

The result sensitively depends on A, the amplitude of the density field. When A is equal to a few, the combination of thermal support and tidal forces stabilizes most perturbations. Only fluctuations denser than about ten times the background density can eventually collapse. This is in good agreement with simulations such as those performed by Girichidis et al. (2011), which show that clouds with r⁻² density profiles are not prone to fragmentation.

5.4 Step 4: final mass from the integrated probability of the formation of a nearby first Larson core

In the previous section we obtained the critical density contrast, η_c, for forming a second self-gravitating core near a point mass M_L as a function of the distance. The question now is under which circumstances would such fluctuations be encountered, or more precisely, what is the probability for a given perturbation amplitude to occur? To answer this question, we need to know the fluctuation spectrum within a collapsing clump around an accreting point mass. Since there is no available statistics, this last step is certainly the least quantitative. To obtain some estimate, we now use the typical density fluctuations produced by supersonic turbulence to determine the characteristic distance at which a perturbation with sufficient amplitude could be found and a second core could form. The mass within this radius is expected to be protected from fragmentation and therefore should be accreted.

We use the simple assumption of a lognormal PDF (e.g., Vázquez-Semadeni 1994; Federrath et al. 2008; Hennebelle & Falgarone 2012), that is, the probability of finding perturbations of η at distance is (26)

where δ = log(1 + η) and σ is the width of the density PDF. The density fluctuation is related to the local turbulent Mach number such that , with b = 0.5. To estimate the local Mach number, we assume that the turbulent Mach number is proportional to the that obtained from the infall velocity. This implies that the collapse is able to amplify existing perturbations until local energy equipartition is reached. Given the mass of the central core and the envelope, the turbulent Mach number therefore is (27)

where the second equivalence is obtained using Eq. (22) and the third using Eq. (21). The coefficient, ϵ ≤ 1, represents the amount of released gravitational energy that goes into the turbulent fluctuations, rather than the infalling component. As estimated in Paper I, A ≃ 10, leading to . We refer to this assumption as model I. We integrate the mass that exceeds the critical density at all r_p, using the density of the fluctuation , the local background multiplied by the relative density fluctuations. Dividing by M_L, we obtain the typical number of self-gravitating cores with a mass at least equal to M_L contained within r_p (28)

where we have used Eq. (21) to normalize the expression. With the definition of , the mass contained within a sphere of radius is equal to aM_L. When , the probability of getting another fragment with mass M_L is equal to 1. The mass enclosed within this radius is accreted by the central object since this mass is protected from fragmentation. Thus the mass of the final object is certainly higher than .

Figure 11 shows M_F as a function of A. The typical values of A ~ 10–20 and ϵ ~ 0.5 give a final stellar mass of ~4 M_L when we integrate up to (black lines), ~8–9 or 20–30 M_L is obtained when we consider (blue lines) or (cyan lines). Up to which radius should be integrated, that is to say, how far the protostar is able to accrete, is a difficult question. Typically, we can estimate that about six fragments (one for each direction) should prevent any further accretion, but on the otherhand, these fragments will attract mass that should not be attributed to the central object. Therefore, three fragments sounds like a reasonable number. In Appendix F we explore another model for the Mach number dependence, which leads to similar conclusions.

We recall that the mass is integrated without considering whether the masses are spatially connect. The value of is thus a lower limit, and so is thederived . Although the model we presented is not very accurate, it nevertheless suggests that the nonlinear factor of ≃ 10 inferred from the simulations is entirely reasonable. Further work is certainly required here.

Fig. 10 Critical value of the perturbation amplitude, ηc (upper panel) and corresponding normalized size of the perturbation (lower panel) as a function of the normalized distance . 3D gravitational energy is used, and thermal support is included for various density profiles. Fluctuations have to be at least several times the local density to be unstable. The low-amplitude perturbations are stabilized by tidal forces or thermal support.

5.5 Tidal protection in simulations

To verify our model of tidal protection that inhibits nearby core formation, we examined the distance between the sinks in our simulations. For each sink, the distance of formation position was calculated for all younger sinks. In Fig. 12 we plot the distances of the first and tenth closest neighbors against the instantaneous sink mass at the moment of neighbor formation for run C10h+. The distances of the first few closest sinks are similar and almost identical in most cases (which suggests that our model needs further refinement). The formation of a nearby sink is generally inhibited within the distance of several hundreds of AU for established sinks.

Fig. 11 Final mass MF, in units of ML, plotted against A, the amplitude of the density envelope (Eq. (10)), and for various ϵ, the efficiency at which infall triggers turbulent fluctuations. The critical radius is equivalent to the amount of protected envelope mass in units of ML, and thus . Black lines show the expected mass when the envelope is truncated at , and MF is expected to be larger than 4 for typical values of A. Blue thin lines show the same result, while truncated at , and the expected MF is increased to ~8–9. Cyan lines show results for .

Fig. 12 Distance of nearby sink formation plotted against sink mass. The distance is calculated at the moment of the neighboring sink formation. The distances from the each sink to its first and tenth (small and large dots interlinked by a line) closest neighbors are plotted against its mass at the moment of the formation of the neighbors. Around an existing sink particle, any neighbor can hardly form within the distance of a few hundreds of AU.

6 Discussions

6.1 Necessary developments

In our simulations, the physics of the first Larson core sets the peak of the stellar distribution, and we infer M_peak ≃ 10 M_Larson. The mass of the first Larson core at the point where it undergoes second collapse is typically on the order of a few 0.01 M_⊙ (Masunaga et al. 1998; Vaytet et al. 2012; Vaytet & Haugbølle 2017), which is very close to the observed peak of the IMF (Bastian et al. 2010). The most recent radiative simulations of core collapse by Vaytet & Haugbølle (2017) showed that the mass of the first Larson core is ~ 0.02 M_⊙. One of the difficulties encountered by our calculations is that it is not possible to properly describe the second collapse andthe disappearance of the first Larson core as it is being accreted. This causes some uncertainties in the value of the first Larson core that should be used. A smoothed two-slope polytropic eos was used in this study to represent the results of radiative transfer balancing in dense gas. Detailed simulations with complete radiation calculations would give more accurate results, although the numerical resolution is a major obstacle. On the other hand, varying the density at which the sink particles are introduced (runs C10f and C10fd in Fig. 5) does not change the position of the peak, although it reduces the number of objects. This may indicate that the exact way in which the first core is accreted onto the protostar may be not too critical in determining the peak position.

Finally, we reiterate that the physics used in the present paper and its companion may be oversimplified in other aspects since it does not include the magnetic field and heating by accretion luminosity. Both are known to reduce fragmentation (Krumholz et al. 2007a; Bate 2009b; Commerçon et al. 2011; Hennebelle et al. 2011; Myers et al. 2013) and may drastically modify the fragmentation we obtain.

6.2 What sets the peak of the IMF: is there a cosmic conspiracy?

Unlike the prediction by simple gravo-turbulent fragmentation models (e.g., Hennebelle & Chabrier 2008, 2013; Hopkins 2012), which posit that the IMF peak depends on the turbulent Jeans mass, the IMF characteristic mass in the high-density regime has a lower limit imposed by the local thermal balancing and the formation of the first Larson core. The traditional view that the typical fragmentation mass is linked to the Jeans mass is no longer valid in these conditions because self-gravity triggers the development of a power-law density PDF. As discussed in Paper I, the existence of a peak may depend on the environment, however. In somewhat diffuse clouds, or more precisely, when the local thermal support dominates the turbulent support, the stellar distribution is flat (dN∕d log M ∝ M⁰) and the mass of the first Larson core, M_Larson, may truncate the distribution rather than leading to a peak.

One important issue, however, is the initial conditions and the exact way the matter is assembled. Although we verified in Paper I (see the appendix) that the initial fluctuations are not crucial, setting up a dense reservoir without preexisting cores may be a significant oversimplification. In other words, the initial conditions that should be used to describe the ISM are not well understood yet, and we caution that the conditions used in this paper may introduce some biases. In particular, as stressed in Paper I, the high-mass part of the mass spectra tends to present a power law that is too shallow compared to the canonical Salpeter distribution. If dense clusters turn out to be assembled from material where cores are already developed, the preexisting cores may play a stronger role and the peak of the IMF could be inherited from the peak of the core mass function (CMF), as proposed theoretically in Padoan et al. (1997), Hennebelle & Chabrier (2009), Hennebelle (2012), Hopkins (2012), Lee & Hennebelle (2016a), and Lee et al. (2017), and suggested by simulations (Gong & Ostriker 2015; Hennebelle 2018). This picture is suggested by the observations of the CMF, which resembles the IMF (Motte et al. 1998; Alves et al. 2007; André et al. 2010; Könyves et al. 2015), and it could be the dominant mode in some environments. It is remarkable that the peak of the CMF (particularly if an efficiency of ≃ 30–50% is taken into account) and the peak of the stellar distribution imposed by the physics of the first Larson core, as discussed in this paper, lead to values that are similar within a factor of 2–3. This coincidence may be for a large part responsible of the apparent universality of the IMF.

7 Conclusions

We have performed a series of numerical simulations describing the collapse of a 1000 M_⊙ cloud. These numerical experiments explicitly ignored magnetic field and radiative transfer at this stage. The transition from the isothermal to the adiabatic phase is of great importance, and several effective equations of state were explored. This includes varying γ, the effective adiabatic exponent, and ρ_ad, the density at which the transition occurs. We paid great attention to the numerical spatial resolution and its influences on the resulting stellar mass spectrum, in particular, its peak position. The influence of the initial conditions were studied in Paper I.

We found that the isothermal simulations show no sign of convergence as the peak position shifts toward lower masses when the spatial resolution improves. On the other hand, when a polytrope-like eos is used (typically with γ > 4∕3), the peak position becomes independent of the numerical resolution if high enough. Using simple 1D hydrostatic models, we inferred the expected mass of the first Larson core that we compared to the measured peak position of the stellar distribution, for which we compared our simulations and various published results. We found a clear correlation between the two, with the peak position typically about ten times higher than the mass of the first Larson core.

We presented an analytical model that may account for this factor of roughly ten. It is based on the idea that in the vicinity of a point mass, such as a first Larson core, tidal forces are important and tend to quench further gravitational fragmentation by shearing out the density fluctuations. Therefore the mass located in the neighborhood of such a point mass tends to be accreted, which increases the mass of the central objects in consequence. As discussed in Paper I, this effect is true only when the gas is initially dense and turbulent, generating many fluctuations that eventually form stars close one to the other.

Since the physics of the first Larson core is not expected to significantly vary from place to place, this mechanism may provide a robust explanation for the apparent universality of the IMF peak. We emphasize, however, that an IMF-like distribution is produced only if the cloud is sufficiently massive and cold.

Acknowledgements

We thank the anonymous referee for comments that have improved the manuscript. We thank Gilles Chabrier for a critical reading of the manuscript. We thank Michel Rieutord for providing the ESTER code and helping us to run it to explore the influence of the rotation on the mass of a polytrope. This work was granted access to HPC resources of CINES under the allocation x2014047023 made by GENCI (Grand Equipement National 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). Y.-N. Lee acknowledges the financial support of the UnivEarthS Labex program at Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02).

Appendix A Sink accretion scheme

The thresholding scheme for sink accretion by Bleuler & Teyssier (2014) was used in this study. Within the sink accretion radius, four times the smallest cell size, a fraction c_acc of the mass exceeding the density threshold n_sink is accreted onto the sink, and by default, c_acc = 0.75. When the mass from a cell is accreted by multiple sinks, the algorithm treats the sinks in decreasing order in terms of mass to ensure that more massive sinks accrete more mass. We performed a simulation with all parameters identical to C10h, except that c_acc = 0.1. As shown in Fig. A.1, the stellar mass distribution does not differ significantly from that of run C10h (Fig. 2), ensuring that the accretion scheme does not have a strong effect on the final mass of the stars.

Fig. A.1 Stellar distributions of runs with a smoothed two-slope polytropic eos with γ = 5∕3, nad = 1010 cm−3, and cacc = 0.1 instead of 0.75 as in all the other runs. The peak of the distribution is not varied compared to run C10h.

Appendix B Perturbation mass

The mass of the perturbation is equal to (B.1)

where the normalized quantities u = δr∕r_p and u_p = δr_p∕r_p.

Appendix C Tidal forces in the simulations

In the vicinity of an accreting point mass, two contributions have to be considered. First the r⁻² gravitational force induced by the central object, and second the r⁻¹ force associated with the r⁻² density profile envelope. Both have a positive eigenvalue, λ₃, and two negative ones, λ₁ = λ₂. It can be shown that for the r⁻² gravitational field, λ₃∕λ₁ = −2 while for the second one λ₃∕λ₁ = −1.

In order to assess the importance of tidal forces in our simulations, we calculated the gravitational stress tensor. We first averaged λ₃ and λ₁ in density bins, and then plotted the absolute value of their ratio. Figure C.1 shows the result from run C10h at 6.1 and 9.3 kyr. Low densities correspond to the outer region of the cloud that is dominated by the tidal field of the global density distribution of the cloud and by the central cluster. At high densities (> n_sink = 10¹⁰ cm⁻³), corresponding to regions in the vicinity of sink particles, the destructive tidal force is strong, such that no core formation is allowed nearby. This illustrates the dominating effect of the gravitational field generated by the sink particle.

Fig. C.1 Tidal forces in run C10h. The absolute value of the ratio between the highest and lowest eigenvalues of the gravitational stress tensor λ3∕λ1 is plotted against local density at 6.1 and 9.3 kyr. At very low density, corresponding to the outer region of the cloud, the gravitational field is dominated by the global density distribution of the cloud and by the central cluster, and thus − 2 < λ3∕λ1<-1. The tidal forces become less destructive at slightly higher densities. At extremely high density (> nsink = 1010 cm−3), corresponding to regions in the vicinity of existing sinks, the tidal forces become locally very destructive and prevent nearby sink formation.

Appendix D 1D calculation of the gravitational energy

The expression of the 1D gravitational energy is (D.1)

Appendix E 3D calculation of gravitational energy

We define ξ the angle between r_e and δr, and θ the angle between δr and − r_p. Projected onto the direction of δr, the gravitational fields have the values (E.1) (E.2) (E.3)

where the subscripts a and c signify the value at the point considered and the value at the center of the perturbation. The local density is (E.4)

We have the trigonometric relations

and we deduce . The contribution to the gravity at the center of the perturbation from the central star and from the envelop behave in the same way, and the normalized integration gives (E.7)

Considering the integration of Eq. (12) separately for each gravitation term without subtracting the mean gravity field, we obtain (E.8) (E.9) (E.10)

leading to (E.11) (E.12)

where u = δr∕r_p, and u_p = δr_p∕r_p. The condition for a second core forming near the initial central core is that E_g (r_p, δr_p) ≤ 0 and that M(r_p, δr_p, η) = M_L. SubstitutingEqs. (E.9), (E.11), and (E.12) and (17) into Eq. (12) reduces the latter to Eq. (25).

Appendix F Density fluctuations: model II

Since the density fluctuations are important and they control the fragmentation around the accreting objects, we explored a second model (model II) to determine the robustness of our conclusions. We assumed that the turbulence follows a Larson relation inside the cloud and is conserved during collapse (this assumption may not be accurate as the collapse likely generates turbulent fluctuations), we have (F.1)

where v_t, r₀, v_t,cloud, r_cloud, and are the turbulent velocity at r_p, the radius that corresponds to the mass enclosed inside radius r_p before the local mass concentration (at cloud average density), the turbulent velocity at cloud scale, the radius of the cloud, and the Mach number of the cloud. The mass M(r_p) contained inside r_p takes the envelope and the central mass into account and can be expressed with the normalized form . This leads to (F.2)

The local lognormal dispersion is therefore (F.3)

The initial Mach number of our canonical run is ~20, and M_L∕M_cloud ~ 10⁻²∕10³ = 10⁻⁵. These values typically yield –3. An important difference between the two assumptions made in models I and II is that the second estimate depends on the initial conditions, while the first is entirely local.

We performed numerical integrations for several values of . Figure F.1 shows the resulting , without considering the thermal support. We defined the critical radius at which one self-gravitating mass M_L is expected, as marked with circles in the figure (model I is also shown for comparison).

The resulting is displayed in Fig. F.2 and should be compared to Fig. 11. When we also consider the thermal energy of the perturbation, the protected radius is increased since more support is present. In Fig. F.2 we also show results with thermal support using A = 2, 5, 10, and 20. The effect of thermal support is more prominent when A is small, imposing a lower limit on below which a star prevents its surrounding gas form fragmenting and will grow to become massive. The physical interpretation is that when the envelope density is low, stronger fluctuations are needed to create multiple stars, otherwise one star will dominate and prevent any further formation once it is formed. This is in good agreement with the result of Girichidis et al. (2011), who reported that clumps with an r⁻² profile and weak turbulence do not fragment. For A ~ 10, as measured in Paper I (left panels of Fig. 4), the condition gives at least M_F ~ 10M_L for any lower than ~2. Figure 8 shows that this is in good agreement with the value inferred from the simulations in the sense that the model predicts that the accreting envelope, protected from gravitational fragmentation, is several times the mass of the first Larson core.

Fig. F.1 Probable number of self-gravitating cores found within radius . Upper panel:Model I plotted with ϵ = 1 and for several values of the density envelope amplitude, A. Lower panel: Model II without thermal support for several values of the global cloud Mach number, . The circles mark the values where one self-gravitating core is found. In the range of our interest, that is, where equals a few, there is an almost linear dependance between and .

Fig. F.2 Final mass MF in units of ML for model II. MF is plotted against the level of turbulence for cases with A = 2, 5, 10, and 20, and without thermal support (A →∞, magenta). Model II also shows that MF is expected to be larger than 3 when the envelope is truncated at , and > 5 when (blue thin lines).

References

All Tables

All Figures

Fig. 1 Stellar distribution of isothermal runs with varied resolutions of 17, 8, and 4 AU in colored histograms, with time and total accreted mass shown in the legend. Lognormal fits are presented with dashed thin curves. Power laws of M−3∕4 and M−1 are shown for comparison with black and gray lines, respectively. The characteristic mass of the stellar distribution shifts to lower values with increasing resolution.

In the text

	Fig. 2 Stellar distribution of runs with a smoothed two-slope polytropic eos with γ = 5∕3 and nad = 1010 cm−3, with varying resolutions of 17, 8, 4, and 2 AU. The characteristic mass of the stellar distribution is ≃ 0.1 M⊙, regardless of numerical resolution.
In the text

	Fig. 3 Stellar distributions of runs with a smoothed two-slope polytropic eos with γ = 5∕3, lowered nad at 3 × 109 (C3h–) and 109 cm−3 (C1h–, C1h– –). The characteristic mass of the stellar distribution increases with decreasing nad. At nad = 109 cm−3, the stellar distribution becomes top-heavy with improved resolution from 17 to 8 AU.
In the text

	Fig. 4 Stellar distributions of runs with a smoothed two-slope polytropic eos with γ = 4∕3, nad = 1010 cm−3 (C10s–) and 109 cm−3 (C1s–). The characteristic mass of the stellar distribution increases by a factor ~ 3 when nad is lowered by a factor ten.
In the text

	Fig. 5 Stellar distributions of runs with a full eos with nad = 1010 cm−3 (C10f, C10fd) and 109 cm−3 (C1f). The transitions to γ = 5∕3, 7∕5, and 1.1occur at nad, 30 nad, and 30 000 nad. Using a higher sink density threshold (C10fd) results in fewer low-mass stars but does not change the peak of the stellar distribution. Decreasing nad by a factor ten increases the peak mass by ~10.
In the text

Fig. 6 First Larson cores for various equations of state. The hydrostatic equilibrium equations are integrated with n0 = 100 nad as boundary condition. We present cases with γ = 5∕3 (blue), γ = 4∕3 (cyan), full eos (red), and isothermal (black). The legends correspond to the labels of the simulations, with solid lines for nad = 1010 cm−3 and dashed lines for 109 cm−3. We also show the case γ = 7∕5 with nad = 2.45 × 1010 cm−3 by Bate et al. (2003, green). The values of n0 are chosenfor illustrative purpose, and the actual values reached in simulations depend on the resolution (highest density resolved, see Fig. 7). The circle indicates the radius at which the density reaches 107   cm −3, a reference ambient density that truncates the core. As long as the density profile decreases steeply enough, the choice of this truncation density has no strong effect on the resulting mass.

In the text

Fig. 7 Properties of the first Larson core plotted against the central density n0. The hydrostatic equilibrium Eq. (7) is solved with various nad and γ values, and truncated at the ambient density n = 107 cm−3. Three quantities are shown: the radius at which the density decreases to n = 107   cm−3, the local power law exponent α = −rdn∕dr∕n, and the integrated mass. Top left: γ = 5∕3; top right: γ = 4∕3; bottom left: full eos. Bottom right: eos taken from Bate et al. (2003, solid), Bate (2005, dashed), and Bonnell et al. (2011, dotted). Each time the isothermal case is plotted in black for comparison. A core is well defined only when α ≫ 2, and thus giving a characteristic mass of the first Larson core. When a core is well defined, its mass is almost constant regardless of n0 when γ = 4∕3 or full eos is used. In contrast, the core mass increases with increasing n0 when γ = 5∕3, which explains the peak mass shift between runs C1h– and C1h– –.

In the text

Fig. 8 Peak mass of the sink mass spectrum plotted against the mass of the first Larson core calculated with the eos we used in the simulations. The peak mass is read from the stellar distribution, with the bin size as error. The error bars from literature results are larger as they have lower statistics and higher fluctuations in the spectrum. The mass of the first Larson core is read from Fig. 7, with error bars indicating the mass range with the probable n0 range. Whenwe examine this by eye, there is a correlation Mpeak ~ 5MLarson. We also performed a simple linear regression, without considering the errors, to determine the scaling factor, and this led to .

In the text

	Fig. 9 Critical value of ηc (blue dashed) and corresponding (black solid) as functions of the normalized distance . The gravitational potential energy is calculated in 1D.
In the text

Fig. 10 Critical value of the perturbation amplitude, ηc (upper panel) and corresponding normalized size of the perturbation (lower panel) as a function of the normalized distance . 3D gravitational energy is used, and thermal support is included for various density profiles. Fluctuations have to be at least several times the local density to be unstable. The low-amplitude perturbations are stabilized by tidal forces or thermal support.

In the text

Fig. 11 Final mass MF, in units of ML, plotted against A, the amplitude of the density envelope (Eq. (10)), and for various ϵ, the efficiency at which infall triggers turbulent fluctuations. The critical radius is equivalent to the amount of protected envelope mass in units of ML, and thus . Black lines show the expected mass when the envelope is truncated at , and MF is expected to be larger than 4 for typical values of A. Blue thin lines show the same result, while truncated at , and the expected MF is increased to ~8–9. Cyan lines show results for .

In the text

Fig. 12 Distance of nearby sink formation plotted against sink mass. The distance is calculated at the moment of the neighboring sink formation. The distances from the each sink to its first and tenth (small and large dots interlinked by a line) closest neighbors are plotted against its mass at the moment of the formation of the neighbors. Around an existing sink particle, any neighbor can hardly form within the distance of a few hundreds of AU.

In the text

	Fig. A.1 Stellar distributions of runs with a smoothed two-slope polytropic eos with γ = 5∕3, nad = 1010 cm−3, and cacc = 0.1 instead of 0.75 as in all the other runs. The peak of the distribution is not varied compared to run C10h.
In the text

Fig. C.1 Tidal forces in run C10h. The absolute value of the ratio between the highest and lowest eigenvalues of the gravitational stress tensor λ3∕λ1 is plotted against local density at 6.1 and 9.3 kyr. At very low density, corresponding to the outer region of the cloud, the gravitational field is dominated by the global density distribution of the cloud and by the central cluster, and thus − 2 < λ3∕λ1<-1. The tidal forces become less destructive at slightly higher densities. At extremely high density (> nsink = 1010 cm−3), corresponding to regions in the vicinity of existing sinks, the tidal forces become locally very destructive and prevent nearby sink formation.

In the text

Fig. F.1 Probable number of self-gravitating cores found within radius . Upper panel:Model I plotted with ϵ = 1 and for several values of the density envelope amplitude, A. Lower panel: Model II without thermal support for several values of the global cloud Mach number, . The circles mark the values where one self-gravitating core is found. In the range of our interest, that is, where equals a few, there is an almost linear dependance between and .

In the text

	Fig. F.2 Final mass MF in units of ML for model II. MF is plotted against the level of turbulence for cases with A = 2, 5, 10, and 20, and without thermal support (A →∞, magenta). Model II also shows that MF is expected to be larger than 3 when the envelope is truncated at , and > 5 when (blue thin lines).
In the text