Simulating Star Clusters Across Cosmic Time: I. Initial Mass Function, Star Formation Rates and Efficiencies

TL;DR

This study uses advanced simulations to explore star formation in molecular clouds, reproducing observed stellar initial mass functions and star formation laws, and identifying conditions for globular cluster formation.

Contribution

It introduces detailed radiation-magneto-hydrodynamic simulations that replicate observed stellar mass distributions and star formation efficiencies across different cloud densities and metallicities.

Findings

Reproduces Salpeter slope of the initial mass function.

Identifies star formation laws depending on gas density.

Suggests conditions for globular cluster progenitors.

Abstract

We present radiation-magneto-hydrodynamic simulations of star formation in self-gravitating, turbulent molecular clouds, modeling the formation of individual massive stars, including their UV radiation feedback. The set of simulations have cloud masses between $m_{\rm gas}=10^3$~M$_\odot$ to $3 \times 10^5$~M$_\odot$ and gas densities typical of clouds in the local universe ($\overline n_{\rm gas} \sim 1.8\times 10^2$~cm$^{-3}$) and 10$\times$ and 100$\times$ denser, expected to exist in high-redshift galaxies. The main results are: {\it i}) The observed Salpeter power-law slope and normalisation of the stellar initial mass function at the high-mass end can be reproduced if we assume that each star-forming gas clump (sink particle) fragments into stars producing on average a maximum stellar mass about $40\%$ of the mass of the sink particle, while the remaining $60\%$ is distributed into smaller mass stars. Assuming that the sinks fragment according to a power-law mass function flatter than Salpeter, with log-slope $0.8$, satisfy this empirical prescription. {\it ii}) The star formation law that best describes our set of simulation is $d\rho_*/dt \propto \rho_{gas}^{1.5}$ if $\overline n_{gas}<n_{cri}\approx 10^3$~cm$^{-3}$, and $d\rho_*/dt \propto \rho_{\rm gas}^{2.5}$ otherwise. The duration of the star formation episode is roughly $6$ cloud's sound crossing times (with $c_s=10$~km/s). {\it iii}) The total star formation efficiency in the cloud is $f_*=2\% (m_{\rm gas}/10^4~M_\odot)^{0.4}(1+\overline n_{\rm gas}/n_{\rm cri})^{0.91}$, for gas at solar metallicity, while for metallicity $Z<0.1$~Z$_\odot$, based on our limited sample, $f_*$ is reduced by a factor of $\sim 5$. {\it iv)} The most compact and massive clouds appear to form globular cluster progenitors, in the sense that star clusters remain gravitationally bound after the gas has been expelled.

Tables

Table 1. Table 1: Initial conditions of our 16 simulations.

	$m_{\mathrm{gas}}\left(M_{\odot }\right)$ d	$1.0\times {10}^3$	$3.2\times {10}^3$	$1.0\times {10}^4$	$3.2\times {10}^4$	$1.0\times {10}^5$	$3.2\times {10}^5$
	Cloud Name e		XS-F	S-F	M-F	L-F	XL-F
${\overline{n}}_{\mathrm{gas}}$a $=1.8\times {10}^2$cm-3	$r_{\mathrm{gas}}$ (pc) f		5.0	7.3	11	16	23
	$\Sigma$ ($M_{\odot }$ pc-2) x		41	61	89	131	193
	$v_{\mathrm{esc}}$ (km/s) x		2.3	3.4	5.1	7.4	11
	$\Delta x_{\min }$ (AU) g		500	730	1100	1600	2300
$t_{\mathrm{ff}}$b $=4.4$Myr	$n_{\mathrm{sink}}$ (cm-3) h		$1.2\times {10}^7$	$5.6\times {10}^6$	$2.6\times {10}^6$	$1.2\times {10}^6$	$5.6\times {10}^5$
	$M_J(M_{\odot }$) i		0.3	0.4	0.6	0.9	1.3
$l_{\max }$c $=15$	$ℳ$ j		4.6	6.8	10	15	22
	$t_{\mathrm{cr}}$ (Myr) k		0.5	0.7	1.1	1.5	2.3
	$Z/Z_{\odot }$ l		1	1	1	1	1
	Cloud Name		XS-C	S-C	M-C	L-C, L-C-lm, L-C-xlm m
${\overline{n}}_{\mathrm{gas}}=1.8\times {10}^3$cm-3	$r_{\mathrm{gas}}$		2.3	3.4	5.0	7.3
	$\Sigma$		193	283	415	609
	$v_{\mathrm{esc}}$		3.4	5.1	7.4	11
	$\Delta x_{\min }$		460	680	1000	1500
$t_{\mathrm{ff}}=1.4$Myr	$n_{\mathrm{sink}}$		$1.4\times {10}^7$	$6.5\times {10}^6$	$3.0\times {10}^6$	$1.4\times {10}^6$
	$M_J$		0.3	0.4	0.6	0.8
$l_{\max }=14$	$ℳ$		6.8	10	15	22
	$t_{\mathrm{cr}}$		0.23	0.33	0.5	0.7
	$Z$		1	1	1	1, $1/10$, $1/40$
	Cloud Name	XXS-VC	XS-VC	S-VC	M-VC	L-VC
${\overline{n}}_{\mathrm{gas}}=1.8\times {10}^4$cm-3	$r_{\mathrm{gas}}$	0.7	1.1	1.6	2.3	3.4
	$\Sigma$	609	894	1312	1925	2827
	$v_{\mathrm{esc}}$	3.4	5.1	7.4	11	16
	$\Delta x_{\min }$	150	220	320	460	680
$t_{\mathrm{ff}}=0.44$Myr	$n_{\mathrm{sink}}$	$1.4\times {10}^8$	$6.5\times {10}^7$	$3.0\times {10}^7$	$1.4\times {10}^7$	$6.5\times {10}^6$
	$M_J$	0.08	0.12	0.17	0.26	0.38
$l_{\max }$ $=14$	$ℳ$	7	10	15	22	32
	$t_{\mathrm{cr}}$	0.07	0.10	0.15	0.23	0.33
	$Z$	1	1	1	1	1

Table 2. Table 2: A collection of results.

Cloud name	$m_{\mathrm{gas}}\left(M_{\odot }\right)$	${\overline{n}}_{\mathrm{gas}}\left({\mathrm{cm}}^{-3}\right)$	$\Sigma \left(M_{\odot }{\mathrm{pc}}^{-2}\right)$	$Z\left(Z_{\odot }\right)$	$m_{\mathrm{cl}}\left(M_{\odot }\right)$ a	TSFE (%) b	SFRff c	IMF slope d	$n_{\mathrm{SN}}$ e
XS-F	$3.2\text{×}{10}^{03}$	$1.8\text{×}{10}^{02}$	$41$	1	$3.8\text{×}{10}^{02}$	$12.1$	$0.18$	${1.0}_{-0.3}^{+0.4}$	$2$
S-F	$1.0\text{×}{10}^{04}$	$1.8\text{×}{10}^{02}$	$61$	1	$5.1\text{×}{10}^{02}$	$5.1$	$0.062$	${1.3}_{-0.3}^{+0.3}$	$2$
M-F	$3.2\text{×}{10}^{04}$	$1.8\text{×}{10}^{02}$	$89$	1	$1.4\text{×}{10}^{03}$	$4.3$	$0.042$	${1.1}_{-0.2}^{+0.2}$	$12$
L-F	$1.0\text{×}{10}^{05}$	$1.8\text{×}{10}^{02}$	$131$	1	$5.7\text{×}{10}^{03}$	$5.7$	$0.053$	${1.2}_{-0.1}^{+0.2}$	$38$
XL-F	$3.2\text{×}{10}^{05}$	$1.8\text{×}{10}^{02}$	$193$	1	$2.5\text{×}{10}^{04}$	$7.8$	$0.043$	${1.1}_{-0.1}^{+0.1}$	$142$
XS-C	$3.2\text{×}{10}^{03}$	$1.8\text{×}{10}^{03}$	$193$	1	$1.0\text{×}{10}^{02}$	$3.3$	$0.033$	${0.5}_{-0.0}^{+0.8}$	$0$
S-C	$1.0\text{×}{10}^{04}$	$1.8\text{×}{10}^{03}$	$283$	1	$5.3\text{×}{10}^{02}$	$5.3$	$0.052$	${1.6}_{-0.3}^{+0.1}$	$1$
M-C	$3.2\text{×}{10}^{04}$	$1.8\text{×}{10}^{03}$	$415$	1	$3.0\text{×}{10}^{03}$	$9.4$	$0.047$	${1.2}_{-0.2}^{+0.2}$	$5$
L-C	$1.0\text{×}{10}^{05}$	$1.8\text{×}{10}^{03}$	$609$	1	$1.4\text{×}{10}^{04}$	$13.7$	$0.099$	${1.2}_{-0.1}^{+0.1}$	$47$
L-C-lm	$1.0\text{×}{10}^{05}$	$1.8\text{×}{10}^{03}$	$609$	1/10	$3.4\text{×}{10}^{03}$	$3.4$	$0.021$	${1.2}_{-0.2}^{+0.2}$	$5$
L-C-xlm	$1.0\text{×}{10}^{05}$	$1.8\text{×}{10}^{03}$	$609$	1/40	$3.3\text{×}{10}^{03}$	$3.3$	$0.025$	${1.0}_{-0.2}^{+0.2}$	$5$
XXS-VC	$1.0\text{×}{10}^{03}$	$1.8\text{×}{10}^{04}$	$609$	1	$9.8\text{×}{10}^{01}$	$9.8$	$0.099$	${0.5}_{-0.0}^{+0.9}$	$0$
XS-VC	$3.2\text{×}{10}^{03}$	$1.8\text{×}{10}^{04}$	$894$	1	$5.1\text{×}{10}^{02}$	$16.1$	$0.2$	${1.0}_{-0.2}^{+0.3}$	$0$
S-VC	$1.0\text{×}{10}^{04}$	$1.8\text{×}{10}^{04}$	$1312$	1	$3.2\text{×}{10}^{03}$	$32.2$	$0.31$	${1.5}_{-0.2}^{+0.2}$	$0$
M-VC	$3.2\text{×}{10}^{04}$	$1.8\text{×}{10}^{04}$	$1925$	1	$1.5\text{×}{10}^{04}$	$46.6$	$0.25$	${1.4}_{-0.1}^{+0.1}$	$0$
L-VC	$1.0\text{×}{10}^{05}$	$1.8\text{×}{10}^{04}$	$2827$	1	$2.7\text{×}{10}^{04}$	$27.4$		${1.3}_{-0.1}^{+0.1}$	$0$

