<figure><img src="image_1.png" /></figure>

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

Patrick Hennebelle¹, Yueh-Ning Lee²,³, and Gilles Chabrier⁴,⁵

¹Laboratoire AIM, Paris-Saclay CEA/IRFU/SAP—CNRS—Université Paris Diderot, F-91191 Gif-sur-Yvette Cedex, France; patrick.hennebelle@cea.fr

²Department of Earth Sciences, National Taiwan Normal University, 88 Sec. 4, Ting-Chou Road, Taipei 116, Taiwan

³Institut de Physique du Globe de Paris Sorbonne Paris Cité, Université Paris Diderot UMR 7154 CNRS, F-75005 Paris, France

⁴École Normale Supérieure de Lyon CRAL, UMR CNRS 5574, F-69364, Lyon Cedex 07, France

⁵School of Physics, University of Exeter, Exeter, EX4 4QL, UK

Received 2019 May 25; revised 2019 July 26; accepted 2019 August 19; published 2019 September 27

## Abstract

The stellar initial mass function plays a critical role in the history of our universe. We propose a theory that is based solely on local processes, namely the dust opacity limit, the tidal forces, and the properties of the collapsing gas envelope. The idea is that the final mass of the central object is determined by the location of the nearest fragments, which accrete the gas located farther away, preventing it from falling onto the central object. To estimate the relevant statistics in the neighborhood of an accreting protostar, we perform high-resolution numerical simulations. We also use these simulations to further test the idea that fragmentation in the vicinity of an existing protostar is a determinant in setting the peak of the stellar spectrum. We develop an analytical model, which is based on a statistical counting of the turbulent density fluctuations, generated during the collapse, that have a mass at least equal to the mass of the first hydrostatic core, and sufficiently important to supersede tidal and pressure forces to be self-gravitating. The analytical mass function presents a peak located at roughly 10 times the mass of the first hydrostatic core, in good agreement with the numerical simulations. Since the physical processes involved are all local, occurring at scales of a few 100 au or below, and do not depend on the gas distribution at large scale and global properties such as the mean Jeans mass, the mass spectrum is expected to be relatively universal.

**Key words:** gravitation – hydrodynamics – ISM: clouds – methods: numerical – stars: formation – turbulence

## 1. Introduction

The question of the origin of the mass distribution of stars, the initial mass function (IMF), is a long-standing issue (e.g., Salpeter 1955; Kroupa 2001; Chabrier 2003; Bastian et al. 2010; Offner et al. 2014), and several authors have attempted to provide explanations either using analytical modeling of the gas fragmentation (Padoan et al. 1997; Inutsuka 2001; Hennebelle & Chabrier 2008), numerical simulations of a fragmenting cloud (Bonnell et al. 2011; Girichidis et al. 2011; Ballesteros-Paredes et al. 2015), or a statistical description of accretion (Basu & Jones 2004; Maschberger et al. 2014; Basu et al. 2015). In general, these models have been reasonably successful in reproducing the high-mass tail of the IMF and obtained power-law mass spectra with a slope close to the Salpeter value. What is truly intriguing, however, is the lack of variations of the IMF given the large amount of variations in the star-forming environments and particularly the peak of the distribution, around $0.3 M_\odot$, which constitutes a characteristic mass for the stars. Traditional explanations based on turbulence and Jeans mass (Padoan et al. 1997; Hennebelle & Chabrier 2008; Hopkins 2012), gas cooling and variation of the equation of state (Jappsen et al. 2005), or radiative feedback (Bate 2009; Guszejnov et al. 2016; Krumholz et al. 2016) are all facing difficulty in explaining the lack of variations of the peak position and require the scale dependence of some quantities, such as the Mach number and the density or the radiative feedback and the opacities, to cancel out.

Recently, Lee & Hennebelle (2018a, 2018b; hereafter Papers I and II) ran a series of simulations of a collapsing $1000 M_\odot$ clump, changing the initial conditions, namely the initial density and Mach numbers, as well as the magnetic field (Lee & Hennebelle 2019), by orders of magnitude and found

that the peak of the stellar mass spectrum is insensitive to these variations in a large range of parameters. Moreover, they found it to be about 5–10 times the mass, $M_\mathrm{L}$, of the first hydrostatic Larson core (FHSC; Larson 1969; Masunaga et al. 1998; Vaytet et al. 2013; Vaytet & Haugbølle 2017), which is the hydrostatic core that forms when the dust becomes opaque to its own radiation. To establish this, the effective equation of state, which describes the transition from the isothermal to the adiabatic phase, has been varied, and it has been found that the peak of the IMF occurs at a mass that is proportional to $M_\mathrm{L}$. The advantage of this scheme is that it naturally explains the weak variations of the IMF, at least for its peak. The latter is a mere consequence of the physics of the FHSC and its surrounding collapsing envelope being nearly independent of the large-scale conditions.

In Paper II, it was suggested that the origin of the factor of 5–10 is due to further accretion from the envelope. Eventually it is limited by the presence of nearby fragments that shield, and compete for, further accretion. The presence of the nearby fragments is regulated on one hand by the tidal forces induced by the Larson core and the gas envelope, and on the other hand by the density fluctuations arising within its surrounding accreting envelope. Qualitatively, the mechanism we propose can be summarized as follows. Given that the FHSC is adiabatic, it is clear that fragmentation stops at least up to the point where the second collapse restarts (during which in principle tight binaries could possibly form; see, e.g., Bate 1998; Machida et al. 2008). Since $M_\mathrm{L}$ is the minimum mass to collapse, it seems natural that the peak of the IMF should be larger than this value. The infalling gas accumulates and piles up until a mass of at least $M_\mathrm{L}$ has been accumulated, and it appears unavoidable that further accretion will proceed.
Yet it seems plausible that the hydrostatic core is surrounded by a mass of several times $M_L$ and that one needs to go beyond this mass to find a new hydrostatic core.

This mechanism, which is largely based on (1) a drastic change of the effective equation of state and (2) the influence of other objects, puts together two ideas that have been proposed earlier. First is the role played by thermodynamics (Jappsen et al. 2005), although it emphasizes that the important change of the effective adiabatic equation-of-state exponent is between below and above 4/3 (Jappsen et al. 2005; consider the change between 0.7 and 1.1). Second is the role played by the fragmentation-induced starvation envisioned by Peters et al. (2011), in which the fragmentation around massive stars limits their growth. It presents also similarities to the competitive accretion scenario advocated by Bonnell et al. (2001) and Bonnell & Bate (2006), though originally this scenario has been proposed to explain the slope of the IMF and the formation of the most massive stars rather than the peak of the stellar distribution.

