How First Hydrostatic Cores, Tidal Forces, and Gravoturbulent Fluctuations Set the Characteristic Mass of Stars

The simulations were run with the adaptive mesh refinement (AMR) magnetohydrodynamics code RAMSES (Teyssier 2002; Fromang et al. 2006) using a base grid of (256)³ cells. Six to seven AMR levels have been added, leading in the first case to an effective resolution of 16,384 and a spatial resolution of about 4 au. The Jeans length is resolved with at least 10 points.

2.1.1. Initial Conditions

2.1.2. Equation of State and Sink Particles

As found in Paper II, the effective equation of state (EOS) is playing a fundamental role in setting the peak of the IMF. We adopt the following expression:

$$\begin{eqnarray}&&T={T}_{0}\left\{1+\displaystyle \frac{{\left(n/{n}_{\mathrm{ad}}\right)}^{({\gamma }_{1}-1)}}{1+{\left(n/{n}_{\mathrm{ad},2}\right)}^{({\gamma }_{1}-{\gamma }_{2})}}\right\},\end{eqnarray} \tag{ 1 }$$

where T₀ = 10 K, ${n}_{\mathrm{ad}}$ = 10¹⁰ cm⁻³, ${n}_{\mathrm{ad},2}$ = 3 × 10¹¹ cm⁻³, γ₁ = 5/3, and γ₂ = 7/5. This EOS mimics the thermal behavior of the gas when it becomes nonadiabatic. It is well known that an EOS with an adiabatic exponent larger than 4/3 leads to a hydrostatic first core (see Appendix A for orders of magnitude). As discussed in Paper II with this EOS, the mass of the FHSC is about 0.02–0.03 M_⊙. Let us stress that its radius is about 15 au, which is larger than the 5 au usually assumed. This makes the FHSC easier to resolve.

We used the sink particle algorithm from Bleuler & Teyssier (2014; see Paper II for a brief description). Sink particles are formed at the highest refinement level at the peak of clumps whose density is larger than n = 10¹¹ cm⁻³. The sinks are introduced at a density n ≥ 10¹² cm⁻³. In addition, the clumps must be thermally and virially unstable, and they should also be contracting. The accretion onto the sinks occurs within a sphere of radius 4dx, that is, 4 times the smallest resolution element in the simulation. At each time step, 10% of the gas within this sphere and above the density at which sinks are introduced is retrieved from the computational cells and assigned to the sink. In Paper II, the value at which the sinks are formed and accrete and the fraction of the dense gas that is transferred to the sinks have been varied. It has been found that they do not have a strong effect on the resulting mass spectrum. The influence of numerical resolution, to which the sink radius is proportional, has also been checked, and it has been found that it has no significant influence on the peak position.

To illustrate the simulation results, the column density is portrayed in Figure 1. Left is run 1, and right is run 2. The sink particles have been overlaid and are represented by the dark dots. In both runs, the dense gas is organized in dense filaments. The distributions of the filaments are, however, different. In run 1, they tend to be more numerous and more uniformly distributed than in run 2, which is organized in a few massive objects in which most of the star formation is taking place. The simulations have been run respectively for a time of 0.013 Myr and 1.04 Myr (which corresponds to 0.04 Myr after gravity is switched on). In the two cases, this is about one freefall time of the dense gas.

Zoom In Zoom Out Reset image size

Figure 2 displays the mass distribution of the sink particles for run 1 (left) and run 2 (right) at several time steps that correspond to various total accreted masses. In run 1, a marked peak at ≃0.2 M_⊙ is inferred. At higher mass, the distribution is compatible with a power law M^−0.75 as found in Paper I. In run 2, the peak is less clearly pronounced because there are more low-mass objects than in run 1, and the corresponding distribution tends to be flatter. There is, nevertheless, also a peak located at 0.1–0.2 M_⊙. At higher mass, the distribution again is broadly compatible with a power law M^−0.75.

Zoom In Zoom Out Reset image size

This is in good agreement with the results obtained in Paper I and II.

At the end of the calculations, the fraction of gas that has been accreted onto the stars is about 17% for run 1 and 10% for run 2. The total number of sinks formed is respectively 762 and 723. While these numbers are sufficient to produce reasonable statistics, the question arises as to whether the mass spectrum may vary at a later time. However, Figure 2 reveals that the shape of the mass spectrum, and particularly the peak, does not vary very significantly between the times where 30 M_⊙ and 170 M_⊙ has been accreted. In Paper I, we were able to go up to the point where 300 M_⊙ of gas has been accreted. Therefore, we expect this integration time to be sufficient to estimate the mass spectrum. Another related issue is the stellar feedback, which should become important once enough mass has been converted into stars, and may even disperse the clouds before all of the gas is turned into stars.

We now turn to the analysis of the sink neighbor distribution. Our goal here is to investigate to what extent the presence of sinks in the immediate vicinity may influence the peak of the stellar mass function as proposed in Paper II. To address this question, for each sink we have calculated the number of neighbors, n_neigh, located below a series of specific distances and for various sink ages. The statistics are obtained by analyzing the sink configurations during the entire simulation. In practice, we consider all sinks formed in the simulations, and at each time step we calculate the number of neighbors lying within a series of specific distances. To define the number of neighbors, we simply select all time steps for which the corresponding objects have an age below a specific value and take the maximum number of neighbors found. In practice, we do not find much fluctuations in the number of neighbors. Once its maximum value is reached, it rarely decreases. As it is more common to draw the statistics of the number of systems (that is, the number of groups of neighboring stars) and also because this is what the analytical model developed below produces, we mainly present the statistics of the number of systems of size n_sys = n_neigh + 1. The number N_sys of systems with n_sys objects is related to the number N_neigh of sinks with n_neigh numbers through N_sys = N_neigh/(n_neigh + 1). Note that the systems we are interested in may not be necessarily gravitationally bound (although it is the case for most of them), because we specifically address the question of the influence of fragmentation onto the sink mean mass.