Equations

$${\alpha }_{\mathrm{vir}}=\frac{5{\sigma }_{3D}^{2}R}{3GM}\approx 0.4,$$

$${\lambda }_{\mathrm{J}}=c_{\mathrm{s}}\sqrt{\frac{\pi }{G\rho }}>N_{\mathrm{sink}}\Delta x.$$

$${\rho }_{\mathrm{J}}=\frac{\pi c_{\mathrm{s}}^2}{G{N}^2\Delta x^2}.$$

$$M_{\mathrm{J}}=\frac{4\pi }{3}{\rho }_{\mathrm{J}}{\left(\frac{{\lambda }_{\mathrm{J}}}{2}\right)}^3=0.55M_{\odot }\left(\frac{\Delta x_{\min }}{1000\mathrm{AU}}\right),$$

$$q(m_i)=V(0.3m_i)\left(\frac{S_{\mathrm{cl}}({\Sigma }_i{m}_i)}{{\Sigma }_{i}V(0.3m_i)}\right),$$

$$\ln \mathcal{L}$$

$$=N_0\left[\ln \Gamma -\ln (m_{\min }^{-\Gamma }-m_{\max }^{-\Gamma })\right]-\ln 10\Gamma \sum _{\mathrm{i}}{\mu }_{\mathrm{i}}.$$

$$M_{\max }/{\mathrm{M}}_{\odot }\approx 205m_4^{0.66},$$

$$m_{\mathrm{cl}}=200\text{leavevmode}{\mathrm{M}}_{\odot }\cdot {\left(\frac{m_{\mathrm{gas}}}{10^4\text{leavevmode}{\mathrm{M}}_{\odot }}\right)}^{1.4}{\left(1+\frac{{\overline{n}}_{\mathrm{gas}}}{n_{\mathrm{cri}}}\right)}^{0.91}+m_{\mathrm{fl}},$$

$$f_{\ast ,\mathrm{tot}}=2.0\%{\left(\frac{m_{\mathrm{gas}}}{10^4\text{leavevmode}{M}_{\odot }}\right)}^{0.4}{\left(1+\frac{{\overline{n}}_{\mathrm{gas}}}{n_{\mathrm{cri}}}\right)}^{0.91}.$$

$$f_{\mathrm{F}}(\tau )=\frac{f_0}{e^{-(\tau -{\tau }_0)/\Delta \tau }+1}.$$

$${\mathrm{SFR}}_{\mathrm{ff}}{|}_{\max }\approx 1.1\%{\left(\frac{m_{\mathrm{gas}}}{10^4{M}_{\odot }}\right)}^{0.36}{\left(1+\frac{{\overline{n}}_{\mathrm{gas}}}{n_{\mathrm{cri}}}\right)}^{{\alpha }_{\mathrm{f}}}$$

$$f_{\ast ,\mathrm{tot}}\approx f_0=\frac{\mathrm{d}{f}_{\mathrm{F}}}{\mathrm{d}\tau }{|}_{\max }\times \Delta {\tau }_{\mathrm{SF}}.$$

$$\frac{\mathrm{d}{\rho }_{\ast }}{\mathrm{d}t}=ϵ{\left(\frac{m_{\mathrm{gas}}}{10^4{M}_{\odot }}\right)}^{0.36}{\left(1+\frac{{\overline{\rho }}_{\mathrm{gas}}}{{\rho }_{\mathrm{cri}}}\right)}^{1.0}\frac{{\overline{\rho }}_{\mathrm{gas}}}{t_{\mathrm{ff}}}\propto ({\overline{\rho }}_{\mathrm{gas}}{)}^{2.5},$$

$$\text{if }{\overline{n}}_{\mathrm{gas}}>n_{\mathrm{cri}}\approx 10^3\text{leavevmode}{\mathrm{cm}}^{-3}$$

$$N_{\mathrm{ref}}\sqrt{f_{\mathrm{c}}}<N_{\mathrm{sink}}<2N_{\mathrm{ref}}\sqrt{f_{\mathrm{c}}},$$

Full text

††pubyear: 2019††pagerange: Simulating Star Clusters Across Cosmic Time: I. Initial Mass Function, Star Formation Rates and Efficiencies–4

1 Introduction

Stars formation in galaxies is a complex and only partially understood astrophysical phenomenon. It is difficult to formulate a general theory in part because of the wide range of scales and of physical processes involved. From an observational point of view, quantifying star formation efficiency (SFE) in nearby molecular clouds has been the focus of much recent research (e.g., Lada et al., 2010; Heiderman et al., 2010; Gutermuth et al., 2011). A power-law relationship between the gas surface density of galaxies and their star formation rate (SFR) was first proposed by Schmidt (1959) and later tested by large, multi-galaxy data (Kennicutt, 1998). This relationship has been widely used in cosmological simulations of galaxy formation. However, on sub-galactic scales the dispersion of star formation rates for a given gas surface density of H i is large, and other parameters such as gas metallicity (Bolatto et al., 2011; Krumholz, 2013) and stellar surface density (Leroy et al., 2008) appear to become important. As the resolution of surveys improved, numerous studies have shown that star formation on kpc scales is more strongly correlated with H2 surface density (e.g. Krumholz, 2014) rather than atomic gas. Therefore, modern cosmological simulations of galaxy formation aim at reproducing the molecular phase of the interstellar medium (ISM) and adopt an empirical sub-grid recipe for star formation within partially resolved molecular clouds of the form $\dot{\rho}_{\mathrm{*}}\propto\rho_{\mathrm{H_{\mathrm{2}}}}^{n}$, where typically $n=1$ or $1.5$. The molecular phase of the ISM is treated in simulations using different prescriptions: Robertson & Kravtsov (2008) pre-computed a grid of models from a photo-chemistry code, Gnedin et al. (2009) directly solved the formation and dissociation equations for H2, but with and increased formation rate to model unresolved clumping, and Kuhlen et al. (2012) used an analytic model to estimate the equilibrium H2 abundance. The sub-grid recipe (with grid maximum resolution typically between few parsecs to few kpc) is calibrated to reproduce observational data in galaxies at $z=0$.

However, the conditions in the ISM of high-redshift galaxies are likely different to those found in the present day. Krumholz et al. (2012) argue that the SFR is in fact correlated to the local free-fall time set by the gas density, not to the column density. Simulations show that densities and pressures of star-forming regions in high-redshift galaxies are much higher than in today’s ISM (e.g., Ricotti, 2002; Wise et al., 2014; Ricotti, 2016). Using Adaptive Mesh Refinement (AMR) simulations of the first stars and galaxies with parsec-resolution, Ricotti et al. (2016) found that compact molecular clouds in primordial galaxies can either form gravitationally bound star clusters that resemble the progenitors of today’s globular clusters (GC)111The compact bound stellar objects found in the simulations are actually not only globular clusters progenitors, but also ultra-compact dwarfs and dwarf-globular transition objects, depending on whether the stellar cluster forms at the centre of the halo, in the disk’s spiral arms, or even in satellite minihalos., or the clusters may disperse and fill up a large fraction of the dark matter halo of primordial dwarf galaxies. In this second case the stars would appear as spheroids 20-200 pc in radius, dark matter dominated and with very low surface brightness. These objects would be identified today as “ultra-faint” dwarf galaxies observed in the Local Group (e.g., Willman et al., 2005; Zucker et al., 2006a, b; Belokurov et al., 2007; Walsh et al., 2007; Majewski et al., 2007; Martin et al., 2009). Star formation in compact star clusters appears to be especially important, even perhaps the dominant mode of star formation at high-redshift. Thus, in order to make progress in understanding the formation of the first dwarf galaxies and the sources of reionisation, it is important to focus on understanding the small-scale physics of this process, which is poorly resolved in cosmological simulations.

Most numerical work on star formation in molecular clouds focuses on star formation in the local universe, aiming at explaining observed young star forming regions. In this paper we analyse the results of a large grid of simulations of realistic molecular clouds with initial conditions chosen to reproduce not only local molecular clouds but also clouds that form in higher density and pressure environments, typical of star formation in high redshift galaxies. We vary the masses of the clouds, their compactness (central density), and in few cases explore the effect of changing the gas metallicity and therefore the gas cooling function.

