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

# What Is the Role of Stellar Radiative Feedback in Setting the Stellar Mass Spectrum?

Patrick Hennebelle¹, Benoît Commerçon², 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

²École normale supérieure de Lyon CRAL, UMR CNRS 5574, F-69364, Lyon Cedex 07, France

³Department of Earth Sciences, National Taiwan Normal University, 88, Section 4, Ting-Chou Road, Taipei 11677, Taiwan

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

Received 2020 July 26; revised 2020 October 4; accepted 2020 October 7; published 2020 December 4

## Abstract

In spite of decades of theoretical efforts, the physical origin of the stellar initial mass function (IMF) is still debated. Particularly crucial is the question of what sets the peak of the distribution. To investigate this issue, we perform high-resolution numerical simulations with radiative feedback exploring, in particular, the role of the stellar and accretion luminosities. We also perform simulations with a simple effective equation of state (EOS), and we investigate 1000 solar-mass clumps having, respectively, 0.1 and 0.4 pc of initial radii. We found that most runs, both with radiative transfer or an EOS, present similar mass spectra with a peak broadly located around 0.3–0.5 $M_\odot$ and a power-law-like mass distribution at higher masses. However, when accretion luminosity is accounted for, the resulting mass spectrum of the most compact clump tends to be moderately top-heavy. The effect remains limited for the less compact one, which overall remains colder. Our results support the idea that rather than the radiative stellar feedback, this is the transition from the isothermal to the adiabatic regime, which occurs at a gas density of about $10^{10} \text{ cm}^{-3}$, that is responsible for setting the peak of the IMF. This stems from (i) the fact that extremely compact clumps for which the accretion luminosity has a significant influence are very rare and (ii) the luminosity problem, which indicates that the effective accretion luminosity is likely weaker than expected.

*Unified Astronomy Thesaurus concepts: Initial mass function (796); Star formation (1569); Collapsing clouds (267); Stellar feedback (1602); Stellar accretion (1578); Radiative transfer (1335); Hydrodynamical simulations (767)*

## 1. Introduction

In trying to understand the origin of the mass distribution of stars, the initial mass function (IMF; Salpeter 1955; Kroupa 2001; Chabrier 2003; Bastian et al. 2010; Offner et al. 2014; Lee et al. 2020) is a fundamental issue to unravel the history of the universe. In particular, the fact that the IMF seems, at first sight, to be universal, that is to say, weakly varies from one environment to another, remains a puzzle, although some recent variations have been claimed (e.g., Cappellari et al. 2012; Chabrier et al. 2014; Schneider et al. 2018). Various theories have been proposed to explain the origin of the IMF. This includes a correspondence between the core mass function and the IMF essentially through analytical modeling (Inutsuka 2001; Padoan et al. 1997; Hennebelle & Chabrier 2008; Hopkins 2013), numerical simulations of a fragmenting cloud using sink particles to represent the stars (Bonnell et al. 2011; Girichidis et al. 2011; Ballesteros-Paredes et al. 2015; Guszejnov et al. 2020; Padoan et al. 2020), or analytical statistical description of stellar accretion (Basu & Jones 2004; Basu et al. 2015). In general, the high-mass tail of the IMF is reasonably reproduced in these models, although the physical reasons invoked are different. In these models, the universality of the slope relies on the invoked scale-free processes, gravity, and/or turbulence.<fnref n="5" /> The question of the peak appears, however, to be more complicated because most theories are based on the Jeans mass, which depends on the gas density and temperature,

and thus inferring a characteristic mass, say, around $0.3 M_\odot$, that does not vary significantly with the physical conditions is a challenge. Most proposed explanations consist in identifying mechanisms that could result in a weak dependence of the effective Jeans mass on gas density (Lee & Hennebelle 2018b; Guszejnov et al. 2020). For instance, Hennebelle & Chabrier (2008), Hennebelle (2012), and Lee & Hennebelle (2016b) proposed that there is a compensation between the density and Mach number variations, Jappsen et al. (2005) argued that the change of the effective equation of state (EOS) at a density of about $10^5 \text{ cm}^{-3}$ makes the corresponding Jeans mass play a dominant role, while Bate (2009), Krumholz et al. (2016), and Guszejnov et al. (2016) proposed that radiative feedback heats up the gas at very high density (Krumholz et al. 2007), i.e., $10^{8-10} \text{ cm}^{-3}$, setting up again a Jeans mass that weakly depends on density, for instance.

A somewhat different explanation has recently been proposed by Lee & Hennebelle (2018b) and Hennebelle et al. (2019), who argue that the peak of the IMF is directly linked to the mass, $M_L$, of the first hydrostatic Larson core (FHSC Larson 1969; Masunaga et al. 1998; Vaytet et al. 2013; Vaytet & Haugbølle 2017; Bhandare et al. 2018, 2020), which is the hydrostatic core that forms when the dust becomes opaque to radiation. The mass of FHSC is about $M_L \simeq 0.03 M_\odot$, which is about 10 times below the peak of the IMF. However, performing high resolution of collapsing 1000 $M_\odot$ clumps, these authors infer that the peak of the IMF is about 5–10 $M_L$. This is due to further accretion from the envelope onto the FHSC, which is eventually halted when new fragments form (Hennebelle et al. 2019). Because, in particular, of the stabilizing effect of the tidal forces (Lee & Hennebelle 2018b; Colman & Teyssier 2020), the immediate

<footnote n="5">⁵ Strictly speaking, turbulence and gravity are not entirely scale-free. However, in the fully nonlinear regime, there is usually a broad range of scales over which a high Reynolds turbulent flow is self-similar. The same is true for a self-gravitating fluid that has developed a power-law density PDF and for which there is not a unique Jeans length.</footnote>
neighborhood of the FHSC is stable against gravitational instability, and finding another FHSC requires to go at a distance L such that the mass enclosed in the sphere of radius L is about 5–10 $M_L$. Very importantly, changing the initial conditions of the initial clumps, initial density, Mach number, or magnetic field (Lee & Hennebelle 2019) by orders of magnitude is found to leave the peak of the stellar mass spectrum almost unaffected.

So far the simulations performed to investigate the FHSC-based theory have been using an effective barotropic EOS aimed at mimicking the thermal behavior of the gas at densities above $10^{10} \text{ cm}^{-3}$. While the approach was rather useful to establish and test these ideas, radiative transfer calculations are mandatory for a more realistic treatment. In particular, since the origin of the FHSC is the high optical depth that makes the gas adiabatic, it is important to test the theory in this context. Several attempts have been made to study the IMF using radiative transfer calculations. Urban et al. (2010) performed SPH calculations and introduce the sink particles at a density of about $10^8 \text{ cm}^{-3}$; they include radiative feedback onto the sink particles, which includes both the stellar and accretion luminosity. They found that radiative transfer calculations are quite different from the isothermal ones, in particular, the stars are much more massive when radiation is considered. Bate (2009) also performed high-resolution SPH calculations but introduced the sink particles at very high density, i.e., $n > 10^{19} \text{ cm}^{-3}$. This includes the optically thick regime, which occurs at a density $n > 10^{10} \text{ cm}^{-3}$, but the simulations do not add any stellar feedback onto the sink particles. By doing this, the IMF presents a peak at about 0.3 $M_\odot$ and a mass spectrum at high mass that is clearly flatter than the Salpeter's exponent of 1.3 (in $dN/d\log m$). Krumholz et al. (2012) performed adaptive mesh refinement (AMR) calculations with a resolution of 20–40 au. The sinks are introduced when the Jeans conditions get violated, that is to say, when the mesh size is larger than 1/10 of the local Jeans length, and they consider only objects more massive than 0.05 $M_\odot$ as being stars, the smaller ones being allowed to merge. Both intrinsic luminosity and accretion luminosity are added to the sinks. By doing so, they obtain mass spectra that are almost flat, that is, $dN/d\log M \propto M^0$, when winds are not considered, while in the presence of winds, which allow the radiation to escape, the mass spectra present a peak around 0.3 $M_\odot$ and a power law, $dN/d\log M \propto M^{-\alpha}$, with $\alpha \simeq 0.5–1$. Recently Mathew & Federrath (2020) performed simulations with a spatial resolution of 200 au and compared runs that either use a polytropic EOS or take into account the stellar heating assuming either spherical symmetry or a polar distribution. They infer stellar mass spectra that peaked at about 2 $M_\odot$ and found that when heating is included more massive stars would form.

In the present paper, we want to explore further the influence of the radiative feedback in establishing the stellar mass spectrum during the collapse of a massive clump. In particular, we stress that so far studies that do consider the accretion luminosity (e.g., Urban et al. 2010; Krumholz et al. 2012; Mathew & Federrath 2020) have been performed at relatively coarse resolution. The lack of resolution can be particularly severe in this context because it may result in overestimating the mass at which the stellar distribution peaks (e.g., Ntormousi & Hennebelle 2019). This, in turn, implies that too much mass is contained in massive stars, and since they exert a strong

feedback onto their environment, this may lead to overestimating the importance of radiative feedback.

To determine the impact of various contributions, we present a set of both barotropic and radiative transfer calculations, taking into account the different contributions of the feedback luminosity and for different types of initial conditions. By performing these various runs, we can in particular distinguish between the influence of the optically thick and hydrostatic phase (the FHSC) and the heating of the collapsing envelope by radiation. The plan of the paper is as follows. In Section 2, we present the numerical setup and the various assumptions done to perform the two types of simulations. In Section 3, we present and discuss our results regarding the temperature distribution through the clouds. Section 4 is devoted to the stellar mass spectra and how they depend on the radiative feedback and on the initial conditions. Section 6 concludes the paper.

## 2. Numerical Simulations

### 2.1. Numerical Methods and Setup

All simulations were run with the AMR magnetohydrodynamics (MHD) code RAMSES (Teyssier 2002; Fromang et al. 2006), though in this work magnetic field is not considered. Two types of simulations have been carried out.