The results are portrayed in Figure 3. The top panels show the cumulative histogram of system numbers (N_sys, solid lines) and sink numbers (N_neigh, dashed lines) of various ages with (up to) n_neigh = n_sys − 1 neighbors at a distance below 200 au. As can be seen at early time (t ≤ 100 yr), more than one-half of the sinks have zero to two neighbors below 200 au, and about three-quarters of the systems are constituted by a single object. Interestingly, a fast evolution is visible between t = 100 yr and t = 1000 yr, especially for n_neigh < 2. Indeed, we see that the distance between the curves of age 100 and 1000 yr is larger than what it is between 2000 and 3000 yr, for example. This indicates that fragmentation in the vicinity of young sinks is commonly occurring. A significant number of objects are surrounded by a relatively large number of neighbors, n_neigh ≥ 10, but the corresponding system number is small. This corresponds to a restricted number of regions where the clustering is particularly important.

Zoom In Zoom Out Reset image size

The middle panels display the cumulative histogram of the number of systems with n_neigh = n_sys − 1 objects for several distances and at time 10² yr (dotted line) and at time 10³ yr (solid line). As expected, the number of neighbors increases with the distance, but we clearly see that this increase is not very fast and is clearly sublinear with distance.

To quantify this better, the bottom panels show the mean system size $\langle {n}_{\mathrm{sys}}\rangle$ as a function of distance. The red curves are for t ≤ 10² yr, while blue ones are for t ≤ 10³ yr. The black curves correspond to the analytical model and will be discussed later. By comparing the red and blue curves of run 1 (left panel), it is very clear that $\langle {n}_{\mathrm{sys}}\rangle$ below 100 au is ≤2 at t = 100 yr and ≃2 at t = 1000 yr. This is similar for run 2. At 1000 au, there is more difference between run 1 and run 2, with about four to five objects for the former and less than or equal to three for the latter.

To demonstrate that the fragmentation close to existing objects is critical in setting the peak of the IMF, we seek a clear test. We cannot increase the formation of sinks without changing the EOS, which affects the mass of the FHSC and the peak of the stellar mass function as shown in Paper II. We can, however, prevent the formation of new sinks below a given distance from an existing sink by enforcing it in the simulation. We have therefore performed two runs, called respectively run 1a and run 1b, where sink formation is prevented below respectively 140 au and 280 au from any existing sink.

The resulting mass spectra are portrayed in Figure 4. Clearly, the sink distributions shift to higher masses as the distance from an existing sink, below which new sinks cannot form, increases. The solid lines are for a total mass of sinks equal to 150 M_⊙, while dashed ones are for 50 M_⊙. Both the peak and the mass of the most massive sinks shift by nearly a factor of 2. In run 1b, there are about two times fewer sinks, and the sinks are about two times more massive than in run 1. These numbers are entirely compatible with the neighbor distributions displayed in Figure 3, which reveal that one or two objects on average sit at a distance below 100–200 au from an existing sink.

Zoom In Zoom Out Reset image size

These simulations demonstrate that indeed the close fragmentation that limits the accretion onto most of the objects is a determinant for their masses. By preventing fragmentation at a distance of 100–200 au from existing sinks, the peak of the stellar mass function shifts by a factor of 2 because there is more gas available for accretion. In the presence of more fragments, less gas is available and less-massive stars form.

As fragmentation around existing objects appears to be determinant, we need to characterize the physical conditions that prevail around newly formed sinks. This is necessary to construct a quantitative model for explaining this fragmentation. To perform this estimate, we have selected all sinks at the age of 1000 yr within 10 snapshots regularly spaced in time. Since we want to limit the influence that sinks may have on the gas dynamics, we select the objects that have at most one neighbor within 200 au. We focus on three quantities that are of particular importance for the analytical model developed below: the mean density, the Mach number, and the probability density function (PDF) of the density fluctuations. For each selected object, we have computed the mean density, radial velocity, and rms velocity (from which the mean radial velocity has been subtracted) in concentric spherical shells with logarithmic spacing. We have also computed the PDF of the density fluctuations (that is, the density divided by the mean density within the shell) within two spheres of radius 240 au and 800 au. We have then summed all of the results to compute the mean values.

Figure 5 displays the results. The top panels show the mean density (red curve) as well as the standard deviation from it (blue shaded area). The blue line is the singular isothermal sphere, n_SIS = ρ_SIS/m_p = c_s²/(2πGr²)/m_p. The mean density is proportional to r⁻² and about 10 and 7 times n_SIS for run 1 and run 2, respectively.

Zoom In Zoom Out Reset image size

The Mach number is shown in the middle panels. For run 1, it goes from about 5 at 10 au to 10 at 300 au, while for run 2 it is more on the order of 4–6. In both cases, there are considerable fluctuations, which go from zero to about 2 times the mean value. As discussed below, such Mach number values are indeed expected from energy equipartition.

The bottom panels portray the PDF of the density fluctuations. Note that the natural logarithm is used to compare with the usual log-normal distribution (Hennebelle & Falgarone 2012). The upper and lower blue curves represent the fluctuation PDF within 800 au and 240 au, respectively. The red curve is a log-normal distribution with a width ${\sigma }^{2}=\mathrm{ln}(1+{b}^{2}{{ \mathcal M }}^{2})$ for ${ \mathcal M }=6$ and b = 0.8. It is not a fit, but as can be seen, it matches reasonably well with the distributions, except for the high-density power-law tails that are reminiscent of the gravitational infall that is taking place. For run 1, these values are entirely reasonable because the Mach number goes from 6 to 10 (the b parameter is obviously degenerated, so b = 0.6 with ${ \mathcal M }=8$ is equivalent). It is compatible with the idea that many compressible modes (Federrath et al. 2008) are generated by the infall. Interestingly, the shapes of the four PDFs are all similar, and the value b = 0.8 with ${ \mathcal M }=6$ may sound a bit large for run 2. This is likely a consequence of the gravoturbulence that dominates in these collapsing regions. Indeed, the power-law tails, which are of gravitational origin, indicate that the fluctuations with respect to the mean density (which we recall is itself proportional to ρ_SIS) are themselves collapsing.

We now turn to the analysis of the self-gravitating gas structures that lead to the formation of stars and that we name "cores," although we warn that the structures we identify may be too dense and also not sufficiently isolated to correspond to most of the observed cores (Ward-Thompson et al. 2007). The underlying picture is, however, the same, having a coherent gas structure that is dominated by their gravity. The motivation to perform such an analysis is to understand how fragmentation occurs in the vicinity of protostars; we need to understand how self-gravitating perturbations develop. In particular, as our analytical model below uses the virial theorem, we want to verify that the cores are virialized, but we also want to estimate the importance of surface terms that appear in the virial theorem.