The motivation for this paper is twofold. The first goal is to deepen our understanding of the physics of star formation in high-pressure environments to justify and inform the sub-grid star formation recipe used in cosmological simulations. A closely related important question in Near Field Cosmology is: how does the formation of self-gravitating bound star-clusters relate to the star formation efficiency, compactness, mass and gas metallicity of molecular clouds found in cosmological simulations? We will only touch on this questions in the present paper, but more detailed work will be presented in a followup paper.

The second goal is to estimate the escape fraction of H i ionising radiation from molecular clouds as a function of cloud compactness and mass. This is the first necessary step for a realistic estimate of the escape fraction from galaxies. Ricotti (2002) have shown that, if a non-negligible fraction of today’s GCs formed at $z>6$ with $\langle f_{\rm esc}\rangle$ $\sim 1$, their progenitors would be a dominant source of ionising radiation during reionisation. Katz & Ricotti (2014) presented arguments in support of significant fraction of today’s old GCs forming before the epoch of reionisation. However, although it is naively expected, it has not been shown with numerical simulations that $\langle f_{\rm esc}\rangle$ from GC progenitors forming in compact molecular clouds is higher than $\langle f_{\rm esc}\rangle$ in more diffuse clouds. The answer to this question and the contribution of compact star clusters to reionisation will be presented in a separate companion paper.

This paper is organised as follows. In Section 1.1 we present a brief review of the current status of numerical simulations of star cluster formation. In Section 2 we provide an overview of our numerical methods, including details on the initial conditions of our simulations and the recipes for formation of sink particles and feedback. In Section 3 we present some results from the analysis of our large set of simulations with emphasis on the stellar initial mass function (IMF) and the star formation rate (SFR) and efficiency (SFE). A summary and conclusions are presented in Section 4.

1.1 The IMF and SFE of Molecular Clouds

Simulations of molecular cloud dynamics are valuable tools in understanding the conditions in the ISM. Typically these simulations adopt idealised initial conditions similar to those in observed clouds: a gas cloud $\sim 1-10$ pc in size supported against gravity by a turbulent velocity field such that the initial virial ratio, i.e. the ratio of the kinetic energy to the potential energy of the cloud, is $\lesssim 0.5$. One model involves injecting turbulence into a volume of gas in the initial conditions and allowing it to decay over time. This can be done by either using smoothed particle hydrodynamics (SPH, e.g. Klessen, 2001; Bonnell et al., 2006) or grid-based methods (e.g., Gammie et al., 2003). Another model involves adding turbulence continuously over time, simulating the effect of momentum injection from outside flows or energy from massive stars inside the cloud (Vazquez-Semadeni et al., 1997; Ballesteros-Paredes et al., 2006; Padoan et al., 2007). Many of these models adopted an isothermal equation of state, while others have included self-consistent cooling and heating functions (e.g. Koyama & Inutsuka, 2004; Audit & Hennebelle, 2005) and molecular chemistry (e.g. Glover et al., 2010).

The fragmentation of molecular clouds into stars is a long-standing problem. Observational studies (Salpeter, 1955; Kroupa, 2002; Chabrier, 2005) have found that the masses of stars follow a “Initial Mass Function” (IMF) with a power law $\mathop{}\!\mathrm{d}N/\mathop{}\!\mathrm{d}\log M\propto M^{-\Gamma}$ at the high-mass end (Salpeter (1955) calculate $\Gamma\approx 1.35$). Various theoretical models have been constructed to explain this (Padoan & Nordlund, 2002; Mac Low & Klessen, 2004; Hennebelle & Chabrier, 2008; Hopkins, 2012) based on gravoturbulent fragmentation of the host cloud. Radiative stellar feedback has been invoked to explain the precise shape of the IMF, using both simulations (e.g., Bate, 2009) and analytic models (e.g., Guszejnov & Hopkins, 2016). Early pioneering simulations of cluster formation approached the problem of producing a well-defined IMF (e.g., Bate et al., 2003; Bate & Bonnell, 2005; Klessen et al., 2008; Offner et al., 2008), but were often limited in terms of statistics or resolution. More recent work, with increasing computing power, provided more reliable statistics and IMF distributions (e.g., Bonnell et al., 2003; Bate, 2009; Bonnell et al., 2011; Girichidis et al., 2011; Krumholz et al., 2011; Bate, 2012; Ballesteros-Paredes et al., 2015).

Simulations attempting to capture the stellar IMF require a high dynamic range to resolve both brown dwarfs and OB stars. Most recently, Bate (2019) resolve in detail the mass spectrum of brown dwarfs while only producing stars of up to 3 ${\rm M}_{\odot}$, finding that low metallicities do not produce observable differences in the stellar IMF, while increasing fragmentation. Gavagnin et al. (2017) have lower mass resolution but capture more massive stars that emit significant quantities of ionising radiation, arguing that this alters the high mass end of the IMF. In the absence of radiation and cooling, Lee & Hennebelle (2018a) and Lee & Hennebelle (2018b) study the early formation of protostellar Larson cores, and find that the choice of equation of state (eos) has a strong influence on the peak of the IMF. In general, these works are relatively successful at reproducing not only the IMF but also stellar multiplicity and separation.

Previous authors have included ideal MHD in their simulations (Myers et al., 2013; Krumholz et al., 2016; Cunningham et al., 2018). However, since these authors only form stars up to $\sim$20 ${\rm M}_{\odot}$, they neglect ionising radiation. Non-ideal MHD effects, while challenging to include in resolving the IMF for reasons of computational cost, appear to affect the dynamics of protostar formation on small scales (Masson et al., 2016; Vaytet et al., 2018). The physics that shapes the IMF is complex, and a full treatment that covers non-ideal MHD, both low and high energy radiation, chemistry and the full mass range of stars remains difficult with modern computational resources.

As well as the stellar IMF, an important consideration is how many stars are formed out of a given mass of gas, or the Star Formation Efficiency (SFE). The efficiency of conversion of gas into stars is typically much lower than 100% since energetic processes from massive stars are able to disperse the cloud in which a star cluster forms before all of the gas collapses into protostars. These processes are widely termed “feedback” (see review by Dale et al., 2015). Recent work favours ionising radiation as the main driver of molecular cloud dispersal (Dale et al., 2005; Gritschneder et al., 2009; Peters et al., 2010; Walch et al., 2012; Dale et al., 2012), as opposed to other effects such as stellar winds (Dale et al., 2014), although Howard et al. (2016) find that UV photoionisation has little effect on the initial evolution of the SFE.

The relationship between gas properties and the SFE is a matter of ongoing study. Lada et al. (2010) and Heiderman et al. (2010) argue for a constant ratio between gas above a certain surface or column density, although this is still subject to discussion (Gutermuth et al., 2011; Hony et al., 2015). There is no clear theoretical link between the SFE and projected column density, although Clark & Glover (2014) argue that there may be a link between the observed column density and the local density around the star. Geen et al. (2017) reproduce the SFE observed by Lada et al. (2010), although they find that the result is likely to be dependent on the average density of the neutral gas in the cloud. These simulations produce similar results to the simulations of Colin et al. (2013). Geen et al. (2018) finds that SFE can change by up to a factor of 4 by varying the initial velocity field of the cloud and the stellar IMF, although relationships can be found between the early cloud state and the final SFE. Semi-analytic models by Vázquez-Semadeni et al. (2018) also find considerable scatter in the SFE.

2 Numerical Simulations and Methods

We conduct our numerical simulations using the AMR radiative magneto-hydrodynamical code ramses (Teyssier, 2002; Bleuler & Teyssier, 2014). Radiative transfer is implemented using a first-order moment method with M1 closure described in Rosdahl et al. (2013).

Kim et al. (2017) demonstrates that M1 closure method is inaccurate near sources only in regions where the flux is about an order of magnitude smaller than the mean value (due to shielding), while it agrees with adaptive ray-tracing methods (e.g., Wise et al., 2014; Hartley & Ricotti, 2016) both at larger distances from individual sources and on global scales. M1 closure, however, is significantly more computationally efficient than ray-tracing methods.

The ionising photons interact with neutral gas and we track the ionisation state and cooling/heating processed of hydrogen and helium (see Geen et al., 2017, for details). Our simulations include magnetic fields in the initial conditions, but we do not include the chemistry of molecular species (i.e., formation/dissociation). 3-D ‘zoom-in’ simulations of the chemical evolution of molecular clouds suggest that, for gas at solar metallicity, the cloud is almost fully molecular with H2 fractions around 0.9 in the later stages of transition to dense molecular phase (Seifried et al., 2017).

We simulate a set of isolated and turbulent molecular clouds that collapse due to their own gravity. We explore a grid of simulations varying the initial gas mass and compactness (i.e., the core density) of the clouds. In our simulations, dense proto-stellar cores collapsing below the resolution limit of the simulations produce sink particles. These sinks may represent single stars or multiple stars or even clusters of stars if the resolution is not sufficiently high. However, in all our simulations we aim at reproducing a realistic high-mass end of the stellar IMF and therefore realistic feedback from individual massive stars. To accomplish this goal, sink particles emit hydrogen and helium ionising photons according to their mass as described in § 2.3. The gas is ionised and heated by massive stars, producing over-pressurised bubbles that blow out the gas they encounter. In our simulations low mass stars and proto-stellar cores do not produce any feedback. In this work we do not include mechanical feedback from supernova (SN) explosions and from stellar winds and we also neglect the effect of radiation pressure from infrared radiation. However, with the exception of the two most massive clouds in the set of simulations representing today’s molecular clouds (the lowest density set), all the simulations stop forming stars before the explosion of the first SN. Therefore, neglecting SN feedback is well justified in these cases.

