11institutetext:  AIM, CEA, CNRS, Université Paris-Saclay, Université Paris
Diderot, Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France, 22institutetext:
INAF-IAPS, via del Fosso del Cavaliere 100, I-00133 Roma Italy,
33institutetext: Jodrell Bank Centre for Astrophysics , Department of Physics
and Astronomy, The University of Manchester, Manchester M13 9PL, United
Kingdom 44institutetext: Physikalisches Institut, University of Cologne,
Zülpicher Str. 77, 50937 Köln, Germany 55institutetext: INAF - Osservatorio
Astronomico di Cagliari, Via della Scienza 5, I-09047 Selargius (CA), Italy
66institutetext: Instituto de Radioastronomía y Astrofísica, Universidad
Nacional Autónoma de México, Apdo. Postal 3-72, 58089 Morelia, Michoacán,
México 77institutetext: Universität Heidelberg, Zentrum für Astronomie,
Institut für theoretische Astrophysik Albert-Ueberle-Str. 2, 69120,
Heidelberg, Germany

# How magnetic field and stellar radiative feedback influences the collapse
and the stellar mass spectrum of a massive star forming clump

Patrick Hennebelle How magnetic field and stellar radiative feedback
influences the collapse and the stellar mass spectrum of a massive star
forming clumpHow magnetic field and stellar radiative feedback influences the
collapse and the stellar mass spectrum of a massive star forming clump Ugo
Lebreuilly How magnetic field and stellar radiative feedback influences the
collapse and the stellar mass spectrum of a massive star forming clumpHow
magnetic field and stellar radiative feedback influences the collapse and the
stellar mass spectrum of a massive star forming clump Tine Colman How magnetic
field and stellar radiative feedback influences the collapse and the stellar
mass spectrum of a massive star forming clumpHow magnetic field and stellar
radiative feedback influences the collapse and the stellar mass spectrum of a
massive star forming clump Davide Elia How magnetic field and stellar
radiative feedback influences the collapse and the stellar mass spectrum of a
massive star forming clumpHow magnetic field and stellar radiative feedback
influences the collapse and the stellar mass spectrum of a massive star
forming clump Gary Fuller How magnetic field and stellar radiative feedback
influences the collapse and the stellar mass spectrum of a massive star
forming clumpHow magnetic field and stellar radiative feedback influences the
collapse and the stellar mass spectrum of a massive star forming clumpHow
magnetic field and stellar radiative feedback influences the collapse and the
stellar mass spectrum of a massive star forming clumpHow magnetic field and
stellar radiative feedback influences the collapse and the stellar mass
spectrum of a massive star forming clump Silvia Leurini How magnetic field and
stellar radiative feedback influences the collapse and the stellar mass
spectrum of a massive star forming clumpHow magnetic field and stellar
radiative feedback influences the collapse and the stellar mass spectrum of a
massive star forming clump Thomas Nony How magnetic field and stellar
radiative feedback influences the collapse and the stellar mass spectrum of a
massive star forming clumpHow magnetic field and stellar radiative feedback
influences the collapse and the stellar mass spectrum of a massive star
forming clump Eugenio Schisano How magnetic field and stellar radiative
feedback influences the collapse and the stellar mass spectrum of a massive
star forming clumpHow magnetic field and stellar radiative feedback influences
the collapse and the stellar mass spectrum of a massive star forming clump
Juan D. Soler How magnetic field and stellar radiative feedback influences the
collapse and the stellar mass spectrum of a massive star forming clumpHow
magnetic field and stellar radiative feedback influences the collapse and the
stellar mass spectrum of a massive star forming clump Alessio Traficante How
magnetic field and stellar radiative feedback influences the collapse and the
stellar mass spectrum of a massive star forming clumpHow magnetic field and
stellar radiative feedback influences the collapse and the stellar mass
spectrum of a massive star forming clump Ralf S. Klessen How magnetic field
and stellar radiative feedback influences the collapse and the stellar mass
spectrum of a massive star forming clumpHow magnetic field and stellar
radiative feedback influences the collapse and the stellar mass spectrum of a
massive star forming clump Sergio Molinari How magnetic field and stellar
radiative feedback influences the collapse and the stellar mass spectrum of a
massive star forming clumpHow magnetic field and stellar radiative feedback
influences the collapse and the stellar mass spectrum of a massive star
forming clump Leonardo Testi How magnetic field and stellar radiative feedback
influences the collapse and the stellar mass spectrum of a massive star
forming clumpHow magnetic field and stellar radiative feedback influences the
collapse and the stellar mass spectrum of a massive star forming clump

###### Abstract

Context. In spite of decades of theoretical efforts, the physical origin of
the stellar initial mass function (IMF) is still debated.

Aims.  We aim at understanding the influence of various physical processes
such as radiative stellar feedback, magnetic field and non-ideal magneto-
hydrodynamics on the IMF.

Methods. We present a series of numerical simulations of collapsing 1000 M⊙
clumps taking into account radiative feedback and magnetic field with spatial
resolution down to 1 AU. Both ideal and non-ideal MHD runs are performed and
various radiative feedback efficiencies are considered. We also develop
analytical models that we confront to the numerical results.

Results. The sum of the luminosities produced by the stars in the calculations
is computed and it compares well with the bolometric luminosities reported in
observations of massive star forming clumps. The temperatures, velocities and
densities are also found to be in good agreement with recent observations. The
stellar mass spectrum inferred for the simulations is, generally speaking, not
strictly universal and in particular varies with magnetic intensity. It is
also influenced by the choice of the radiative feedback efficiency. In all
simulations, a sharp drop in the stellar distribution is found at about
$M_{min}\simeq$ 0.1 M⊙, which is likely a consequence of the adiabatic
behaviour induced by dust opacities at high densities. As a consequence, when
the combination of magnetic and thermal support is not too large, the mass
distribution presents a peak located at 0.3-0.5 M⊙. When magnetic and thermal
support are large, the mass distribution is better described by a plateau,
i.e. $dN/d\log M\propto M^{-\Gamma}$, $\Gamma\simeq 0$. At higher masses the
mass distributions drop following power-law behaviours until a maximum mass
$M_{max}$ whose value increases with field intensity and radiative feedback
efficiency. Between $M_{min}$ and $M_{max}$ the distributions inferred from
the simulations agree well with an analytical model inferred from gravo-
turbulent theory. Due to the density PDF $\propto\rho^{-3/2}$ relevant for
collapsing clouds, values on the order of $\Gamma\simeq 3/4$ are inferred both
analytically and numerically. More precisely, after 150 $M_{\odot}$ of gas
have been accreted, the most massive star has a mass of about 8 M⊙ when
magnetic field is significant, and 3 M⊙ only when both radiative feedback
efficiency and magnetic field are low, respectively.

Conclusions. When both magnetic field and radiative feedback are taken into
account, they are found to have a significant influence on the stellar mass
spectrum. In particular both reduce fragmentation and lead to the formation of
more massive stars.

###### Key Words.:

ISM: clouds – ISM: structure – Turbulence – gravity – Stars: formation

##  1 Introduction

Star formation is a topic of fundamental importance in astrophysics. In
particular the mass distribution of stars, described by the initial mass
function (IMF Salpeter 1955; Kroupa 2001; Chabrier 2003; Bastian et al. 2010;
Offner et al. 2014; Lee et al. 2020) , plays a crucial role in setting the
abundances of heavy elements and regulating stellar feedback, which in turn
play major roles in the formation and evolution of galaxies and the
interstellar medium. In the efforts to find a complete description of the IMF,
it is sometimes overlooked that observed stellar masses span more than three
orders of magnitudes: from 0.1 M⊙ to more than 100 M⊙. This facts likely
implies the existence of several regimes of dominant physical processes and
star formation conditions. Clearly this problem requires a long standing
community effort and during the last decades several teams have conducted
systematic investigations with the help of numerical simulations, introducing
progressively more and more physical processes with increasingly higher
numerical resolution.

The first attempts to obtain stellar mass spectra from numerical simulations
in isothermal, self-gravitating, supersonic turbulent flows have been made by
Klessen (2001) and Bate et al. (2003). Together with several high-resolution
studies performed by various authors (e.g. Girichidis et al. 2011a; Bonnell et
al. 2011a; Ballesteros-Paredes et al. 2015; Lee & Hennebelle 2018a), they find
stellar mass spectra that present similarities with the observationally
inferred mass spectra. In particular, at high masses the distributions are
compatible with powerlaws, i.e. $dN/dlogM\propto M^{-\Gamma}$, although in
many runs, values of $\Gamma=3/4$ to 1, seemingly shallower than the canonical
$\Gamma\simeq$1.3 value inferred by Salpeter (1955) have been obtained (see
the discussion in Lee & Hennebelle 2018a). The inferred distributions also
present a peak, which however, when the simulations are strictly isothermal,
is due to limited spatial resolution. A robust, numerically converged peak is
obtained when an effective equation of state with an adiabatic index larger
than $4/3$ is taken into account (Lee & Hennebelle 2018b).

The influence of the magnetic field on the stellar mass spectrum has been
investigated by Haugbølle et al. (2018), Lee & Hennebelle (2019), Guszejnov et
al. (2020) perfoming high spatial resolution simulations with various
magnetisations. The resulting mass spectra have been found to be similar to
those inferred from simulations without magnetic field. In particular
Guszejnov et al. (2020) stress that magnetic field cannot provide a
characteristic mass that may explain the peak of the IMF and that thermal
processes have to be considered.

Several attempts have been made to study the IMF using radiative transfer
calculations. Urban et al. (2010) considered radiative feedback, i.e. stellar
and accretion luminosity, by adding them onto the sink particles. They
concluded that isothermal and radiative transfer calculations are
significantly different, in particular the stars are much more massive in
simulations with radiative feedback. Bate (2009) performed high resolution
calculations, introducing the sink particles at very high density the released
gravitational energy, i.e. $n>10^{19}$ cm-3 but stellar feedback onto the sink
particles is not explicitly included, which makes it much weaker than it
should. Krumholz et al. (2012) performed adaptive mesh refinement calculations
with a resolution of 20-40 AU. Both stellar and accretion luminosity are added
to the sinks. A relatively flat mass spectrum that is to say such that
$\Gamma\simeq 0$ is inferred when winds are not considered while in the
presence of stellar winds the mass spectra present a peak around 0.3
$M_{\odot}$ and a power-law with $\Gamma\simeq 0.5-1$. Likely enough when
winds are present, the radiation escape along the cavities and the heating is
reduced. Mathew & Federrath (2020) presented simulations with a spatial
resolution of 200 AU and perform calculations which use either a polytropic
equation of state or heating from stars. They found that when heating is
included more massive stars would form. Hennebelle et al. (2020b) conducted
adaptive mesh simulations with a spatial resolution of 4 AU and down to 1 AU.
Both stellar and accretion luminosity are treated, with various efficiencies,
$f_{acc}$, ranging from 0 to 50$\%$, as well as two sets of initial
conditions, namely a very compact and more standard clumps have been
considered. For the most compact clumps and when $f_{acc}$ is high, a flat
mass spectrum develops. Otherwise all runs present mass spectra with a peak
around 0.3-0.5 M⊙ and a powerlaw at higher masses, even when radiative
feedback is not considered, i.e. $f_{acc}=0$, and when a barotropic equation
of state is used instead. High efficiency radiative feedback runs however tend
to present a broader distribution, both at the low mass and high mass end,
with high mass stars up to 2 - 3 times more massive than in the barotropic and
low feedback efficiency runs.