The first type of simulation uses radiative transfer using the flux diffusion method and the gray approximation as described in Commerçon et al. (2011, 2014). At high density, the EOS is the one given by Saumon & Chabrier (1992) and Saumon et al. (1995) which models the thermal properties of a gas containing the species H2, H, H+, He, He+, and He2+ (the He mass concentration is 0.27). The opacities are as described in Vaytet et al. (2013). For the range of temperature and densities covered in this work, the opacities are essentially the ones calculated in Semenov et al. (2003).

The second type of simulations employs an effective EOS and no radiative transfer. The prescription is the same as the one used in Hennebelle et al. (2019),

$$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 \text{ K}$, $n_{\text{ad},2} = 30 n_{\text{ad}}$, $\gamma_1 = 5/3$, and $\gamma_2 = 7/5$. This EOS mimics the thermal behavior of the gas when it becomes non-isothermal. Two values of $n_{\text{ad}}$ have been explored, namely, $(n_{\text{ad}})_1 = 4 \times 10^{10} \text{ cm}^{-3}$ and $(n_{\text{ad}})_2 = 1.2 \times 10^{11} \text{ cm}^{-3}$. The latter one has been chosen because it appears to be closer to the actual temperature of the radiative transfer calculations.

The boundary conditions used in this work are periodic, and we simulate a spherical cloud whose radius is four times lower than the computational domain size. All simulations were run on a base grid of $2^8$, and typically seven to eight (up to 10) AMR levels have been added, leading to a total number of 15 or 16 AMR levels. The Jeans length is resolved with at least 10 points. In Appendix B less and more resolved runs are presented to investigate the issue of numerical convergence and test the influence of sink particle numerical parameters.

### 2.2. Sink Particles

We used the sink particle algorithm of Bleuler & Teyssier (2014). Sink particles are formed at the highest refinement level at the peak of clumps whose density is larger than $n \ge n_{\text{acc}}/10$,
while the sinks are introduced at a density $n = n_{acc}$. Only clumps that satisfy a series of criteria indicating sufficient gravitational boundness (see Bleuler & Teyssier 2014) may lead to sink formation. The value of $n = n_{acc}$ is equal to either $n_{acc} = 10^{12} \text{ cm}^{-3}$ or $n_{acc} = 10^{13} \text{ cm}^{-3}$ and is discussed below. Typically, the value of $n_{acc}$ is chosen such that a computational cell having a density equal to $n_{acc}$ contains a mass that is about 1%–2% of the mass of the first hydrostatic core, i.e., about $M_L = 0.03 M_\odot$. At each time step, 10% of the gas mass that is located within the sinks and has a density larger than $n_{acc}$ is removed from the grid and transferred to the sink. Lee & Hennebelle (2018b) and Hennebelle et al. (2020) have modified this value and conclude that it does not affect significantly the accretion rate onto the stars (essentially because the accretion is controlled by the larger scales), while, on the other hand, it may affect the disk that forms around the sink. The sinks are not allowed to merge. We stress that, according to us, it is necessary to describe sufficiently the FHSC in order to get the peak of the IMF. Introducing sink particles when the Jeans criterion is not satisfied, for instance, is a viable approach in the isothermal phase only and is not suited to the physics of the FHSC, which is essentially adiabatic.

## 2.3. Stellar Feedback and Accretion Luminosity

An important aspect is the radiative feedback emitted by the sink particles. Two types of contributions have to be taken into account. The first is the accretion luminosity, which is given by

$$L_{\text{acc}} = \frac{f_{\text{acc}} GM_* \dot{M}}{R_*} \quad (2)$$

When assuming that all of the accretion gravitational energy is radiated away, we have $f_{acc} \simeq 1$. This has been shown to be the dominant source of heating at early times and to have important consequences on the cloud (e.g., Krumholz et al. 2007; Offner et al. 2009). The second is the intrinsic luminosity of the protostars (mimicked by the sinks), $L_*$.

The difficulty for star-forming calculations is that the protostars are embedded and still heavily accreting. In our calculations, we use the radius and stellar luminosity given by Kuiper & Yorke (2013; see also Hosokawa & Omukai 2009), which have developed models that take into account the accretion. As will be seen later, both effects can have significant influence on the outcome of the calculations. There are, however, serious uncertainties here. These radiations are emitted at the stellar surface, i.e., at a few solar radii. To what extent it is accurate to introduce them isotropically at a scale of a few astronomical units is highly uncertain. Indeed, Krumholz et al. (2012) concluded that in the presence of a wind cavity, much of the radiation may escape, and this would limit the impact of the radiation heating. Performing collapse up to stellar densities, Bate (2010) found that the accretion luminosity may halt the accretion onto the young protostar and drive a jet. Therefore, investigating the exact consequences of accretion luminosity requires detailed small-scale calculations that are not reachable in calculations aimed at getting the mass distribution of stars. More generally, it is now firmly established that the luminosity of protostars is significantly below the expected values (Kenyon & Hartmann 1995; Evans et al. 2009; Offner & McKee 2011; Stamatellos et al. 2011) and presents a considerable scatter, which has been

A fundamental question is, What is the value of the parameter $f_{acc}$ in Equation (2), that is to say, what is the value of the effective accretion luminosity? Based on an accretion burst frequency estimate, Offner & McKee (2011), following Dunham et al. (2010), estimated $f_{acc} \simeq 0.23$. We note that if some of the accretion energy is radiated away through mechanical processes, such as jets, or if the emission of the accretion luminosity is isotropic and concentrated in the direction of the jets, for instance, $f_{acc}$ could be further reduced.

Since it is likely that a substantial fraction of the accretion energy is radiated in short bursts, a key aspect lies in the comparison between the cooling and the dynamical time. The former can be estimated as

$$\tau_{\text{cool}} \simeq \frac{E_{\text{therm}}}{\partial_{R}F_{R}} \simeq \frac{k_{\text{B}}Tnr^{2}\kappa\rho}{caT^{4}}, \quad (3)$$

and the latter is simply the freefall time $\tau_{ff} = \sqrt{3\pi/(32G\rho)}$, where all expressions have their usual meaning, $\kappa$ is the opacity, $c$ is the speed of light, $a$ is the radiative constant, $r$ is the radius, $k_B$ is the Boltzmann constant, $F_R$ is the radiative flux, and $T$ is the temperature. Computing the ratio of these two times for physical conditions corresponding to the left panel of Figure 6, we found that $\tau_{\text{cool}}/\tau_{ff}$ is typically between $10^{-2}$ and $10^{-6}$; therefore, the gas temperature adjusts instantaneously to the source luminosity, and thus if it is experiencing short bursts, as inferred from observations, the gas dynamics is not affected. For this reason we have performed runs in which the effective accretion luminosity has been multiplied by a factor $f_{acc} = 0.5, 0.1$, and even 0 (that is to say, $L_{\text{acc}}$ is ignored).

Another parameter that needs to be determined is when, i.e., for which mass of the sink, to start injecting the accretion luminosity onto the sink particle, and we choose to do so when the sink has a mass of about $2 M_L$, i.e., $0.07 M_\odot$. The reason is that, due to the limited spatial resolution, when the sink is introduced, the protostar is not formed yet. Since the size of the sink particles is not very different from the radius of the FHSC, it seems reasonable to assume that the protostar is formed only when the sink reaches a mass equal at least to $M_L$. Since in the time delay more gas falls into the sink, and since it is not very clear what is the minimum mass of the protostar for which the accretion luminosity can be described by Equation (2), in particular, because the accretion rate is measured at a scale of the sink particle, we choose to start the accretion luminosity at $2 M_L$.

## 2.4. Initial Conditions and Runs Performed

We consider spherical clouds in which turbulence has been added and is freely decaying. The velocity perturbations present a power spectrum that is equal to $11/3$ and aim at reproducing a standard turbulent flow while the phases are random. Note that we do not start by running the code without gravity, as Lee & Hennebelle (2018a) found that it makes little difference at least for the cases they explored. To get relevant initial conditions, we looked at distributions of observed star-forming clumps, such as the ones of the ATLASGAL (Urquhart et al. 2014) and Hi-GAL
<table><caption>**Table 1**
Summary of the Runs Performed</caption><thead><tr><th>Name</th><th>R<sub>c</sub> (pc)</th><th>M</th><th>Imax</th><th>dx (au)</th><th>n<sub>ad</sub> (cm<sup>-3</sup>)</th><th>n<sub>acc</sub> (cm<sup>-3</sup>)</th><th>Stellar Lum</th><th>f<sub>acc</sub></th></tr></thead><tbody><tr><td>COMP-ACLUM</td><td>0.1</td><td>10</td><td>15</td><td>2.3</td><td>N/A</td><td>10<sup>13</sup></td><td>yes</td><td>0.5</td></tr><tr><td>COMP-LOWACLUM</td><td>0.1</td><td>10</td><td>15</td><td>2.3</td><td>N/A</td><td>10<sup>13</sup></td><td>yes</td><td>0.1</td></tr><tr><td>COMP-NOACLUM</td><td>0.1</td><td>10</td><td>15</td><td>2.3</td><td>N/A</td><td>10<sup>13</sup></td><td>yes</td><td>0</td></tr><tr><td>COMP-NOFEED</td><td>0.1</td><td>10</td><td>15</td><td>2.3</td><td>N/A</td><td>10<sup>13</sup></td><td>no</td><td>0</td></tr><tr><td>COMP-bar1</td><td>0.1</td><td>10</td><td>15</td><td>2.3</td><td>1.2 10<sup>11</sup></td><td>10<sup>13</sup></td><td>N/A</td><td>N/A</td></tr><tr><td>COMP-bar2</td><td>0.1</td><td>10</td><td>15</td><td>2.3</td><td>4 10<sup>10</sup></td><td>10<sup>13</sup></td><td>N/A</td><td>N/A</td></tr><tr><td>STAN-ACLUM</td><td>0.4</td><td>5</td><td>16</td><td>4.6</td><td>N/A</td><td>10<sup>12</sup></td><td>yes</td><td>0.5</td></tr><tr><td>STAN-LOWACLUM</td><td>0.4</td><td>5</td><td>16</td><td>4.6</td><td>N/A</td><td>10<sup>12</sup></td><td>yes</td><td>0.1</td></tr><tr><td>STAN-NOFEED</td><td>0.4</td><td>5</td><td>16</td><td>4.6</td><td>N/A</td><td>10<sup>12</sup></td><td>no</td><td>0</td></tr><tr><td>STANHMAC-ACLUM</td><td>0.4</td><td>10</td><td>16</td><td>4.6</td><td>N/A</td><td>10<sup>12</sup></td><td>yes</td><td>0.5</td></tr><tr><td>COMP-NOACLUMlrls</td><td>0.1</td><td>10</td><td>14</td><td>4.6</td><td>N/A</td><td>10<sup>12</sup></td><td>yes</td><td>0</td></tr><tr><td>COMP-NOACLUMls</td><td>0.1</td><td>10</td><td>15</td><td>2.3</td><td>N/A</td><td>10<sup>12</sup></td><td>yes</td><td>0</td></tr><tr><td>COMP-NOACLUMhr</td><td>0.1</td><td>10</td><td>16</td><td>1.15</td><td>N/A</td><td>10<sup>13</sup></td><td>yes</td><td>0</td></tr><tr><td>STAN-ACLUMhrhs</td><td>0.4</td><td>5</td><td>17</td><td>2.3</td><td>N/A</td><td>10<sup>13</sup></td><td>yes</td><td>0.5</td></tr><tr><td>STAN-ACLUMvhrhs</td><td>0.4</td><td>5</td><td>18</td><td>1.15</td><td>N/A</td><td>10<sup>13</sup></td><td>yes</td><td>0.5</td></tr></tbody></table>