We find that in all simulations star formation has ceased before $\sim 5-6\leavevmode\nobreak\ t_{\rm ff}$, which is the typical time it takes for feedback to act. We stop simulations no earlier than this point. A few simulations are continued beyond this time. This does not have an affect on the IMF since the mass function of sink particles does not change after star formation has ceased.

For the simulations in the set in which SN explosions should occur while star formation is ongoing, in order to compensate for this missing feedback, we do not shut down UV radiation feedback from massive stars after the time the star should have exploded as SN. In the following sections we provide some more details on the simulations set up.

2.1 Initial Conditions

We run a grid of 14 simulations of clouds with a range of central densities and initial gas masses. We also run some additional simulations varying the initial gas metallicity and therefore the gas cooling function. The magnetic field strength in the initial conditions is set such that $v_{\mathrm{a}}=0.2\leavevmode\nobreak\ \sigma_{\mathrm{3D}}$, where $v_{\mathrm{a}}$ is Alfven wave velocity and $\sigma_{\mathrm{3D}}$ is the turbulence velocity dispersion. This $\nu_{\mathrm{a}}$ is $\sim 2$ times smaller than that measured in a group of molecular clouds by Crutcher (2012) who finds $\nu_{\mathrm{a}}\approx 0.5\sigma_{\mathrm{3D}}$.

The clouds have initially a spherically symmetric structure with density profile of a non-singular isothermal sphere with core density $n_{\mathrm{\rm c}}$. The cloud extends out to $r_{\mathrm{\rm gas}}=3r_{\mathrm{\rm c}}$, where $r_{\mathrm{\rm c}}$ is the core radius. Beyond $r_{\mathrm{\rm gas}}$ the cloud is embedded in a uniform density envelope that extends to $6r_{\mathrm{\rm c}}$ with a density $0.01\leavevmode\nobreak\ n_{\mathrm{\rm c}}$. Outside of the envelope the number density is constant at $1$ cm*-3*. The box length $L_{\mathrm{box}}$ is set to $48r_{\mathrm{\rm c}}$ in each simulation. The initial value of the (isothermal) sound speed of the cloud is set to $c_{\mathrm{\rm s}}=0.24$ km/s, while the envelope and background densities are in pressure equilibrium.

The initial density profile is perturbed with a turbulent velocity field, analogously to the set up used in Geen et al. (2017). The initial turbulence of the clouds follows a Kolmogorov power spectrum with random phases and has an amplitude such that the cloud is approximately in virial equilibrium. All simulations have the same set of random phases. The initial cloud virial ratio

[TABLE]

is kept constant in all the simulations. Therefore the ratio $t_{\mathrm{ff}}/t_{\mathrm{turb}}$, where $t_{\mathrm{turb}}\equiv R/\sigma_{\mathrm{3D}}$ is kept constant in all the simulations. However, the sound crossing time $t_{\mathrm{cr}}\equiv R/c_{\mathrm{s}}$, where troughout this paper we assume $c_{\mathrm{s}}=10$ km/s, is not constant. The virial parameter, $\alpha_{\mathrm{\rm vir}}$, is small enough to ensure collapse and fragmentation, but sufficiently large to prevent a rapid radial collapse of the cloud. Before allowing any star formation in the cloud we evolve these idealised initial conditions for $\sim 3t_{\mathrm{\rm ff}}$, so that the turbulent velocities develop into density perturbations and the initial conditions relax into a quasi-equilibrium turbulent medium. If we do not allow the initial conditions to relax before forming stars, the stars form mostly near the centre of the cloud during the transient relaxation phase.

A detailed list of the parameters in our simulations is shown in Table 1. The clouds in Table 1 are labelled with letters of two or three parts. The first part is either ‘XXS’ (extra-extra-small), ‘XS’ (extra-small), ‘S’ (small), ‘M’ (medium), ‘L’ (large), or ‘XL’ (extra-large), representing various initial gas masses of ${10}^{3}$, $3.16\text{\times}{10}^{3}$, ${10}^{4}$, $3.15\text{\times}{10}^{4}$, ${10}^{5}$, and $3.16\text{\times}{10}^{5}$ M*⊙, respectively. The second part is either ‘F’ (fiducial, which are the most similar to clouds in the solar neighbourhood), ‘C’(compact), or ‘VC’ (very compact), in order of increasing initial mean gas density. The mean particle number density of the cloud, $\overline{n}_{\mathrm{gas}}=\overline{\rho}_{\mathrm{gas}}/(\mu m_{\mathrm{p}})$, where $\mu=1.4$ is the mean molecular weight of the atomic gas, increases by a factor of ten between each set of simulations from $1.8\text{\times}{10}^{2}$ cm-3* to $1.8\text{\times}{10}^{4}$ cm*-3*. The ‘L-C’ setup has two more simulations that have a third part in the name, ‘lm’ and ‘xlm’, representing ‘low-metallicity’ and ‘extra-low-metallicity’. A comparison of our setups with the literature is shown in Figure 1.

2.2 Resolution and Sink Formation

We use a Cartesian grid with an octree structure with cells that we subdivide into $2^{3}$ child cells as the simulation evolves (“adaptive refinement”). Our starting refinement level is $\ell_{\mathrm{\rm min}}=7$ (corresponding to $\Delta x=L_{\mathrm{\rm box}}/2^{7})$ and maximum level of refinement is $\ell_{\mathrm{\rm max}}=15$ for runs with the lowest mean density and $\ell_{\mathrm{\rm max}}=14$ for all the other runs. The resolution is therefore $\Delta x_{\mathrm{\rm min}}=L_{\mathrm{\rm box}}/2^{15}$ for the "fiducial" runs, which corresponds to resolutions between $500$ AU and $2300$ AU. The "compact" clouds have resolution between $460$ AU and $1500$ AU and the "very compact" clouds between $150$ AU and $680$ AU.

In order to resolve the Jeans length with $N$ grid cells it is required that

[TABLE]

From Equation (2), the Jeans length is resolved with at least $N_{\mathrm{sink}}$ grid points if $\rho<\rho_{\mathrm{J}}$, where

[TABLE]

In our simulations we enforce the refinement criterion that the Jeans length is resolved with at least $N_{\mathrm{\rm ref}}=10$ cells. Hence, when the local density goes up and reaches a point where $\lambda_{\mathrm{\rm J}}$ becomes smaller than $N_{\mathrm{\rm ref}}\Delta x$, each cell is refined individually into eight new children cells. This refinement condition is always true, up to the maximum refinement level (when $\Delta x=\Delta x_{\mathrm{\rm min}}$). When the gas density exceeds $\rho_{\mathrm{\rm J}}^{\rm max}=\rho_{\mathrm{\rm J}}(\Delta x=\Delta x_{\mathrm{\rm min}},N)$ at the maximum refinement level, we cannot continue to resolve the Jeans length with at least $N_{\mathrm{\rm ref}}$ cells. We therefore create sink particles to trace material above these densities. We set $\rho_{\mathrm{\rm sink}}=\rho_{\mathrm{\rm J}}(\Delta x=\Delta x_{\mathrm{\rm min}},N=N_{\mathrm{\rm sink}})$ as critical density threshold to form sink particles. Sink particles are created on the fly using a peak detection algorithm (see Bleuler & Teyssier, 2014, for details on sink particle formation in RAMSES). We first detect density clumps above a density threshold $f_{\mathrm{\rm c}}\rho_{\mathrm{\rm sink}}$, with $f_{\mathrm{\rm c}}=0.1$. Then, the algorithm performs a peak density check, a collapsing check ($\nabla\cdot v=0$), and virial check before forming a sink particle.

In order to avoid numerical fragmentation it is usually suggested that $N_{\mathrm{\rm sink}}\geq 4$ (Truelove et al., 1997). In our simulations we adopt $N_{\mathrm{\rm sink}}=5$ for reasons detailed Appendix A. With an initial sound speed $c_{\mathrm{\rm s}}=0.24$ km/s the Jeans mass at the sink density threshold is

[TABLE]

which results in $M_{\mathrm{\rm J}}\sim 0.08M_{\mathrm{\odot}}$–$0.8M_{\mathrm{\odot}}$ for the compact and very compact clouds and $\sim 0.3M_{\mathrm{\odot}}$–$1.3M_{\mathrm{\odot}}$ for the fiducial clouds.

The sink particles are then treated like point masses and accrete gas based on the mechanism described as ‘threshold accretion’ in Bleuler & Teyssier (2014). The dynamics of the sink particles takes into account gravitational force from gas and stars and it is evolved using a leap-frog integration scheme. The effect of gas dynamical friction is not included.

2.3 Feedback and Properties of UV source

In our simulations, ionising UV photons are emitted from sink particles from the time they form to the end of the simulation. Massive stars have lifetime of few Myrs, shorter than the duration of some of our simulations, and they may explode as SNe during the simulation. Since we are not implementing SNe feedback, we keep the stars emitting radiation after their death to compensate for the lack of SNe in the attempt of avoiding underestimating feedback effects. While SN explosions produce a significant amount of mechanical energy (typically $10^{51}$ egs), the energy associated with ionising radiation from massive stars integrated through their main-sequence lifetime is comparable (or larger for more massive stars) and this feedback starts acting earlier than SN feedback. For an O-star, more than half of the radiation is emitted in hydrogen ionising photons. Typically $\sim 10\%$ of a star’s hydrogen is burned in the nuclear fusion process, with an energy efficiency of $\sim 0.7\%$. Thus the amount of energy radiated by a massive star during its life time is $\sim 2\times 10^{-3}M_{\mathrm{*}}$, or $\sim 4\times 10^{52}$ ergs for a 20 M*⊙* star.