In the present paper we pursue the investigation of the origin of the stellar
mass spectrum within a massive star forming clump. In particular, we focus on
the role that magnetic field, in conjunction with radiative feedback may have.
A number of studies performed calculations with both magnetic field and
radiative feedback although most of the time, without predicting mass
spectrum. As revealed by previous work (Peters et al. 2010, 2011; Commerçon et
al. 2011a; Myers et al. 2013) both these physical processes significantly
influence the collapse and star formation, particularly by reducing the
fragmentation. Moreover, their joint effect is not a mere superposition. These
studies, however, did not present sufficient statistics to draw conclusions
regarding the stellar mass spectrum. A stellar mass spectrum has been obtained
by Li et al. (2018) where a magnetized and radiative calculation is performed
with a spatial resolution of about 30 AU. These statistics need to be expanded
and various initial conditions must be systematically explored. To do so we
perform high resolution simulations of massive star forming clumps where both
magnetic field and radiative feedback are accounted for. To get a good
description of the small scales which are mandatory to describe the formation
of low mass stars, we employ an adaptive mesh reffinement with a spatial
resolution down to 1 AU. As the magnetic intensity is likely varying from
clump to clump and not many constraints from observations are available yet,
we explore three magnetisations. Also the radiative feedback efficiency is
subject to large uncertainties, so we consider two different values.
Importantly, we also perform a simulation in which non-ideal MHD effects,
namely ambipolar diffusion (Mestel & Spitzer 1956), are explicitly taken into
account. We stress that these runs are the first for which both magnetic field
and radiative feedback are taken into account, while considering a
configuration which leads to sufficient statistics and spatial resolution to
provide a reliable stellar distribution in the range 0.1 to 10 M⊙.

The paper is structured as follows. The second section presents the equations
that are being solved, the relevant physical processes as well as the
numerical methods used to solve these equations. It also presents the initial
conditions and describes the various runs presented in the paper. In the third
section, we look at the evolution of the clump during its collapse and
investigate the effect of the magnetisation and radiative feedback. The global
properties such as the total accreted mass and radiated energy, the
temperature, magnetic field and mass distribution are studied. An analytical
model, which is presented in an appendix is developed to understand the
temperature distribution in the simulations. Comparisons are made with
observational results. The fourth section presents the stellar mass spectrum
obtained in the simulations. They are quantitatively compared with an
analytical model which gives more insight into the effect of the different
physical processes and is also presented in an appendix. In the fifth section
a discussion is given while the sixth section concludes the paper.

Figure 1: Top row shows the total accreted mass (top-left panel) and the total
accretion rate (top-right panel) as a function of time for the various runs.
Bottom-left panel portrays the total luminosity while bottom-right panel shows
the same quantity divided by $L_{glob}=0.5\times
GM_{*,tot}\dot{M}_{*,tot}/(2R_{\odot})$.  Figure 2: Column density at $t=0.11$
Myr for the six runs. The dark circles represent the sink particles aiming at
describing the stars.
HYDROf05MHD10f05HYDROf05MHD10f05HYDROf05MHD10f05HYDROf05MHD10f05 Figure 3:
Radial profiles at several timesteps in HYDROf05 and MHD10f05, respectively.
The black dashed line in the density panel shows the singular isothermal
sphere density. HYDROf01HYDROf05NMHD10f05MHD100f05MHD10f01MHD10f05 Figure 4:
Mass distribution as a function of density at several timesteps for the six
runs. Section 4.2 describes the corresponding physical regimes. The dotted
line shows a powerlaw behaviour $M\propto n^{-1/2}$ as expected for a density
PDF $\propto n^{-3/2}$ (see Sect. 4.2.3).
HYDROf01HYDROf05NMHD10f05MHD100f05MHD10f01MHD10f05 Figure 5: Mass-weighted
temperature in density intervals as a function of density at six time steps
for the six runs. The dotted lines show a powerlaw behaviour of $T\propto
n^{1/2}$ and $T\propto n^{0.3}$, respectively.
NMHD10f05MHD100f05MHD10f01MHD10f05 Figure 6: Magnetic field as a function of
density at several timesteps for the four magnetized runs. The dotted line
shows a powerlaw behaviour $B\propto n^{1/2}$.

##  2 Numerical simulations

###  2.1 Equations, numerical methods and setup

In this paper we solve the equations of the radiative magneto-hydrodynamics.
All the radiative quantities are estimated in the co-moving frame and assuming
the grey approximation; that is to say the radiative energies are integrated
over the entire frequency spectrum (e.g. Commerçon et al. 2011a). The
equations are

$\begin{array}[]{llll}\partial_{t}\rho+\nabla\cdot\left[\rho\textbf{u}\right]&=0,\\\
\partial_{t}\rho\textbf{u}+\nabla\cdot\left[\rho\textbf{u}\otimes\textbf{u}+P\mathbb{I}\right]&=-\lambda\nabla
E_{\mathrm{r}}+\textbf{F}_{\mathrm{L}}-\rho\nabla\phi,\\\
\partial_{t}E_{\mathrm{T}}+\nabla\cdot\left[\textbf{u}\left(E_{\mathrm{T}}+P_{\mathrm{}}+\frac{B^{2}}{8\pi}\right)\right.&\\\
\hskip 30.0pt\left.-{1\over
4\pi}\left(\textbf{u}\cdot\textbf{B}\right)\textbf{B}-\textbf{E}_{\mathrm{AD}}\times\textbf{B}\right]&=-\mathbb{P}_{\mathrm{r}}\nabla:\textbf{u}-\lambda\textbf{u}\nabla
E_{\mathrm{r}}\\\ &\hskip
9.0pt+\nabla\cdot\left(\frac{c\lambda}{\rho\kappa_{\mathrm{R}}}\nabla
E_{\mathrm{r}}\right)+S_{\star}\\\ &\hskip
9.0pt-\rho\textbf{u}\cdot\nabla\phi,\\\
\partial_{t}E_{\mathrm{r}}+\nabla\cdot\left[\textbf{u}E_{\mathrm{r}}\right]&=-\mathbb{P}_{\mathrm{r}}\nabla:\textbf{u}+\nabla\cdot\left(\frac{c\lambda}{\rho\kappa_{\mathrm{R}}}\nabla
E_{\mathrm{r}}\right)\\\ &\hskip 9.0pt+\kappa_{\mathrm{P}}\rho
c(a_{\mathrm{R}}T^{4}-E_{\mathrm{r}})\\\ &\hskip 9.0pt+S_{\star},\\\
\partial_{t}\textbf{B}-\nabla\times\left[\textbf{u}\times\textbf{B}+\textbf{E}_{\mathrm{AD}}\right]&=0,\\\
\nabla\cdot\textbf{B}&=0,\\\ \Delta\phi&=4\pi G\rho,\end{array}$ (1)

where $\rho$ is the material density, u is the velocity, $P$ the thermal
pressure, $\lambda$ is the radiative flux limiter (Minerbo 1978),
$E_{\mathrm{r}}$ is the radiative energy,
$\textbf{F}_{\mathrm{L}}=1/(4\pi)(\nabla\times\textbf{B})\times\textbf{B}$ is
the Lorentz force, $\phi$ is the gravitational potential, $E_{\mathrm{T}}$ the
total energy $E_{\mathrm{T}}=\rho\epsilon+1/2\rho
u^{2}+B^{2}/(8\pi)+E_{\mathrm{r}}$ ($\epsilon$ is the gas specific internal
energy), B is the magnetic field, $\textbf{E}_{\mathrm{AD}}$ is the ambipolar
electromotor field (EMF), $\kappa_{\mathrm{P}}$ is the Planck mean opacity,
$\kappa_{\mathrm{R}}$ is the Rosseland mean opacity, $\mathbb{P}_{\mathrm{r}}$
is the radiation pressure, $S_{\star}$ the luminosity source, and $T$ is the
gas temperature. The ambipolar EMF is given by

$\textbf{E}_{\mathrm{AD}}=\frac{\eta_{\mathrm{AD}}}{B^{2}}\left[\left(\nabla\times\textbf{B}\right)\times\textbf{B}\right]\times\textbf{B},$
(2)

where $\eta_{\mathrm{AD}}$ is the ambipolar diffusion resistivity, calculated
as a function of the density, temperature, and magnetic field amplitude.

The numerical method is overall very similar to the one used in Hennebelle et
al. (2020a). The simulations were performed with the adaptive mesh refinement
(AMR) magnetohydrodynamics (MHD) code RAMSES (Teyssier 2002; Fromang et al.
2006). When non-ideal MHD, i.e. ambipolar diffusion, is included the scheme is
the one described in Masson et al. (2012) and used in previous studies (Masson
et al. 2016; Hennebelle et al. 2020b; Mignon-Risse et al. 2021; Commerçon et
al. 2021; Lebreuilly et al. 2021). The resistivities are the ones calculated
in Marchand et al. (2016).

In all simulations presented here, radiative transfer is accounted for using
the flux limited diffusion method assuming grey approximation (see Commerçon
et al. 2011b, 2014). The flux limited diffusion method is known to present
some restrictions for instance it does not treat shadows well due to its
isotropic nature. More accurate methods such as the M1 method (González et al.
2007), the hybrid method (Kuiper et al. 2010; Mignon-Risse et al. 2020) or the
VETTAM method (Menon et al. 2022) have been developed and deal significantly
better with anisotropic radiative transfer (see also Jaura et al. 2018; Peter
et al. 2022). However, they tend to be more costly than the flux limited
diffusion method employed in this work, and are often limited in their current
implementations. This certainly represents line of future improvements.

At high density, the equation of state is taken from Saumon & Chabrier (1992)
and Saumon et al. (1995) which takes into account 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 temperatures and densities covered in this
work, the opacities are the ones calculated in Semenov et al. (2003).

The boundary conditions are periodic. The cloud is initially spherical and has
a radius four times lower than the computational domain size. All simulations
were run on a regular grid of $256^{3}$ computing cells and 10 AMR levels have
been further added during the course of the calculation leading to a total
number of 18 AMR levels. The resolution criterion is the Jeans length and it
is resolved with at least 10 points. In the present paper, the issue of
numerical resolution is not further discussed and we refer to the appendix of
Hennebelle et al. (2020a) for an investigation of the impact of numerical
resolution.

###  2.2 Sink particles and stellar feedback

The sink particle algorithm is described in Bleuler & Teyssier (2014). Sink
particles are formed at the highest refinement level at the peak of clumps
whose maximum density is larger than $n_{\rm acc}$. The sink particles are
created if the parent clump has a density $n>n_{\rm acc}$ and if it is
sufficiently gravitationally bound (see Bleuler & Teyssier 2014). The value of
$n_{\rm acc}$ is equal to $10^{13}$ cm-3. With this value of $n_{\rm acc}$,
the computational cells having a density equal to $n_{\rm acc}$ possess a mass
of roughly 1-2$\%$ of the mass of the first hydrostatic core,
$M_{L}=0.03M_{\odot}$. At each time step, 10$\%$ of the gas mass inside the
sink’s accretion radius and with a density above $n_{\rm acc}$ is retrieved
from the grid and accreted by the sink. The sinks are not allowed to merge.
The impact of changing the value of $n_{\rm acc}$ has been discussed in
Hennebelle et al. (2020a). It has been found that both the spatial resolution
and the value of $n_{\rm acc}$ may influence the peak of the stellar
distribution. However once the first hydrostatic core is sufficiently
resolved, this should not be the case.

Sink particles are also a source of radiation due to the stellar luminosity
and gas accretion. The accretion luminosity is given by

$\displaystyle L_{\rm acc}={f_{\rm acc}GM_{*}\dot{M}\over R_{*}}.$ (3)

where $M_{*}$ and $R_{*}$ are respectivelly the star’s mass and radius while
$\dot{M}$ is the accretion rate. If all the kinetic energy of the infalling
gas was radiated away, we would have $f_{\rm acc}\simeq 1$. The accretion
luminsosity has been shown to be the dominant source of gas heating at early
time and has important effects on the surrounding gas (e.g. Krumholz et al.
2007; Offner et al. 2009). The stellar luminosity of the protostars, $L_{*}$
and $R_{*}$ are taken from Kuiper & Yorke (2013) (see also Hosokawa & Omukai
2009). As discussed in Hennebelle et al. (2020a), the value of $f_{\rm acc}$
that should be used is not clearly established. In particular, the radiation
is emitted at very small scale, i.e. few stellar radii and is not expected to
propagate uniformly because of the highly anisotropic density distribution
(e.g. Krumholz et al. 2012). As in Hennebelle et al. (2020a), we perform
simulations in which we use an effective accretion luminosity and explore the
values $f_{\rm acc}=0.1$ and 0.5. By considering an effective luminosity
smaller than the estimated total luminosity, we envisage that the rest of the
energy either escape preferentially along the cavities open by winds and jets
or is converted into jet or a wind kinetic energy. This is obviously an
important source of uncertainties which requires further investigations.