To identify the cores, we proceed as in Hennebelle (2018) and Ntormousi & Hennebelle (2019), using the HOP algorithm (Eisenstein & Hut 1998) with a density threshold of 10⁸ cm⁻³. As the algorithm captures many density fluctuations that are close to the density threshold, we keep only those for which the gravitational energy dominates over the thermal one. Figure 6 presents various statistics inferred from this sample.

Zoom In Zoom Out Reset image size

The top panel displays the core mass spectrum. Interestingly, it peaks at about (1–2) × 10⁻² M_⊙, that is, close to M_L. At large masses, it drops relatively steeply. This is, therefore, in good agreement with the idea that cores with a mass close to M_L form an object (a sink in our case) with a mass that is about M_L and then further keep accreting the surrounding gas until the coherence with it is lost.

Before estimating the various contributions of the virial theorem, it is worth remembering that a virialized core is such that

$$\begin{eqnarray}&&{E}_{{\rm{g}}}+2{E}_{\mathrm{th}}+2{E}_{\mathrm{kin}}-{\mathrm{VP}}_{\mathrm{ext}}-{\mathrm{VP}}_{\mathrm{ram}}=0,\end{eqnarray} \tag{ 2 }$$

where the various terms are the gravitational energy and the thermal and kinetic energies, while VP_ext and VP_ram stand for the integrated thermal and ram pressure in the outskirt of the cores. As detailed in Appendix C, their exact expressions are

$$\begin{eqnarray}\begin{array}{rcl}{\mathrm{VP}}_{\mathrm{ext}} & = & \oint {P}_{\mathrm{th}}\,{\boldsymbol{r}}\cdot {\boldsymbol{dS}},\\ {\mathrm{VP}}_{\mathrm{ram}} & = & \oint \rho \,({\boldsymbol{v}}\cdot r)({\boldsymbol{v}}\ \cdot \ {\boldsymbol{dS}}).\end{array}\end{eqnarray} \tag{ 3 }$$

The exact values of these two terms are not straightforward to estimate. Dimensionally they are similar to the thermal and kinetic energy, respectively, but their exact values require a surface integration. Particularly important are the ratios VP_ext/2E_th and VP_ram/2E_kin. As discussed in Appendix C, it is expected that for small cores the contrast between the bulk and the surface thermal and kinetic energies is likely small, while it becomes large for massive cores. Dib et al. (2007) have estimated the importance of the virial surface terms in the context of collapsing cores and found them to be comparable with the volume terms. The top right panel of Figure 6 shows the mean value of (−E_grav + VP_tot)/(2E_kin + 2E_th) in logarithmic bins of mass. The selected cores appear to be remarkably virialized, therefore stressing the importance of the surface terms. It is therefore important to estimate them quantitatively to develop analytical models.

The middle panels of Figure 6 portray the mean values of VP_ext/2E_th and VP_ram/2E_kin per bin of core mass. It is found that both of them are almost zero for M ≥ 0.1 M_⊙. For M < 10⁻² M_⊙ ≃ M_L, it is found that VP_ext/2E_th ≥ 0.5 while VP_ram ≃ 2E_kin. Between these two regimes, both ratios decrease linearly with log M. As discussed in Appendix C, however, it is likely the case that the Mach number within the cores is more relevant for inferring the dependence of their ratios. Therefore, the bottom panels display the same quantities as a function of the Mach numbers. We see that roughly linear dependencies are inferred. At small Mach numbers, the ratios are the highest with values near unity, while for large Mach numbers, the surface terms become negligible. The following expressions inferred from these plots will be used:

$$\begin{eqnarray}\begin{array}{rcl}{\mathrm{VP}}_{\mathrm{ext}} & \simeq & {E}_{\mathrm{th}}\times (0.6\mbox{--}0.1{ \mathcal M }),\\ {\mathrm{VP}}_{\mathrm{ram}} & \simeq & {E}_{\mathrm{kin}}\times (1.3\mbox{--}0.18{ \mathcal M }).\end{array}\end{eqnarray} \tag{ 4 }$$

In this paper, we attempt to construct a "minimum" model that characterizes the essence of the physics responsible for determining the peak of the IMF. In this model, star formation proceeds chronologically throughout the following fundamental stages:

• (1)  
  At the very early stages of star formation, large-scale turbulence generates overdense fluctuations with a power-law mass spectrum. Some of these fluctuations, which can be seen as initial mass reservoirs, become gravitationally unstable, decoupling themselves from the surrounding low-density medium.
• (2)  
  In order to lead eventually to the formation of a protostar, the above density fluctuations must have a minimum mass equal to the mass of the FHSC, M_L, which is typically of the order of a few 10⁻² M_⊙ (Larson 1969; Masunaga et al. 1998; Vaytet et al. 2013; Vaytet & Haugbølle 2017). As discussed in Appendix A, the main reason is that cores with M < M_L are not massive enough to trigger the second collapse without radiating away a large fraction of their thermal energy. However, it takes a long time for objects of mass below M_L to cool (see the discussion in Appendix A).
• (3)  
  Within the envelope, turbulence-induced density fluctuations are likely to be generated on top of the average r⁻² profile, leading to the formation of new collapsing fragments, within the same reservoir, providing their own mass is also at least M_L. These new objects, which we call shielding fragments, accrete the gas in their surrounding and thus prevent further gas accretion onto the central object. Another effect that likely plays a role in the protocluster is that the newly formed object will act differently on the central object and its gas reservoir because the latter dissipates kinetic energy but not the former. This would push the star away from its initial mass reservoir and certainly reduces, or even possibly stops in the most extreme case, the accretion onto the central object. Appendix B explores this effect and suggests that this indeed is happening.
• (4)  
  The total number of such fragments within the collapsing envelope is thus crucial in setting up the final mass of the central object. It is where tidal forces induced by the central object and the surrounding envelope play a crucial role by shearing apart density fluctuations. Therefore, at a given distance r from the central core, only density fluctuations above a certain threshold accounting for this tidal shear contribution can become unstable and lead to the formation of a nearby core of mass ≥ M_L.