For each simulation we estimate the total hydrogen-ionising photon emission rate at a given time as $S_{\mathrm{\rm cl}}(m_{\mathrm{cl}})=8.96\times 10^{46}\,\textrm{s}^{-1}\,(m_{\mathrm{cl}}/M_{\mathrm{\odot}})$ (see Geen et al., 2017), where $m_{\mathrm{cl}}$ is the total mass of the sink particles. This is calculated by Monte Carlo sampling a stellar population as described in Geen et al. (2016) (See Sec. B). The fraction of the total hydrogen ionising photon emission rate attributed to each sink particle is based on the following relation

[TABLE]

where $V(m)$ is the hydrogen-ionising photon emission rate from a star with mass $m$, using the fits from Vacca et al. (1996).

The factor $0.3$ is an empirical factor to account for the scaling between the masses of sink particles and those of massive stars, necessary because we do not fully resolve the fragmentation of sink particles into proto-stars. See Section 3.1 for further discussion. The correction factor $X\equiv S_{\mathrm{\rm cl}}(\Sigma_{\mathrm{i}}m_{\mathrm{i}})/\Sigma_{\mathrm{i}}V(0.3m_{\mathrm{i}})$ is very close to unity in most simulation in which we resolve massive stars and it is introduced only to prevent overproducing ionising radiation in case massive stars are poorly resolved. In all simulations we also impose $X\leq 1$.

For the fiducial clouds we also include He0 and He+ ionising photons with total emission rates being $S_{\mathrm{\rm He^{0}}}=(1.178\times 10^{46}\,\textrm{s}^{-1})(M_{\mathrm{*}}/M_{\mathrm{\odot}})$ and $S_{\mathrm{\rm He^{+}}}=(2.422\times 10^{43}\,\textrm{s}^{-1})(M_{\mathrm{*}}/M_{\mathrm{\odot}})$. These rates are calculated using the same method as for hydrogen ionising photons described above using a Kroupa IMF (Kroupa, 2002). For the luminosity of individual stars we use Schaerer (2002) fitting for $Q$(He0) and $Q$(He+) with extrapolations above 150 $M_{\mathrm{\odot}}$. This is contrasted by the model of Gavagnin et al. (2017), who assume blackbody spectra for each star. See Appendix B for details.

Various authors have concluded that UV photoionisation is typically the most important process in regulating star formation on a cloud scale. Haworth et al. (2015) find that additional processes beyond hydrogen photoionisation have a correcting factor of 10% at best. Radiation pressure mainly becomes important at very high surface densities, which principally affects smaller scales than the ones studied here - see Crocker et al. (2018) for idealised conditions and Kim et al. (2018) for simulations with self-consistent star formation feedback. Dale et al. (2014) further find that winds have a minimal effect on the star formation efficiency of molecular clouds. We are thus justified in our choice to focus on UV photoionisation feedback in this work, but discuss cases where this may not be sufficient later in the paper.

2.4 Cooling

We use the radiative cooling function described in Geen et al. (2016). The cooling in neutral gas is based on the prescription in Audit & Hennebelle (2005), which includes cooling from carbon, oxygen and dust grains as well as the effect of the ambient UV background in the ISM. For collisionally ionized gas at temperatures $>10^{4}$ K we use Sutherland & Dopita (1993) cooling function. The out-of-equilibrium cooling of photoionised hydrogen and helium is treated as described in Rosdahl et al. (2013). Out-of-equilibrium cooling of photoionised metals is treated with a piecewise fit to the cooling curve given in Ferland (2003). We assume a uniform metallicity as listed in Table 1. For most simulations, this is solar metallicity, though we perform some simulations at sub-solar metallicity. We do not implement out-of-equilibrium molecular chemistry.

3 Results

In this section we present and discuss the results of our simulations. In § 3.1, we study the mass function of cores (sink particles) and the IMF. In § 3.2 we focus on the star formation efficiency and in § 3.3, on the star formation rate.

A representative sample of snapshots from our simulation set is shown in Figures 2 and 3. These figures show the time evolution of the density-weighted projections of the gas density and temperature as well as the position of the sink particles along line of sight for three simulations with different mean clund densities (fiducial on the left column, compact, middle column, and very-compact on the right column), and a cloud mass of $3.2\times 10^{4}$ M*⊙*.

We observe clearly that the star formation efficiency increases with increasing cloud density, and the stellar cluster that is formed at the end of the simulation remains more compact and self-gravitating for the densest cloud. The effect of radiative feedback from massive stars is also clearly visible in the density and temperature projections. H ii regions break out of the dense filaments destroying them and reducing the overall mass in dense gas in which stars are formed.

3.1 Stellar Initial Mass Function

3.1.1 Cores Fragmentation and Initial Mass Function

Maps in the continuum of cluster regions and larger areas in star-forming systems allow to construct a core mass function (CMF), that is the mass function of high-density gas concentrations (starless cores) with mass sufficiently large to be identified given the resolution of the observations (e.g., Motte et al., 1998). Observations of a well-resolved CMF in the Pipe nebula show a striking similarity to the stellar IMF, but shifted to higher masses by a factor of a few, which suggests that the IMF is the direct product of the CMF with a roughly constant core-to-star conversion efficiency $\sim 30\%$ (e.g., Matzner & McKee, 2000; Alves et al., 2007)

Previous works on star formation in molecular clouds which adopted sink particles (like in the present paper) have investigated the mapping between the masses of pre-stellar cores at the time they become self-gravitating and the final masses of the stars that form within them (e.g., Padoan et al., 2001; Smith et al., 2009). For instance, Smith et al. (2009), using SPH simulations, find that at early times the relationship between stellar masses and the parent cores can be reproduced within a modest statistical dispersion with the star being about one-third of the parent core mass.

We find results in agreement with these previous studies. Figure 4 (and Figure 7) show the stellar mass function for our grid of simulations obtained assuming that the sink particles are about a factor of 2.5 more massive than the corresponding massive stars they produce. This means that the number of stars of a given mass ($>1$ M*⊙*) is given by a Kroupa IMF for a star cluster with total mass equal to the total mass of the sinks. In Figure 7 the IMF is shown at the end of the simulation when star formation has stopped, while in Figure 4 we also show the time evolution of the IMF, together with the Kroupa IMF for a cluster with total mass equal to the total mass in sink particles (dashed lines). This means that we are assuming nearly 100% efficiency of star formation in the cores (i.e., the cores fragment into stars) or a lower efficiency but the gas expelled by feedback is later transformed into low-mass stars. However, we think that this second model is less physically motivated. We can see that the shifted sink mass function (SMF) matches the Kroupa IMF at the high-mass end, both in terms of the slope and normalisation at any time during the formation of the star cluster. The figure suggests that the birth of stars in a cluster follows the same random sampling of the universal mass function throughout the star-formation process. In other words, it appears that there is not a bias toward formation of high mass-stars or low-mass stars during the early times when the cluster is in the formation process.

As discussed above the sink particles can be interpreted as pre-stellar cores, and each sink particle converts $\sim 40\%$ of its mass to a single massive star, with the rest of the mass ending up in low-mass stars, filling up the lower mass end of the IMF. The flattening or cut-off of the IMF at the low-mass end observed in our simulations is likely due to insufficient spatial and mass resolution to capture the formation of low mass cores or the fragmentation of more massive pre-stellar cores. To further clarify, we note that the interpretation that only 40% of the core mass is converted into a star and the remaining 60% is returned to the gas phase, never to participate in star formation, would not produce the correct normalisation of the IMF. This is because in this scenario, although the stellar masses are $\sim 40\%$ of the core masses, the total stellar mass and the star formation efficiency would be reduced by a factor of $\sim 2.5$, lowering the expected number of massive stars below the value found in the simulations.

So far we have assumed that all the gas in the cores fragments into stars with $\eta=100\%$ efficiency. In this case we find a conversion factor $\varepsilon=0.4$ between the CMF and the IMF (such that the normalisation of the IMF agrees with the observed one). However, it is possible to match the observed IMF also in models in which $\eta<1$. Note that in this case the TSFE shown in all our plots should be re-scaled by a factor $\eta$. In models with $\eta<1$, the conversion factor $\varepsilon$ that matches the mass-normalised empirical IMF, is $\varepsilon=0.4\eta^{1/\Gamma}$. For $\eta=1,0.69$ and $0.4$, $\varepsilon=0.4,0.3$ and $0.2$, respectively.

The results of this section justify our assumptions to model radiative feedback in Section 2.3 and it further implies that our simulations are self-consistently treating the formation of individual massive stars, their feedback effects, and can be used reliably to estimate of the escape fraction of hydrogen and helium ionising photons.

3.1.2 Resolution Studies