We start considering accretion and stellar luminosities when the sink has a
mass of about 2 $M_{L}$, i.e. 0.07 M⊙. The reason is that due to the limited
spatial resolution, when the sink is introduced the protostar is not truly
formed yet. Since the size of the sink particles is not very different from
the radius of the first hydrostatic core, it seems reasonable to assume that
the protostar is formed only when the sink reaches a mass equal to a few
$M_{L}$. Note that although reasonable, this assumption clearly requires
further investigation. For instance, Bhandare et al. (2020) who have performed
two-dimensional simulations of the second Larson cores, i.e. the young
protostar, found that they grow with time far beyond the solar radii. This
clearly suggests that at least some of the accretion energy is not fully
radiated away but somehow stored in the star for some time. Indeed, the
accretion shock at the edge of the second Larson core is subcritical (Vaytet
et al. 2013) meaning that most of the accretion energy is advected inside the
protostar and not immediately radiated away. This constitutes an important
source of uncertainty for calculations such as the ones performed in this
work.

###  2.3 Initial conditions and runs performed

Our initial conditions consist in spherical clouds in which a turbulent
velocity field has been added. The velocity field has a classical Kolmogorov
power-spectrum equal to 11/3 with random phases. A fully self-consistent
approach would require to also set up the density and magnetic field
fluctuations. This is however not an easy task. In practice it requires
running a large scale simulations and zooming-in or at least performing a
preliminary run without self-gravity (see for instance Lane et al. 2022). Note
that Lee & Hennebelle (2018a) have compared various approaches including
starting from a previous phase during which the simulation is run without
self-gravity and starting directly from a prescribed turbulent field as it is
done here. They found very similar results. This suggests that, at least in
the context of collapsing clumps, the choice of the initial turbulent field
may not be so important, probably because as the collapse proceeds, the
fluctuations evolve and the initial perturbations are largely forgotten.

The clump we consider has a mass of $10^{3}M_{\odot}$ and an initial radius of
0.4 pc corresponding to a uniform density of about 8$\times 10^{4}$ cm-3
initially. Observationally, this corresponds to relatively standard massive
star forming clumps (e.g. Urquhart et al. 2014; Elia et al. 2017, 2021; Lin et
al. 2022). With an initial temperature of 10 K, the ratio of the thermal over
gravitational energy is about 0.008. The clump density leads to a freefall
time of about 110 kyr. The initial value of the Mach number is equal to 7
leading to a turbulent over gravitational energy ratio of about 0.4, i.e. the
clumps are close to be initially virialised.

We set up the simulations with a uniform initial magnetic field through the
cloud and intercloud medium. We considered two initial mean-field strengths,
with mass-to-flux ratios, $\mu$, of respectivelly 10 and 100 (corresponding to
about 100 $\mu$G and 10 $\mu$G respectivelly). These values are motivated by
the observations of $\mu$ on the order of a few in dense cores (e.g. Crutcher
2012; Myers & Basu 2021). This selection also aims to account for the broad
dispersions in the $\mu$ values of the massive clumps identified in the 1-kpc
scale simulation presented in Hennebelle (2018). In these MHD simulations,
which have a spatial resolution down to 400 AU, it has been found that the
mass-to-flux ratio, which presents a broad dispersion, is indeed on the order
of a few for Solar mass cores but is lower for the more massive ones. More
precisely, Fig. 5 of Hennebelle (2018) shows the distribution of self-
gravitating objects with density larger than $10^{4}$ cm-3. It shows a clear
trend that the mass-to-flux ratio, $\mu$, is increasing with mass in spite of
a broad distribution. The most massive clumps displayed has a mass of about
100 M⊙ and thus one needs to extrapolate to get a hint on 1000 M⊙ clumps.
Based on this figure, we would expect that the typical $\mu$ of a 1000 M⊙
clumps is certainly larger than 10. This may at first sight be surprising
because low mass cores have been observed to present values of $\mu$ on the
order of a few (Pattle et al. 2022). It should however be remembered that the
mass-to-flux ratio is the ratio of a volume over a surface weighted quantity.
Thus considering objects with similar mean density, the mass-to-flux is
expected to increase with the object size.

Another fundamental aspect is the physics of the magnetic field evolution.
Whereas many studies have assumed ideal MHD, it is however clear that this is
a poor approximation at high density in star forming regions (e.g. Zhao et al.
2020). We performed one run with ambipolar diffusion and a magnetisation
corresponding to $\mu=10$. Due to the small time steps induced by the second
order derivative in the ambipolar diffusion operator, this numerical
simulation is quite challenging and has required about 500,000 cpu hours.

Table 1 summarizes the various runs performed.

Name |  $R_{c}$ (pc) | $\mu$ | $\mathcal{M}$ | $lmax$ |  $dx$ (AU) | $f_{\rm acc}$ | NMHD
---|---|---|---|---|---|---|---
NMHD10f05 | 0.4 | 10 | 7 | 18 | 1.15 | 0.5 | yes
MHD10f05 | 0.4 | 10 | 7 | 18 | 1.15 | 0.5 | no
MHD10f01 | 0.4 | 10 | 7 | 18 | 1.15 | 0.1 | no
MHD100f05 | 0.4 | 100 | 7 | 18 | 1.15 | 0.5 | no
HYDROf05 | 0.4 | $\infty$ | 7 | 18 | 1.15 | 0.5 | no
HYDROf01 | 0.4 | $\infty$ | 7 | 18 | 1.15 | 0.1 | no
Table 1: Summary of the runs performed. $R_{c}$ is the initial clump radius.
$\mu$ is the mass-to-flux over critical mass-to-flux ratio. $\mathcal{M}$ is
the initial clump Mach number. $lmax$ is the maximum level of grid used and
$dx$ corresponds to the maximum resolution which is equal to 1.15 AU in these
runs. $f_{\rm acc}$ gives the fraction of the accretion luminosity which is
taken into account in the calculation. NMHD tells whether non-ideal MHD is
accounted for.

Two ideal MHD runs (MHD10f05 and MHD100f05) and one purely hydrodynamical runs
(HYDROf05) allow us to investigate the role of the magnetic field during the
collapse notably on the initial mass function. Whereas the non-ideal MHD run
(NMHD10f05) is the most realistic simulation and to our knowledge is the first
simulation of a massive star forming clump which includes both radiative
feedback and non-ideal MHD. These four runs are complemented by two runs with
a lower accretion luminosity parameter, $f_{acc}=0.1$, namely MHD10f01 and
HYDROf01 which allow to assess and discuss the influence of radiative feedback
on our results. All simulations are carried out until at least 150 $M_{\odot}$
of gas have been accreted onto the sink particles, except run NMHD10f05 for
which the final mass accreted by the sinks is 100 $M_{\odot}$. This is
because, as explained above, this simulation is more computationally
demanding.

##  3 General clump description

In this section, we look at the global evolution and final properties of the
collapsing clump as a whole. We start by describing the general morphology,
before proceeding to discuss the star formation and luminosity. We then study
the gas density and temperature distribution inside the clump.

###  3.1 Total accreted mass

Figure 1 displays the total mass accreted (top-left panel) by sink particles,
$M_{tot,*}$ as a function of time for the six simulations. The total accretion
rate is also plotted (top-right panel). Since it is a heavily fluctuating
quantity, the latter is calculated by averaging its instantaneous value over
100 uniformly spaced time intervals. As indicated above, all simulations are
run until about 150 $M_{\odot}$ have been accreted, corresponding to a star-
formation efficiency of 15$\%$. There are two exceptions, NMHDf05 for which at
the end of the simulation 100 $M_{\odot}$ have been turned into the sinks, and
MHD10f05 for which the final total mass of sinks is equal to 200 $M_{\odot}$.
Two groups of simulations are easily distinguished. On one hand the
hydrodynamical simulations and the low magnetized one, MHD100f05 and on the
other hand the more magnetized ones, i.e. MHD10f01, MHD10f05 and NMHDf05. As
expected, due to magnetic support, the latter group of simulations collapses a
bit more slowly. Two points are worth mentioning, first the simulations with
$f_{acc}=0.1$ and with $f_{acc}=0.5$ behave very similarly showing that in
spite of strong heating, radiation does not significantly alter the large
scale dynamics. Similar conclusion is also reached for the ambipolar
diffusion. This is because i) thermal support is rather weak and an increase
of temperature even by a factor of several does not make thermal support
sufficiently strong to provide a significant support at the clump scale; and
also because ii) magnetic field is only significantly modified by non-ideal
MHD processes, at high density (say $n>10^{7}$ cm-3). Interestingly, we see
that after a fast increase the accretion rate, $\dot{M}_{tot,*}$, reaches
values, which are nearly identical for all simulations and equal to about
$10^{-2}$M⊙ yr-1. This is because the accretion rate is controlled by the
largest scale, here the clump, which is globally collapsing. The magnetic
intensity considered here is too weak to significantly modify this global
dynamics.

###  3.2 Total luminosity

Bottom-left panel of Fig. 1 portrays the total luminosities,
$\sum(L_{*}+L_{acc})$ of the sink particles (see also Fig. 9). As for the
accretion rate, after an increase which takes about 0.02 Myr, it reaches, in
the case with $f_{acc}=0.5$, a plateau at about $1-2\times 10^{5}$ L⊙. In the
case $f=0.1$, the total luminosity is 10-20 times lower for run HYDROf01 and 3
times lower for run MHD10f01. Note that Fig. 1 shows that the clumps spend
about 0.02 Myr in the protostellar phase with a luminosity several times lower
than the peak values.

It is interesting to compare these values with observations. Although what is
observationally available is the bolometric luminosities rather than the total
source luminosities, they are obviously related and can be compared. In
particular, it is expected that the radiation emitted by the star in the
visible domain is quickly absorbed and reemitted in the infrared by the dust.
For instance Fig. 12 of Elia et al. (2017) displays the bolometric
luminosities as a function of mass for a sample of clumps. As can be seen our
values agree well with the luminosity distribution of the protostellar 1000 M⊙
clumps, which range from 102 to few 105 L⊙ although correspond to their upper
values when $f_{acc}=0.5$. Note, however, that the bolometric luminosities
calculated in Elia et al. (2017) correspond to wavelengths longer than 20$\mu
m$. Therefore these values represent themselves lower limits of the real
luminosities.

Lin et al. (2022) present detailed observations for several massive star
forming clumps with comparable mass and radii (see their Table 6).
Luminosities of a few $10^{5}$ L⊙ are also reported.

Since the mass of the clump is 1000 $M_{\odot}$, this means that once the
luminosity is about 105 L⊙, the luminosity per solar mass is about $10-100$ L⊙
/ M⊙. Again this is in good agreement with the values seen in Fig. 13 of Elia
et al. (2017).

In order to define a reference with which the luminosities can be compared, we
define a quantity $L_{glob}=0.5\times GM_{*,tot}\dot{M}_{*,tot}/(2R_{\odot})$,
which would correspond to the accretion luminosity of an object of mass
$M_{*,tot}$ and radius 2 $R_{\odot}$, accreting at a rate $\dot{M}_{*,tot}$
with an efficiency $f_{acc}=0.5$. The ratio $\sum L_{*}/L_{glob}$ is expected
to be smaller than 1 because the luminosity is a non-linear quantity, which
decreases with the number of stars. It gives a sense of how fragmented is the
clump and how efficiently is the gravitational energy converted into
radiation. Bottom-right panel shows $\sum L_{*}/L_{glob}$ for the six runs as
a function of time. As expected, this quantity decreases with time and, at
later time, it reaches values as small as $10^{-3}$ and as high as $2\times
10^{-2}$, depending on the run. Interestingly, there is a clear trend for the
magnetized runs to have values 1.5-2 times larger than their hydrodynamical
counterpart. As we show later this is because magnetic field tends to reduce
fragmentation, therefore building more massive stars which present higher
luminosities.

###  3.3 General morphology

Figure 2 portrays the column density of the whole clump at time $t=0.11$ Myr,
which as seen from Fig. 1, corresponds to a time where approximately 100-120
$M_{\odot}$ have been accreted. The dark circles show the sink particles,
which represent individual stars. The six simulations present a similar
pattern. A complex network of intervowen and interconnected filaments have
formed and three of them appear to be a little more prominent. Their length is
approximately $\simeq 0.3$ pc and is comparable to the whole clump size. The
three main filaments intersect, forming a hub located approximately at
$y=0.75$ pc and $z=0.65$ pc. The stars are represented by the dark circles,
are mainly, though not exclusively, located in the hub and in the filaments.
Let us recall that density filaments are naturally produced both by MHD
turbulence (Hennebelle 2013; Federrath 2016; Xu et al. 2019), shocks (Abe et
al. 2021), and by gravity (e.g. Smith et al. 2014; Abe et al. 2021), which for
different reasons tend both to amplify anisotropies.