Note. Imax is the maximum level of grid used. For the COMP runs Imax = 15 corresponds to about 2.3 au of resolution, while for the STAN runs Imax = 16 corresponds to a resolution of about 4.6 au. n_ad is the density at which the gas becomes adiabatic. Stellar luminosity indicates whether it is taken into account, and f'acc gives the fraction of the accretion luminosity that is taken into account in the calculation. N/A stands for not applicable.

surveys (Elia et al. 2017). In both surveys the clump mass spans a range that typically goes from $100 \sim 10^4 M_\odot$, with a few clumps that have lower or larger values. The radius has been found to depend on the mass; typically one has $M \propto R^2$ (see Lee & Hennebelle 2016a, 2016b, for an explanation of this relation), but for a given mass, there is a spread in radius. Typically for 1000 solar-mass clumps, the observed radius goes from 0.1 to 1 pc. We stress that the final galactic IMF should definitely be obtained by summing the stellar mass distribution of a clump distribution that reflects these observations (see Lee et al. 2017).

In this work, we consider clumps having a uniform density initially of mass $10^3 M_\odot$. They have initially a radius of either about 0.1 pc, corresponding to an initial density of $5 \times 10^6 cm^{-3}$, or 0.4 pc, corresponding to a density of about $8 \times 10^4 cm^{-3}$. Observationally, this seems to correspond to a very compact star-forming clump and to a standard one. Below we refer to the first type of initial condition as COMP (for compact) and to the second as STAN (for standard). Indeed, observations of massive star-forming clumps found that a radius of 0.4 pc is typical for a clump of $10^3 M_\odot$, while a radius of 0.1 pc corresponds to more extreme clouds (e.g., Urquhart et al. 2014; Elia et al. 2017). The COMP runs have a freefall time of about 14 kyr, while for the STAN ones it is about 110 kyr. The initial temperature is equal to 10 K, and the ratio of thermal over gravitational energy is about 0.002 for the COMP cases and 0.008 for the STAN ones. The initial value of the Mach number is 10 for the COMP runs and 5 for the STAN runs, which leads to the same turbulent over gravitational energy ratio. To investigate the influence of the initial Mach number, a STAN run with $M = 10$ is also performed (STANHMAC). This means that the turbulent over gravitational energy ratio is about 0.2 for all clumps except for STANHMAC, for which it is 0.8. Let us remind that Lee & Hennebelle (2018a), exploring the influence of the initial velocity dispersion on the stellar mass spectrum, concluded that its initial amplitude has a modest influence as long as the cloud is bound and that it is not too small (Lee & Hennebelle 2018a infer that it should be larger than $\sim 0.1$).

To understand the impact of the radiative feedback processes, we perform various runs that include none of them (NOFEED), only the stellar luminosity (NOACLUM), or both the stellar and the accretion luminosity (ACLUM). Moreover, as already mentioned, since the actual value of the effective accretion luminosity that must be used is unclear, we perform runs for which the accretion luminosity is $GM_*dM/dt/R_*$ divided by 2 and by 10 (LOWACLUM).

For comparisons with the radiative transfer calculations, two runs with a barotropic EOS are also done with two different values of the parameter $n_{ad}$.

Note that since the STAN-type clouds are four times more spatially extended than the COMP-type ones, it has not been possible to run the STAN simulations with the same spatial resolution as the COMP ones except for two runs (STAN-ACLUMhrhs and STAN-ACLUMvhrhs, which are employed to investigate the issue of numerical convergence). The COMP simulations have a nominal spatial resolution of about 2 au, while the STAN ones have 4 au. Because of this difference, the value of $n_{acc}$ in STAN runs is chosen to be $10^{12} cm^{-3}$, while it is equal to $10^{13} cm^{-3}$ for the COMP runs.

Finally, to test the influence of numerical parameters, we have also performed runs that have both lower and higher spatial resolutions, as well as runs that investigate the influence of $n_{acc}$. These runs are discussed in Appendix B.

Table 1 summarizes the various runs performed.

## 3. General Description

### 3.1. Illustrations for a Specific Case

To illustrate the simulation results, large-scale and small-scale images are displayed. Figure 1 portrays two snapshots of runs STAN-ACLUMhrhs at early time after only a few sink particles (black dots) have formed and later when a significant fraction of the gas has turned into stars. The stars tend to form in dense filaments and are strongly clustered. A small-scale view is given in Figure 2, where the column density around the
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Column density at two snapshots for run STAN-ACLUMhrhs. The black dots represent the sink particles. In the first and second snapshots about 2 and 120 $M_{\odot}$ of gas have been accreted onto the sink particles.</figcaption></figure>

second sink particle that has formed in the simulation is displayed at three snapshots. On the second and third snapshots, the clustering is also clear. The object distribution is clearly hierarchical. A disk-like structure is seen in the second snapshot, and two objects appear to have formed as a consequence of disk fragmentation. Two objects have formed slightly farther away, and at least one of them is surrounded by a disk. Two more objects can be seen at time 0.104 Myr that are already decoupled from their gas reservoir. As expected, more objects weakly correlated to the gas appear in the last snapshot. It is worth reminding here that since the magnetic field is not included here, the disks are likely too large and fragment too easily (see, e.g., Wurster & Li 2018; Hennebelle et al. 2020; Zhao et al. 2020).

Temperature cuts through the yz plane are also given (right column). As can be seen, the central region around the sink presents temperatures that peak around 1500 K, which is typical of the temperatures reached at scales smaller than a few astronomical units in a collapsing protostar. At a distance of a few hundreds of astronomical units, the temperature drops below or becomes comparable to 100 K. Overall the temperature distribution remains less structured than the

column density, for instance. As time goes on, the mean temperature increases, which is due to the larger number of stars that have formed as quantified below.

## 3.2. Accretion as a Function of Time

As accreted mass onto sink/star is of primordial importance in these simulations, Figure 3 displays the total sink mass, $M_{*,tot}$, as a function of time for the COMP and STAN-type runs. As the initial densities in the two series of runs differ by almost two orders of magnitude (a factor of 64), the freefall time and therefore the accretion times differ by a factor of about 8–10. As can be seen, the accretion luminosity has only a modest influence on the global accretion, except initially for COMP-type runs, where we see that it almost stops accretion for a brief period of time before the first solar mass of gas has been accreted. As time goes on and after a few solar masses of gas is accreted, the mass ratio is roughly a factor of 2, and this ratio keeps decreasing with time. The difference between NOFEED and ACLUM runs for the more diffuse clumps (STAN-type) is weaker, with mass differences of only a few tens of percent. These curves constitute a first indication that the accretion luminosity is playing some role during the collapse of a massive clump without changing drastically the final result. This is extensively discussed below.

Also plotted is STANHM-ACLUM, which we remind has a Mach number of 10 initially, compared to 5 for STAN-ACLUM. Clearly, the higher turbulence modifies the accretion history. Star formation starts a bit earlier, and this is because some velocity fluctuations help compress the gas locally. However, globally the accretion rate is a bit lower (by 20%–30%), and this is because turbulence exerts some support on the cloud at large scales.

# 4. Temperature Distributions

The gas temperature is strongly influenced by the radiative feedback, and here we investigate its distribution within the clumps.

## 4.1. Temperature-Density Histogram

To get insight on the physical conditions that prevail within the simulated clouds, we now present in Figures 4 and 5 the mean temperature as a function of gas density for several time steps (which are more easily referenced by their accreted mass $M_{*,tot}$). This is obtained by simply computing the mass-weighted temperature in all density intervals. While the information carried by the mean temperature is incomplete, it is relatively simple, which facilitates the comparisons between runs. Bidimensional histograms, which contain much more detailed information, are given in Appendix A. The dotted lines visible in Figures 4 and 5 represent the analytical expression stated by Equation (1).

In all runs, we observe a change between an isothermal and a non-isothermal, adiabatic-like regime around $10^8-10^9 \text{ cm}^{-3}$. Note that the adiabatic regime appears to be only poorly described by the analytic functions. This clearly is a consequence of the heating that results from the emitted radiation. However, the discrepancy is amplified partly by the averaging procedure and partly due to insufficient resolution (see Appendix A). The nonmonotonic behaviors, in particular the bump located around $10^{10} \text{ cm}^{-3}$, are due to the averaging
<figure><img src="image_3.png" /></figure>

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

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

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

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

<figure><img src="image_8.png" /><figcaption>**Figure 2.** Column density (left) and temperature cuts (right) at three snapshots for run STAN-ACLUMhrhs around one of the sink particles.</figcaption></figure>

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