We conduct four extra runs at lower resolution to evaluate the numerical convergence of our simulations. All other simulations have the highest resolution we could afford computationally, and increasing the resolution is unfeasible for the relatively large cloud masses considered in this study. More massive sink particles and fewer low-mass sink particles form in simulations with lower resolution (see Figure 5). The mean mass of the IMF, represented by the vertical dashed line in the figure, increases by a factor $\sim 2-3$ when the spatial resolution halves. As the resolution increases, while there is no significant change in the total mass in sinks, the mean mass of sinks decreases, suggesting that some of the sinks fragment into smaller sub-clumps. A model in which the cores form stars with $\sim 30\%$ efficiency and, as we increase the resolution, additional small mass cores form from unused gas at the low-mass end of the CMF, is instead less consistent with our results for two reasons. i) Allowing more of the diffuse gas to form low mass cores, would produce, in some simulations, a total core formation efficiency above unity (in simulations that have f_* > 0.3-0.4). ii) Figure 5 shows that the core formation efficiency, which is $\sim f_{\mathrm{*}}$, is close to being converged. Thus, by increasing the resolution we do not add new cores from the gas, otherwise we would observe an increase of $f_{\mathrm{*}}$, which instead slightly decreases with increasing resolution. Increasing the resolution simply changes the CMF, but the total mass in cores remains nearly the same. For these reasons, we find that the cores-fragmentation model is the most likely, although we cannot rule out alternative scenarios.

3.1.3 Monte-Carlo Numerical Experiments for Fragmentation

To demonstrate more convincingly that our interpretation is robust, we perform a simple numerical experiment. We assume that each pre-stellar core (sink particle) fragments into smaller sub-units with a power-law mass function (MF) with index $\Gamma$, and with limits on the fragment masses between 0.1 M*⊙* and the sink mass. We draw randomly from this distribution until the total mass of the fragments equals the sink mass. We repeat this procedure for all the sinks. Such sampling is done 20 times and the average of the bins is taken. Figure 6 shows the resulting mass function for the XL-F cloud, but we obtain similar results for all the simulations. The mass function we obtain by fragmenting the sinks is shown as the solid histogram (blue for $\Gamma=1.35$ and orange for $\Gamma=0.8$). The original SMF is shown by the green solid histogram and the shifted mass function is shown by the dashed histogram. The black solid curves show the Chabrier IMF for a cluster with total mass equal to the total mass in sink particles, which are in very good agreement with the mass function of the fragmented pre-stellar cores assuming $\Gamma=0.8$.

This sampling method does not produce a modal mass for the IMF. To address this, we tried another sampling method. However, if we assume that the lower mass limit in the sampling is set to max$(0.01\leavevmode\nobreak\ {\rm M}_{\odot},0.01m_{\mathrm{sink}})$, instead of $0.1\leavevmode\nobreak\ {\rm M}_{\odot}$, the resulting IMF has a similar shape to the SMF but peaks at a mass 100 times smaller, resulting in a model mass of the IMF $\sim 0.1\leavevmode\nobreak\ {\rm M}_{\odot}$ (see bottom panel in Figure 6).

In summary, the fragmentation of the pre-stellar cores into numerous small mass stars, a process which is not captured in our simulations due to limited resolution, explains the deficit of stars with mass below $\sim 1\leavevmode\nobreak\ {\rm M}_{\odot}$ in our simulations with respect to the number expected assuming a Chabrier IMF.

3.1.4 High-mass Slope of the IMF

In this section we quantify more rigorously the slope of the IMF. In Figure 7 the IMF is shown at the end of the simulations when star formation has stopped. The green lines show the best fit power-law at the high-mass end of the IMF using Bayesian inferences as explained below. We do not notice any significant relationship between the high-mass end slope of the IMF and the mass or compactness of the cloud. We notice, however, a flattening of the IMF at 1-10 M*⊙* in very compact clouds of high-mass (M$>10^{3}$ M*⊙*), which is instead not observed in the fiducial and massive clouds. Since these clouds have the highest star formation efficiency and the strongest radiative feedback, a speculative interpretation would be that we are observing the effect of photo-evaporation of small fragments. When a proto-star is exposed to the ionising flux of a new-born OB star, the disk mass decreases rapidly with time. This may regulate the mass accretion rate through the disk and therefore to the star.

Here are some details of the Bayesian inference of the IMF slope. We assume a power-law slope mass distribution with general form $\mathrm{d}N/\mathrm{d}\log m=Am^{-\Gamma}$ where $A$ is a constant of normalisation and $m_{\mathrm{\rm min}}<m<m_{\mathrm{\rm max}}$. When the total number of stars is $N_{\mathrm{0}}$, this constant becomes $A=\ln{10}\,N_{\mathrm{0}}\leavevmode\nobreak\ \Gamma/(m_{\mathrm{\rm min}}^{-\Gamma}-m_{\mathrm{\rm max}}^{-\Gamma})$. The likelihood is proportional to the distribution function, $\mathcal{N}(\mu_{\mathrm{i}}|\Gamma)=Am_{\mathrm{i}}^{-\Gamma}$, where $\mu_{\mathrm{i}}\equiv\log(m_{\mathrm{i}})$. To find the most likely $\Gamma$ we calculate the value of $\Gamma$ that maximises the log of the likelihood:

[TABLE]

Model-independent constants are removed from this equation. We do the un-binned fitting only to stars with masses above a critical value. This value is chosen somewhat arbitrarily as the point at which the IMF starts to deviate from the Kroupa IMF. In each panel the best fit line is shown as a segment between the critical mass and the maximum stellar mass along with the slope $\Gamma$ and 1-$\sigma$ errors. The errors are calculated as the $16\%$ and the $84\%$ points of the cumulative likelihood for $\Gamma$ between $\Gamma=0.5$ to $1.8$. The fitted value of $\Gamma$ has a dependence on the critical minimum mass for the points included in the fit, but we find that the values of $\Gamma$ agree with a Kroupa IMF within the 1-$\sigma$ errors in most cases. Here, we adopt a critical mass such that particles above this mass account for $70\%$ of the total mass in stars.

3.1.5 Maximum Stellar Mass in the Cluster

Figure 8 shows the maximum stellar mass ($M_{\mathrm{\rm max}}$) as a function of the mass of the star cluster. The relationship between $M_{\mathrm{\rm max}}$, obtained by multiplying the maximum sink mass by 0.4, and the stellar cluster mass, $m_{\mathrm{cl}}$, is tight. The best fit power-law is

[TABLE]

where $m_{\mathrm{4}}=m_{\mathrm{cl}}/10^{4}\leavevmode\nobreak\ {\rm M}_{\odot}$, valid when $m_{\mathrm{cl}}\gtrsim 100\leavevmode\nobreak\ {\rm M}_{\odot}$. The relationship is well correlated, with a coefficient of determination R${}^{2}=0.93$. A power-law relationship between the maximum stellar mass and the cluster mass is consistent with observations, although the observed power-law slope is 0.45 (Larson, 1969), which is slightly flatter than the value found in our simulations. However, the slope we find is in good agreement with numerical studies of star formation in clusters using SPH codes (Bonnell et al., 2003, 2004). We also neglect smaller-scale feedback from protostellar outflows that can reduce the final mass of stars. In addition, it should be kept in mind that the maximum stellar mass here is defined as 0.4 the maximum sink mass, therefore it is possible that the fragmentation of the largest sinks may produce stellar masses systematically smaller than $M_{\mathrm{\rm max}}$.

Given the uncertainty due to Poisson statistical fluctuations, the SMF appears to be consistent with power-law all the way to the mass bin that is expected to have $\sim 1$ particle in it (the horizontal dashed-dotted lines in Figure 7). Hence we do not have strong evidence for a high-mass truncation of the CMF. We conclude that the CMF, as represented by the SMF, does not have a fundamental upper mass limit below $\sim 1000\leavevmode\nobreak\ {\rm M}_{\odot}$ (the maximum sink mass in all simulations). Since our simulations have the same initial turbulence field and we have only one random realisation for each set of parameters (mass, and density of the cloud), we are not able to address the question of whether the maximum stellar mass in a cluster is determined by physical (Kroupa & Weidner, 2003) or statistical effects (e.g. Fumagalli et al., 2011). In addition, we use an empirical relationship between sinks mass and massive stars, rather than resolving the fragmentation of sinks into massive stars using a physical model. This also prevents us from drawing robust conclusions about this open question.

3.2 Star Formation Efficiency

We define star formation efficiency (SFE, or $f_{\mathrm{*}}$) in our simulated clouds as the fraction of the initial gas mass that is converted into sink particles. Figure 9 shows the SFE as a function of time in units of the free-fall time $t_{\mathrm{\rm ff}}$ (shown at the top-right of each panel), for the simulations in Table 4. The top panel refers to the fiducial clouds, the middle panel to the compact clouds and the bottom panel to the very compact clouds. Lines in each panel refer to different cloud masses as explained by the simulation IDs in the legend. The vertical lines mark the time of the explosion of the first two SNe in the simulation, where the lifetimes of stars are given by Schaller et al. (1992) fitting functions. As discussed before we do not include mechanical feedback from SNe, but star formation has already stopped or it is mostly terminated before the explosion of the first SN in all simulation but XL-F, i.e. the fiducial run with mass $m_{\mathrm{gas}}=3.2\times 10^{5}$ M*⊙*.

When time is measured in units of the free-fall time, the shape of the SFE curves are qualitatively similar: the SFE increases rapidly with time and peaks at $t\approx 2-3t_{\mathrm{\rm ff}}$. Generally the total SFE at the end of the simulations increases with increasing cloud mass and with increasing cloud compactness.

This is shown more clearly in Figure 10. The top panel in Figure 10 shows the stellar mass of the cluster $m_{\mathrm{\rm cl}}$ as a function of the cloud gas mass for the 3 set of simulations with different compactness (as shown in the legend). The smaller open circle with the label $Z=1/40$ Z*⊙*, shows a compact cloud simulation but with lower gas metallicity (see Section 3.4). The dot-dashed line shows SFE$=100\%$, while the dashed lines are fits to the simulation results with the following function:

[TABLE]

where $n_{\mathrm{\rm cri}}\approx 10^{3}\leavevmode\nobreak\ \leavevmode\nobreak\ {\rm cm}^{-3}$ is the critical density and $m_{\mathrm{fl}}$ is the mass floor. The dashed lines show the fit assuming $m_{\mathrm{fl}}=0$, while the dotted line has $m_{\mathrm{fl}}=10$ M*⊙*. Equation (8) is a good fit to the points when excluding the 3 lowest mass simulations for the fiducial run (shown as smaller sized open squares). The motivation for excluding these 3 simulations from the fits is explained below.

The open symbols show star cluster that become dynamically unbound (i.e., open star clusters), while the solid symbols show star cluster that at the end of the simulations, after most of the gas has been used up for star formation or expelled, remain gravitationally bound (i.e., globular cluster progenitors).

The star symbols show the results of simulations by Jones & Bate (2018) for clouds with mass $m_{\mathrm{gas}}=500$ ${\rm M}_{\odot}$ and for mean densities $\overline{n}_{\mathrm{gas}}=$3\text{\times}{10}^{2}$$, $3\text{\times}{10}^{3}$, and $3\text{\times}{10}^{4}$ cm*-3*, from bottom to top, respectively. These densities are slightly different from the mean densities in our fiducial, compact and very compact simulations, thus we show as diamonds the corresponding points obtained using our fitting formula in Equation (8) with $m_{\mathrm{fl}}=10$ M*⊙. These simulations do not include feedback by massive stars being very small mass clouds in which the most massive star that forms has is $<10$ M⊙. However, the resolution of these simulations is higher than our simulations and, contrary to our simulations, feedback by IR radiation is included. In addition, these simulation are run using an SPH code. It is interesting to note that despite the different codes and physics included, the results are consistent with the extrapolation of our fitting formulae to low mass clouds if we assume a minimum mass floor for the star cluster mass of $\sim 10$ M⊙*.

The bottom panel in Figure 10 is the same as the top panel but shows the total star formation efficiency $f_{\mathrm{*,tot}}\equiv m_{\mathrm{\rm cl}}/m_{\mathrm{\rm gas}}$ and the best fit:

[TABLE]

The solid horizontal line at TSFE $\sim 15\%$ roughly separates star clusters that become globular cluster progenitors ($f_{\mathrm{*}}>15\%$) from open star clusters ($f_{\mathrm{*}}<15\%)$. This separation is based on the dynamical state of the cluster at the end of the simulations, but a more detailed analysis of the dynamics of the stellar cluster will be the subject of a followup study.

Let’s now address the reason why we excluded the 3 lower mass fiducial simulations from our analysis. We observe that the star cluster mass in these simulations does not obey a simple power-law relationship with the initial gas mass of the molecular cloud. The discrepancy does not appear to be a convergence issue due to insufficient resolution, as confirmed by the lower-resolution simulations (shown as lighter colour small squares), but rather lack of the necessary physics for self-regulation feedback. This can be understood inspecting Figure 8 which relates the mass of massive stars to the cloud gas mass. The low TSFE of the diffuse clouds in combination with the small cloud gas mass produces stellar masses below $10^{2}$ M*⊙, which corresponds to a maximum stellar mass $M_{\mathrm{max}}<10$ M⊙. Such stars do not produce significant quantities of ionising UV radiation, therefore the cloud can continue to form stars. This is due to our neglecting feedback mechanisms from lower mass stars. This requirement for stars that produce ionising radiation to disperse the cloud leads to a minimum cluster mass floor $m_{\mathrm{\rm cl}}\sim 300$ M⊙, much larger than the $\sim 10$ M⊙* floor which is a good fit to the simulations of Jones & Bate (2018).

This large mass floor is not evident in the compact and very compact clouds: if it exists, it must be at masses $m_{\mathrm{\rm cl}}<100$ M*⊙*. The reason for this apparent inconsistency is not fully understood, but it appears to be related to the smaller ratio of the crossing to free-fall time for the fiducial cloud when compared to the more compact clouds (Sec. 3.3). We offer the following hypothesis: Inspecting the middle and bottom panels in Fig. 9, we observe a longer delay for onset of star formation in the small mass clouds for the compact and very compact runs, which is not observed in the fiducial runs. This can be understood in terms of the necessary number of crossing times required by the supersonic turbulence to create dense clumps for star formation (with $n>n_{\mathrm{\rm sink}}$). In the fiducial cloud this enhancement of the density due to supersonic turbulence is faster when compared to the free-fall timescale, hence the steeper rise of $f_{\mathrm{*}}$ as a function of time. When feedback from massive stars is absent due to random sampling of a small mass stellar cluster, this rapid increase of $f_{\mathrm{*}}$ can lead to significant overshooting of star formation above the threshold expected from self-regulation. This overshooting does not happen, or is milder, for more compact clouds in which $f_{\mathrm{*}}$ increases with $t/t_{\mathrm{ff}}$ more slowly. The existence of a minimum cluster mass floor, however, should eventually become evident also in more compact clouds when decreasing further the initial cloud masses.

We observe a power-law relation between the mass of the cloud and the mass of the star cluster. Howard et al. (2018) find that the stellar mass of the most massive cluster that forms from a molecular cloud has a power-law dependence on the mass of the cloud with an exponent of $0.78$. In our work, this relation, taking all sink particles as the cluster, has an exponent of 1.4 (Equation 8). By multiplying it with the exponent of the $M_{\mathrm{max}}$-$m_{\mathrm{cl}}$ relation, 0.66 (Equation 7), we get an exponent of 0.92. Similar to Howard et al. (2018), our work suggests that young massive star clusters are natural extensions of low-mass cluster formation.

Our results (Figure 10, or Table 4) are in good agreement with Kim et al. (2018), who find that the TSFE depends primarily on the initial gas surface density, such that the TSFE increases from $4\%$ to $51\%$ as $\Sigma$ increases from 13 to 1300 ${\rm M}_{\odot}\leavevmode\nobreak\ {\rm pc}^{-2}$.

To summarise, we believe that the increase in TSFE observed for the fiducial simulations with masses $m_{\mathrm{\rm gas}}\leq 10^{4}$ M*⊙* is unphysical, meaning that it is due to missing feedback processes in our simulations. When the most massive star has mass $M<10$ M*⊙*, IR radiation feedback or proto-stellar jets feedback should be included in the simulation. In all the other simulations UV feedback by massive stars is likely the dominant feedback at play; therefore these simulations incorporate the relevant physics for the formation of realistic star clusters.

3.3 Star formation law in molecular clouds

Next, we ask the question of what is the physical interpretation of the empirical relationship we derived for the star formation efficiency as a function of cloud mass and compactness. To answer this question we first fit the SFE $f_{\mathrm{*}}(\tau)$ with an analytic function, where $\tau\equiv t/t_{\rm ff}$, in order to minimise the stochastic noise of the simulations. The $f_{\mathrm{*}}(\tau)$ has a shape that can be fit by an $\arctan$ function or the Fermi function:

[TABLE]

Both fits give similar results for the purpose of interpreting $f_{\mathrm{*}}(\tau)$. In Figure 11 we show the fit to $f_{\mathrm{*}}(\tau)$ using the Fermi function $f_{\mathrm{F}}$ (orange solid curves) and its time derivative (blue curves), or the dimensionless SFR per free-fall time, SFR${}_{\mathrm{ff}}\equiv\mathrm{d}f_{\mathrm{*}}/\mathrm{d}\tau\approx\mathrm{d}f_{\mathrm{F}}/\mathrm{d}\tau$. The fits are a good approximations to the data points from the simulations (solid points), except for a few clouds where $f_{\mathrm{*}}(\tau)$ has a pit near the end of the star formation process.

The value of the peak of SFRff has a weak dependence on the cloud mass (see top panel in Figure 12) and a stronger dependence on the cloud mean density. We fit the SFR${}_{\mathrm{ff}}|_{\mathrm{\rm max}}$ with a power-law similar to Eq. (9):

[TABLE]

where $\alpha_{\mathrm{f}}\approx 1.0$ and $n_{\mathrm{cri}}$ is the same critical density as in Eq. (8)222The value of $\alpha_{\mathrm{f}}$ is somewhat correlated with $n_{\mathrm{\rm cri}}$. We sample a sequence of $n_{\mathrm{\rm cri}}$ for which we obtain a good fit and find that for $n_{\mathrm{\rm cri}}$ in the range $\sim 400$ – $1600\leavevmode\nobreak\ {\rm cm}^{-3}$, the corresponding $\alpha_{\mathrm{\rm f}}$ is in the range $0.85\--1.1$.. The duration of the star formation burst in units of $t_{\rm ff}$, $\Delta\tau_{\mathrm{SF}}$, is proportional to the width of the SFRff shown as the blue lines in Fig. 11. The function $\mathrm{d}f_{\mathrm{F}}/\mathrm{d}\tau$ has a peak value $f_{\mathrm{0}}/4\Delta\tau$ and a full-width half-maximum $3.526\Delta\tau$. We define $\Delta\tau_{\mathrm{SF}}\equiv 4\Delta{\tau}$ so that

[TABLE]