• (5)  
  Therefore, the final mass of a collapsing core is about the mass contained within a typical accretion radius r_acc in the envelope, which is the radius containing the characteristic number n_f of unstable fluctuations necessary to prevent accretion of gas located farther away than r_acc.

Whereas the first item is the basic building block of the gravoturbulent theories (Padoan et al. 1997; Hennebelle & Chabrier 2008; Hopkins 2012), steps 2–5 highlight the importance of further fragmentation within the collapsing prestellar core envelope and the key role played by tidal interactions in limiting the final mass of the central object, as revealed by recent numerical simulations of a collapsing 1000 M_⊙ clump (Paper II).

Quantifying precisely this series of processes is difficult because they are all complex and nonlinear. The number of fragments, n_f, necessary to halt accretion for r ≥ r_acc, in particular, is not straightforward to estimate, as it is a highly nonlinear phenomenon that involves N-body interactions. We thus start with a set of simplifying approximations, among which we assume that n_f is typically equal to 2–5, as suggested by the numerical simulation results. In reality, as shown in Figure 3, there is a distribution of n_f values, which may vary from place to place with the local conditions. However, as will be seen later, the exact value of n_f has a modest influence on the result.

3.2.1. Mean Envelope Density

The density distribution in the collapsing envelope around the accreting object is assumed to follow that of a collapsing sphere (Larson 1969; Shu 1977; Murray & Chang 2015):

$$\begin{eqnarray}&&{\rho }_{\mathrm{env}}(r)=\displaystyle \frac{{{Ac}}_{{\rm{s}}}^{2}}{2\pi G}\displaystyle \frac{1}{{r}^{2}},\end{eqnarray} \tag{ 5 }$$

where A defines the amplitude of the density profile and c_s denotes the speed of sound. As shown above (see Figure 5), numerical simulations suggest A ∼ 5–10. For convenience, we define r_e,L as M_L = 2A(${c}_{{\rm{s}}}^{2}/G$)r_e,L. That is to say, r_e,L is the radius for which the mass within the envelope is equal to M_L. We stress that it is not the radius of the Larson core itself, which is typically on the order of 5–10 au. For typical values, A = 10 leads to r_e,L ≃ 20 au. This yields

$$\begin{eqnarray}&&{M}_{\mathrm{env}}(r)=\int {\rho }_{\mathrm{env}}(r)4\pi {r}^{2}{dr}=2A\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{G}r={M}_{{\rm{L}}}\widetilde{r},\end{eqnarray} \tag{ 6 }$$

where $\widetilde{r}$ = r/r_e,L.

3.2.2. Density Fluctuations within the Envelope

As mentioned above, the envelope is dominated by turbulent motions, leading to density fluctuations δ_ρ on top of the mean density ρ_env(r), given, in reasonable approximation, by a log-normal distribution

$$\begin{eqnarray}&&P(\widetilde{r},{\delta }_{\rho }(\widetilde{r}))=\displaystyle \frac{1}{\sqrt{2\pi }\sigma (\widetilde{r})}\exp \left[\displaystyle \frac{-{\left({\delta }_{\rho }(\widetilde{r})+{\sigma }^{2}(\widetilde{r})/2\right)}^{2}}{2{\sigma }^{2}(\widetilde{r})}\right],\end{eqnarray} \tag{ 7 }$$

with δ_ρ($\widetilde{r}$) = ln (ρ($\widetilde{r}$)/ρ_env($\widetilde{r}$)) = ln (1 + η($\widetilde{r}$)), where η($\widetilde{r}$) = ρ($\widetilde{r}$)/ρ_env($\widetilde{r}$) − 1, and σ defines the width of the density PDF. The density fluctuation is related to the local turbulent Mach number as ${\sigma }^{2}(\widetilde{r})=\mathrm{ln}\left(1+{b}^{2}{ \mathcal M }{\left(\widetilde{r}\right)}^{2}\right)$ (e.g., Hennebelle & Falgarone 2012). To estimate the local Mach number ${ \mathcal M }(\widetilde{r})$, we assume that the energy in turbulent motions is a fraction  ≤ 1 of the gravitational energy:

$$\begin{eqnarray}&&{ \mathcal M }(\widetilde{r})=\displaystyle \frac{1}{{c}_{{\rm{s}}}}\sqrt{\epsilon \displaystyle \frac{G[{M}_{{\rm{L}}}+{M}_{\mathrm{env}}(r)]}{r}}=\sqrt{2A\left(1+\displaystyle \frac{1}{\widetilde{r}}\right)\epsilon }.\end{eqnarray} \tag{ 8 }$$

As seen from Figure 5, which displays the PDF of the density fluctuations, δρ/ρ, around sink particles, the log-normal PDF with b = 0.8 and ${ \mathcal M }\simeq 6$ is a good approximation. This corresponds to  ≃ 1. Note that the Mach number distribution seen in Figure 5 shows a tendency for ${ \mathcal M }$ to increase with r instead of decreasing. This is obviously a consequence of the large-scale environment that is not accounted for in our simple estimate.

Let us stress that the typical scales discussed here are a few hundred astronomical units, which are difficult to probe observationally because of the lack of reliable tracers. Supersonic collapsing motions are clearly observed toward the center of at least some cores (di Francesco et al. 2007; Ward-Thompson et al. 2007), and it may be difficult to separate the infall from the turbulence. Also, the regions under investigation correspond to massive star-forming clumps (Lee & Hennebelle 2016; Traficante et al. 2018) rather than low-mass cores.

3.2.3. Instability Threshold

As discussed above and calculated in Paper II, only density fluctuations exceeding a certain threshold, η_crit ($\widetilde{r}$), eventually collapse. This is because the self-gravity of the perturbation must supersede the tidal forces and the local thermal and turbulent support. Writing the virial theorem for a density perturbation, located at r_p, of size $\delta {r}_{p}$ and of overdensity η_p = ρ_p/ρ_e − 1, we get

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{vir}}({r}_{p},\delta {r}_{p},{\eta }_{p}) & = & {E}_{{\rm{g}}}({r}_{p},\delta {r}_{p},{\eta }_{p})+2{E}_{\mathrm{ther}}\\ & & +\ 2{E}_{\mathrm{kin}}-{\mathrm{VP}}_{\mathrm{ext}}-{\mathrm{VP}}_{\mathrm{ram}}.\end{array}\end{eqnarray} \tag{ 9 }$$