<figure><img src="image_10.png" /><figcaption>**Figure 3.** Accreted mass as a function of time for the various runs.</figcaption></figure>
<figure><img src="image_11.png" /></figure>

COMP-NOFEED

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

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

COMP-NOACLUM

<figure><img src="image_14.png" /><figcaption>**Figure 4.** Mass-weighted temperature in density intervals as a function of density at five time steps (corresponding to five accreted masses, $M_{*,tot}$) for the four COMP runs. In runs COMP-ACLUM and COMP-LOWACLUM, 50% and 10%, respectively, of the accretion luminosity is taken into account. In run COMP-NOFEED the stellar feedback is not taken into account, while only the stellar luminosity is considered in COMP-NOACLUM.</figcaption></figure>

procedure and to the presence of high-temperature gas as revealed in Figure 10.

As expected, the temperatures in run COMP-LOWACLUM are lower than the ones of run COMP-ACLUM, typically by a factor on the order of 1.5–2. Anticipating the analytical development made below, this is expected since the temperature typically varies like $L_{acc}^{1/3-1/4}$ and the luminosities of the two simulations differ by a factor of 5.

The comparison with run COMP-NOACLUM (which we remind takes into account the stellar luminosity but not the accretion one) shows that in a first phase ($M_{*,tot} < 100 M_{\odot}$) the temperature in run COMP-NOACLUM remains typically 3–4 times below the temperature of the runs that take the accretion luminosity into account. However, at later times, when more gas has been turned into stars, several stars more massive than a few solar masses formed, and the stellar luminosity leads to temperatures that are roughly only a factor of 2 below the ones of run COMP-LOWACLUM.

Finally, the bulk temperatures of run COMP-NOFEED remain low, typically around 20–30 K, even when $M_{*,tot}$ is larger than 200–300 $M_{\odot}$.

Overall the temperatures of the series of STAN-type runs are up to three times lower. For instance, in run STAN-ACLUM the temperature at low density is on the order of 50 K for $M_{*,tot} = 330 M_{\odot}$, while for run STAN-LOWACLUM it is roughly 30 K. This clearly is because (i) the cloud is more extended, so the distances from the sources are more important in STAN-type runs than in the COMP-type ones, and (ii) the

accretion rate is lower for the former than for the latter. This is quantified in the next section.

As for the COMP-NOFEED run, the temperature of the STAN-NOFEED run is significantly lower than when the accretion luminosity is taken into account; the largest temperature obtained at low density is only about 20 K.

To explore further how initial conditions influence the temperature distribution, the bottom right panel of Figure 5 presents run STANHmach-ACLUM, which initially has a Mach number of 10 instead of 5 for run STAN-ACLUM. The accreted mass remains below 200 $M_{\odot}$ because the clump is marginally bound. Comparing the temperatures of runs STAN-ACLUM and STANHmach-ACLUM when the same amount of mass has been accreted, we see that the temperatures are slightly lower for run STANHmach-ACLUM. This is because the accretion rate is lower in this run than in run STAN-ACLUM.

## 4.2. Analytical Developments: Predicting the Temperature Distribution

As temperature distribution plays an important role in the clump evolution with regard to both its fragmentation and its chemical composition, we provide here analytical estimates. More specifically, we will estimate here the clump mass per units of accreted mass, which lies above a certain temperature threshold chosen to be 100 K. The calculation entails several steps. First, we estimate the temperature profile of an envelope (with a density profile assumed to be $\propto r^{-2}$) around a source
<figure><img src="image_15.png" /></figure>

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

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

<figure><img src="image_18.png" /><figcaption>**Figure 5.** Same as Figure 4, but for the STAN runs. Clearly, when accretion luminosity is taken into account, temperatures are significantly lower than for the COMP runs. This is due to lower accretion rate but also larger distance, on average, from the source. The run STANHMich has a Mach number equal to 10 initially instead of 5 and is globally accreting significantly less than the runs for which it is equal to 5 initially.</figcaption></figure>

that is emitting a flux $f_{\text{acc}} GM_*dM/dt/R_*$. Second, we obtain the accretion rate for a source of mass $M_*$. Third, we choose (based on former studies) the mass spectrum of the stars. Finally, we perform an integration over the mass spectrum to get the heated mass of gas per units of accreted mass.

### 4.2.1. Temperature Distribution around a Single Source

We consider a spherically symmetric clump with a central source of mass $M_*$ accreting at a rate $\dot{M}_*$. The gas density is further assumed to be

$$ \rho(r) = \frac{\delta_{\rho} C_{s,0}^2}{2\pi Gr^2}, \qquad (4) $$

where $\delta_{\rho}$ is a dimensionless factor that typically is equal to 10–200 as discussed in Appendix C. We further assume that gas, dust, and radiation have the same temperature and are all stationary. A single radiation frequency is considered, and since the medium is optically thick, we have

$$ -4\pi r^2 \frac{c}{3\kappa(T)\rho(r)} \partial_r(aT^4) = f_{\text{acc}} \frac{GM_*\dot{M}_*}{R_*}. \qquad (5) $$

Following Semenov et al. (2003), we can distinguish two regimes of temperature,

$$ \begin{aligned} \kappa(T) &\approx 5 \text{ cm}^2 \text{ g}^{-1} \text{ for } T > T_{\text{crit}} \approx 100 \text{ K}, \\ \kappa(T) &\approx 5 \text{ cm}^2 \text{ g}^{-1} \left(\frac{T}{T_{\text{crit}}}\right)^{\alpha} \text{ for } T < T_{\text{crit}}, \end{aligned} \qquad (6) $$

where $\alpha$ is typically between 1 and 2. In this work we adopted $\alpha = 1.5$. Combining Equations (4)–(6), we get

$$ \begin{aligned} T(r) &= \left( T_{\text{crit}}^4 + K \left( \frac{1}{r^3} - \frac{1}{r_{\text{crit}}^3} \right) \right)^{1/4} \quad \text{for } T > T_{\text{crit}}, \\ T(r) &= \left( T_{\infty}^{4-\alpha} + K T_{\text{crit}}^{-\alpha} \frac{4-\alpha}{4} \frac{1}{r^3} \right)^{1/(4-\alpha)} \quad \text{for } T < T_{\text{crit}}, \end{aligned} \qquad (7) $$

where $T_{\infty}$ is the temperature at infinity and $T_{\infty} = T_0 = 10$ K initially,

$$ K = B_{\text{rad}} f_{\text{acc}} \delta_{\rho} C_{s,0}^2 M_* \dot{M}_*, \qquad (8) $$

$$ B_{\text{rad}} = \frac{3\kappa}{24\pi^2 R_* a c}, \qquad (9) $$

and $r_{\text{crit}}$ is the radius at which $T = T_{\text{crit}}$ and is given by

$$ r_{\text{crit}}^3 = K \frac{4 - \alpha}{4} \frac{T_{\text{crit}}^{-\alpha}}{T_{\text{crit}}^{4-\alpha} - T_{\infty}^{4-\alpha}} \approx K \frac{4 - \alpha}{4} T_{\text{crit}}^{-4}. \qquad (10) $$

The left panel of Figure 6 displays $T(r)$ for three cases corresponding roughly to a low-mass protostar ($M_* = 0.1 M_\odot$, $dM/dt = 10^{-5} M_\odot \text{ yr}^{-1}$, and $\delta_{\rho} = 10$; red line), a protostar with intermediate mass ($1 M_\odot$, $10^{-4} M_\odot \text{ yr}^{-1}$, $\delta_{\rho} = 30$; black line), and a more massive protostar ($5 M_\odot$, $10^{-3} M_\odot \text{ yr}^{-1}$, $\delta_{\rho} = 100$; blue line).
<figure><img src="image_19.png" /><figcaption>**Figure 6.** Left panel: temperature stated by Equation (7) for three cases corresponding to three values of $M_*$, $dM/dt$, and $\delta_\rho$, namely, 0.1 $M_\odot$, $10^{-5} M_\odot$ yr$^{-1}$, 10 (red line); 1 $M_\odot$, $10^{-4} M_\odot$ yr$^{-1}$, 30 (black line); and 5 $M_\odot$, $10^{-3} M_\odot$ yr$^{-1}$, 100 (blue line). Right panel: mass above 100 K per mass of gas turned into stars for the various runs. The dotted lines correspond to the analytical expression stated by Equation (20), with values of $R_c$, $f_{acc}$, and $\delta_\rho$ that correspond to each specific run. While the values of $R_c$ and $f_{acc}$ are specified, $\delta_\rho$ (stated by Equation (4)) is measured in the simulations. As can be seen, the agreement between the analytical expression and values inferred from the simulation is entirely reasonable.</figcaption></figure>

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

### 4.2.2. Accretion Rate

To estimate the accretion rate on each star, we proceed like in Lee & Hennebelle (2018a), who have estimated it to be $\dot{M} \simeq M/\tau_{ff}$, where M is the mass of the reservoir from which the star is building its mass and $\tau_{ff}$ is the associated freefall time, $\sqrt{3\pi/(32G\rho)}$. To get the accretion reservoir, for simplicity we assume that its mass is nearly one of the star (for instance, jets that may change this efficiency are not considered in the present work), meaning that all mass losses are neglected, while the reservoir radius is determined by the virial theorem, which leads to

$$M = \frac{\pi^{5/2} \left[ (C_s)^2 + (\sigma_c^2/3)(R/R_c)^{2\eta} \right]^{3/2}}{6 \sqrt{G^3 \rho}}, \quad (11)$$

where $\sigma_c$ is the velocity dispersion at the cloud scale, $R_c$ is the clump radius, and $\eta \simeq 0.5$ is the exponent through which the velocity dispersion varies with spatial scale, $\sigma \propto R^\eta$ as expected for a turbulent fluid. This leads to

$$\tau_{ff} = \sqrt{\frac{3}{2}} \pi^{-1/4} \frac{R}{((C_s)^2 + (\sigma_c^2/3)(R/R_c)^{2\eta})^{1/2}}. \quad (12)$$

Combining Equations (11) and (12), we have the freefall time and therefore the accretion rate as a function of the mass.