Inspecting Fig. 11 we see that $\Delta\tau_{\mathrm{SF}}$ increases with the cloud mass, and appears to be proportional to the dimensionless sound crossing time of the cloud. Here we define the sound crossing time, $t_{\mathrm{\rm cr}}$, as the ratio of the time it takes for a sound wave with $c_{\mathrm{s}}=10$ km/s to cross the cloud radius. Similarly to the dimensionless $\Delta\tau_{\mathrm{\rm SF}}$, we define $\tau_{\mathrm{cr}}\equiv t_{\mathrm{\rm cr}}/t_{\rm ff}$, where the free-fall time is defined at the cloud’s mean density. We find that $\Delta\tau_{\mathrm{\rm SF}}/\tau_{\mathrm{\rm cr}}=\Delta t_{\mathrm{\rm SF}}/t_{\mathrm{\rm cr}}\approx 6$ (the horizontal line in the bottom panel of Fig. 12). This results makes physical sense because the feedback mechanism stops star formation by creating over-pressured H ii regions which require a constant number of crossing times to expel the gas.

Since $t_{\mathrm{\rm cr}}\propto r_{\mathrm{\rm gas}}\propto(m_{\mathrm{gas}}/\overline{n})^{1/3}$, we have $\Delta\tau_{\mathrm{SF}}\propto t_{\mathrm{cr}}/t_{\rm ff}\propto m_{\mathrm{gas}}^{1/3}\overline{n}^{1/6}$. From Equation (12) we derive $f_{\mathrm{*,\rm tot}}\propto m_{\mathrm{gas}}^{\;\;0.69}\overline{n}_{\mathrm{gas}}^{\;\;0.17}(1+\overline{n}_{\mathrm{gas}}/n_{\mathrm{\rm cri}})^{1.0}$, which is in good agreement with Eq. (9) for $\overline{n}>n_{\mathrm{\rm cri}}$. The agreement can be improved further by considering a more accurate fit to $\tau_{\mathrm{SF}}/\tau_{\mathrm{cr}}$ rather than assuming a constant value $\sim 6$. Namely, considering the weak dependence of the star formation timescale on the cloud mass and density: $\Delta\tau_{\mathrm{SF}}/\tau_{\mathrm{cr}}\propto m_{\mathrm{gas}}^{\;\;-0.3}\overline{n}_{\mathrm{gas}}^{\;\;-0.2}$.

From the analysis and interpretation of these results we can thus derive a star formation law in molecular clouds that can be used as a more accurate sub-grid recipe in cosmological simulations that resolve the molecular cloud phase. Assuming a constant mean volume for the cloud we have $f_{\mathrm{*}}\equiv m_{\mathrm{*}}/m_{\mathrm{\rm gas}}\approx\rho_{\mathrm{*}}/\rho_{\mathrm{\rm gas}}$. Therefore, assuming $\rho_{\mathrm{\rm gas}}={\rm const}$ (i.e., assuming $f_{\mathrm{*}}\ll 1$) during the episode of star formation, which has a duration $\Delta t_{\mathrm{SF}}$, we have SFR${}_{\mathrm{ff}}|_{\mathrm{max}}\equiv df_{\mathrm{*}}/d\tau|_{\mathrm{\rm max}}\approx d\rho_{\mathrm{*}}/dt|_{\mathrm{\rm max}}(t_{\mathrm{\rm ff}}/\rho_{\mathrm{\rm gas}})$, which implies:

[TABLE]

with $\epsilon=1.1\%$ for solar metallicity and $\epsilon=0.36\%$ for $Z<0.1$ Z*⊙* (see § 3.4). A star formation law $d\rho_{\mathrm{*}}/dt\propto\rho_{\mathrm{\rm gas}}^{n}$ with $n=1$ or $n=1.5$ is most often used as a sub-grid star formation recipe in cosmological simulations. Therefore we suggest that a steeper power-law index $n\sim 2.5$ is a better description of the star formation rate at densities typical of molecular clouds in high-redshift galaxies. This theoretical result can, in principle, be tested against observations of young stellar clusters in our galaxy.

Krumholz et al. (2012) suggests that star formation law is universal in which the star formation rate is $\sim 1.5\%$ of the molecular gas mass per local free-fall time. Eq. 11 results in SFR${}_{\mathrm{ff}}\approx 1\%-2\%$ at $\overline{n}_{\mathrm{gas}}\lesssim 10^{3}$\mathrm{c}\mathrm{m}^{-3}$$, in agreement with this work for local molecular clouds. However, Krumholz et al. (2012) finds this universal value also for high-redshift galaxies but averaged over the whole galaxy. We find that SFRff can be as large as $\sim 10\%$ for more compact clouds typical of high-redshift galaxies, and/or more massive clouds (see also the top panel of Fig. 12). A direct comparison to Krumholz’s results is not trivial for the galaxy as a whole, as it depends on modelling the multi-phase ISM of high-z galaxies.

3.4 Effects of Lowering the Gas Metallicity

The set of compact and very compact molecular clouds we have analysed are meant to represent clouds typical of the ISM in dwarf galaxies forming at high-redshift. However, we also know that the gas metallicity in these dwarf galaxies is less than solar. In order to keep the parameter study consistent we have not changed the gas metallicity in the compact and very compact clouds, but in this section we briefly test the influence of gas metallicity [Fe/H] on the star formation rate and IMF. In our simulations, changing the gas metallicty affects the cooling of the gas (see Section 2.4).

Figure 13 is the same as Figure 7 but for the LC clouds with metallicity $Z=1\leavevmode\nobreak\ Z_{\mathrm{\odot}}$, $0.1\leavevmode\nobreak\ Z_{\mathrm{\odot}}$, and $0.025\leavevmode\nobreak\ Z_{\mathrm{\odot}}$. The shape of the IMF is not affected by the gas metallicity. Only the normalisation of the IMF is influenced because of the lower SFE in the low-metallicity simulations. This is in agreement result of previous theoretical works (e.g. Myers et al., 2011; Bate, 2014).

Lower metallicity translates into lower cooling rates, which should result in lower efficiency of star formation. Figure 14 shows $f_{\mathrm{*}}$ as a function of time for the large compact cloud (LC) with intermediate (0.1 $Z_{\mathrm{\odot}}$) and low (0.025 $Z_{\mathrm{\odot}}$) metallicity. The effect of lowering the metallicity by a factor of ten, from $Z=1$ $Z_{\mathrm{\odot}}$ to $Z=0.1$ $Z_{\mathrm{\odot}}$ is to lower $f_{\mathrm{*}}$ at the end of the simulation by roughly a factor of 5. But lowering further the metallicity from $Z=0.1$ $Z_{\mathrm{\odot}}$ to $Z=0.025$ $Z_{\mathrm{\odot}}$ does not change $f_{\mathrm{*}}$, suggesting that $f_{\mathrm{*}}$ decreases almost linearly with the metallicity from solar to $Z=0.2$ $Z_{\mathrm{\odot}}$, but this effect saturates when further lowering the metallicity. The SFE decreases mainly because the peak SFR decreases by roughly a factor of 3 with decreasing metallicity, while the duration of the star formation episode is nearly unchanged (see small circles in Fig. 12).

In order to better understand what is causing a decrease of the SFR at lower metallicity, we have analysed the density and temperature structure of these two simulations. We found that lowering the metallicity causes the temperature and the thermal pressure inside H ii regions to increase by roughly a factor of $3$, as shown in Figure 15. This result is in agreement with observations and theoretical models of H ii regions. The strength of feedback, due to the increase of thermal pressure inside the H ii regions, is therefore stronger at lower metallicity, resulting in a lower star formation efficiency. This result on the effect of the gas metallicity goes in the opposite direction of what found by Howard et al. (2018). In their work, lowering the metallicity of the gas cloud reduces the opacity of the gas to radiation and results in higher gas accretion which leads to an increase of the total star formation efficiency. However, this can be understood because in their simulations the dominant feedback mechanism is IR radiation pressure while, contrary to our work, UV feedback does not play a major role. However, their simulations describe more massive clouds and have much lower resolution than the simulations in our work.

4 Summary and Conclusions

In this paper, the first of a series, we present a large set of radiation-magneto-hydrodynamic simulations of star formation in self-gravitating, turbulent molecular clouds. The initial conditions for the clouds are isothermal spheres initially close to virial equilibrium, being supported by turbulent motions.

We model the formation of individual massive stars, replacing self-gravitating clumps that are collapsing below the resolution of the simulations with sink particles, which represent individual massive stars, therefore including their UV radiation feedback self-consistently. We consider a grid of simulations varying the cloud masses between $m_{\mathrm{\rm gas}}=10^{3}$ M*⊙* to $3\times 10^{5}$ M*⊙. Depending on the cloud mass, we resolve scales between 200 AU to 2000 AU. In addition, we consider three compactness for the molecular clouds. The fiducial clouds have gas mean number densities typical of those observed in the local universe ($\overline{n}_{\mathrm{\rm gas}}=1.8\times 10^{2}$ cm-3*). Compact ($\overline{n}_{\mathrm{\rm gas}}=1.8\times 10^{3}$ cm*-3*) and very compact ($\overline{n}_{\mathrm{\rm gas}}=1.8\times 10^{4}$ cm*-3*) clouds represent clouds expected to exist in high-redshift galaxies. We also partially explore varying the gas metallicity. Our goal is to run a realistic set of simulations of formation of star clusters in molecular clouds to understand the physics of star formation across cosmic time: from conditions typical of present-day ISM to the the higher-pressure environments found in the ISM of higher redshift galaxies.

In this paper we focus on understanding the IMF, the SFR and SFE as a function of the cloud mass and compactness. We derive a star formation law valid at densities typical of high-redshift molecular clouds that will help to justify and inform the sub-grid star formation recipe used in cosmological simulations.