In the present paper, we propose a model that attempts to quantify this picture and to estimate the mass at which the peak of the mass spectrum occurs. The model consists of counting the density fluctuations in the vicinity of an existing protostar that are sufficiently high to be gravitationally unstable, given a protostar mass that is equal to at least $M_L$. From the probability of finding a certain number of fragments, we can infer the typical mass that is unshielded and therefore available for accretion into the central object. To constrain and test the model, we have performed new simulations (1) to check that the peak is still located at about 5–10 × $M_L$, where $M_L$ is the mass of the FHSC, when different initial conditions are explored, (2) to measure the statistics of the flow in the vicinity of existing sink particles, and (3) to perform a direct test of the importance of the fragmentation at a few hundreds of astronomical units from an existing protostar in setting the peak of the stellar mass function.

The plan of the paper is as follows. In the second section we present the numerical simulations, which use two different setups. We perform several statistical estimates, including numbers of sink neighbors, Mach numbers and density distributions, and core properties. In the third section, we detail the assumptions of the model and the mathematical formalism. We also compare the predictions of the model with the simulation results. The fourth section concludes the paper. It is complemented by a long appendix in which the physical properties of the FHSC are discussed and simple orders of magnitude are obtained, as well as a discussion on the surface term of the virial theorem.

## 2. Numerical Simulations

### 2.1. Numerical Methods and Setup

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

In order to verify that our results are not too dependent on a particular choice of initial conditions, we have performed two types of runs.

Run 1 considers a spherical cloud with uniform density initially, in which turbulence has been added (the fluctuations have random phases and a velocity power spectrum that follows the classical power law with a $-11/3$ exponent) and is freely decaying. The initial conditions consist of a cloud of $10^3 M_\odot$, which initially has a radius of about 0.1 pc corresponding to an initial density of $5 \times 10^6 \text{ cm}^{-3}$. The initial Mach number is about 10. The number of AMR levels being used is 14, leading to an effective resolution of about 4 au. As explained below, this run is complemented by two more runs that are identical except that the formation of new sink particles is prevented 140 au (run 1a) and 280 au (run 1b) around already-existing sink particles.

Run 2 consists of a periodic box in which turbulent forcing is applied for wavenumbers $k = 1–3$. The size of the box is 0.5 pc, the total mass is also equal to $10^3 M_\odot$, and this corresponds to a density of about $10^5 \text{ cm}^{-3}$. The turbulence is driven for 1 Myr before gravity is applied. By this time, the velocity dispersion is about 2 km s⁻¹, which corresponds to a Mach number of about 10. The number of AMR levels being used is 15, leading to an effective resolution of about 4 au. Note that we have also performed a run with a density 8 times smaller (and a computing box 2 times larger, i.e., of 1 pc) and inferred very similar results. Therefore, this run is not presented here.

### 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:

$$T = T_0 \left\{ 1 + \frac{(n/n_{\text{ad}})^{(\gamma_1 - 1)}}{1 + (n/n_{\text{ad,2}})^{(\gamma_1 - \gamma_2)}} \right\}, \quad (1)$$

where $T_0 = 10$ K, $n_{\text{ad}} = 10^{10} \text{ cm}^{-3}$, $n_{\text{ad,2}} = 3 \times 10^{11} \text{ cm}^{-3}$, $\gamma_1 = 5/3$, and $\gamma_2 = 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_\odot$. 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^{11} \text{ cm}^{-3}$. The sinks are introduced at a density $n \ge 10^{12} \text{ cm}^{-3}$. 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
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Column density of run 1 (left) and run 2 (right). The top panels show most of the whole cloud, while the bottom panels represent zoomed-in regions. The black dots represent the sink particles aiming to model the stars.</figcaption></figure>

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.

## 2.2. Qualitative Description

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.

## 2.3. Sink Mass Distributions

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 $\simeq 0.2 M_\odot$ 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_\odot$. At higher mass, the distribution again is broadly compatible with a power law $M^{-0.75}$.

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_\odot$ and $170 M_\odot$ has been accreted. In Paper I, we were able to go up to the point where $300 M_\odot$ 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.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Mass spectrum at various times for run 1 (left) and run 2 (right). The total mass accreted by the sink particles (expressed in solar mass) is indicated in the legend. The distribution peaks at about 0.2 $M_\odot$ for run 1 and 0.1 $M_\odot$ for run 2.</figcaption></figure>

<figure><img src="image_7.png" /></figure>