Beyond these general similarities, significant differences between the six
simulations are clearly visible. First of all, we see that radiative feedback
has a clear influence on the cloud evolution and its fragmentation (Krumholz
et al. 2007; Hennebelle et al. 2020a). For instance there are more sinks in
run HYDROf01 than in run HYDROf05 and in run MHD10f01 than in MHD10f05 (this
will be further quantified in § 5.1). This is a clear consequence of less
heating when $f=0.1$ than when $f=0.5$. The impact on the gas structure
appears to remain more limited. Second of all, clearly magnetic field reduces
significantly the numbers of sinks. This is particularly obvious by comparing
run HYDROf05 with run MHD10f05 as well as run HYDROf01 with run MHD10f01. It
is also clear that sink particles tend to form in higher column density
regions in the magnetized runs. This clearly is a consequence of the support
provided by the magnetic field which efficiently stabilizes the gas
particularly when its column density is not too high. This happens obviously,
only if magnetic field is strong enough. Indeed in run MHD100f05 which has an
initial magnetic field 10 times lower than run MHD10f05, the sink distribution
is very similar to run HYDROf05. Interestingly, the sink distribution in run
NMHDf05 is comparable to run MHD10f05 except near the high column densities
areas where a small excess of sinks is sometimes visible. This stems from the
fact that ambipolar diffusion is efficient only at small scales and at high
density.

##  4 Gas and magnetic field distribution

###  4.1 Density, velocity and temperature profiles

Figure 3 displays radial profiles of various density-weighted quantities for
runs MHD10f05 (left) and HYDROf05 (right) and at several timesteps. For the
sake of conciseness only 2 simulations are being displayed and discussed here.
The adopted center is the position of the most massive sink particle, which is
located in the hub at $y=0.75$ pc and $z=0.65$ pc.

The first row displays the gas temperature that we remind is initially uniform
and equal to 10 K. As time goes on, temperature increases by 2 to 3 orders of
magnitudes in the center and about one order of magnitude in the clump’s outer
part. The temperature profile in the clump inner part broadly behaves as
$r^{-1}$ while it is almost flat in the clump outer part. We further note the
presence of many temperature peaks associated to sink particles distributed
through the clouds. We also stress that the temperature is clearly larger by a
factor $\simeq 1.5$ in run MHD10f05 than in run HYDROf05. We will come back on
this particular point later but this effect is similar to what has been
reported by Commerçon et al. (2011a), where radiative MHD calculations where
also performed and higher temperatures have been reported in the MHD case.
This is a consequence of the non-linearity of the accretion luminosity
proportional to $\dot{M}M$. By reducing fragmentation and extracting angular
momentum, magnetic field increases both $M$ and $\dot{M}$ leading to higher
accretion luminosity. These temperatures appear to be in good agreement with
the ones presented in Fig. 10 of Lin et al. (2022). For instance, the
temperature at few 0.01 pc is about 100-200 K while at 0.1 pc it is typically
50-70 K.

Second row shows the density profiles. The straight line represents the
density of the singular isothermal sphere (SIS), i.e. $\rho=c_{s,0}^{2}/(2\pi
Gr^{2})$, where $c_{s,0}$ is the sound speed taken here equal to 0.2 km s-1.
As collapse proceeds, the density increases from outside-in and after roughly
0.1 Myr, it presents a powerlaw-like shape close to, but slightly shallower
than, $r^{-2}$. As we see the values evolve with times and also slightly
depend on the radius. We see however than it is nearly 2 orders of magnitude
denser than the SIS, which is an expected consequence of the low initial
thermal energy and the compactness of the cloud. The density in run HYDROf05
is slightly lower than in run MHD10f05, which is a consequence of the magnetic
support. Let us remember that $r^{-2}$ density profile is the expected density
structure of a spherical collapsing cloud (e.g. Larson 1969; Shu 1977) as,
together with a uniform radial velocity, it leads to a roughly constant
accretion rate through the cloud (Li 2018; Gómez et al. 2021). When turbulence
is included, it is however common to find profiles slightly shallower, for
instance $\rho\propto r^{-1.5}$ is often reported (e.g. Murray & Chang 2015;
Li et al. 2018), though in the present case this value seems a little too
shallow.

The radial velocity through the cloud is presented in the third row. At 0.1
Myr, a constant radial velocity of $\simeq 1.5-2$ km s-1 appears to reasonably
represent the cloud radial velocity for radius between 0.03-0.3 pc. The radial
velocity increases towards the cloud inner part where it reaches $\simeq$ 10
km s-1. The parallel velocity, which represent both turbulent and rotation
(i.e. the non-radial component) is displayed in the fourth row. Due to the
chosen initial conditions, it is of the order of $\simeq 2$ km s-1 in the
cloud outer part. As the collapse proceeds and due to the increase of $v_{r}$,
turbulence is further amplified toward the cloud center (Hennebelle 2021) and
this behaviour explains the density profile being shallower than $r^{-2}$.
These velocity values are in good agreement with the values presented in Fig.
22 of Lin et al. (2022). For instance at 0.1 pc, values of about 3 km s-1 are
reported.

###  4.2 Mass distribution

The mass distribution, which we remind is equivalent to the mass weighted
density PDF, is displayed in Fig. 4 for the six simulations at several
timesteps. The distribution contains several features and going from low to
high densities four domains can be identified: the interclump medium, the
clump outer part, the collapsing envelopes and the high density material. We
stress that the last two do not correspond to a single physical region but
rather develop around each individual collapse center.

####  4.2.1 The interclump medium

At low density the mass distribution presents a roughly lognormal shape which
peaks at about 500 cm-3 (e.g. Vázquez-Semadeni 1994; Federrath et al. 2008;
Kritsuk et al. 2011) and remains stationary through time. It is due to the
development of turbulence in the cloud outerpart. The latter has formed by the
turbulent-driven expansion of the cloud external layer. Clearly it contains a
small amount of mass.

####  4.2.2 The clump outer part

At higher density, i.e. $\rho\simeq 10^{5}-10^{6}$ cm-3, a second peak of the
mass distribution located at the cloud initial mean density, is visible. It
contains most of the mass of the cloud and shifts toward higher densities as
collapse proceeds. Meanwhile as expected the mass it contains, declines over
time. Overall the mass distribution of this density range is similar for the
six simulations. We can nevertheless note that the peak is a bit broader for
the two hydro runs, than for the more magnetized runs MHD10f01, MHD10f05 and
NMHDf05. This is likely a consequence of the magnetic field which is known to
reduce the turbulent dispersion of the density distribution (e.g. Molina et
al. 2012).

####  4.2.3 The collapsing envelopes

At densities higher than its peak value, the mass distribution is better
described by a powerlaw behaviour up to densities of $10^{9}$ and even 1010
cm-3. This part of the mass distribution corresponds to the $n\propto
r^{-\alpha}$, $\alpha\simeq 2$, envelope discussed in Fig. 3.

Let us remember that there is a simple correspondence between $\alpha$ and the
index of the mass distribution. Let $d{\mathcal{N}}$ be the number of fluid
particles located between radius $r$ and $r+dr$. We have
$d{\mathcal{N}}\propto r^{2}dr$. But since $n\propto r^{-\alpha}$, we have
$d{\mathcal{N}}/d\log n\propto n^{-3/\alpha}$ and the mass weighted density
PDF is

$\displaystyle n{d\mathcal{N}\over d\log n}\propto n^{-3/\alpha+1}.$ (4)

For $\alpha\simeq 2$, we thus find that $n{d\mathcal{N}}/d\log n\propto
n^{-1/2}$, which indeed is close to the observed bevaviour of the mass
distribution between $10^{7}$ and $10^{9}-10^{-10}$ cm-3 as shown by a
comparison with the dotted lines.

Several aspects are worth noticing. First at early stages (black and red
curves), the mass distributions evolve with time. More mass is gradually
accumulated at high densities as collapse proceeds. Once the $n\propto r^{-2}$
envelope is fully developed, the mass distribution is stationary. This is all
consistent with the stationarity observed in Fig. 3 illustrating that the
accretion rate remains broadly constant with time.

####  4.2.4 The high density material

At density larger than $10^{9}-10^{10}$ cm-3, the mass distribution becomes
flatter, meaning that mass is pilling up. This is a consequence of rotational
and thermal supports. Indeed protoplanetary disks form (see Lebreuilly et al.
2021, for a description of disks in similar simulations). Clearly the amount
of mass significantly varies with magnetisation and it is several times higher
in the hydro runs than in the significantly magnetized ones (MHD10). This is a
clear consequence of magnetic braking, which by extracting angular momentum
leads to smaller and less massive disks.

###  4.3 Temperature vs density distributions

Figure 5 shows the mean temperature as a function of density in the six
simulations. In each density bin, the mean temperature is simply the mass
weighted temperature. The overall behaviour is as suggested by the temperature
profiles shown in § 4.1.

The temperature associated to the high density material is typically larger
than $\simeq 300$ K and reaches values of few thousands K. As expected the
temperature increases with $f$.

For the lower density material (i.e. $n<10^{9}$ cm-3), we see first that the
temperature decreases roughly as $T\propto n^{-0.3-0.5}$ (as indicated by the
dotted line) and then at density of about $\simeq 10^{7}$ cm-3, it reaches a
plateau and remains constant, $T=T_{ext}$, at lower densities. Depending of
the runs and the time the temperatures vary between 10 and up to $\simeq$ 30
K.

To interpret these temperatures, we developed a simple spherical model which
is presented in § A. Although we see from Fig. 2 that the clouds are not
spherical and that the sources are not clustered in the center as assumed in
our model, this nevertheless allows us to get a deeper understanding of these
temperatures. The inferred powerlaw behaviours are as described by Eq. (A) and
Eq. (15). More precisely, Eq. (A) combined with Eq. (8) predicts that for
$T>100$ K, $T\propto n^{0.5}$ while for $T<100$ K, $T\propto n^{0.3}$.

To quantitatively estimate the values of $T_{\rm ext}$, we use Eq. (16)

$\displaystyle T_{ext}=33.5K\;\left({\tau_{0}\over
0.5}\right)^{3/(4+2\alpha)}\left({L_{tot}\over
10^{5}L_{\odot}}\right)^{1/(4+2\alpha)}\left({\delta\over
100}\right)^{-1/(2+\alpha)}$ (5)

where we remind that $\delta$ is as defined by Eq. (8) and $\tau_{0}$ is the
optical depth at which the radiation is free streaming.

From Fig. 1 and Fig. 5, we see that when $L_{tot}\simeq 10^{5}$ L⊙,
$T_{ext}\simeq 30$ K, whereas when $L_{tot}\simeq 10^{4}$ L⊙ $T_{ext}\simeq
20$ K, which is close to what Eq. (5) predicts. Looking at Fig. 5 of Elia et
al. (2017), we see that 20-30 K corresponds to the temperature of the warmest
star forming clumps, which agrees well with the relatively high luminosities
that we inferred. The HiGAL-based temperature is the average termperature of
the cold dust in a clump. They are derived from 160-to-500 (and 870, 1100,
when available) $\mu m$ grey-body fit, so that probed temperatures cannot be
higher than that. Since the mass in the outer part of the clump dominates,
this cold component corresponds to that of the outer layers and of most of the
volume of the clump. This should therefore broadly correspond to what
$T_{ext}$ is.

###  4.4 Magnetic field distributions