In this expression, the gravitational energy is given by

$$\begin{eqnarray}&&{E}_{{\rm{g}}}({r}_{p},\delta {r}_{p},{\eta }_{p})={\int }_{{V}_{p}}({\rho }_{{\rm{e}}}+{\rho }_{p})\,({{\boldsymbol{g}}}_{{\rm{L}}}+{{\boldsymbol{g}}}_{{\rm{e}}}+{{\boldsymbol{g}}}_{p})\cdot \delta {\boldsymbol{r}}{dV},\end{eqnarray} \tag{ 10 }$$

where ρ_e and ${\rho }_{p}$ are the density of the envelope and the perturbation, respectively, and ${{\boldsymbol{g}}}_{{\rm{L}}}$, ${{\boldsymbol{g}}}_{{\rm{e}}}$, and ${{\boldsymbol{g}}}_{p}$ are the gravitational fields due to the central Larson core, the envelope, and the perturbation itself. Their expressions are given in Paper II.

To compute the thermal energy, it is necessary to know the temperature. For the sake of simplicity, the calculations presented below assume that the gas remains isothermal. This assumption is discussed in Appendix D.

The kinetic energy in Equation (9) requires knowledge of the velocity dispersion within the perturbation. We assume that the usual scaling laws also apply here so that the velocity dispersion, δv, within a fluctuation of size δr is proportional to δr^0.5. Thus,

$$\begin{eqnarray}&&{E}_{\mathrm{kin}}=0.5\,{M}_{p}\delta {v}^{2}\simeq 0.5{\alpha }_{\mathrm{turb}}{M}_{p}{c}_{s}^{2}{{ \mathcal M }}^{2}(r)\times \displaystyle \frac{\delta r}{r},\end{eqnarray} \tag{ 11 }$$

where M_p is the mass of the perturbation and α_turb is a coefficient of the order of a few.

Finally, the effect of the external confining pressure at the perturbation boundary must also be considered in Equation (9). As discussed above (see Figure 6) and in the appendix, estimating this quantity is difficult. The numerical estimate shown in Figure 6 suggests that in the range M_L < M < 10 M_L we have the expressions stated by Equation (4).

The top panel of Figure 7 displays the critical density fluctuation, η_crit($\widetilde{r}$) = exp(δ_p,crit) − 1, and the perturbation relative size, that is, u = δr_p/r_p for various masses of the perturbation. The model parameters correspond to the ones measured for run 1: A ≃ 10 and ${ \mathcal M }\simeq 8$ (see Figure 5). As can be seen, typical unstable density fluctuations are about 20–40 times the local mean density. This high value stems from the high turbulence that is adopted here. As expected, bigger masses require density fluctuations that are less intense but bigger in size than smaller ones.

Zoom In Zoom Out Reset image size

3.2.4. Mean Number of Self-gravitating Fluctuations of Mass ≥ M_L

We now estimate the mean number of self-gravitating fluctuations of mass at least equal to M_L as a function of the distance, r, from the central fragment. Since the density threshold δ_c varies with the radial distance r, the probability of finding unstable fluctuations must be performed in concentric shells inside which this probability is uniform. The first step is to estimate the number of relevant fluctuations within the shells. Two types of fluctuations must be considered.

First, we have the ones whose mass is larger than M_L while their density, η_p,crit, and radius, δr_p, are the results of the virial equilibrium expressed by Equation (9). As the volume they occupy is δV = 4π/$3\delta {r}_{p}^{3}$, a shell of radius r and thickness dr contains up to 3(r/δr_p)³dr/r perturbations of this size. The number of unstable fluctuations in the shell is equal to the probability of finding a density fluctuation in the desired range of density times the number of fluctuations of this size. Within the shell, the mean number of self-gravitating fluctuations of mass ≥M_L is therefore given by

$$\begin{eqnarray}\begin{array}{rcl}{{dN}}_{1}(r) & = & {\int }_{{\delta }_{\min }}^{{\delta }_{\max }}\displaystyle \frac{3}{\sqrt{2\pi }{\sigma }_{0}}\exp \left[\displaystyle \frac{-{\left(\delta +{\sigma }_{0}^{2}/2\right)}^{2}}{2{\sigma }_{0}^{2}}\right]\\ & & \times \ {\left(\displaystyle \frac{r}{\delta {r}_{p}}\right)}^{3}\displaystyle \frac{{dr}}{r}d\delta ,\end{array}\end{eqnarray} \tag{ 12 }$$

where δ_max = δ_p,crit is the density corresponding to a fluctuation of mass M_L (Figure 7), and δ_min is the density for which δr_p/r = 0.5; that is, the size of the fluctuation is one-half the radius r that we adopt as an upper threshold. We have tested other values and found that they do not drastically affect the result.

The second type of fluctuation that must be considered is those whose mass is equal to M_L but the density is larger than δ_p,crit. The size of these fluctuations is given by $\delta {r}_{p,\mathrm{crit}}\times {\left(\tfrac{{\rho }_{p,\mathrm{crit}}}{{\rho }_{p}}\right)}^{1/3}$. The density term takes into account the fact that denser fluctuations have a smaller size because we integrate over fluctuations of mass M_L. The number of self-gravitating fluctuations with mass equal to M_L is thus

$$\begin{eqnarray}\begin{array}{rcl}{{dN}}_{2}(r) & = & {\displaystyle \int }_{{\delta }_{p,\mathrm{crit}}}^{\infty }\displaystyle \frac{3}{\sqrt{2\pi }{\sigma }_{0}}\exp \left[\displaystyle \frac{-{\left(\delta +{\sigma }_{0}^{2}/2\right)}^{2}}{2{\sigma }_{0}^{2}}\right]\\ & & \times \ {\left(\displaystyle \frac{r}{\delta {r}_{p,\mathrm{crit}}}\right)}^{3}\displaystyle \frac{{\rho }_{p}}{{\rho }_{p,\mathrm{crit}}}\displaystyle \frac{{dr}}{r}d\delta .\end{array}\end{eqnarray} \tag{ 13 }$$

Within a sphere of radius r, the number of self-gravitating fluctuations of mass larger than or equal to M_L is thus