To get a further physical hint, it is worth simplifying this expression, which can be achieved by neglecting the sound speed with respect to the turbulent dispersion in Equations (11) and (12). This leads to

$$\begin{aligned} \tau_{ff} &= \frac{3}{\sqrt{2}} \pi^{-1/4} \sigma_c^{-1} R \left(\frac{R}{R_c}\right)^{-1/2} \\ &= K_{ff} G^{1/4} \sigma_c^{-3/2} R_c^{3/4} M^{1/4}, \end{aligned} \quad (13)$$

where we have assumed $\eta = 0.5$ and where $K_{ff} = 3\pi^{-1/4}(6\pi)^{1/3}/\sqrt{2}$. This expression, which will be used later in the final expression of the heated gas mass per mass of stars, $f_M$, implies that only the contribution of the most massive stars will be accurately represented. However, by comparing the two estimates inferred from Equations (12) and (13), we found that this is a valid approximation.

### 4.2.3. Source Distribution

To get the temperature distribution inside the clouds, we need to know the source distribution. We assume that the mass spectrum is given by

$$\mathcal{N} = \frac{dN}{d \log M} = A M^{-\beta}. \quad (14)$$

This mass spectrum applies between a minimum and maximum mass, $M_{min}$ and $M_{max}$. While the former is typically equal to 0.3 $M_\odot$, which corresponds to the peak of the IMF, the latter increases with the total accreted mass, $M_{tot,*}$, and is equal to a few solar masses. Obviously this is a simplification since one should sum over the full mass spectrum. However, as seen below, the dependence on $M_{min}$ is quite shallow, and its exact value is not really consequential. We have

$$\int_{M_{\min}}^{M_{\max}} M \frac{dN}{d \log M} = \frac{A}{-\beta + 1} (M_{\max}^{-\beta+1} - M_{\min}^{-\beta+1}) = M_{\text{tot},*}, \quad (15)$$

which leads to

$$A = \frac{(-\beta + 1)}{(M_{\max}^{-\beta+1} - M_{\min}^{-\beta+1})} M_{\text{tot},*}. \quad (16)$$

### 4.2.4. Heated Mass Fraction

The mass enclosed in the sphere of radius $r_{crit}$ is given by

$$m_{\text{crit}} = \int_{0}^{r_{\text{crit}}} 4\pi r^2 \delta_{\rho} \frac{C_{s,0}^2}{2\pi Gr^2} dr \approx \frac{2\delta_{\rho} C_{s,0}^2}{G} r_{\text{crit}}. \quad (17)$$

Thus, the total mass heated above $T_{crit}$ can be estimated as

$$\begin{aligned} M_{\text{crit}} &= \int_{m_{\min}}^{m_{\max}} m_{\text{crit}} \mathcal{N} d \log M \\ &= \frac{2 f_{\text{acc}}^{1/3} \delta_{\rho}^{4/3} C_{s,0}^{8/3}}{G T_{\text{crit}}^{4/3}} A B_{\text{rad}}^{1/3} \left(\frac{4-\alpha}{4}\right)^{1/3} \\ &\quad \times \int_{m_{\min}}^{m_{\max}} \left(\frac{M^2}{\tau_{\text{ff}}(M)}\right)^{1/3} M^{-\beta} d \log M, \end{aligned} \quad (18)$$
which is the expression that we will use below to confront with the simulation results. Note at this stage that a difficulty arises regarding the choice of $\beta$, the exponent of the stellar mass spectrum (as stated by Equation (14)). Most observations found that above a mass of a few solar masses the exponent is close to 1.3, which is the value originally inferred by Salpeter. However, as explained below and in Lee & Hennebelle (2018a), the mass spectra obtained in numerical simulations of massive collapsing clumps tend to be a bit flatter and are more accurately described by an exponent of 3/4–1 (followed by an exponential cutoff at the highest masses). Therefore, for the purpose of comparing with the numerical results, we will use from this point the value of 3/4, although observationally it would be more logical to use the value of 1.3. Fortunately, it makes little difference, as the results with the two values of $\beta$ vary by a few tens of percent, which is far below the expected accuracy of our analytical approach.

As it is useful to obtain a simpler expression, we use the simplified expression for the freefall time as stated by Equation (13). Assuming further that $\eta = 0.5$ and $\beta = 3/4$, we then get

$$f_{M,crit} = \frac{M_{crit}}{M_{tot,*}} = 3^{2/3} K_{ff}^{-1/3} (4 - \alpha)^{1/3} f_{acc}^{1/3} \\ \times \frac{B_{\text{rad}}^{1/3}}{G^{13/12}} \frac{M_{\text{min}}^{-1/6} - M_{\text{max}}^{-1/6}}{M_{\text{max}}^{1/4} - M_{\text{min}}^{1/4}} \frac{\delta_{\rho}^{4/3} C_{s,0}^{8/3}}{T_{\text{crit}}^{4/3}} \sigma_{c}^{1/2} R_{c}^{-1/4}, \quad (19)$$

where $f_{M,crit}$ is the mass of gas having a temperature larger than $T_{\text{crit}}$ per units of accreted mass. As we see, it weakly depends on the minimum and maximum stellar mass present in the sample. It depends on the clump physical conditions through the radius $R_c$, the velocity dispersion $\sigma_c$, and the overdensity $\delta_{\rho}$. The values of the latter vary with the clump parameters and are typically on the order of 10–100. Through $B_{\text{rad}}$ and Equation (9), we see that $f_{M,crit} \propto f_{acc}^{1/3}$. Thus, variations of $f_{acc}$ induce limited changes of $f_{M,crit}$.

If we further assume that $\sigma_c = \sqrt{GM_c/R_c}$, which is close to the chosen initial value that simply reflects energy equipartition, we obtain

$$f_{M,crit} = \frac{M_{crit}}{M_{tot,*}} \approx 10^{-2} \times \delta_{\rho}^{4/3} f_{acc}^{1/3} \left( \frac{M_c}{1000 M_{\odot}} \right)^{1/4} \left( \frac{R_c}{1 \text{ pc}} \right)^{-1/2}. \quad (20)$$

### 4.2.5. Feedback Efficiency: Comparisons between Simulations and Theory

Figure 6 portrays the values of $f_{M,crit}$, i.e., the mass of gas having a temperature higher than $T_{\text{crit}}$ per unit of accreted mass, for the various runs performed as a function of the accreted mass, $M_{\text{tot,*}}$. Solid lines represent the COMP-type runs, while dashed lines display the STAN-type ones. The dotted lines represent the analytical models where the parameters entering in Equation (20), namely, $R_c$ and $f_{acc}$, are taken from Table 1. On the other hand, the parameter $\delta_{\rho}$ is estimated from the simulations (see Appendix C). For COMP-type runs we have $\delta_{\rho} \simeq 30–300$, and $\delta_{\rho} \simeq 10–100$ for the STAN-type ones. To perform the calculations of the models in Figure 6, we have

used the values $\delta_{\rho} = 150$ and 50, respectively. Note that from Figure 15 it is seen that $\delta_{\rho}$ tends to increase with $M_{*}$; the dependence on $M_{\text{max}}$ stated in Equation (20) may actually be underestimated.

As expected for the runs where $f_{acc} \neq 0$ (i.e., ACLUM and LOWACLUM-type runs), $f_M$ is almost independent of $M_{\text{tot,*}}$. The change of slope at 30–100 $M_{\odot}$ for the COMP-ACLUM runs is due to the fact that all the gas is warm in the computational box, and so $M_{\text{crit}}$ does not increase while $M_{\text{tot,*}}$, the mass within the sinks, keeps increasing. The agreement between the model (dotted lines) and the simulations (solid and dashed lines) is reasonable. In particular, the trends are well reproduced.

Typically we have $f_M \simeq 30$ for the COMP-ACLUM run, while $f_M \simeq 3$ for the STAN-ACLUM run. Altogether, a star of mass $M_{*}$ is able to hit above $T_{\text{crit}}$ almost 10 times more gas in the compact cloud than in the diffuse one. This is a consequence of the density, which is 3–4 times larger for the former compared to the latter. As expected, $f_M$ is considerably smaller, about a factor of 10 in runs where there is no accretion luminosity ($f_{acc} = 0$), than in run COMP-ACLUM. We also note that $f_M$ is a decreasing function of $M_{\text{tot,*}}$ for runs for which $f_{acc} = 0$.

We, therefore, conclude that Equation (20), which gives the expression of $f_M$, is accurate within a factor of a few and reproduces the qualitative behavior observed in the simulations.

### 4.3. Qualitative Comparison with Observations

As discussed above, the temperature distribution reflects the gas mass distribution and the evolution of the clumps, i.e., the fraction of gas that has been converted into stars assuming, as mentioned previously, that all the mass of the star is about the mass of the reservoir. Therefore, comparing with observations is not an easy task, as these quantities are generally poorly known. There are also various techniques, such as spectral energy distribution (SED) fitting and molecular spectroscopy, that provide different results depending on which region of the clump is actually probed (see, e.g., Figure 11 of Giannetti et al. 2017).

Both the ATLASGAL (Urquhart et al. 2014) and Hi-GAL surveys (Elia et al. 2017) provide mean temperature distributions. Looking, for instance, at Figure 5 of Elia et al. (2017), we see that the temperatures of protostellar sources are higher than those of the prestellar clumps, indicating internal heating. However, the peak of the distribution is about 13 K, and only a few protostellar sources present temperatures above 30 K. Since the sample contains both massive and compact clumps, comparable to the ones simulated here, this may place constraints on the effective $f_{acc}$, although the SED fitting has been restricted to temperatures below 40 K. Similar numbers are provided by Urquhart et al. (2014) (Figure 10), where the NH$_3$ molecules have been used, although sources with temperatures larger than 45 K have been discarded. Using other molecular tracers (such as CH$_3$OH, for instance), Giannetti et al. (2017) infer temperatures for massive star-forming clumps selected for the TOP100 sample (Csengeri et al. 2016). Temperatures as high as a few hundreds of kelvin are reported, but this may correspond to the inner part of the clumps. Indeed, the temperature spatial distribution is a clue to assess the importance of thermal feedback.