Figure 6 portrays the volume weighted magnetic intensity as a function of gas
density for the 4 magnetized simulations. Overall we see that, at least for
$n$ between $10^{7}$ and $10^{9}$ cm-3, the magnetic field scales with density
broadly as $B\propto n^{1/2}$, a result observed in previous works (see for
instance Hennebelle & Inutsuka 2019, for a review). This is a consequence of
the field amplification induced by field lines dragging by collapsing motion.
Even more simply, this is likely a consequence of energy equipartition. As
seen in Fig. 3, $v^{2}$ depends weakly on $r$ while $n\propto r^{-2}$,
therefore the kinetic energy scales as $r^{-2}$ and thus $B\propto
r^{-1}\propto n^{1/2}$. Interestingly, this implies that the Alfvén velocity,
$V_{a}$, remains roughly constant in this range of density. We see however
that its value is not identical for the four runs. We estimate that for runs
MHD10f05 and MHD10f01, $V_{a}\simeq 1-1.2$ km s-1 while for run MHD100f05,
$V_{a}$ is less than half this value. When non-ideal MHD is treated, the
Alfvén velocity is reduced by tens of percents at $n\simeq 10^{9}$ cm-3.

At lower densities, the behaviour depends on the field intensity. For run
MHD100f05, the dependence of $B$ on $n$, is a bit stiffer. This is expected as
when the field is weak, the clump contraction tends to be spherical in which
case $B\propto n^{2/3}$ (Li et al. 2015). This explains why the magnetic field
at high density in run MHD100f05 is larger than a tens of the $B$ values in
run MHD10f05. Magnetic intensity is more vigorously amplified when it is
weaker.

At high densities, i.e. $n>10^{9}$ cm-3, the magnetic field is further
amplified up to density values on the order of $10^{11}$ cm-3. The highest
magnetic intensities vary from one run to the other. In the most magnetized
runs, MHD10f05 and MHD10f01, it reaches $\simeq$100 G and about one third of
this in run MHD100f05. Run NMHDf05 presents different behaviour. For
$n>10^{10}$ cm-3, the intensity is nearly independent of $n$ and the largest
intensities is about $\simeq$30 G. This behaviour, which has been discussed
previously (e.g. Masson et al. 2016; Wurster & Li 2018) is a consequence of
ambipolar diffusion, that tends to diffuse the field. This implies that the
influence of magnetic field on the high density gas is significantly reduced
compared to ideal MHD runs.

Figure 7: The number of sinks (left panel), $N_{tot}$ and the largest sink
mass (right panel) as a function of the total accreted mass, $M_{tot}$, in
several runs.

##  5 Stellar mass spectrum

###  5.1 Fragmentation and massive stars

Figure 7 portrays the number of sink particles as a function of accreted mass
(left panel) as well as the mass of the most massive star (right panel).

The number of sinks at the end of the simulations is typically between 100 and
300 depending of the runs. As anticipated from the clump images, both magnetic
field and radiative feedback reduce fragmentation. Here we see that the
differences between runs HYDROf05 and MHD10f05 or between HYDROf01 and
MHD10f01 is about a factor of 2, the difference being more pronounced for the
two runs with $f_{acc}=0.5$. On the other hand the differences between runs
MHD10f05 and MHD10f01 is on the order of 50$\%$, showing that whereas
radiative feedback contributes to reduce fragmentation, its effect is
comparatively lower than magnetic field. Indeed, although the initial
magnetisation of run MHD100f05 is quite weak, it nevertheless reduces the
fragmentation by tens of percents compared to run HYDROf05. Interestingly run
NMHDf05, that treats ambipolar diffusion and has the same magnetisation than
run MHD10f05, presents a number of sinks similar to run MHD100f05.

In all runs but HYDROf01, two phases can be distinguished. When $M_{tot,*}$ is
smaller than $\simeq 3$ M⊙ ($\simeq 10$ for run HYDROf05), the number of sinks
increases fast and is nearly proportional to $M_{tot,*}^{m_{tot}}$ with
$m_{tot}\simeq 0.7-1$. Beyond this value, the number of sinks increases much
less rapidly, and typically $m_{tot}\simeq 0.2-0.3$. For instance for run
MHD10f05, the number of sinks has roughly doubled between the time when
$M_{tot,*}=10$ M⊙ and $M_{tot,*}=100$ M⊙. This is most certainly related to
radiative feedback and to the global increase of temperature within the
clumps. The consequence is obvioulsy that the sink particles, build their
masses in this second phase after fast fragmentation has occured.

At the end of the runs, the mass of the most massive star is between $3$ and
$10$ M⊙. The observed trends are in good agreement with the sink numbers. The
mass of the most massive star is higher when magnetic field and radiative
feedback are larger and magnetisation is comparatively slightly more efficient
than radiative feedback in producing massive stars. Two phases of growth can
also be distinguished, typically below and above $M_{tot,*}\simeq 10$ M⊙,
where $M_{max}$ grows respectively slowly and fastly. We observe that
$M_{max}\propto M_{tot,*}^{m_{max}}$ with ${m_{max}}\simeq 0.5$ when
$M_{tot,*}<10$ M⊙ while ${m_{max}}\simeq 1$ otherwise.

Note an important feature of the stellar mass distribution is that in a group
of stars which in total contains about 100-120 $M_{\odot}$, a star more
massive than 8 M⊙ is expected. We see that in our simulations only runs
MHD10f05 and MHD10f01 have reached this value. Runs HYDROf05 and MHD100f05 are
slightly below while run HYDROf01 is almost a factor 3 below. This may
constitute a hint that magnetic field is playing a role regarding the building
of the massive stars, essentially by reducing the cloud fragmentation.

HYDROf01HYDROf05NMHD10f05MHD100f05MHD10f01MHD10f05 Figure 8: Mass spectra at
various times characterised by the total accreted mass, for the six runs and
for three values of the accreted mass. The red dotted lines represent the
analytical model presented in the paper for comparison. The black dotted one
shows for reference a $M^{-1}$ power laws.

###  5.2 The sink mass function

Figure 8 displays the sink mass function, ought to represent the initial mass
function, for the six runs and 3 values of $M_{tot,*}$.

####  5.2.1 Analytical model

Before presenting the stellar mass spectra induced from the simulations, we
discuss an analytical model that will be useful to interpret the results. It
is in essence the model proposed in Hennebelle & Chabrier (2008) in which the
density PDF is the one appropriated to the gravitational collapse and stated
by Eq. (4) as proposed in Lee & Hennebelle (2018a). For the sake of
completeness, it is described in appendix B. Equations (22) and (23) are the
final equations to be used. Let us remember that the model predicts two
asymptotic behaviours. At small mass, when thermal and/or magnetic support
dominates, $\Gamma\rightarrow 0$, while at larger mass, when turbulent support
dominates, $\Gamma\rightarrow 3/4$. The transition between these two regimes
occurs at scales or equivalently masses (see Eq. 22 ) for which
thermal/magnetic and turbulent supports are comparables.

In order to be compared with the numerical simulations, one needs to specify
the values of the sound speed $c_{s}$, of the Alfvén speed, $V_{a}$, and of
the turbulent velocity dispersion, $V_{0}$. All these values can be inferred
from the results presented in § 4. Another important point when comparing
simulations with the analytical model is the normalisation. For this purpose,
we write

$\displaystyle M_{tot,*}$ $\displaystyle=$
$\displaystyle\int_{M_{min}}^{M_{max}}M{\cal N}(M)dM$ (6) $\displaystyle=$
$\displaystyle\int_{\log_{10}(M_{min})}^{\log_{10}(M_{max})}M{\cal
N}(M){M\over\log 10}d\log_{10}M.$

Thus, ${\cal N}_{0}$ as defined by Eq. 23, is determined once $M_{tot,*}$,
$M_{min}$ and $M_{max}$ are specified. The various parameters are reported in
table 2. Since $c_{s}$, $V_{a}$ and $V_{0}$ are all evolving with time and
positions, the reported values are global estimates.

Name |  $R_{c}$ (pc) |  $c_{s}$ (km s-1) |  $V_{a}$ (km s-1) |  $V_{0}$ (km s-1) |  $M_{tot,*}$ (M⊙) |  $M_{max}$ (M⊙) |  $M_{min}$ (M⊙)
---|---|---|---|---|---|---|---
NMHD10f05 | 0.3 | 0.35 | 1 | 3 | 100 | 3 | 0.1
IMHD10f05 | 0.3 | 0.35 | 1 | 3 | 150 | 8 | 0.1
IMHD10f01 | 0.3 | 0.25 | 1 | 3 | 150 | 7 | 0.1
IMHD100f05 | 0.3 | 0.35 | 0.3 | 3 | 150 | 7 | 0.1
HYDROf05 | 0.3 | 0.35 | 0 | 3 | 150 | 7 | 0.1
HYDROf01 | 0.3 | 0.25 | 0 | 3 | 150 | 3 | 0.1
Table 2: Parameters used to confront the stellar initial mass function
inferred from the simulations with the analytical model stated in § B. $R_{c}$
is the clump radius, $c_{s}$ is the typical sound speed, $V_{a}$ the Alfvén
speed, $V_{0}$ the velocity dispersion, $M_{tot,*}$ is the total accreted
mass, $M_{max}$ the mass of the most massive stars formed and $M_{min}$ the
smallest mass for which the comparison is meaningful.

We recall that the model is isothermal in nature. The sound speed may vary for
instance over time but remains uniform within the whole cloud. This has an
important consequence, which is that the model does not predict a minimum
stellar mass. One should however remember that the isothermal assumption
becomes invalid when the density reaches density on the order of $10^{10}$
cm-3 when the gas becomes progressively adiabatic. As discussed in Hennebelle
et al. (2019), the change of thermal behaviour, which leads to the formation
of the first hydrostatic core (Larson 1969), results in a peak/cut-off for the
IMF at typically several times the mass of the first hydrostatic cores,
$M_{L}$ (that we recall is about $M_{L}\simeq$0.03 M⊙). This implies that the
analytical model is valid for masses larger than a few times $M_{L}$ and this
is why we choose $M_{min}=0.1$ M⊙. The values of $M_{max}$ are taken from Fig.
7.

####  5.2.2 The hydrodynamical runs

The two top panels of Fig. 8 show results for run HYDROf05 and HYDROf01. The
mass spectra are similar to those obtained in Hennebelle et al. (2020a) with
slightly different initial conditions and less spatial resolution. Essentially
most of the sinks have their mass between few 10-2 and few M⊙. The
distributions present a plateau that ranges between $\simeq$0.1 and
$\simeq$0.5 M⊙. A relatively sharp drop occurs around 0.1 M⊙ and we get a
small number of objects at lower mass, particularly in run HYDROf05. At mass
larger than $\simeq$ 0.5 M⊙, the distribution drops following a powerlaw-like
behaviour, whose index cannot be reliably determined due to the lack of
statistics. A tentative M-1 distribution (dotted line) is represented for
comparison. This value is similar to what previous authors have inferred from
simulations (Bonnell et al. 2011b; Girichidis et al. 2011b; Ballesteros-
Paredes et al. 2015; Lee & Hennebelle 2018a, b; Padoan et al. 2020). Overall
we see that there is a good agreement between the analytical model (red dotted
line), and the sink mass distribution, or M$>0.1$ M⊙. We stress that the main
effect of increasing radiative feedback is to broaden the distribution toward
larger masses. From the analytical model, we see that this is compatible with
this being a consequence of the mean cloud temperature increasing due to the
radiative heating.

The peak of the distribution however is barely affected. This confirms, as
claimed in Hennebelle et al. (2020a) that radiative feedback is not
responsible of setting the peak of the IMF. In fact, at early time (total
accreted mass of 50 M⊙), the distribution is clearly peaked toward 0.1-0.2 M⊙
which is several time the mass of the first hydrostatic core. As time goes on,
the mass of the most massive stars increases while the number of low mass
objects remains constant or increases moderately. This is entirely compatible
with the idea that the stars inherite a minimum mass reservoir equal to a few
times the mass of the first hydrostatic core (Hennebelle et al. 2019; Colman &
Teyssier 2020), which is fastly accreted. After this, the stars keep accreting
from their mass reservoir which likely is set by gravo-turbulence (Padoan et
al. 1997; Hennebelle & Chabrier 2008; Hopkins 2012). While this process should
largely be deterministic in nature, it is also likely the case that stochastic
processes modulate this accretion as well (Bonnell et al. 2001; Basu & Jones
2004; Basu et al. 2015).

####  5.2.3 The influence of magnetic field on the stellar mass spectrum