$$\begin{eqnarray}&&N(r)={\int }_{0}^{r}({{dN}}_{1}+{{dN}}_{2}).\end{eqnarray} \tag{ 14 }$$

3.2.5. Mean Mass Calculation

To obtain the mass distribution, we proceed as follows.

The probability of finding n_f fragments inside the sphere of radius r_i = r_acc with at least a fragment in the shell i is

$$\begin{eqnarray}&&P({r}_{\mathrm{acc}},{n}_{f})=\sum _{{n}_{i}=1,{n}_{f};\sum {n}_{j}={n}_{f}}\prod _{j=0,i}{p}_{j}^{{n}_{j}},\end{eqnarray} \tag{ 15 }$$

where ${p}_{j}^{{n}_{j}}$ is the probability of finding n_j fragments in shell j. When dN is small, we simply have ${p}_{j}^{1}$ ≃ dN(r_j) and p_j⁰ ≃ 1 − ${p}_{j}^{1}$. The product in Equation (15) runs over all shells inside (including) shell i, the total number of fragments being $\sum {n}_{i}={n}_{f}$, and the sum is over the whole fragment distribution. For example, the probability of finding two fragments inside r₃ is given by P(r₂, 3) = ${p}_{3}^{2}{p}_{2}^{0}{p}_{1}^{0}$ + ${p}_{3}^{1}{p}_{2}^{1}{p}_{1}^{0}$ + ${p}_{3}^{1}{p}_{2}^{0}{p}_{1}^{1}$ . Note that because we use a large number of shells (typically 90–100, logarithmically spaced in r), the shells are all thin, and the probability of having more than one fragment within the same shell is negligible. We have verified that our result does not depend on the shell number.

The mass M_acc = M_env(r_acc) inside radius r_i = r_acc is given by Equation (6). In combination with Equation (15), this gives the overall mass distribution, dN_acc/dM = P_acc. All of the mass, M_acc, however, is not going to be accreted because some of the fragments lie in shells contained within shell i and thus can accrete a fraction of M_acc. Therefore, we compute the mass that is located at a radius smaller than the fragment position, with the assumption that a fragment located in shell j accretes a fraction 1/n_f of the mass located outside the sphere of radius r_j. That is to say, we assume that all fragments shield 1/n_f of the infalling gas. We call this mass (a fraction of M_acc) that has not been accreted by the other (shielding) fragments and then will eventually be accreted by the central seed the unshielded mass M_unsh. It is straightforward to see that

$$\begin{eqnarray}&&{M}_{\mathrm{unsh}}=\displaystyle \frac{{\sum }_{l=1,{n}_{f}}{m}_{l}}{{n}_{f}},\end{eqnarray} \tag{ 16 }$$

where m_l = M_env(r_l). Note that M_unsh depends on the fragment distribution, not only on the position of the most distant fragment. To get the mass distribution, dN_unsh/dM = P_unsh, we have thus computed the histogram of M_unsh weighted by the probability P(r_acc, n_f).

Clearly this mass is indicative, and in particular the assumption that all fragments shield a fraction 1/n_f of the mass, while reasonable for large samples, is likely to undergo large statistical fluctuations. Let us remember that the goal here is to estimate the peak position, that is, a preferred mass, and not the full mass spectrum.

3.2.6. Summary of Model Parameters

Table 1 summarizes the parameters of the model. We stress that their values are consequences of the physics, notably the dynamics that develop during the collapse. They are estimated through numerical simulations and are fluctuating by a factor of a few.

Table 1.  Summary of the Model Parameters

Name	Nature	Typical Value
ML	Mass of the first Larson core	0.02–0.03 M⊙
A	Density of the envelope	5–10
	Density fluctuations and Mach number	0.3–1.
αturb	Turbulent support	0–1
nf	Number of shielding fragments	2–5

Download table as:  ASCIITypeset image

We focus on model values that are directly inspired from the values measured in run 1. We adopt A = 10, ${ \mathcal M }=8$, and b = 0.6. For the external pressure, we adopt a fit taken from Section 2.7. Let us stress that those are typical values, but as seen from Figure 5, large deviations from these values have been inferred, implying that in reality one should sum over a large number of configurations. Another limit is that Equation (4) is obtained for cores that are already formed and developed, while the fluctuations we are identifying are more core progenitors.

3.3.1. Results for Value Parameters Inferred from the Simulations

The top panel of Figure 8 shows N (blue curve), dN₁ (red curve), and dN₂ (black curve) as a function of $\widetilde{r}$ for a series of logarithmic shells (see Equations (12), (13), and (14) for definitions). As can be seen, the fragments of type 2 dominate (from Figure 7, this is because they are smaller in size and therefore more numerous). The mean number of fragments, N, becomes equal to 1 for r/r_e,L ≃ 8−9. Therefore, qualitatively, assuming that all of the mass within this radius will be accreted onto the central fragment, we see that the final mass of the central fragment is expected to be ≃10 M_L.

Zoom In Zoom Out Reset image size

The corresponding distribution of the unshielded masses as stated by Equation (16) is portrayed in the bottom panel of Figure 8, where the mass distribution is shown for two, three, and four shielding fragments. As can be seen, it peaks at about 10 M_L and weakly depends on the number of shielding fragments, n_f.

Note that the validity of the distribution is not expected to extend over the whole mass spectrum and is more likely limited in the region of the peak. In particular, the high masses are likely determined by the reservoir distribution, as discussed in Paper I.

3.3.2. Comparison between Model and Simulation Results

First, we assume that the core is virialized:

$$\begin{eqnarray}&&\displaystyle \frac{3}{5}\displaystyle \frac{{{GM}}_{{\rm{L}}}{m}_{p}}{R}\simeq { \mathcal H }\times 3\,{k}_{B}T,\end{eqnarray} \tag{ 17 }$$

where m_p is the mass per particle, M_L is the mass of the Larson core, and T and R its temperature and radius. Here, ${ \mathcal H }$ is a parameter of the order of a few, which reflects that the hydrostatic equilibrium equations should be solved instead of mere virial equilibrium. We do so below, but it is very insightful to get a simple and explicit expression valid within a factor of a few. It will be convenient to express M_L as a function of the density, M_L = 4π/3ρR³, so