## 2.4. Sink Neighbor Distributions

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_{\text{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_{\text{sys}} = n_{\text{neigh}} + 1$. The number $N_{\text{sys}}$ of systems with $n_{\text{sys}}$ objects is related to the number $N_{\text{neigh}}$ of sinks with $n_{\text{neigh}}$ numbers through $N_{\text{sys}} = N_{\text{neigh}}/(n_{\text{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_{\text{sys}}$, solid lines) and sink numbers ($N_{\text{neigh}}$, dashed lines) of various ages with (up to) $n_{\text{neigh}} = n_{\text{sys}} - 1$ neighbors at a distance below 200 au. As can be seen at early time ($t \le 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_{\text{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_{\text{neigh}} \ge 10$, but the corresponding system number is small. This corresponds to a restricted number of regions where the clustering is particularly important.

The middle panels display the cumulative histogram of the number of systems with $n_{\text{neigh}} = n_{\text{sys}} - 1$ objects for several distances and at time $10^2$ yr (dotted line) and at time $10^3$ 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_{\text{sys}} \rangle$ as a function of distance. The red curves are for $t \le 10^2$ yr, while blue ones are for $t \le 10^3$ 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_{\text{sys}} \rangle$ below 100 au is $\le 2$ at $t = 100$ yr and $\simeq 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.

## 2.5. A Simple Numerical Experiment

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_\odot$, while dashed ones are for 50 $M_\odot$. 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
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Neighbor statistics. Left is for run 1 and right for run 2. The top panels give the cumulative histogram of the number of systems, $N_{sys}$, (solid lines) of size 200 au that contain $n_{sys} = n_{neigh} + 1$ objects. The dashed lines show the number of sinks ($N_{neigh} = N_{sys}n_{sys}$) having $n_{neigh}$ neighbors located at a distance of less than 200 au at various times after their formation. The middle panels provide the cumulative histogram of the number of systems possessing $n_{neigh} + 1$ objects at various distances and at time $t = 10^2$ yr (dotted line) and at time $t = 10^3$ yr (solid line). The bottom panels display the mean system number, $n_{sys}$, as a function of distance. Red curves are for times shorter than $10^2$ yr, while blue ones are for $t < 10^3$ yr. The black lines show the result of the analytical model inferred in Section 3 for various Mach numbers (see text in Section 3.3). Good agreement between the model and the simulations for the highest Mach numbers measured is obtained. When the external pressure terms are not taken into account, the number of objects is below the numbers inferred in the simulations.</figcaption></figure>

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.

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
<figure><img src="image_9.png" /><figcaption>**Figure 4.** Mass spectra for run 1, run 1a (sinks are forbidden to form below a distance of 140 au from an existing one), and run 1b (sinks are forbidden to form below a distance of 240 au from an existing sink). Solid lines are for a total sink mass of 150 $M_\odot$, while dashed lines stand for 50 $M_\odot$. Clearly, the sink distribution shifts to larger masses from run 1 to run 1a and run 1b. There are altogether 2 times fewer sinks in run 1b than in run 1, and their masses are about 2 times larger. This clearly demonstrates that fragmentation within the inner hundreds of astronomical units has a drastic effect on the sink/stellar mass distribution.</figcaption></figure>

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.

## 2.6. Sink Environment

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} = \rho_{SIS}/m_p = c_s^2/(2\pi Gr^2)/m_p$. The mean density is proportional to $r^{-2}$ and about 10 and 7 times $n_{SIS}$ for run 1 and run 2, respectively.

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 = \ln(1 + b^2 M^2)$ for $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 $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 $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 $\rho_{SIS}$) are themselves collapsing.

## 2.7. Core Statistics

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^8$ cm$^{-3}$. 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.

The top panel displays the core mass spectrum. Interestingly, it peaks at about $(1-2) \times 10^{-2} M_\odot$, 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

$$E_g + 2E_{th} + 2E_{kin} - V_{P_{ext}} - V_{P_{ram}} = 0, \quad (2)$$

where the various terms are the gravitational energy and the thermal and kinetic energies, while $V_{P_{ext}}$ and $V_{P_{ram}}$ stand for the integrated thermal and ram pressure in the outskirt of the
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Gas statistics around the sinks younger than 1000 yr. Left: run 1; right: run 2. The top panels show the mean density (red line) and the mean standard deviation (blue shaded area). For reference, the blue line shows the density of the singular isothermal sphere, $n_{SIS} = \rho_{SIS}/m_p$. We see that on average the density is about 10 × $\rho_{SIS}$ for run 1 and about 7 × $\rho_{SIS}$ for run 2. The middle panels display the mean Mach numbers, which below 300 au are found to be about 5–10 for run 1 and 4–5 for run 2. The bottom panels portray the distribution of the density fluctuations (i.e., with respect to the mean density shown on the top panel) within a sphere of radius 800 au (upper curve) and 250 au (lower curve). For reference, the red curves represent a log-normal distribution with a width $\sigma^2 = \ln(1 + b^2 M^2)$ with $b = 0.8$ and $M = 6$.</figcaption></figure>

cores. As detailed in Appendix C, their exact expressions are

$$ 
\begin{aligned}
VP_{\text{ext}} &= \oint P_{\text{th}} \mathbf{r} \cdot d\mathbf{S}, \\
VP_{\text{ram}} &= \oint \rho (\mathbf{v} \cdot \mathbf{r})(\mathbf{v} \cdot d\mathbf{S}).
\end{aligned} 
\quad (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
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 6.** Core analysis and statistics for run 1. The top left panel displays the core mass spectrum. The top right panel shows the ratio of confining to support terms that appear in the virial theorem, showing that cores are remarkably virialized. The middle left and middle right panels respectively portray the ratio of thermal and kinetic surfaces of bulk terms. As can be seen, surface terms, while negligible for cores of mass larger than 0.1 $M_\odot$, are comparable to the bulk values for cores of mass below $M_\mathrm{L}$. The bottom right and middle panels show the same quantities but plotted as a function of the Mach numbers instead.</figcaption></figure>

$VP_{\mathrm{ext}}/2E_{\mathrm{th}}$ and $VP_{\mathrm{ram}}/2E_{\mathrm{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_{\mathrm{grav}} + VP_{\mathrm{tot}})/(2E_{\mathrm{kin}} + 2E_{\mathrm{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_{\mathrm{ext}}/2E_{\mathrm{th}}$ and $VP_{\mathrm{ram}}/2E_{\mathrm{kin}}$ per bin of core mass. It is found that both of them are almost zero for $M \ge 0.1 M_\odot$. For $M < 10^{-2} M_\odot \simeq M_\mathrm{L}$, it is found that $VP_{\mathrm{ext}}/2E_{\mathrm{th}} \ge 0.5$ while $VP_{\mathrm{ram}} \simeq 2E_{\mathrm{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{aligned}
VP_{\text{ext}} &\approx E_{\text{th}} \times (0.6-0.1M), \\
VP_{\text{ram}} &\approx E_{\text{kin}} \times (1.3-0.18M).
\end{aligned}
\quad (4) $$

### 3. Relation between the Mass of the First Hydrostatic Core and the Mass of the Stars: Analytical Model

#### 3.1. Model Assumptions and Definitions

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^{-2} M_\odot$ (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^{-2}$ 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 $\ge 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_\odot$ 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 \ge 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. Mathematical Description

##### 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):

$$ \rho_{\text{env}}(r) = \frac{A c_s^2}{2\pi G r^2}, \quad (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 \sim 5–10$. For convenience, we define $r_{e,L} = M_L = 2A(c_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} \approx 20$ au. This yields

$$ M_{\text{env}}(r) = \int \rho_{\text{env}}(r) 4\pi r^2 dr = 2A \frac{c_s^2}{G} r = M_L \tilde{r}, \quad (6) $$

where $\tilde{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 $\delta_\rho$ on top of the mean density $\rho_{\text{env}}(r)$, given, in reasonable approximation, by a log-normal distribution

$$ P(\tilde{r}, \delta_{\rho}(\tilde{r})) = \frac{1}{\sqrt{2\pi}\sigma(\tilde{r})} \exp\left[\frac{-(\delta_{\rho}(\tilde{r}) + \sigma^2(\tilde{r})/2)^2}{2\sigma^2(\tilde{r})}\right], \quad (7) $$
with $\delta_p(\tilde{r}) = \ln(\rho(\tilde{r})/\rho_{\text{env}}(\tilde{r})) = \ln(1 + \eta(\tilde{r}))$, where $\eta(\tilde{r}) = \rho(\tilde{r})/\rho_{\text{env}}(\tilde{r}) - 1$, and $\sigma$ defines the width of the density PDF. The density fluctuation is related to the local turbulent Mach number as $\sigma^2(\tilde{r}) = \ln(1 + b^2 M(\tilde{r})^2)$ (e.g., Hennebelle & Falgarone 2012). To estimate the local Mach number $M(\tilde{r})$, we assume that the energy in turbulent motions is a fraction $\epsilon \le 1$ of the gravitational energy:

$$M(\tilde{r}) = \frac{1}{c_s} \sqrt{\epsilon \frac{G[M_L + M_{\text{env}}(r)]}{r}} = \sqrt{2A \left(1 + \frac{1}{\tilde{r}}\right) \epsilon}. \quad (8)$$

As seen from Figure 5, which displays the PDF of the density fluctuations, $\delta\rho/\rho$, around sink particles, the log-normal PDF with $b = 0.8$ and $M \simeq 6$ is a good approximation. This corresponds to $\epsilon \simeq 1$. Note that the Mach number distribution seen in Figure 5 shows a tendency for $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, $\eta_{\text{crit}}(\tilde{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 $\eta_p = \rho_p/\rho_e - 1$, we get

$$E_{\text{vir}}(r_p, \delta r_p, \eta_p) = E_g(r_p, \delta r_p, \eta_p) + 2E_{\text{ther}} + 2E_{\text{kin}} - V P_{\text{ext}} - V P_{\text{ram}}. \quad (9)$$

In this expression, the gravitational energy is given by

$$E_g(r_p, \delta r_p, \eta_p) = \int_{V_p} (\rho_e + \rho_p) (g_L + g_e + g_p) \cdot \delta r dV, \quad (10)$$

where $\rho_e$ and $\rho_p$ are the density of the envelope and the perturbation, respectively, and $g_L$, $g_e$, and $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, $\delta v$, within a fluctuation of size $\delta r$ is proportional to

<figure><img src="image_17.png" /><figcaption>Figure 7. Critical density fluctuation amplitude, $\eta_{\text{crit}} = \exp(\delta_{p,\text{crit}}) - 1$ (top panel), and relative perturbation radius, $u = \delta r_p / r_p$ (bottom panel), for the run 1 case (i.e., $A = 10$, $M = 8$) and various fragment masses. The surface terms (see Equation (3)) are given by Equation (4). Typically, only density fluctuations 20–40 times above the local mean density are unstable.</figcaption></figure>

$\delta r^{0.5}$. Thus,

$$E_{\text{kin}} = 0.5 M_p \delta v^2 \simeq 0.5 \alpha_{\text{turb}} M_p c_s^2 M^2(r) \times \frac{\delta r}{r}, \quad (11)$$

where $M_p$ is the mass of the perturbation and $\alpha_{\text{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, $\eta_{\text{crit}}(\tilde{r}) = \exp(\delta_{p,\text{crit}}) - 1$, and the perturbation relative size, that is, $u = \delta r_p / r_p$ for various masses of the perturbation. The model parameters correspond to the ones measured for run 1: $A \simeq 10$ and $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.

### 3.2.4. Mean Number of Self-gravitating Fluctuations of Mass $\ge 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 $\delta_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, $\eta_{p,crit}$, and radius, $\delta r_p$, are the results of the virial equilibrium expressed by Equation (9). As the volume they occupy is $\delta V = 4\pi/3\delta r_p^3$, a shell of radius $r$ and thickness $dr$ contains up to $3(r/\delta r_p)^3 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 $\ge M_L$ is therefore given by

$$dN_1(r) = \int_{\delta_{\min}}^{\delta_{\max}} \frac{3}{\sqrt{2\pi}\sigma_0} \exp\left[\frac{-(\delta + \sigma_0^2/2)^2}{2\sigma_0^2}\right] \times \left(\frac{r}{\delta r_p}\right)^3 \frac{dr}{r} d\delta, \quad (12)$$

where $\delta_{\max} = \delta_{p,crit}$ is the density corresponding to a fluctuation of mass $M_L$ (Figure 7), and $\delta_{\min}$ is the density for which $\delta 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 $\delta_{p,crit}$. The size of these fluctuations is given by $\delta r_{p,crit} \times \left(\frac{\rho_{p,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

$$dN_2(r) = \int_{\delta_{p,crit}}^{\infty} \frac{3}{\sqrt{2\pi}\sigma_0} \exp\left[\frac{-(\delta + \sigma_0^2/2)^2}{2\sigma_0^2}\right] \times \left(\frac{r}{\delta r_{p,crit}}\right)^3 \frac{\rho_p}{\rho_{p,crit}} \frac{dr}{r} d\delta. \quad (13)$$

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

$$N(r) = \int_{0}^{r} (dN_{1} + dN_{2}). \quad (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

$$P(r_{\text{acc}}, n_f) = \sum_{n_i=1, n_f} \prod_{\sum n_j=n_f} p_j^{n_j}, \quad (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 \simeq dN(r_j)$ and $p_j^0 \simeq 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_3$ is given by $P(r_2, 3) = p_3^2 p_2^0 p_1^0 + p_3^1 p_2^1 p_1^0 + p_3^1 p_2^1 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

<table><caption>**Table 1**
Summary of the Model Parameters</caption><thead><tr><th>Name</th><th>Nature</th><th>Typical Value</th></tr></thead><tbody><tr><td><i>M</i><sub>L</sub></td><td>Mass of the first Larson core</td><td>0.02–0.03 <i>M</i><sub>☉</sub></td></tr><tr><td><i>A</i></td><td>Density of the envelope</td><td>5–10</td></tr><tr><td><i>ε</i></td><td>Density fluctuations and Mach number</td><td>0.3–1.</td></tr><tr><td><i>α</i><sub>turb</sub></td><td>Turbulent support</td><td>0–1</td></tr><tr><td><i>n</i><sub>f</sub></td><td>Number of shielding fragments</td><td>2–5</td></tr></tbody></table>

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_{\text{env}}(r_{\text{acc}})$ inside radius $r_i = r_{\text{acc}}$ is given by Equation (6). In combination with Equation (15), this gives the overall mass distribution, $dN_{\text{acc}}/dM = P_{\text{acc}}$. All of the mass, $M_{\text{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_{\text{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_{\text{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_{\text{unsh}}$. It is straightforward to see that

$$M_{\text{unsh}} = \frac{\sum_{l=1, n_f} m_l}{n_f}, \quad (16)$$

where $m_l = M_{\text{env}}(r_l)$. Note that $M_{\text{unsh}}$ depends on the fragment distribution, not only on the position of the most distant fragment. To get the mass distribution, $dN_{\text{unsh}}/dM = P_{\text{unsh}}$, we have thus computed the histogram of $M_{\text{unsh}}$ weighted by the probability $P(r_{\text{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.

### 3.3. Results and Comparison with Simulations

We focus on model values that are directly inspired from the values measured in run 1. We adopt $A = 10$, $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.
<figure><img src="image_18.png" /><figcaption>**Figure 8.** Top panel: mean number of fragments as a function of distance. Blue curve: mean number of fragments (see Equation (14)) as a function of $\bar{r} = r/r_e$, $L$ for the run 1 case (i.e., $A = 10$, $M = 6$). Red and black are the mean number of type 1 and type 2 fragments, respectively (see Equations (12) and (13)), in logarithmic shells. The mean number of fragments is equal to one at about $r \simeq 9 r_{e,L}$. Bottom panel: mass distribution of $M_{\text{unsh}}$ (see Equation (16)) for $n_f = 2, 3$, and 4 shielding fragments. The distribution peaks at typically $10 M_L$.</figcaption></figure>

### 3.3.1. Results for Value Parameters Inferred from the Simulations

The top panel of Figure 8 shows $N$ (blue curve), $dN_1$ (red curve), and $dN_2$ (black curve) as a function of $\bar{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} \simeq 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 $\simeq 10 M_L$.

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

The distribution of the unshielded mass, which peaks at about $10 M_L$, is in reasonable agreement with the peak of the distribution observed in run 1 and run 2 (see Figure 2). We remind the reader that $M_L$ is about 0.02–0.03 $M_\odot$ with our numerical implementation, which implies that in run 1 and run

2 the ratio between $M_L$ and the observed peak position is about 5–10.

Since the distribution is inherited from the mean fragment number (see Equation (14)), it is worth comparing it with the simulation results, that is, the mean number of objects in systems, $\langle n_{\text{sys}} \rangle$, shown in the bottom panels of Figure 3. In the left panel (run 1), the four black curves represent respectively (from top to bottom) the values $M = 6, 8, 10$, and 10 but with the assumption $\text{VP}_{\text{ext}} = \text{VP}_{\text{ram}} = 0$. In the right panel (run 2), the adopted values are 4, 6, and 8. For run 2, we used $A = 7$ instead of 10 for run 2. The width of the density perturbation distribution is chosen to be the same for both runs. In both cases, a reasonable agreement is obtained for the highest values of the Mach number explored. The model with $\text{VP}_{\text{ext}} = \text{VP}_{\text{ram}} = 0$ leads to too-few fragments.

The most important disagreement is that the model predicts very few fragments below 100 au (i.e., $\langle n_{\text{sys}} \rangle = 1$), while in the simulations this number is about 0.3–1 (i.e., $\langle n_{\text{sys}} \rangle \simeq 1.3–2$ depending on the time and the run). The most likely explanation is that in the simulations fragments have a tendency to migrate, approaching each other because they are gravitationally bound to each other. Indeed, the typical density at 100 au is about $10^9 \text{ cm}^{-3}$, which corresponds to a freefall time of $10^3 \text{ yr}$. Because this time is short and the fragments likely fall toward each other even before a sink particle is introduced, this effect certainly introduces some differences between the simulations and the model. It may, for instance, lead to a somewhat smaller ratio between the peak of the stellar spectrum and $M_L$.

## 4. Conclusion

We have performed numerical simulations with the AMR code RAMSES to investigate the origin of the peak in the stellar mass function, that is, the characteristic mass of stars. We use two completely different setups and initial conditions to verify that the results are not due to a specific choice. We confirm that the peak is typically equal to about 5–10 times the mass of the FHSC.

To verify that the factor 5–10 is due to further accretion that is eventually limited by the fragmentation of the collapsing cloud at a distance of typically 100–200 au, we have performed complementary numerical experiments in which the fragmentation is forbidden up to 140 au and 280 au from existing stars. We found that, as expected, the peak shifts to higher masses. To quantify this fragmentation, we have computed the typical number of neighbors at various distances and ages of the sink particles. We have also measured several properties such as density, Mach number, and the PDF of density fluctuations in the vicinity of the sink particles. And finally we have extracted the cores that eventually lead to new objects and compute the various terms that enter in the virial theorem.

We have then presented an analytical model to derive the position of the peak of the stellar IMF. It is based on the idea that, as the collapse of a star-forming clump is stopped when dust becomes opaque to its own radiation, matter piles up until at least the mass of the FHSC, $M_L$, is accumulated. Because there is gas surrounding the FHSC, it grows by accretion, and this growth is very important, according to us, in setting the characteristic mass of stars. Due to the density fluctuations within the collapsing envelope surrounding this central core, however, new fragments form, accreting material around them and thus limiting the accretion onto the central object. The
induced tidal forces between core and envelope, however, shear apart these density fluctuations and thus, together with thermal and turbulent support, limit further fragmentation in the envelope.

As a result, the distribution of the final central accreted mass presents a peak at about 10 $M_L$, as indeed inferred in numerical simulations. A quantitative comparison between the mean number of neighbors as a function of distance from the stars measured in the simulations and inferred from the analytical model reveals good agreement, suggesting that the analytical model captures the most important physical processes responsible for setting the peak of the stellar mass function.

We thank the anonymous referee for their constructive comments that have improved the paper. This research has received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007–2013 grant Agreement No. 306483).

# Appendix A

## The Physics of the First Hydrostatic Core: Orders of Magnitude and Trends

Because the present theory relies on the first hydrostatic or so-called first Larson core, we first start by discussing orders of magnitudes and trends. We are in particular interested in getting a simple expression for its mass and typical lifetime. While several studies have performed detailed calculations, simple expressions would be useful for inferring simple dependence.

### A.1. Physics of the First Hydrostatic Core: Simplifying Assumptions

First, we assume that the core is virialized:

$$ \frac{3}{5} \frac{GM_L m_p}{R} \approx \mathcal{H} \times 3 k_B T, \qquad (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\pi/3\rho R^3$, so

$$ M_L \approx \mathcal{H}^{3/2} \times 5^{3/2} \left( \frac{k_B T}{G m_p} \right)^{3/2} \left( \frac{3}{4\pi \rho} \right)^{1/2} . \qquad (18) $$

The hydrostatic core ends when the H₂ molecule starts dissociating, and a second collapse occurs at about $T = T_{dis} \approx 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, $\rho_{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, $\rho_{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 $\gamma = 5/3$. Then when the rotational levels of H₂ get excited, $\gamma$ drops, and at about $T_{ex} \approx 150$ K, we have $\gamma = 7/5$. Thus, we have

$$ \rho_{\text{dis}} \approx \rho_{\text{ex}} \left( \frac{T_{\text{dis}}}{T_{\text{ex}}} \right)^{5/2} \approx \rho_{\text{ad}} \left( \frac{T_{\text{ex}}}{T_0} \right)^{3/2} \left( \frac{T_{\text{dis}}}{T_{\text{ex}}} \right)^{5/2}, \qquad (19) $$

where $T_0$ is the mean cloud temperature in the isothermal regime. We now need to estimate the density $\rho_{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 $\rho_{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_{ad}}$, to the energy radiated per units of time, $\mathcal{L}$:

$$ E_{\text{therm}} = \frac{M_J}{m_p} \times \frac{3}{2} k_B T_0, \qquad (20) $$

$$ \mathcal{L} \approx 4\pi \left(\frac{\lambda_J}{2}\right)^2 \frac{4acT_0^3}{3\kappa_R(T_0)\rho_{\text{ad}}} \frac{T_0}{\lambda_J/2}. \qquad (21) $$

This last expression is simply the radiative flux through a sphere of diameter $\lambda_J$ in which the gradient $\partial_r T$ has been replaced by $T/r$, $\kappa_R$ is the Rosseland opacity, $c$ is the speed of light, and $a = 4\sigma/c$, $\sigma$ being the Stefan–Boltzmann constant. The cooling time is simply estimated as

$$ \tau_{\text{cool}} = \frac{E_{\text{therm}}}{\mathcal{L}}, \qquad (22) $$

while the freefall time is

$$ \tau_{\text{ff}} = \sqrt{\frac{3\pi}{32G\rho_{\text{ad}}}}. \qquad (23) $$

Requiring that $\tau_{\text{cool}} \approx \tau_{\text{ff}}$, we get

$$ \rho_{\text{ad}} = \left( \frac{32}{3\pi} G m_p^4 k_B^{-4} a^2 c^2 \right)^{1/3} T_0^{4/3} \kappa_R^{-2/3}. \qquad (24) $$

### A.2. Results, Prediction, and Qualitative Comparisons

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

$$ \begin{aligned} M_L &\approx \mathcal{H}^{3/2} \times \frac{5^{3/2} 3^{2/3}}{2^{11/6} \pi^{1/3}} (G^{10} m_p^{13} k_B^{-13} a^2 c^2)^{-1/6} \\ &\quad \times T_{\text{dis}}^{1/4} T_{\text{ex}}^{1/2} T_0^{1/12} \kappa_R(T_0)^{1/3} \\ &= \mathcal{H}^{3/2} \times 0.08 M_{\odot} \left(\frac{T_0}{10 \text{ K}}\right)^{1/12} \left(\frac{\kappa_R(T_0)}{0.1 \text{ cm}^2 \text{ g}^{-1}}\right)^{1/3}. \end{aligned} \qquad (25) $$

The estimated value of $M_L$ from numerical simulations is about 0.03 $M_\odot$ (Vaytet & Haugbølle 2017), so $\mathcal{H} \approx 0.5$ leads to the correct value. Indeed, integrating the proper hydrostatic equilibrium (see Equation (7) of Paper II) starting from $\rho_{dis}$ and using the piecewise polytropic equation of state discussed in the text, we obtain a mass equal to 0.03 $M_\odot$. The typical value of $\kappa_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_0$. Second, it depends on the opacity $\kappa_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_0$ is unchanged. There is, however, likely another effect that is due to the fact that at low metallicity the temperature $T_0$ 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 $\kappa_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$^{-2}$. According to the dependence stated by Equation (25) and after doing the proper hydrostatic calculation, this leads to an estimate of $M_L \simeq 0.02 M_\odot$ 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_\odot$.

### A.3. Cooling Time of the First Hydrostatic Core

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.

$$E_{\text{therm}} \approx \frac{M_L}{m_p} \times \frac{5}{2} k_B T_{\text{dis}}, \qquad (26)$$

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

$$\mathcal{L} \approx 4\pi R_L^2 \frac{4acT_{\text{dis}}^3}{3\kappa_R(T_{\text{dis}})\rho_{\text{dis}}} \frac{T_{\text{dis}}}{R_L}. \qquad (27)$$

Using again $M_L = 4\pi/3\rho R^3$ to infer the density in Equations (27) and (17) to infer $RT_{\text{dis}}$, we get

$$\begin{aligned} \tau_{\text{cool}} &= \frac{3^2 5^5}{2^7 \pi^2} \frac{k_B^5 T_{\text{dis}}}{acG^4 m_p^5} \frac{\kappa_R(T_{\text{dis}})}{M_L^2} \\ &\approx 5.6 \times 10^4 \text{ yr} - \frac{\kappa_R(T_{\text{dis}})}{10 \text{ g cm}^{-2}} \left( \frac{M_L}{0.03 M_\odot} \right)^{-2}, \end{aligned} \qquad (28)$$

where $\kappa_R(T_{\text{dis}}) = 10 \text{ cm}^2 \text{ g}^{-1}$ is from Semenov et al. (2003). By contrast, the freefall time at a density $\rho_{\text{dis}}$ is about 1 yr. More generally, a freefall time of $5 \times 10^4 \text{ yr}$ corresponds to a density of roughly $10^6 \text{ cm}^{-3}$, which is far below the density of the first hydrostatic core.

The other relevant time is the accretion one, $\tau_{\text{acc}}$, that is to say the time necessary to accumulate a mass equal to $M_L$,

$$\tau_{\text{acc}} = \frac{M_L}{\frac{dM}{dt}} = 300 \text{ yr} \frac{M_L}{0.03 M_\odot} \left( \frac{\frac{dM}{dt}}{10^{-4} M_\odot \text{ yr}^{-1}} \right)^{-1}, \qquad (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^{-5} M_\odot \text{ yr}^{-1}$ (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 ($\tau_{\text{cool}} \propto M^{-2}$). 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$_2$ 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 \simeq 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 $\rho_{\text{ad}}$ (as stated from Equation (24)), is not clear because the fragments have typically a mass of about $10^{-3} M_\odot$, which is too low to induce the dissociation of H$_2$ 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.

## Appendix B

### Sink Velocity as a Function of Age and Its Dependence on the Number of Neighbors

To investigate the dynamical effect that sink particle neighbors may have on each other, we have computed the velocity variation (i.e., the velocity with respect to the sink initial velocity) as a function of sink age. The idea is to determine how fast the velocity of the sink changes with respect to its initial environment and whether the effect depends on the number of neighbors. Figure 9 displays bidimensional histograms for run 1 and run 2 and for the sinks having zero, two to four, or more than four neighbors within 200 au. Clearly we see that the sink velocity variation increases with time, reaching values on the order of 3–10 km s$^{-1}$ in a few kiloyears. This time is comparable to the accretion time of the sink particles. We also see that the velocity variation is stronger when the sinks have more neighbors. This supports the idea that the neighbors limit the accretion into the sinks by two ways. First they screen the gas that would otherwise have fallen into the primary sink, and second they push it away from its initial gas reservoir.
<figure><img src="image_19.png" /><figcaption>Figure 9. Bidimensional histogram of the sink velocity variation (i.e., with respect to the velocity at the sink birth) versus the sink age (i.e., time since its birth). Top: sinks having no neighbor within 200 au. Middle: sinks having a number of neighbors between two and four. Bottom: sinks having more than four neighbors within 200 au. Left: run 1; right: run 2.</figcaption></figure>

<figure><img src="image_20.png" /></figure>

## Appendix C
Estimating the External Pressure in the Virial Theorem

As discussed in Section 3.2.3, the turbulent pressure could have a very strong impact on the result by decreasing drastically the number of fragments. Therefore a more detailed analysis is required, particularly on the role of the external pressure that must also be taken into account.

We remember the virial theorem

$$ \frac{1}{2} \frac{\partial^2 I}{\partial t^2} = -E_g + 2E_{\text{th}} + 2E_{\text{kin}} - \oint P_{\text{th}} \mathbf{r} \cdot d\mathbf{S} - \oint \rho (\mathbf{v} \cdot \mathbf{r})(\mathbf{v} \cdot d\mathbf{S}), \quad (30) $$

where the various terms are the second derivative of the inertia momentum, the gravitational energy, and the thermal and kinetic energies, while the last two terms stand for the thermal and ram pressure in the outskirts of the cores. We remind the
reader that the last term comes from the equality

$$\int \rho \frac{dv}{dt} r dV = \frac{1}{2} \frac{\partial^2 I}{\partial t^2} - 2E_{\text{kin}} + \int \nabla \cdot (\rho(\mathbf{v} \cdot \mathbf{r})\mathbf{v}) dV. \quad (31)$$

Below, we seek expressions for $2E_{\text{ther}}$, $2E_{\text{kin}}$, VP$_{\text{ext}}$, and VP$_{\text{ram}}$. Assuming spherical symmetry, we have

$$VP_{\text{ext}} = \oint P_{\text{th}} r \cdot dS = 3P_{\text{ext}} V_{\text{vol}}, \quad (32)$$

where $V_{\text{vol}}$ is the volume, and

$$VP_{\text{ram}} = \oint \rho (\mathbf{v} \cdot \mathbf{r})(\mathbf{v} \cdot dS) = 3P_{\text{ram}} V_{\text{vol}}, \quad (33)$$

where $P_{\text{ram}} = \rho_{\text{ext}} v_n^2$ is the ram pressure at the clump boundary.

Let us consider a spherical clump of external radius $r_0$, outer density $n_0$, temperature $T_0$, and typical velocity dispersion $\sigma_0$. To estimate the various contributions, we also need to know the density distribution within the clump, which for the sake of simplicity we assume to be $n(r) \propto r^{-\alpha}$, with $\alpha$ ranging typically between 0 and 2. We also assume that $\sigma(r) = \sigma_0(r/r_0)^{0.5}$. The external density is $n_{\text{ext}}$. Typically we expect it to be related to the outer density (i.e., the density at the inner edge of the clump), $n_0$, by a factor on the order of $M^2$, the typical Mach number at the clump scale. We have $n_{\text{ext}} = Cn_0$, where $C$ ranges between 0 and 1. Given the Mach numbers, which at the scale of the clumps are on the order of 1–2, this leads to a contrast, $C$, of 0.25–1. This contrast, however, is along the direction of the confining shock and therefore for a spherical clump is less pronounced.

This leads to

$$2E_{\text{ther}} = \frac{9}{3-\alpha} M_0 C_{s,0}^2, \quad (34)$$

where $M_0 = 4\pi/3r_0^3n_0m_p$, $C_{s,0}^2 = k_bT_0/m_p$,

$$2E_{\text{kin}} = \frac{3M^2}{4-\alpha} M_0 C_{s,0}^2, \quad (35)$$

$$3P_{\text{ext}}V = 3CM_0C_{s,0}^2, \quad (36)$$

and the ram pressure

$$3P_{\text{ram}} V_{\text{vol}} = \oint \rho_0 v_n^2 r ds = 3CM_0 \sigma_{0,\text{sph}}^2 \\ = 3CM_{\text{sph}}^2 M_0 C_{s,0}^2, \quad (37)$$

where $\sigma_{0,\text{sph}}^2 = \oint v_n^2 r ds / \oint r ds$ is the spherical part of the velocity field at the surface of the clump. Thus,

$$2E_{\text{ther}} + 2E_{\text{kin}} - 3P_{\text{ext}}V = 3M_0 C_{s,0}^2 \\ \times \left( \frac{3}{3-\alpha} - C + \frac{M^2}{4-\alpha} - CM_{\text{sph}}^2 \right). \quad (38)$$

Thus we see that for $\alpha = 1$, the effective thermal energy, corrected from the external pressure, should be multiplied by $3/(3 - \alpha) - C$, while the effective kinetic energy should be multiplied by $1 - 3CM_{\text{sph}}^2/M^2$. The external pressure exerted by the velocity at the edge of the clump reduces the turbulent pressure. If the spherical mode represents one-third of the energy and if $C = 1$, then the turbulent pressure cancels out exactly. The exact value of $1 - 3CM_{\text{sph}}^2/M^2$ is hard to infer analytically, but it seems likely that the value $2E_{\text{kin}}$ overestimates the turbulent pressure. Let us stress in particular that if there is no clear contrast between the external density and the clump density, that is, $C \simeq 1$, then the value of 1–3

$CM_{\text{sph}}^2/M^2 \simeq 0$ since $M_{\text{sph}}^2 \simeq M^2/3$ (as is the case for a plane-parallel converging flow, for example).

## Appendix D
Envelope Temperature

Strictly speaking, the isothermal assumption is no longer valid, due to the heating induced by the gravitational contraction and the fact that the optical depth becomes close to or larger than 1. To address this issue, we use the effective equation of state inferred to fit the thermal behavior of dense material (Masunaga et al. 1998; Vaytet et al. 2013; Vaytet & Haugbølle 2017) used in the present simulations:

$$T = T_0[1 + (\rho/\rho_{\text{ad}})^{(\gamma-1)}], \quad (39)$$

where we adopt $T_0 = 10$ K, $\rho_{\text{ad}} = 3.4 \times 10^{-14}$ g cm$^{-3}$, and $\gamma = 5/3$. This expression holds at early time close to the formation of the Larson core. At later time, the temperature increases further, but here we are concerned with the very early fragmentation stages for low-mass stars, which have masses of a few 0.1 $M_\odot$. Moreover, the increase in temperature remains limited for most of the dense gas. When solving Equation (9), with the condition that the mass of the fluctuation is $M_p = M_L$, we find either zero or two solutions, corresponding respectively to $\rho_{\text{min}}$ and $\rho_{\text{max}}$. This stems from the fact that at high density, the thermal pressure always dominates gravity due to the stiff value of $\gamma$. The cases without solution arise when $A$ is high, typically $>20$. In practice, in the relevant regime, we found out that this does not yield significant differences compared to the isothermal case. The reason is that the typical distance between the shielding fragments is 3–10 $r_{e,L}$, which corresponds to a distance of 60–200 au. For $A = 10$, this corresponds to densities $3 \times 10^8-3 \times 10^9$ cm$^{-3}$, at which the gas is nearly isothermal.

## ORCID iDs

Patrick Hennebelle https://orcid.org/0000-0002-0472-7202

Gilles Chabrier https://orcid.org/0000-0002-8342-9149

## References

* Ballesteros-Paredes, J., Hartmann, L. W., Pérez-Goytia, N., & Kuznetsova, A. 2015, MNRAS, 452, 566

* Bastian, N., Covey, K. R., & Meyer, M. R. 2010, ARA&A, 48, 339

* Basu, S., Gil, M., & Auddy, S. 2015, MNRAS, 449, 2413

* Basu, S., & Jones, C. E. 2004, MNRAS, 347, L47

* Bate, M. R. 1998, ApJL, 508, L95

* Bate, M. R. 2009, MNRAS, 392, 1363

* Bate, M. R. 2014, MNRAS, 442, 285

* Bleuler, A., & Teyssier, R. 2014, MNRAS, 445, 4015

* Bonnell, I. A., & Bate, M. R. 2006, MNRAS, 370, 488

* Bonnell, I. A., Bate, M. R., Clarke, C. J., & Pringle, J. E. 2001, MNRAS, 323, 785

* Bonnell, I. A., Smith, R. J., Clark, P. C., & Bate, M. R. 2011, MNRAS, 410, 2339

* Chabrier, G. 2003, PASP, 115, 763

* Dib, S., Kim, J., Vázquez-Semadeni, E., Burkert, A., & Shadmehri, M. 2007, ApJ, 661, 262

* di Francesco, J., Evans, N. J., II, Caselli, P., et al. 2007, Protostars and Planets V (Tucson, AZ: Univ. Arizona Press), 17

* Eisenstein, D. J., & Hut, P. 1998, ApJ, 498, 137

* Federrath, C., Klessen, R. S., & Schmidt, W. 2008, ApJL, 688, L79

* Foster, P. N., & Chevalier, R. A. 1993, ApJ, 416, 303

* Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371
Girichidis, P., Federrath, C., Banerjee, R., & Klessen, R. S. 2011, MNRAS, 413, 2741

Guszejnov, D., Krumholz, M. R., & Hopkins, P. F. 2016, MNRAS, 458, 673

Hennebelle, P. 2018, A&A, 611, A24

Hennebelle, P., & Chabrier, G. 2008, ApJ, 684, 395

Hennebelle, P., & Falgarone, E. 2012, A&ARv, 20, 55

Hennebelle, P., Whitworth, A. P., Gladwin, P. P., & André, P. 2003, MNRAS, 340, 870

Hopkins, P. F. 2012, MNRAS, 423, 2037

Inutsuka, S.-i 2001, ApJL, 559, L149

Jappsen, A.-K., Klessen, R. S., Larson, R. B., Li, Y., & Mac Low, M.-M. 2005, A&A, 435, 611

Kroupa, P. 2001, MNRAS, 322, 231

Krumholz, M. R., Myers, A. T., Klein, R. I., & McKee, C. F. 2016, MNRAS, 460, 3272

Larson, R. B. 1969, MNRAS, 145, 271

Lee, Y.-N., & Hennebelle, P. 2016, A&A, 591, A30

Lee, Y.-N., & Hennebelle, P. 2018a, A&A, 611, A88

Lee, Y.-N., & Hennebelle, P. 2018b, A&A, 611, A89

Lee, Y.-N., & Hennebelle, P. 2019, A&A, 622, A125

Low, C., & Lynden-Bell, D. 1976, MNRAS, 176, 367

Machida, M. N., Tomisaka, K., Matsumoto, T., & Inutsuka, S.-i. 2008, ApJ, 677, 327

Maschberger, T., Bonnell, I. A., Clarke, C. J., & Moraux, E. 2014, MNRAS, 439, 234

Masunaga, H., & Inutsuka, S.-i 1999, ApJ, 510, 822

Masunaga, H., Miyama, S. M., & Inutsuka, S.-i. 1998, ApJ, 495, 346

Murray, N., & Chang, P. 2015, ApJ, 804, 44

Ntormousi, E., & Hennebelle, P. 2019, A&A, 625, A82

Offner, S. S. R., Clark, P. C., Hennebelle, P., et al. 2014, Protostars and Planets VI (Tucson, AZ: Univ. Arizona Press), 53

Padoan, P., Nordlund, A., & Jones, B. J. T. 1997, MNRAS, 288, 145

Peters, T., Banerjee, R., Klessen, R. S., & Mac Low, M.-M. 2011, ApJ, 729, 72

Rees, M. J. 1976, MNRAS, 176, 483

Salpeter, E. E. 1955, ApJ, 121, 161

Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713

Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611

Shu, F. H. 1977, ApJ, 214, 488

Stamer, T., & Inutsuka, S.-i 2018, ApJ, 869, 179

Teyssier, R. 2002, A&A, 385, 337

Traficante, A., Duarte-Cabral, A., Elia, D., et al. 2018, MNRAS, 477, 2220

Vaytet, N., Chabrier, G., Audit, E., et al. 2013, A&A, 557, A90

Vaytet, N., & Haugbølle, T. 2017, A&A, 598, A116

Ward-Thompson, D., André, P., Crutcher, R., et al. 2007, Protostars and Planets V (Tucson, AZ: Univ. Arizona Press), 33