The influence of magnetic field can be seen by comparing on one hand runs
HYDROf05, MHD100f05 and MHD10f05 and on the other hand run HYDROf01 with run
MHD10f01. Clearly, magnetic field has a significant impact on the mass
spectrum, that it tends to broaden towards larger masses. In fact, the low
mass distribution is almost unchanged. Again this provides further
confirmation that radiative feedback has no significant impact on the low mass
end of the stellar initial mass function since as discussed above magnetic
field leads to stronger radiative feedback. This also obviously shows that
magnetic field does not influence the low mass end of mass spectrum in good
agreement with the idea that it is mainly linked to the mass of the
hydrostatic core.

Run MHD10f05 presents a plateau that extends from about 0.1 $M_{\odot}$ to
$\simeq 2M_{\odot}$. It is reminiscent of run A presented in Fig 6 of Lee &
Hennebelle (2018a) and in the run presented in Fig 2 bottom panel of Jones &
Bate (2018). These runs have in common to have a high thermal energy
initially, or equivalently a low Mach number. The analytical model suggests
that when thermal support is high, a collapsing clump would indeed develop a
stellar mass spectrum $dN/d\log M\propto M^{0}$, while when turbulent support
dominates the support of the mass reservoir, $dN/d\log M\propto M^{-3/4}$ is
expected. Likely enough run MHD10f05 falls in the regime where thermal and
magnetic field dominates over turbulence at the scale of the mass reservoirs
and this explains the flat mass spectrum. This is indeed what the good
agreement with the MHD models and the simulations suggests since the broad
plateau (where $dN/d\log M\propto M^{0}$) displayed by the analytical models
is due to combination of a high Alfvén velocity and a high sound speed.

Compared to run MHD10f05, the mass spectum of run MHD10f01 presents a plateau
that is less broad. This is the case both for the numerical and the analytical
models, which are again in good agreement. This clearly is due to the lower
temperatures in run MHD10f01, which compared to run MHD10f05, leads to weaker
thermal support.

####  5.2.4 The impact of ambipolar diffusion

The mass spectrum of run NMHDf05 presents similarities with the one of run
MHD10f05 but also significant differences. Overall it is more similar to the
mass spectrum of run MHD10f01. First of all, unlike run MHD10f05, it does not
present a plateau that extends up to $\simeq$3 M⊙ but rather stops at 1 M⊙ and
the most massive stars are also less massive. This is in good agreement with
the slightly lower magnetic field which is found for run NMHDf05 (see Fig. 6)
than for run MHD10f05. The similarity with run MHD10f01 likely comes from the
total support due to both thermal and magnetic supports are closer because run
MHD10f01 has stronger field but lower temperatures than run NMHDf05.

A more surprising difference comes from the low mass objects. As can be seen
there are more sink particles of masses lower than 0.1 M⊙ in run NMHDf05 than
in the ideal MHD runs but also more than in the hydrodynamical runs. The
reason for this remains to be clarified. The most likely explanation is the
relatively weak magnetic field intensity at density above 1010 cm-3 in run
NMHDf05 compared for instance to run MHD10f05. As seen from Fig. 6, the change
of behaviour is relatively sharp, with $B$ being very comparable in runs
NMHDf05 and MHD10f05 below 1010 cm-3. Thus while in both runs, high densities
may develop due to field support, the field support drops at density above
1010 cm-3 for run NMHDf05 and this may favor fragmentation. This may also be
due to the difference in the disk populations that form in the various runs
and presented in Lebreuilly et al. (2021). The disks formed in non-ideal MHD
runs are intermediate in mass and size between the hydrodynamical disks and
the the ones which form in ideal MHD runs. While the latter are usually very
stable due to the fast growth of a toroidal magnetic component, the former
fragments but since more mass is available in bigger disks, they tend to form
bigger objects than in non ideal MHD disks.

##  6 Discussions

###  6.1 Dependence of the high-mass slope of the stellar mass spectrum

As discussed in the previous section, our numerical results suggest that from
a few solar mass to at least 7-8 M⊙, the stellar distribution presents a power
law behaviour $dN/d\log M\propto M^{-\Gamma}$, with $\Gamma\simeq 3/4$.
Analytically, this behaviour is found when at the scale of the individual mass
reservoir, $i)$ the dominant support against gravity is turbulence and $ii)$
when the density PDF is $\propto\rho^{-3/2}$ which is a consequence of
gravitational collapse. On the other-hand, when the density PDF is close to a
lognormal distribution, we do expect $\Gamma\simeq 1.3$ as discussed in Lee &
Hennebelle (2018a). In essence, the density PDF is a direct estimate of how
the gas mass is distributed amongst densities and therefore controls the
number of density fluctuations at a given density. Typically a log-normal
distribution has less dense gas than a PDF $\propto\rho^{-3/2}$ and therefore
less small mass objects are produced with the former than with the latter. The
transition between the two exponents, $\Gamma=3/4$ and $\Gamma=1.3$, is
expected to occur at the density, $n_{trans}$, which typically connects the
turbulent log-normal PDF to the power law $\propto\rho^{-3/2}$ gravitational
PDF. In the present simulations this occurs around $\simeq 10^{6}$ cm-3.
Combining Eqs. (19) and (22), we can estimate the mass, $M_{trans}$ it
corresponds to

$\displaystyle M_{trans}\simeq 15\;M_{\odot}\left({V_{0}\over 3\;{\rm
km\,s^{-1}}}\right)^{6}\left({R_{c}\over 0.3{\rm
pc}}\right)^{-3}\left({n_{trans}\over 10^{6}\;{\rm cm^{-3}}}\right)^{-2}$
$\displaystyle\simeq 15\;M_{\odot}\left({M\over
10^{3}\;{M_{\odot}}}\right)^{3}\left({R_{c}\over 0.3{\rm
pc}}\right)^{-6}\left({n_{trans}\over 10^{6}\;{\rm cm^{-3}}}\right)^{-2}.$ (7)

where for simplicity, we have assumed $\eta=0.5$ and $V_{0}^{2}\simeq
GM/R_{c}$. Because of the sixth power which appears for $V_{0}$ or $R_{c}$,
the clump radius, the value of $M_{trans}$ is clearly not accurate and likely
can abruptly change from one environment to another. Typically we expect a
fast transition around $R_{c}\simeq 0.3$ pc. It is however illustrative and
shows that for our simulations, at high mass, the mass spectra are expected to
be mostly if not exclusively described by the $\Gamma=3/4$ exponent since our
stellar masses are smaller than 15 M⊙. It also shows that in a less dense and
compact clump, the transition should occur at smaller masses since the value
of $V_{0}$, or equivalently the value of $R_{c}$, should be smaller. While
most of the studies which have started from massive clumps, comparable to the
ones studied here, tend to present $\Gamma$ lower than the canonical Salpeter
exponent (see the discussion in Lee & Hennebelle 2018a), works in which the
IMF is obtained from larger scale clouds studies which have attempted to
obtain the IMF in larger scale simulations with initial conditions that
correspond to more standard giant molecular clouds, generally report $\Gamma$
values that are closer to 1.3. This is the case for instance for the run XL-F
presented in Fig. 4 of He et al. (2019) and the run presented in Fig. 3 of
Padoan et al. (2020) for masses between 10 and 50 M⊙, respectively. This is
also the case for the runs presented in Ntormousi & Hennebelle (2019) and the
core mass function extracted from these simulations (Louvet et al. 2021).

###  6.2 Observationally inferred mass distribution in actively star forming
regions

While it may sound at first surprising not to find $\Gamma\simeq 1.3$, which
is the slope inferred by Salpeter (1955), it should be stressed that recent
observations have been inferring that in some actively star forming regions,
the IMF may indeed be top-heavy (Zhang et al. 2018; Lee et al. 2020). More
precisely, in the Arches cluster Hosek et al. (2019) inferred $\Gamma\simeq
0.8$. On the other hand, recent studies of the core mass function also
obtained within massive star forming regions, have also inferred power law
behaviours with indices $\Gamma\simeq 0.95$ (Motte et al. 2018; Pouteau et al.
2022). As cores are widely assumed to be the progenitors of stars out of which
they build their mass, the inferred $\Gamma$ are compatible with the idea that
the shape of the IMF in massive star forming regions is inherited from the
shape of the CMF, at least at high masses, although eventually it should be
compared with the IMF of the very same region.

While more detailed investigations, including careful comparisons between
simulations and observations must be carried out before firm conclusion can be
drawn, there is a clear suggestion coming from both observations and theories
that systematic variations of the IMF may occur, particularly in very compact
star forming regions.

###  6.3 Limits of the present work and the universality of the IMF

Our work presents several important limits that need to be discussed. Indeed,
one of the conclusion is that the combination of magnetic field and radiative
transfer possibly leads to more variability that what observational inferences
of the IMF may have led to conclude. Admittedly, this question even for our
own Galaxy remains difficult to address, particularly because of the
relatively limited samples that are often available but it seems nevertheless
unavoidable that at least some level of fluctuations should be present (see
for instance the comprehensive discussion provided in Dib 2022).

Determining whether the variations observed in the present work are compatible
with the galactic fluctuations of the IMF, is beyond the scope of the present
paper but it is worth to remind that an important source of variations is due
to the efficiency of the accretion luminosity expressed by the parameter,
$f_{acc}$. Whereas there may be some variability of $f_{acc}$, likely enough
it is not a factor of 5 as we have been exploring here. The other possibly
extreme variations we have considered is magnetic intensity since we have
explored a factor of 10 (and even go to pure hydrodynamical cases). This is
not well constraint yet but a 1000 M⊙ clump is a relatively large ensemble and
it is unclear what are the variations of the magnetisation in the galactic
populations.

Finally, we stress that in this work a possibly important process has been
omitted, namely the protostellar jets. Recently, Guszejnov et al. (2021) have
been exploring their impact in simulations comparable to the ones presented
here (with a resolution of few tens of AU). They concluded that protostellar
jets may be playing a significant role in setting the IMF in particular for
the formation of low mass objects in the presence of a significant initial
magnetic field. Whether this process may help explaining the universality of
the IMF is however not clear yet.

##  7 Conclusions

With the goal of understanding how magnetic field and radiative feedback
influence the collapse and the fragmentation of a massive star forming clump,
we have performed high resolution adaptive mesh calculations with a spatial
resolution down to about 1 AU. Six runs in which 2 radiative feedback
efficiencies, 3 magnetic intensities as well as the impact of non-ideal MHD
are explored. We show that the physical characteristics of the simulated star
forming clumps compare well with various observations. This is for instance
the case for the observational bolometric luminosities that we compared with
the total luminosities of the sink particles produced in the simulations as
well as for the gas temperatures. For the latter, we develop an analytical
model which agrees well with the temperatures inferred from the simulations.

The stellar mass spectra of the six runs are analysed in detail and compared
with an analytical model in which thermal, magnetic and turbulent supports are
playing a major role. Overall the analytical model reproduces well the
numerical mass spectra for masses above $\simeq 0.1$ M⊙. At this mass which
corresponds to a few times the mass of the first hydrostatic core the
underlying gas thermodynamics is nearly adiabatic and specific models should
be considered (e.g. Hennebelle et al. 2019). The combination between
simulations and analytical results allows us to clearly assess the role and
influence of each physical process which are as follows:

  * -

in the density range at which the gas is not adiabatic, the density PDF which
is $\propto\rho^{-3/2}$ is deeply shaping the stellar mass spectrum and leads
to two physical distinct regimes for the mass spectra.

  * -

at masses larger than $\simeq$ 0.1 M⊙, thermal pressure and magnetic field may
lead to a flat mass spectrum, i.e. $dN/d\log M\propto M^{-\Gamma}$ with
$\Gamma\simeq 0$ if they are strong enough compared to turbulence.

  * -

at larger scales, turbulence dominates and may lead to a mass spectrum with
$\Gamma\simeq 3/4$. At even larger scales and lower density, the PDF is
expected to be log-normal in shape and stiffer mass spectra, with larger
$\Gamma$ are expected.

  * -

the transition between the regime with $\Gamma\simeq 0$ and $\Gamma\simeq 3/4$
is not universal and depends on the local physical processes such as thermal
support, magnetic field and Mach number.