$$\begin{eqnarray}&&{M}_{{\rm{L}}}\simeq {{ \mathcal H }}^{3/2}\times {5}^{3/2}{\left(\displaystyle \frac{{k}_{B}T}{{{Gm}}_{p}}\right)}^{3/2}{\left(\displaystyle \frac{3}{4\pi \rho }\right)}^{1/2}.\end{eqnarray} \tag{ 18 }$$

The hydrostatic core ends when the H₂ molecule starts dissociating, and a second collapse occurs at about T = T_dis ≃ 1500 K. Equation (17) provides a relation between M_L and R. To estimate, say, the radius as a function of mass, we need to know at which density, ρ_dis, the second collapse starts or equivalently how the density depends on the temperature. To achieve this, we need to know the entropy of the fluid. Using the well-known fact (that we verify below) that the cooling time of the core is several times larger than its dynamical time, we can assume that above a density, ρ_ad, to be determined, the gas is essentially adiabatic. At the beginning of this adiabatic phase, the rotational levels of H₂ are not excited, and the adiabatic index is γ = 5/3. Then when the rotational levels of H₂ get excited, γ drops, and at about T_ex ≃ 150 K, we have γ = 7/5. Thus, we have

$$\begin{eqnarray}&&{\rho }_{\mathrm{dis}}\simeq {\rho }_{\mathrm{ex}}{\left(\displaystyle \frac{{T}_{\mathrm{dis}}}{{T}_{\mathrm{ex}}}\right)}^{5/2}\simeq {\rho }_{\mathrm{ad}}{\left(\displaystyle \frac{{T}_{\mathrm{ex}}}{{T}_{0}}\right)}^{3/2}{\left(\displaystyle \frac{{T}_{\mathrm{dis}}}{{T}_{\mathrm{ex}}}\right)}^{5/2},\end{eqnarray} \tag{ 19 }$$

where T₀ is the mean cloud temperature in the isothermal regime. We now need to estimate the density ρ_ad at which the gas becomes adiabatic. This estimate has been done by several authors (Low & Lynden-Bell 1976; Rees 1976; Masunaga & Inutsuka 1999) and consist in estimating the cooling and freefall time. To estimate ρ_ad, we must compare the cooling time and the freefall time (Masunaga & Inutsuka 1999). To estimate the former, we simply compute the ratio of the thermal energy contained in a Jeans mass, defined as a sphere of radius ${\lambda }_{J}/2=\sqrt{\pi }{C}_{s}/2/\sqrt{G{\rho }_{\mathrm{ad}}}$, to the energy radiated per units of time, ${ \mathcal L }$:

$$\begin{eqnarray}&&{E}_{\mathrm{therm}}=\displaystyle \frac{{M}_{J}}{{m}_{p}}\times \displaystyle \frac{3}{2}{k}_{B}{T}_{0},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{ \mathcal L }\simeq 4\pi {\left(\displaystyle \frac{{\lambda }_{J}}{2}\right)}^{2}\displaystyle \frac{4{{acT}}_{0}^{3}}{3{\kappa }_{{\rm{R}}}({T}_{0}){\rho }_{\mathrm{ad}}}\displaystyle \frac{{T}_{0}}{{\lambda }_{J}/2}.\end{eqnarray} \tag{ 21 }$$

This last expression is simply the radiative flux through a sphere of diameter λ_J in which the gradient ∂_rT has been replaced by T/r, κ_R is the Rosseland opacity, c is the speed of light, and a = 4σ/c, σ being the Stefan–Boltzmann constant. The cooling time is simply estimated as

$$\begin{eqnarray}&&{\tau }_{\mathrm{cool}}=\displaystyle \frac{{E}_{\mathrm{therm}}}{{ \mathcal L }},\end{eqnarray} \tag{ 22 }$$

while the freefall time is

$$\begin{eqnarray}&&{\tau }_{\mathrm{ff}}=\sqrt{\displaystyle \frac{3\pi }{32G{\rho }_{\mathrm{ad}}}}.\end{eqnarray} \tag{ 23 }$$

Requiring that τ_cool ≃ τ_ff, we get

$$\begin{eqnarray}&&{\rho }_{\mathrm{ad}}={\left(\displaystyle \frac{32}{3\pi }{{Gm}}_{p}^{4}{k}_{B}^{-4}{a}^{2}{c}^{2}\right)}^{1/3}{T}_{0}^{4/3}{\kappa }_{{\rm{R}}}^{-2/3}.\end{eqnarray} \tag{ 24 }$$

Combining Equations (18), (19), and (24), we get an expression for M_L:

$$\begin{eqnarray}\begin{array}{rcl}{M}_{{\rm{L}}} & \simeq & {{ \mathcal H }}^{3/2}\displaystyle \frac{{5}^{3/2}{3}^{2/3}}{{2}^{11/6}{\pi }^{1/3}}{\left({G}^{10}{m}_{p}^{13}{k}_{B}^{-13}{a}^{2}{c}^{2}\right)}^{-1/6}\\ & & \times \ {T}_{\mathrm{dis}}^{1/4}{T}_{\mathrm{ex}}^{1/2}{T}_{0}^{1/12}{\kappa }_{{\rm{R}}}{\left({T}_{0}\right)}^{1/3}\\ & = & \ {{ \mathcal H }}^{3/2}\times 0.08\,{M}_{\odot }{\left(\displaystyle \frac{{T}_{0}}{10{\rm{K}}}\right)}^{1/12}{\left(\displaystyle \frac{{\kappa }_{{\rm{R}}}({T}_{0})}{0.1{\mathrm{cm}}^{2}{{\rm{g}}}^{-1}}\right)}^{1/3}.\end{array}\end{eqnarray} \tag{ 25 }$$

The estimated value of M_L from numerical simulations is about 0.03 M_⊙ (Vaytet & Haugbølle 2017), so ${ \mathcal H }\simeq 0.5$ leads to the correct value. Indeed, integrating the proper hydrostatic equilibrium (see Equation (7) of Paper II) starting from ρ_dis and using the piecewise polytropic equation of state discussed in the text, we obtain a mass equal to 0.03 M_⊙. The typical value of κ_R is taken from Semenov et al. (2003) at 10 K. We recall that the goal of Equation (25) is to make explicit the various dependencies.