In this respect, an interesting set of observations has been undertaken by Ginsburg et al. (2017), who mapped several massive star-forming clumps and inferred temperatures using a rotational diagram of CH$_3$OH (see also Figure 3 of Motte et al. 2018, where temperature above 60–80 K are obtained for massive
<figure><img src="image_21.png" /></figure>

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

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

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

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

<figure><img src="image_26.png" /><figcaption>**Figure 7.** Mass spectra at various times for COMP-type runs. We remind that COMP-bar1 and COMP-bar2 use a barotropic EOS and no radiative transfer. The total mass accreted by the sink particles (expressed in solar mass) is indicated in the legend. All distributions but the ones of COMP-ACLUM peak at about 0.2–0.3 $M_\odot$ and have generally a very similar shape. This clearly demonstrates that up to the time investigated here stellar feedback has a very limited impact on the mass spectrum. What is more important is the accretion luminosity, but it is subject to uncertainties (see text).</figcaption></figure>

clumps at scales of several thousands of astronomical units). The data of Ginsburg et al. (2017) reveal temperatures exceeding 100 K extending up to 5000 au. For instance, Figure 6 of Ginsburg et al. (2017) shows for clump e2 temperatures of 100–200 K at distances larger than $10^4$ au from the center of the source. While the mass in the region around e2 is estimated to be on the order of $10^4 M_\odot$, it contains about 500 $M_\odot$ within the central $10^4$ au. Since it contains massive stars and a total stellar mass higher than $50 M_\odot$ (A Ginsburg 2020, private communication; Goddi et al. 2018), it is broadly comparable to our COMP-type clumps at an age where at least $50 M_\odot$ have been accreted. This may be comparable (within a factor of 2–3) with what has been inferred not only for run COMP-ACLUM but

also for run COMP-LOWACLUM. At this stage, because of the broad uncertainties, it does not seem possible to draw strong conclusions, and this remains a challenge for future studies.

## 5. Mass Spectra

We now turn to the stellar mass spectra formed in these calculations.

### 5.1. Results

Figures 7 and 8 portray the mass distribution of the sink particles, respectively, for the six COMP-type runs and four STAN-type runs listed in Table 1 (10 first runs). To follow the
<figure><img src="image_27.png" /><figcaption>**Figure 8.** Mass spectra at various times for STAN-type runs. All distributions peak at about 0.2–0.3 $M_\odot$ and have generally a very similar shape. Only run STAN-ACLUM, which considers the accretion luminosity, is broader in shape. This shows that for moderately dense clumps, radiative feedback, even when it includes high accretion luminosity, does not strongly modify the shape of the mass spectrum.</figcaption></figure>

evolution, the mass spectra are shown at various time steps, which correspond to various amounts of accreted mass. The mass spectra are complemented by Figure 9, which portrays for six of the runs the number of objects formed and the mass of the most massive object as a function of $M_{\text{tot,\*}}$, the total accreted mass.

Let us start with run COMP-NOFEED, which we remind has no stellar feedback and no accretion luminosity. The mass spectra present a clear peak around $0.3 M_\odot$ at early times, i.e., when less than $100 M_\odot$ have been accreted. At later times, the peak broadens and shifts toward $\simeq 0.5 M_\odot$. The high-mass part presents a power-law-like shape with an exponent around $\simeq 1$. This behavior is very similar to several mass spectra published in the literature using either a barotropic EOS (e.g., Bate et al. 2003; Lee & Hennebelle 2018a)<fnref n="6" /> or radiative transfer calculations but no accretion luminosity (e.g., Bate 2009, 2012). The run COMP-NOACLUM (which considers stellar feedback) presents a very similar behavior, although at late times the power-law behavior for the high mass is better defined. Interestingly, we see that the run COMP-bar1, which has an EOS that broadly reproduces the density–temperature relation of run COMP-NOFEED at the transition point between the isothermal and adiabatic regimes

(that is, for $n \simeq 10^{10}-10^{11} \text{ cm}^{-3}$), presents mass spectra that are very similar, with a peak located at $0.3 M_\odot$ and a power-law-like mass spectrum at high mass. In run COMP-bar2, for which the transition from isothermal to adiabatic occurs at slightly lower density, the peak occurs at slightly larger mass. Indeed, if the gas becomes adiabatic at lower density, then more gas piles up before a density of $10^{13} \text{ cm}^{-3}$ is reached and the sink particle is being introduced. Physically, this would correspond to a more massive FHSC. Therefore, the sink has more gas to accrete, and the peak of the IMF is shifted toward larger masses.

The inclusion of the accretion luminosity (run COMP-ACLUM) leads to significant differences. We recall that in these runs $f_{\text{acc}} = 0.5$ for the accretion luminosity. First, at early times when $M_{\text{tot,\*}} = 50 M_\odot$, the peak of the distribution is located at about $0.07 M_\odot$. The reason is that this is precisely the mass at which the accretion luminosity starts being applied. This choice, which simply corresponds to roughly two times the mass of the first hydrostatic core, is somehow arbitrary, and therefore this feature remains rather uncertain. The accretion luminosity feedback is so strong that when it is taken into account, further accretion is abruptly stopped for a moment. This is also very clearly seen in Figure 9, which reveals (top panel) that when $M_{\text{tot,\*}} \simeq 3 M_\odot$ the mass of the most massive object remains up to a factor of 5 below its value in run COMP-NOFEED. Correspondingly, the number of objects for a given value of $M_{\text{tot,\*}}$ is initially the largest in run COMP-ACLUM. As time goes on, the sinks eventually grow in mass, and their

<footnote n="6">⁶ Note that most published calculations, which use a barotropic EOS, generally have mass spectra that peak at mass smaller than $0.3 M_\odot$. This is simply a consequence of the chosen EOS. Some of the simulations presented in Lee & Hennebelle (2018b) and in this work use a barotropic EOS but nevertheless present a peak of the mass spectrum near $0.3 M_\odot$.</footnote>
<figure><img src="image_28.png" /><figcaption>**Figure 9.** Largest sink mass (top panel) and number of sinks, $N_{\text{tot}}$ (bottom panel), as a function of the total accreted mass, $M_{\text{tot}}$, in several runs.</figcaption></figure>

numbers, as well as the largest sink mass, become similar to the other runs (COMP-NOFEED and COMP-LOWACLUM) when $M_{\text{tot},*} \simeq 30 M_{\odot}$. The peak at $0.07 M_{\odot}$ is not visible any more when $M_{\text{tot},*} = 300 M_{\odot}$. Instead, a weak peak located around $0.5 M_{\odot}$ has developed. The mass spectrum tends, however, to be nearly flat between 0.07 and $3 M_{\odot}$. Flat mass spectra (in $dN/d\log M$) are typical of clumps for which the thermal support is very high, so that at the scale of the mass reservoir of an individual sink/star, it dominates over the turbulent dispersion (see the discussion in Lee & Hennebelle 2018a). Indeed, Figure 4 clearly shows that the thermal support is very large and up to 20 times its value in an isothermal run. Let us stress that a similar effect is also probably present in the calculations of Krumholz et al. (2012), where a flat mass spectrum is also inferred.

When accretion luminosity is present but weaker (COMP-ACLUMLOW, which we remind has $f_{\text{acc}} = 0.1$), the peak at $0.07 M_{\odot}$ is much less pronounced, and it quickly disappears. A pronounced peak located at $0.3-0.5 M_{\odot}$ appears, and when $M_{\text{tot},*} = 300 M_{\odot}$, a relatively clear mass spectrum, whose slope is around -1 to 1.2, has developed.

The STAN-NOFEED runs (Figure 4) present mass spectra that resemble the COMP-NOFEED ones, with a peak around $0.5 M_{\odot}$. This is in good agreement with what has been concluded in Lee & Hennebelle (2018a), that beyond a certain threshold of initial

density the peak of the mass spectrum does not depend on the initial density, illustrating in particular that it is not related to the initial Jeans mass (see also Lee & Hennebelle 2019).

The influence of the accretion luminosity on the mass spectrum is significantly less pronounced than for the COMP-type runs. The main difference between runs STAN-ACLUM and STAN-NOFEED is that the mass spectrum is a bit broader for run STAN-ACLUM. There are more objects of mass $0.07 M_{\odot}$ for STAN-ACLUM than for STAN-NOFEED. This is certainly due to the sudden introduction of the accretion luminosity at $0.07 M_{\odot}$. This does not lead, however, to a pronounced peak as for run COMP-ACLUM. This is likely because the accretion rate is much lower in the STAN-type runs than for the COMP-type ones. As expected, the impact of the accretion luminosity is even smaller for runs STAN-LOW-ACLUM, which presents mass spectra that are very similar to the STAN-NOFEED ones. The moderate influence of the accretion luminosity is also visible in Figure 9, where it is seen that both the most massive star and the number of stars are very similar between all STAN-type runs. The most significant difference is found at late times, where the largest mass is about $10 M_{\odot}$, roughly 2 times larger for run STAN-ACLUM than for run STAN-NOFEED, for instance. Clearly, the increase of temperature favors the growth of existing protostars by reducing the amount of fragmentation.

Finally, it is interesting to see that run STANHMACH-ACLUM presents mass spectra that peak at smaller mass than run STAN-ACLUM and remains slightly narrower. This is likely because the higher turbulence makes the accretion rate lower and the clump radius larger. Both effects reduce the influence of the accretion luminosity.

## 5.2. Interpretation and Discussion

By performing numerical simulations of a massive collapsing clump with a barotropic EOS and comparing with an analytical model, Lee & Hennebelle (2018a) have identified two regimes resulting in two different mass spectra. If the thermal support is high, that is, if the initial density is low, or equivalently if the temperature is high, the mass spectrum tends to be flat, that is, $dN/\log M \propto M^0$ (see run A of Lee & Hennebelle 2018a and runs presented in bottom panels of Figure 2 of Jones & Bate 2018). When, on the other hand, turbulent support dominates, the mass spectrum presents a peak at about 10 times the mass of the FHSC and a power law at high mass that is found to be $dN/\log M \propto M^{-3/4}$ (see runs B, C, and D of Lee & Hennebelle 2018a). Given the similarity found between the mass spectra obtained with radiative transfer (except for run COMP-ACLUM) and those obtained with the barotropic EOS, it seems likely that the physical interpretation developed in Lee & Hennebelle (2018a, 2018b) for the peak of the IMF, namely, the impact of tidal effects, remains valid when radiative transfer is considered.