Generally speaking, we find that the main effect of magnetic field and
radiative transfer is to reduce the total number of fragments and to increase
the mass of the most massive stars. These latter have been found to increase
with the magnetic intensity and the radiation feedback efficiency. For
instance, in the present work we found that for the hydrodynamical simulation
with the lowest efficiency, the most massive star produced after 150 M⊙ have
been accreted, is about 3 M⊙. With a higher radiative feedback efficiency or a
sufficiently strong initial field, stars of masses 7-8 M⊙ are produced. We
therefore conclude that whereas magnetic field and radiative feedback may not
be essential to explain the peak or the various slope values of the IMF, they
may be essential to reproduce the exact shape (like the transition between the
various regimes), the level of fragmentation i.e. the number of stars formed,
and the mass of the most massive stars.

###### Acknowledgements.

We thank the anonymous referee for a useful report. This work was granted
access to HPC resources of CINES and CCRT under the allocation x2014047023
made by GENCI (Grand Equipement National de Calcul Intensif). This research
has received funding from the European Research Council synergy grant ECOGAL
(Grant : 855130). G.A.F also acknowledges support from the Collaborative
Research Centre 956, funded by the Deutsche Forschungsgemeinschaft (DFG)
project ID 184018867.

## References

  * Abe et al. (2021) Abe, D., Inoue, T., Inutsuka, S.-i., & Matsumoto, T. 2021, ApJ, 916, 83
  * Ballesteros-Paredes et al. (2015) Ballesteros-Paredes, J., Hartmann, L. W., Pérez-Goytia, N., & Kuznetsova, A. 2015, MNRAS, 452, 566
  * Bastian et al. (2010) Bastian, N., Covey, K. R., & Meyer, M. R. 2010, ARA&A, 48, 339
  * Basu et al. (2015) Basu, S., Gil, M., & Auddy, S. 2015, MNRAS, 449, 2413
  * Basu & Jones (2004) Basu, S. & Jones, C. E. 2004, MNRAS, 347, L47
  * Bate (2009) Bate, M. R. 2009, MNRAS, 392, 1363
  * Bate et al. (2003) Bate, M. R., Bonnell, I. A., & Bromm, V. 2003, MNRAS, 339, 577
  * Bhandare et al. (2020) Bhandare, A., Kuiper, R., Henning, T., et al. 2020, A&A, 638, A86
  * Bleuler & Teyssier (2014) Bleuler, A. & Teyssier, R. 2014, MNRAS, 445, 4015
  * Bonnell et al. (2001) Bonnell, I. A., Bate, M. R., Clarke, C. J., & Pringle, J. E. 2001, MNRAS, 323, 785
  * Bonnell et al. (2011a) Bonnell, I. A., Smith, R. J., Clark, P. C., & Bate, M. R. 2011a, MNRAS, 410, 2339
  * Bonnell et al. (2011b) Bonnell, I. A., Smith, R. J., Clark, P. C., & Bate, M. R. 2011b, MNRAS, 410, 2339
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Colman & Teyssier (2020) Colman, T. & Teyssier, R. 2020, MNRAS, 492, 4727
  * Commerçon et al. (2014) Commerçon, B., Debout, V., & Teyssier, R. 2014, A&A, 563, A11
  * Commerçon et al. (2021) Commerçon, B., González, M., Mignon-Risse, R., Hennebelle, P., & Vaytet, N. 2021, arXiv e-prints, arXiv:2109.10580
  * Commerçon et al. (2011a) Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011a, A&A, 529, A35
  * Commerçon et al. (2011b) Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011b, A&A, 529, A35
  * Crutcher (2012) Crutcher, R. M. 2012, ARA&A, 50, 29
  * Dib (2022) Dib, S. 2022, arXiv e-prints, arXiv:2204.09064
  * Elia et al. (2021) Elia, D., Merello, M., Molinari, S., et al. 2021, MNRAS, 504, 2742
  * Elia et al. (2017) Elia, D., Molinari, S., Schisano, E., et al. 2017, MNRAS, 471, 100
  * Federrath (2016) Federrath, C. 2016, MNRAS, 457, 375
  * Federrath et al. (2008) Federrath, C., Klessen, R. S., & Schmidt, W. 2008, ApJ, 688, L79
  * Fromang et al. (2006) Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371
  * Girichidis et al. (2011a) Girichidis, P., Federrath, C., Banerjee, R., & Klessen, R. S. 2011a, MNRAS, 413, 2741
  * Girichidis et al. (2011b) Girichidis, P., Federrath, C., Banerjee, R., & Klessen, R. S. 2011b, MNRAS, 413, 2741
  * Gómez et al. (2021) Gómez, G. C., Vázquez-Semadeni, E., & Palau, A. 2021, MNRAS, 502, 4963
  * González et al. (2007) González, M., Audit, E., & Huynh, P. 2007, A&A, 464, 429
  * Guszejnov et al. (2020) Guszejnov, D., Grudić, M. Y., Hopkins, P. F., Offner, S. S. R., & Faucher-Giguère, C.-A. 2020, MNRAS, 496, 5072
  * Guszejnov et al. (2021) Guszejnov, D., Grudić, M. Y., Hopkins, P. F., Offner, S. S. R., & Faucher-Giguère, C.-A. 2021, MNRAS, 502, 3646
  * Haugbølle et al. (2018) Haugbølle, T., Padoan, P., & Nordlund, Å. 2018, ApJ, 854, 35
  * He et al. (2019) He, C.-C., Ricotti, M., & Geen, S. 2019, MNRAS, 489, 1880
  * Hennebelle (2013) Hennebelle, P. 2013, A&A, 556, A153
  * Hennebelle (2018) Hennebelle, P. 2018, A&A, 611, A24
  * Hennebelle (2021) Hennebelle, P. 2021, A&A, 655, A3
  * Hennebelle & Chabrier (2008) Hennebelle, P. & Chabrier, G. 2008, ApJ, 684, 395
  * Hennebelle et al. (2020a) Hennebelle, P., Commerçon, B., Lee, Y.-N., & Chabrier, G. 2020a, ApJ, 904, 194
  * Hennebelle et al. (2020b) Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020b, A&A, 635, A67
  * Hennebelle & Inutsuka (2019) Hennebelle, P. & Inutsuka, S.-i. 2019, Frontiers in Astronomy and Space Sciences, 6, 5
  * Hennebelle et al. (2019) Hennebelle, P., Lee, Y.-N., & Chabrier, G. 2019, ApJ, 883, 140
  * Hopkins (2012) Hopkins, P. F. 2012, MNRAS, 423, 2037
  * Hosek et al. (2019) Hosek, Matthew W., J., Lu, J. R., Anderson, J., et al. 2019, ApJ, 870, 44
  * Hosokawa & Omukai (2009) Hosokawa, T. & Omukai, K. 2009, ApJ, 691, 823
  * Jaura et al. (2018) Jaura, O., Glover, S. C. O., Klessen, R. S., & Paardekooper, J. P. 2018, MNRAS, 475, 2822
  * Jones & Bate (2018) Jones, M. O. & Bate, M. R. 2018, MNRAS, 478, 2650
  * Klessen (2001) Klessen, R. S. 2001, ApJ, 556, 837
  * Kritsuk et al. (2011) Kritsuk, A. G., Norman, M. L., & Wagner, R. 2011, ApJ, 727, L20
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
  * Krumholz et al. (2007) Krumholz, M. R., Klein, R. I., & McKee, C. F. 2007, ApJ, 656, 959
  * Krumholz et al. (2012) Krumholz, M. R., Klein, R. I., & McKee, C. F. 2012, ApJ, 754, 71
  * Kuiper et al. (2010) Kuiper, R., Klahr, H., Dullemond, C., Kley, W., & Henning, T. 2010, A&A, 511, A81
  * Kuiper & Yorke (2013) Kuiper, R. & Yorke, H. W. 2013, ApJ, 772, 61
  * Lane et al. (2022) Lane, H. B., Grudić, M. Y., Guszejnov, D., et al. 2022, MNRAS, 510, 4767
  * Larson (1969) Larson, R. B. 1969, MNRAS, 145, 271
  * Lebreuilly et al. (2021) Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2021, ApJ, 917, L10
  * Lee & Hennebelle (2018a) Lee, Y.-N. & Hennebelle, P. 2018a, A&A, 611, A88
  * Lee & Hennebelle (2018b) Lee, Y.-N. & Hennebelle, P. 2018b, A&A, 611, A89
  * Lee & Hennebelle (2019) Lee, Y.-N. & Hennebelle, P. 2019, A&A, 622, A125
  * Lee et al. (2020) Lee, Y.-N., Offner, S. S. R., Hennebelle, P., et al. 2020, Space Sci. Rev., 216, 70
  * Li (2018) Li, G.-X. 2018, MNRAS, 477, 4951
  * Li et al. (2018) Li, P. S., Klein, R. I., & McKee, C. F. 2018, MNRAS, 473, 4220
  * Li et al. (2015) Li, P. S., McKee, C. F., & Klein, R. I. 2015, MNRAS, 452, 2500
  * Lin et al. (2022) Lin, Y., Wyrowski, F., Liu, H. B., et al. 2022, A&A, 658, A128
  * Louvet et al. (2021) Louvet, F., Hennebelle, P., Men’shchikov, A., et al. 2021, A&A, 653, A157
  * Marchand et al. (2016) Marchand, P., Masson, J., Chabrier, G., et al. 2016, A&A, 592, A18
  * Masson et al. (2016) Masson, J., Chabrier, G., Hennebelle, P., Vaytet, N., & Commerçon, B. 2016, A&A, 587, A32
  * Masson et al. (2012) Masson, J., Teyssier, R., Mulet-Marquis, C., Hennebelle, P., & Chabrier, G. 2012, ApJS, 201, 24
  * Mathew & Federrath (2020) Mathew, S. S. & Federrath, C. 2020, MNRAS[arXiv:2007.01875]
  * Menon et al. (2022) Menon, S. H., Federrath, C., Krumholz, M. R., et al. 2022, MNRAS, 512, 401
  * Mestel & Spitzer (1956) Mestel, L. & Spitzer, L., J. 1956, MNRAS, 116, 503
  * Mignon-Risse et al. (2020) Mignon-Risse, R., González, M., Commerçon, B., & Rosdahl, J. 2020, A&A, 635, A42
  * Mignon-Risse et al. (2021) Mignon-Risse, R., González, M., Commerçon, B., & Rosdahl, J. 2021, A&A, 652, A69
  * Minerbo (1978) Minerbo, G. N. 1978, J. Quant. Spec. Radiat. Transf., 20, 541
  * Molina et al. (2012) Molina, F. Z., Glover, S. C. O., Federrath, C., & Klessen, R. S. 2012, MNRAS, 423, 2680
  * Motte et al. (2018) Motte, F., Nony, T., Louvet, F., et al. 2018, Nature Astronomy, 2, 478
  * Murray & Chang (2015) Murray, N. & Chang, P. 2015, ApJ, 804, 44
  * Myers et al. (2013) Myers, A. T., McKee, C. F., Cunningham, A. J., Klein, R. I., & Krumholz, M. R. 2013, ApJ, 766, 97
  * Myers & Basu (2021) Myers, P. C. & Basu, S. 2021, ApJ, 917, 35
  * Ntormousi & Hennebelle (2019) Ntormousi, E. & Hennebelle, P. 2019, A&A, 625, A82
  * Offner et al. (2014) Offner, S. S. R., Clark, P. C., Hennebelle, P., et al. 2014, Protostars and Planets VI, 53
  * Offner et al. (2009) Offner, S. S. R., Klein, R. I., McKee, C. F., & Krumholz, M. R. 2009, ApJ, 703, 131
  * Padoan et al. (1997) Padoan, P., Nordlund, A., & Jones, B. J. T. 1997, MNRAS, 288, 145
  * Padoan et al. (2020) Padoan, P., Pan, L., Juvela, M., Haugbølle, T., & Nordlund, Å. 2020, ApJ, 900, 82
  * Pattle et al. (2022) Pattle, K., Fissel, L., Tahani, M., Liu, T., & Ntormousi, E. 2022, arXiv e-prints, arXiv:2203.11179
  * Peter et al. (2022) Peter, T., Klessen, R. S., Kanschat, G., Glover, S. C. O., & Bastian, P. 2022, arXiv e-prints, arXiv:2207.12848
  * Peters et al. (2011) Peters, T., Banerjee, R., Klessen, R. S., & Mac Low, M.-M. 2011, ApJ, 729, 72
  * Peters et al. (2010) Peters, T., Klessen, R. S., Mac Low, M.-M., & Banerjee, R. 2010, ApJ, 725, 134
  * Pouteau et al. (2022) Pouteau, Y., Motte, F., Nony, T., et al. 2022, arXiv e-prints, arXiv:2203.03276
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161
  * Saumon & Chabrier (1992) Saumon, D. & Chabrier, G. 1992, Phys. Rev. A, 46, 2084
  * Saumon et al. (1995) Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713
  * Semenov et al. (2003) Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611
  * Shu (1977) Shu, F. H. 1977, ApJ, 214, 488
  * Smith et al. (2014) Smith, R. J., Glover, S. C. O., & Klessen, R. S. 2014, MNRAS, 445, 2900
  * Teyssier (2002) Teyssier, R. 2002, A&A, 385, 337
  * Urban et al. (2010) Urban, A., Martel, H., & Evans, Neal J., I. 2010, ApJ, 710, 1343
  * Urquhart et al. (2014) Urquhart, J. S., Moore, T. J. T., Csengeri, T., et al. 2014, MNRAS, 443, 1555
  * Vaytet et al. (2013) Vaytet, N., Chabrier, G., Audit, E., et al. 2013, A&A, 557, A90
  * Vázquez-Semadeni (1994) Vázquez-Semadeni, E. 1994, ApJ, 423, 681
  * Wurster & Li (2018) Wurster, J. & Li, Z.-Y. 2018, Frontiers in Astronomy and Space Sciences, 5, 39
  * Xu et al. (2019) Xu, S., Ji, S., & Lazarian, A. 2019, ApJ, 878, 157
  * Zhang et al. (2018) Zhang, Z.-Y., Romano, D., Ivison, R. J., Papadopoulos, P. P., & Matteucci, F. 2018, Nature, 558, 260
  * Zhao et al. (2020) Zhao, B., Tomida, K., Hennebelle, P., et al. 2020, Space Sci. Rev., 216, 43