Equation (25) makes two important predictions. First, it does not depend significantly on T₀. Second, it depends on the opacity κ_R, and therefore on the metallicity, Z, to the power of one-third, which is rather shallow. This implies that one would not expect an important dependence of the peak of the IMF on the metallicity. Typically, changing Z by a factor of 10 implies a shift in the M_L of about 2 if T₀ is unchanged. There is, however, likely another effect that is due to the fact that at low metallicity the temperature T₀ is higher because the cooling is lower. Yet the opacity increases with the temperature. This effect is particularly visible in the simulations presented by Bate (2014), where four runs with four different metallicities are presented with Z going from 0.01 to 3. The opacities are proportional to Z, making a large variation of κ_R. However, the author found that the background temperature varies from 10 K to 70 K at Z = 0.01. Looking at his Figure 1, we see that the opacity at 10 K and Z = 1 and the one at 70 K and Z = 0.01 are nearly identical and equal 0.02–0.03 g cm⁻². According to the dependence stated by Equation (25) and after doing the proper hydrostatic calculation, this leads to an estimate of M_L ≃ 0.02 M_⊙ and therefore to a peak of the IMF around 10 times this value (as argued in Paper II and below). This is clearly not incompatible with Figure 6 of Bate (2014), where we see that in all cases there is a peak located at a few 0.1 M_⊙.

A central aspect of the physics of the first hydrostatic core is that its lifetime is significantly longer than the dynamical time. Here we provide an estimate and a discussion of the implication.

To get the cooling time, we again use the thermal energy and the radiative flux. However, we now take the typical values that correspond to the inner part of the core, say within 1 au or so.

$$\begin{eqnarray}&&{E}_{\mathrm{therm}}\simeq \displaystyle \frac{{M}_{{\rm{L}}}}{{m}_{p}}\times \displaystyle \frac{5}{2}{k}_{B}{T}_{\mathrm{dis}},\end{eqnarray} \tag{ 26 }$$

where the factor 5 is due to the H₂ molecule being diatomic because in this regime rotational levels are excited (see, e.g., Saumon et al. 1995 for a complete treatment):

$$\begin{eqnarray}&&{ \mathcal L }\simeq 4\pi {R}_{{\rm{L}}}^{2}\displaystyle \frac{4{{acT}}_{\mathrm{dis}}^{3}}{3{\kappa }_{{\rm{R}}}({T}_{\mathrm{dis}}){\rho }_{\mathrm{dis}}}\displaystyle \frac{{T}_{\mathrm{dis}}}{{R}_{{\rm{L}}}}.\end{eqnarray} \tag{ 27 }$$

Using again M_L = 4π/3ρR³ to infer the density in Equations (27) and (17) to infer RT_dis, we get

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{\mathrm{cool}} & = & \displaystyle \frac{{3}^{2}{5}^{5}}{{2}^{7}{\pi }^{2}}\displaystyle \frac{{k}_{B}^{5}{T}_{\mathrm{dis}}}{{{acG}}^{4}\,{m}_{p}^{5}}\displaystyle \frac{{\kappa }_{{\rm{R}}}({T}_{\mathrm{dis}})}{{M}_{{\rm{L}}}^{2}}\\ & \simeq & \ 5.6\times {10}^{4}\,\mathrm{yr}\displaystyle \frac{{\kappa }_{{\rm{R}}}({T}_{\mathrm{dis}})}{10\,{\rm{g}}\ {\mathrm{cm}}^{-2}}{\left(\displaystyle \frac{{M}_{{\rm{L}}}}{0.03{M}_{\odot }}\right)}^{-2},\end{array}\end{eqnarray} \tag{ 28 }$$

where κ_R(T_dis) = 10 cm² g⁻¹ is from Semenov et al. (2003). By contrast, the freefall time at a density ρ_dis is about 1 yr. More generally, a freefall time of 5 × 10⁴ yr corresponds to a density of roughly 10⁶ cm⁻³, which is far below the density of the first hydrostatic core.

The other relevant time is the accretion one, τ_acc, that is to say the time necessary to accumulate a mass equal to M_L,

$$\begin{eqnarray}&&{\tau }_{\mathrm{acc}}=\displaystyle \frac{{M}_{{\rm{L}}}}{\tfrac{{dM}}{{dt}}}=300\,\mathrm{yr}\displaystyle \frac{{M}_{{\rm{L}}}}{0.03\,{M}_{\odot }}{\left(\displaystyle \frac{\tfrac{{dM}}{{dt}}}{{10}^{-4}{M}_{\odot }{\mathrm{yr}}^{-1}}\right)}^{-1},\end{eqnarray} \tag{ 29 }$$

which is very close to the estimated lifetime of the Larson core (Vaytet & Haugbølle 2017) of 100–1000 yr because the typical accretion rate inferred in collapse calculations is a few 10⁻⁵ M_⊙ yr⁻¹ (Foster & Chevalier 1993; Hennebelle et al. 2003; Vaytet & Haugbølle 2017).

This implies that the dynamics of the gas that enters this core is halted. The core stays a period of time that is long with respect to the dynamical timescale, which itself is typical of a self-gravitating object of this density. Moreover, the cooling time is even longer for an object with a smaller mass (τ_cool ∝ M⁻²). Therefore, although in practice objects with a mass smaller than M_L could cool down and eventually form a star or a brown dwarf, in practice it is very difficult if the object is not isolated. There is plenty of time for the object to accrete and acquire a mass of at least M_L, which would trigger H₂ dissociation and a second collapse.

Equations (28) and (29) predict a change of behavior of the lifetime of the first hydrostatic core at a mass M ≃ M_L. This transition is clearly seen in the recent work by Stamer & Inutsuka (2018), who presented a series of 1D calculations (see their Figure 6).

We stress in particular that the minimum mass for fragmentation (Low & Lynden-Bell 1976; Rees 1976; Masunaga & Inutsuka 1999), estimated to be the Jeans mass at the density ρ_ad (as stated from Equation (24)), is not clear because the fragments have typically a mass of about 10⁻³ M_⊙, which is too low to induce the dissociation of H₂ and form a protostar. In most cases, these fragments will simply grow in mass and form a first hydrostatic core or will merge with an existing one or a protostar. Only if the fragments are very isolated will they have the opportunity to cool.