The different behavior obtained for run COMP-ACLUM, for which the mass spectrum is nearly flat for mass between $\simeq 0.3$ and $\simeq 3 M_{\odot}$, is most likely a consequence of the high thermal support induced by the large temperatures due to the high accretion luminosities. As mentioned earlier, however, the *effective* value of $f_{\text{acc}}$ for real protostars is not well established. Indeed, the observed luminosities of protostars are much fainter than the values corresponding to such accretion luminosities (the luminosity problem; Evans et al. 2009). A plausible explanation to the luminosity problem is episodic accretion
(see, e.g., Baraffe et al. 2009, 2012), but this issue is not entirely settled yet. If this indeed happens, it means that the heating of the gas through accretion luminosity may be intermittent as well, and so the *effective* value of $f_{\text{acc}}$, i.e., the value that the clump is experiencing most of the time, may be on average small. Interestingly, Stamatellos et al. (2011) developed a model of intermittent accretion to study disk fragmentation. By comparing simulations with and without accretion luminosity, they concluded that the simulations with episodic accretion resemble the ones without accretion luminosity. While episodic accretion received strong support from observations (Evans et al. 2009) and theory (Baraffe et al. 2012), it is currently unclear whether the proposed mechanism of gravitational instability (Vorobyov & Basu 2010) is sufficiently universal. Another possible effect that could lead to small effective $f_{\text{acc}}$ has been proposed by Krumholz et al. (2012), who claimed that most of the radiation can escape through the wind cavities. It is therefore unclear whether the impact of the accretion luminosity is as found in run COMP-ACLUM. It should also be stressed that the physical conditions corresponding to COMP-type runs are not typical of most star-forming clumps. In this respect, STAN-type runs are more typical. Since there the accretion luminosity appears to have an impact that is more limited, in particular regarding the peak, it seems that, in most star-forming clumps of the Milky Way, the accretion luminosity is not drastically influencing the stellar spectrum.

## 6. Conclusions

We have performed a series of numerical simulations with a spatial resolution of a few astronomical units (and up to 1 au for convergence runs), to investigate the influence of radiative feedback on the mass distribution of stars that form during the collapse of a 1000 $M_{\odot}$ clump. We also performed two runs in which a barotropic EOS is employed. Two types of initial conditions have been explored: one corresponds to a very compact clump initially (with a radius of 0.1 pc), while the other is more typical of Milky Way star-forming clumps (radius 0.4 pc initially).

We found that as long as accretion luminosity is not considered, the stellar mass spectra that form in the various runs, both with radiative transfer and with an effective EOS, present strong similarities and resemble the observed IMF. This suggests that in this case radiative transfer is not fundamental in setting the IMF, even though the gas temperature of the dense star-forming gas is 3–10 times higher than when an EOS is employed.

When accretion luminosity is included, its impact depends on the initial conditions. For the case of the very compact clump, the mass spectrum is initially strongly peaked at small mass. As time goes on, it progressively becomes flat between 0.1 and $3 M_{\odot}$. For the case of the less compact clump (0.4 pc of radii), the effect remains more limited, particularly at late times. The peak is located at roughly the same value as in the absence of accretion luminosity. An interesting effect, though, is that there are few more small and massive objects. While the latter are a consequence of the strong heating induced by the accretion luminosity which prevent fragmentation, the former

are also a consequence of this heating, which, when a low-mass object has just been created, prevent further accretion. We stress that the exact history of accretion luminosity, in particular the mass at which it actually starts, is largely unknown, and therefore the increase of low-mass objects we found must be regarded with care.

We conclude that the accretion luminosity, if its effective value is equal to or at least comparable with the gravitational energy released at the surface of the star, is expected in most galactic star formation clumps to have an impact that mainly consists in producing both smaller and more massive objects than what would have been formed otherwise. In particular, it is likely that in most circumstances the peak of the distribution is a consequence of the change of the effective EOS that is responsible for the first hydrostatic core, rather than due to the feedback heating of the collapsing envelope. It seems, however, that feedback heating leads to the formation of more massive stars (Krumholz et al. 2007; Urban et al. 2010), which otherwise appear to be seldom.

P.H. warmly thanks Timea Csengeri, Neal Evans, Adam Ginsburg, Fabien Louvet, Anaëlle Maury, Sergio Molinari, Frédérique Motte, and Alessio Traficante for enlighting discussions on the temperature interpretations of massive star-forming cores and the issue of accretion luminosity. We thank the anonymous referee for their constructive comments that have improved the paper. P.H. acknowledges financial support from the European Research Council (ERC) via the ERC Synergy Grant ECOGAL (grant 855130).

## Appendix A
### Bidimensional Density–Temperature Histograms

To get more hints on the physical conditions within our modeled clumps, we present here bidimensional temperature–density histograms for the two most extreme runs of the COMP type, namely, COMP-ACLUM and COMP-NOFEED, as well as for the most resolved one, COMP-NOACLUMhr, at three different snapshots. The results are displayed in Figure 10. For completeness, we also show bidimensional temperature–density histograms for runs STAN-ACLUM, STAN-ACLUMhrhs, and STAN-ACLUM-vhrhs in Figure 11.

Bidimensional histograms contain more information than the mean temperature as a function of density presented in Figures 4 and 5. The bidimensional histograms reveal that there is a significant spread in temperature, even regarding the gas within a narrow density range. However, the bulk of the gas tends to lay in regions that appear better defined with significantly weaker dispersion.

By the time of the first snapshot, a few solar masses have been accreted. In the three runs, the temperature distributions are similar; as expected, there is a clear transition between the isothermal and adiabatic regimes. The transition itself is reasonably well described by the barotropic EOS calculations. At higher density, the temperature is clearly higher than the EOS values. However, this is essentially a consequence of insufficient resolution (which we recall is ≈2 au or 1 au for COMP-NOACLUMhr). Indeed, the temperature distribution is
<figure><img src="image_29.png" /><figcaption>**Figure 10.** Bidimensional temperature–density histograms as a function of three time steps for run COMP-ACLUM ($R = 0.1$ pc initially and accretion luminosity is taken into account) and run COMP-NOFEED (no stellar feedback). The left and middle columns are for early times when a few tens of solar masses have been accreted. The right column is for later times when about 200 solar masses have been accreted. The two curves represent the two EOSs as stated by Equation (1).</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>**Figure 11.** Bidimensional temperature–density histograms as a function of three time steps for runs STAN-ACLUM ($R = 0.4$ pc initially and accretion luminosity is taken into account), STAN-ACLUMhrhs, and STAN-ACLUMvhrhs. Comparison between runs at similar times allows us to see the influence of numerical resolution at high density. The two curves represent the two EOSs as stated by Equation (1).</figcaption></figure>

closer to the analytical expression in run COMP-NOACLUMhr, which has more resolution.

## Appendix B
Dependence on Numerical Parameters

To investigate the issue of numerical convergence, several runs have been performed as indicated at the end of Table 1. First, we studied COMP-type runs with no accretion luminosity, and then STAN-type ones taking the accretion luminosity into account.

The corresponding mass spectra for COMP-type runs are depicted in Figure 12. The top left panel reproduces for convenience the result of run COMP-NOACLUM. The bottom left panel presents the run COMP-NOACLUMrls, which has a resolution of 4.6 au, compared to 2.3 au for run COMP-NOACLUM; also sinks get introduced at $n_{\text{acc}} = 10^{12} \text{ cm}^{-3}$ instead of $n_{\text{acc}} = 10^{13}$ for run COMP-NOACLUM. Thus, the mass $dx^3n_{\text{acc}}$, i.e., the mass contained in the finest computational cells, is nearly the same in both runs. Clearly, the peak is located nearly at the same position. There are nevertheless

some differences between the two runs. COMP-NOACLUM has about three times more small objects ($M_* < 0.03 M_\odot$) and roughly three times more big ones ($M_* > 3 M_\odot$). The discrepancy between runs COMP-NOACLUM and COMP-NOACLUMls is even larger. Both runs have the same spatial resolution, but $n_{\text{acc}}$ is 10 times lower in COMP-NOACLUMls, meaning that the sinks are introduced more easily. The peak in COMP-NOACLUMls is located at about $0.1 M_\odot$ instead of $0.3 M_\odot$, and the number of small objects is even larger. Therefore, we see that the mass spectrum is influenced by both resolution and sink threshold. Run COMP-NOACLUMhr explores the influence of further numerical resolution but the same $n_{\text{acc}}$ as COMP-ACLUM. It shows that the mass spectrum shifts toward smaller masses but by a factor of less than 2. In particular, while run COMP-NOACLUMhr has the same value of $dx^3n_{\text{acc}}$ as run COMP-NOACLUMls, we find that the peak does not shift to a much smaller value, as is the case for run COMP-NOACLUMls. We interpret this as being due to the fact that at density $n_{\text{acc}} = 10^{13} \text{ cm}^{-3}$ the gas is nearly adiabatic. This shows that although complete numerical convergence may have not been completely reached, run COMP-NOACLUMhr
<figure><img src="image_31.png" /><figcaption>COMP-NOACLUMIrs</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>COMP-NOACLUMhr</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>COMP-NOACLUMIs</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>**Figure 12.** Mass spectra for the COMP-type runs and no accretion luminosity at various times and for various numerical resolution and sink parameters (see Table 1).</figcaption></figure>

is probably approaching it. Indeed, contrary to the isothermal regime, where the number of fragments increases with resolution, this is not the case in the adiabatic one.