##  Appendix A A simple analytical model for the temperature

To better understand the temperature profils through the clump, we make use of
the simple model discussed in Hennebelle et al. (2020a), improving on various
aspects. The model assumes that the cloud is spherically symmetric and that
all sources are located in the clump center. As seen in § 4.1, the radial
profil of the density field in collapsing envelopes is given by

$\displaystyle\rho(r)={\delta_{\rho}c_{s,0}^{2}\over 2\pi Gr^{2}},$ (8)

where $\delta_{\rho}$ is a dimensionless factor which typically is equal to
$\simeq$30-100. We assume that gas and dust have the same temperature and are
stationary. The grey body approximation that is being used, leads when the
medium is optically thick to

$\displaystyle-4\pi r^{2}{c\over
3\kappa(T)\rho(r)}\partial_{r}(aT^{4})=\sum(L_{*}+L_{acc}).$ (9)

where $a$ is the radiation density constant. In this expression we assume that
all the emitting sources are located in the clump center. The opacity
temperature dependences (e.g. Semenov et al. 2003), suggest that we can
distinguish two regimes of temperature,

$\displaystyle\kappa(T)=\kappa_{0}\simeq 5\,{\rm cm}^{2}{\rm g}^{-1}\;{\rm
for}\;T>T_{\rm crit}\simeq 100\,{\rm K},$ (10)
$\displaystyle\kappa(T)\simeq\kappa_{0}\left({T\over T_{\rm
crit}}\right)^{\alpha}\;{\rm for}\;T<T_{\rm crit}.$

where $\alpha$ is typically between 1 and 2. In this work we adopted
$\alpha=1.5$. Combining Eqs. (8), (9) and (10), we get

$\displaystyle T(r)=\left(T_{\rm crit}^{4}+K\left({1\over r^{3}}-{1\over
r_{\rm crit}^{3}}\right)\right)^{1/4}\;{\rm for}\;T>T_{\rm crit},$
$\displaystyle T(r)=\left(KT_{\rm crit}^{-\alpha}{4-\alpha\over 4}{1\over
r^{3}}\right)^{1/(4-\alpha)}\;{\rm for}\;T<T_{\rm crit}\;{\rm and}\;r<r_{\rm
ext},$ (11) $\displaystyle K=f_{\rm
acc}\delta_{\rho}{3\kappa_{0}c_{s,0}^{2}\over
24\pi^{2}acG}\sum(L_{*}+L_{acc}),$ (12)

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

$\displaystyle r_{\rm crit}^{3}=K{4-\alpha\over 4}{T_{\rm crit}^{-\alpha}\over
T_{\rm crit}^{4-\alpha}-T_{\rm ext}^{4-\alpha}}.$ (13)

Finally, $T_{ext}$ is the temperature where the optical depth is about 1 and
$r_{ext}$ the corresponding radii. We therefore have $\kappa(T_{\rm
ext})\rho(r_{\rm ext})r_{\rm ext}\simeq\tau_{0}$ where $\tau_{0}$ should be on
the order of 1. Combining Eqs. (8), (10), (A), we obtain for $r_{\rm ext}$

$\displaystyle r_{\rm ext}=\left({\kappa_{0}\delta c_{s,0}^{2}\over 2\pi
G\tau_{0}}\right)^{(4-\alpha)/(4+2\alpha)}\left({K\over T_{\rm
crit}^{4}}{4-\alpha\over 4}\right)^{\alpha/(4+2\alpha)}.$ (14)

At this point, the radiative flux becomes simply equal to the term $cE_{R}$
and Eq. (9) becomes invalid. Under the assumption that the temperature remains
the one of a blackbody, it then remains constant at larger radii and thus

$\displaystyle T(r)=T_{\rm ext}=\left(KT_{\rm crit}^{-\alpha}{4-\alpha\over
4}{1\over r_{\rm ext}^{3}}\right)^{1/(4-\alpha)}\;\rm{for}\;r>r_{\rm ext}$
(15)

The expression for $T_{\rm ext}$ is obtained by continuity at $r_{\rm ext}$.
By combining Eq. (15) and Eq. (14), we find that

$\displaystyle T_{\rm ext}=\left(\pi(1-{\alpha\over 4}){L_{tot}\over
ac}{T_{crit}^{2\alpha}G^{2}\over\kappa_{0}^{2}\delta^{2}c_{s}^{4}}\tau_{0}^{3}\right)^{1/(4+2\alpha)}$
(16)

##  Appendix B Analytical model of the mass spectrum

For completeness, we describe here the analytical model developed in
Hennebelle & Chabrier (2008) and Lee & Hennebelle (2018a) that we use in the
paper to interpret the numerical results.

It is based on the equality of mass of the density fluctuation which are
unstable at scale $R$ (left-hand term) and the mass that ends up into the
structures, i.e. the stars:

$\displaystyle{M_{\rm tot}(R)\over
V_{\mathrm{c}}}=\int\limits^{\infty}_{\delta_{R}^{\mathrm{c}}}\overline{\rho}\exp(\delta){\cal
P}_{R}(\delta)d\delta=\int\limits_{0}^{M_{R}^{\mathrm{c}}}M^{\prime}{\cal
N}\\!(M^{\prime})P(R,\\!M^{\prime})\,dM^{\prime},$ (17)

where $\delta=\ln(\rho/\bar{\rho})$, ${\cal P}_{R}$ is the density PDF,
$P(R,M)$ is the probability of finding a self-gravitating clump of mass
$M^{\prime}$ embedded into a self-gravitating clump of mass $M_{R}$ unstable
at scale $R$. It is assumed to be 1.

Taking the derivative with respect to $R$, we get

$\displaystyle{\cal N}(M_{R}^{c})={\overline{\rho}\over M_{R}^{c}}{dR\over
dM_{R}^{c}}\,\left(-{d\delta_{R}^{c}\over dR}\exp(\delta_{R}^{c}){\cal
P}_{R}(\delta_{R}^{c})\right).$ (18)

The mass of the density fluctuations is given by

$\displaystyle M=C_{m}\rho R^{3},\leavevmode\nobreak\ \text{where typically
$C_{m}=4\pi/3$.}$ (19)

Here we assume that the density PDF is given by

$\displaystyle{\cal P}_{R}(\rho)={\cal
P}_{0}\left({\rho\over\rho_{0}}\right)^{-1.5},$ (20)

The gravitational instability criterion for a clump of mass $M$ at scale $R$
is

$\displaystyle M>M_{J}=a_{J}{\Bigl{[}c_{s}^{2}+{V_{a}^{2}\over
6}+{V_{0}^{2}\over 3}\left({R\over R_{c}}\right)^{2\eta}\Bigr{]}^{3\over
2}\over\sqrt{G^{3}\overline{\rho}\exp(\delta)}},$ (21)

where $c_{s}$ is the sound speed, $V_{a}$ the Alfvén speed, $V_{0}$ the rms
velocity dispersion at the cloud scale, $R_{c}$ is the cloud radius and $\eta$
an exponent to describe the turbulent scale dependence. Typically
$\eta=0.3-0.5$ and in this work the value $\eta=0.5$ is assumed for
simplicity. Equation (21) is the standard Jeans mass expression in which the
support is assumed to be as suggested by the virial theorem. Note that the
surface terms are not taken into account, they would typically modify this
expression by a factor of 2. Taking the standard definition of the Jeans mass,
the mass enclosed in a sphere of diameter equal to the Jeans length, we get
$a_{J}=\pi^{5/2}/6$. With Eq. (19), this implies

$\displaystyle M_{R}^{c}={a_{J}^{2\over 3}C_{m}^{1\over 3}\over
G}\left(c_{s}^{2}R+{V_{a}^{2}\over 6}R+{V_{0}^{2}\over 3\,G}\left({R\over
R_{c}}\right)^{2\eta}R\right),$ (22)

where $M_{R}^{c}$ is the critical mass at scale $R$.

With Eq. (20), Eq. (18) leads to

$\displaystyle{\cal N}(M_{R}^{c})={\cal N}_{0}\left({R\over
M_{R}^{c}}\right)^{3/2}{dR\over dM_{R}^{c}}\,\left(-{1\over
M_{R}^{c}}{dM_{R}^{c}\over dR}+{3\over R}\right).$ (23)

Knowing the cloud physical conditions, $c_{s}$, $V_{a}$, $V_{0}$, together
with Eq. (22), Eq. (23) allow to predict the stellar mass spectrum. The
normalisation coefficient ${\cal N}_{0}$ is determined by specifying the total
mass within stars.

It is useful to see that

$\displaystyle M\rightarrow 0\Leftrightarrow{\cal N}\rightarrow
M^{-1}\Leftrightarrow{dN\over d\log M}\rightarrow M^{0}.$ (24)

In this limit, the mass reservoir is thermally supported and the mass spectra
present a plateau, i.e. ${dN\over d\log M}\propto M^{0}$.

On the other hand, in the limit

$\displaystyle M\rightarrow\infty\Leftrightarrow{\cal N}\rightarrow
M^{3/(4\eta+2)-5/2}\Leftrightarrow{dN\over d\log M}\rightarrow
M^{3/(4\eta+2)-3/2}.$ (25)

As revealed by Eq. (22), in this limit the mass reservoir is dominated by the
turbulent dispersion. For $\eta=0.5$, the mass spectrum is ${dN\over d\log
M}\propto M^{-3/4}$.

We recall that at small masses, the asymptotic behaviour will eventually break
down when the gas becomes adiabatic due to the dust opacity and the formation
of the first hydrostatic core, while at large masses, the assumption of the
density PDF being $\propto\rho^{-3/2}$, is eventually invalid (typically it
eventually turns into a log-normal distribution). Therefore while useful,
these asymptotic behaviours must be handled with care.

##  Appendix C Accretion and stellar luminosities

To get a better understanding of the origins of the luminosities, we
investigate the stellar and accretion luminosities separately. The two panels
of Fig. 9 show the sum of the stellar luminosities (left panel) and the sum of
the accretion luminosities (right panel). In a first phase, up to time $\simeq
0.1$ Myr, the accretion luminosity largely dominates. Then as stars of few
solar masses have formed, the stellar luminosities increase steeply and then
reach values comparable to the accretion luminosities.

Figure 9: This figure complements Fig. 1. Left panel displays the total
stellar luminosities and right one the total accretion luminosity.