Figure 13 displays the mass spectra for STAN-type runs. The top panel reproduces run STAN-ACLUM to ease the comparison process. The middle panel presents run STAN-ACLUMhrhs, which has two times more resolution and a value of $n_{\text{acc}}$ that is 10 times higher, leading to roughly the same value of $dx^3n_{\text{acc}}$ when the sink particles get introduced. As can be seen, the agreement is only moderate. There is a tendency for run STAN-ACLUMhrhs to have more massive objects by a factor of about $\simeq 2$. Run STAN-ACLUMvhrhs has a resolution of 1.15 au, which is two times higher than for run STAN-ACLUMhrhs. Both runs have the same $n_{\text{acc}}$. We see that overall the two distributions, without being identical, are nevertheless similar.

Let us stress that in Lee & Hennebelle (2018a, 2018b) systematic investigations of numerical convergence and dependency on $n_{\text{acc}}$ have been performed, and it has been concluded that while numerical convergence could not be achieved when an isothermal equation was used, convergence was achieved when a barotropic one was used. That is to say, when the adiabatic exponent becomes larger than 4/3 above a certain density, numerical convergence was obtained. Regarding the value of $n_{\text{acc}}$, it has been found that if the EOS has an exponent that at very high density is close enough to 4/3 (for instance, 7/5 but not 5/3), then the value of $n_{\text{acc}}$ is not *too* consequential. The situation in the present paper appears to be more difficult, as we found dependence both on numerical resolution and on $n_{\text{acc}}$. The reason is probably that while in barotropic calculations the EOS is imposed irrespectively of the resolution, the fully radiative calculation runs are meant as describing self-consistently the thermal state of the
STAN-ACLUM

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

STAN-ACLUMhrhs

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

STAN-ACLUMvhrhs

<figure><img src="image_37.png" /><figcaption>**Figure 13.** Mass spectra for the STAN-type runs with accretion luminosity at various times and for three numerical resolution values and two values of $n_{\text{acc}}$ (see Table 1).</figcaption></figure>

gas. In particular, the structure of the first hydrostatic core, which is argued to set the peak of the IMF (Lee & Hennebelle 2018b; Hennebelle et al. 2019), is expected to be self-consistently calculated. However, the size of the first hydrostatic cores is about 5 au. It is therefore relatively unsurprising that simulations with only a few astronomical units of resolution have not reached full numerical convergence yet. The similarities between, on one hand, runs COMP-NOACLUM and COMP-NOACLUMhr and, on the other hand, runs STAN-ACLUM, STAN-ACLUMhrhs, and STAN-ACLUMvhrhs also suggest that our simulations are probably approaching convergence.

## Appendix C: The Mean Density Profile

To get an estimate of the $\delta$ parameter that appears in Equation (20), we proceed like in Hennebelle et al. (2019), i.e., we measure the mean density in concentric shells around sink particles. We then compute the mean density value and the standard deviation (shaded area). The result is displayed in Figure 14. The density field is typically $\propto r^{-2}$ and in the case of COMP-ACLUM is about 100 times above the density of the singular isothermal sphere (blue line). For run STAN-ACLUM, it is more on the order of 20–30, which is roughly 4 times lower than for COMP-ACLUM. This is likely a direct consequence of their respective initial radii, which precisely differ by a factor of 4.

As the singular isothermal sphere density is $\propto C_s^2$, it is worth investigating the density profile of run COMP-NOACLUM since its temperatures are factors of 3–5 lower than the ones of COMP-ACLUM. The result is displayed in the bottom left panel, which reveals that the density distribution of run COMP-NOACLUM is very comparable to the one of COMP-ACLUM. The reason is that this is likely the turbulence that is playing here the role of an effective sound speed (Murray & Chang 2015).

Finally, it is also worth displaying the density distribution of run COMP-ACLUMls, as we saw in Appendix B that this run (which we remind uses $n_{\text{acc}} = 10^{12} \text{ cm}^{-3}$) has many more sink particles than run COMP-ACLUM. The density field is typically a factor of nearly 3–5 below that of run COMP-ACLUM. This is probably a consequence of the fact that there are more numerous sinks, meaning that the density field around a given object is partially accreted by the numerous neighbors. Also, the sinks are about three times less massive on average, and as seen in Figure 15, the density around an object increases with its mass.

To get a better estimate of the parameter $\delta_{\rho}$, we have plotted for each sink and at various time steps the mean value of $\delta_{\rho}$, obtained by taking the mean of the density divided by the singular isothermal sphere density, as a function of the sink mass. The results are displayed in Figure 15 for runs COMP-ACLUM and STAN-ACLUM. The two distributions span almost an order of magnitude in the range of 20–300 and 10–100, respectively. There is a trend for $\delta_{\rho}$ to increase with the sink mass.
<figure><img src="image_38.png" /><figcaption>**Figure 14.** Mean density profile around sink particles (red lines) and standard deviation (shaded area). The blue line is the density profile of the singular isothermal sphere.</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>**Figure 15.** Distribution of $\delta_{\rho}$ (ratio between density and singular isothermal sphere density around sink particles). Each point corresponds to a sink particle.</figcaption></figure>

ORCID iDs

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

Baraffe, I., Chabrier, G., & Gallardo, J. 2009, ApJ, 702, L27

Baraffe, I., Elbakyan, V. G., Vorobyov, E. I., & Chabrier, G. 2017, A&A, 597, A19

Baraffe, I., Vorobyov, E., & Chabrier, G. 2012, ApJ, 756, 118

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. 2009, MNRAS, 392, 1363

Bate, M. R. 2010, MNRAS, 404, L79

Bate, M. R. 2012, MNRAS, 419, 3115

Bate, M. R., Bonnell, I. A., & Bromm, V. 2003, MNRAS, 339, 577

Bhandare, A., Kuiper, R., Henning, T., et al. 2018, A&A, 618, A95

Bhandare, A., Kuiper, R., Henning, T., et al. 2020, A&A, 638, A86

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

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

Cappellari, M., McDermid, R. M., Alatalo, K., et al. 2012, Natur, 484, 485
Chabrier, G. 2003, PASP, 115, 763

Chabrier, G., Hennebelle, P., & Charlot, S. 2014, ApJ, 796, 75

Colman, T., & Teyssier, R. 2020, MNRAS, 492, 4727

Commerçon, B., Debout, V., & Teyssier, R. 2014, A&A, 563, A11

Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011, A&A, 529, A35

Csengeri, T., Leurini, S., Wyrowski, F., et al. 2016, A&A, 586, A149

Dunham, M. M., Evans, N. J. I., Terebey, S., Dullemond, C. P., & Young, C. H. 2010, ApJ, 710, 470

Elia, D., Molinari, S., Schisano, E., et al. 2017, MNRAS, 471, 100

Evans, N. J. I., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321

Fischer, W. J., Safron, E., & Megeath, S. T. 2019, ApJ, 872, 183

Frimann, S., Jørgensen, J. K., Padoan, P., & Haugbølle, T. 2016, A&A, 587, A60

Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371

Giannetti, A., Leurini, S., Wyrowski, F., et al. 2017, A&A, 603, A33

Ginsburg, A., Goddi, C., Kruijssen, J. M. D., et al. 2017, ApJ, 842, 92

Girichidis, P., Federrath, C., Banerjee, R., & Klessen, R. S. 2011, MNRAS, 413, 2741

Goddi, C., Ginsburg, A., Maud, L., Zhang, Q., & Zapata, L. 2018, arXiv:1805.05364

Guszejnov, D., Grudić, M. Y., Hopkins, P. F., Offner, S. S. R., & Faucher-Giguère, C.-A. 2020, MNRAS, 496, 5072

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

Hennebelle, P. 2012, A&A, 545, A147

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

Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020, A&A, 635, A67

Hennebelle, P., Lee, Y.-N., & Chabrier, G. 2019, ApJ, 883, 140

Hopkins, P. F. 2013, MNRAS, 433, 170

Hosokawa, T., & Omukai, K. 2009, ApJ, 691, 823

Hsieh, T.-H., Murillo, N. M., Belloche, A., et al. 2018, ApJ, 854, 15

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

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

Jones, M. O., & Bate, M. R. 2018, MNRAS, 478, 2650

Kenyon, S. J., & Hartmann, L. 1995, ApJS, 101, 117

Kroupa, P. 2001, MNRAS, 322, 231

Krumholz, M. R., Klein, R. I., & McKee, C. F. 2007, ApJ, 656, 959

Krumholz, M. R., Klein, R. I., & McKee, C. F. 2012, ApJ, 754, 71

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

Kuiper, R., & Yorke, H. W. 2013, ApJ, 772, 61

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

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

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

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

Lee, Y.-N., Hennebelle, P., & Chabrier, G. 2017, arXiv:1709.01446

Lee, Y.-N., Offner, S. S. R., Hennebelle, P., et al. 2020, SSRv, 216, 70

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

Mathew, S. S., & Federrath, C. 2020, MNRAS, 496, 5201

Motte, F., Nony, T., Louvet, F., et al. 2018, NatAs, 2, 478

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, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 53

Offner, S. S. R., Klein, R. I., McKee, C. F., & Krumholz, M. R. 2009, ApJ, 703, 131

Offner, S. S. R., & McKee, C. F. 2011, ApJ, 736, 53

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

Padoan, P., Pan, L., Juvela, M., Haugbølle, T., & Nordlund, Å. 2020, ApJ, 900, 82

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

Saumon, D., & Chabrier, G. 1992, PhRvA, 46, 2084

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

Schneider, F. R. N., Sana, H., Evans, C. J., et al. 2018, Sci, 359, 69

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

Stamatellos, D., Whitworth, A. P., & Hubber, D. A. 2011, ApJ, 730, 32

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

Urban, A., Martel, H., & Evans, N. J. I. 2010, ApJ, 710, 1343

Urquhart, J. S., Moore, T. J. T., Csengeri, T., et al. 2014, MNRAS, 443, 1555

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

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

Vorobyov, E. I., & Basu, S. 2010, ApJ, 719, 1896

Wurster, J., & Li, Z.-Y. 2018, FrASS, 5, 39

Zhao, B., Tomida, K., Hennebelle, P., et al. 2020, SSRv, 216, 43