Self-consistent proto-globular cluster formation in cosmological simulations of high-redshift galaxies

TL;DR

This study uses high-resolution cosmological simulations to explore how bound star clusters, similar to ancient globular clusters, form in high-redshift galaxies, revealing the conditions and efficiencies involved.

Contribution

It demonstrates for the first time that bound clusters form self-consistently in cosmological simulations across various high-redshift galaxy environments.

Findings

Bound clusters form in high-pressure, turbulent environments.

Cluster mass function follows a power-law of dN/dM~M^-2.

Cluster formation efficiency is consistent across galaxy masses.

Abstract

We report the formation of bound star clusters in a sample of high-resolution cosmological zoom-in simulations of z>5 galaxies from the FIRE project. We find that bound clusters preferentially form in high-pressure clouds with gas surface densities over 10^4 Msun pc^-2, where the cloud-scale star formation efficiency is near unity and young stars born in these regions are gravitationally bound at birth. These high-pressure clouds are compressed by feedback-driven winds and/or collisions of smaller clouds/gas streams in highly gas-rich, turbulent environments. The newly formed clusters follow a power-law mass function of dN/dM~M^-2. The cluster formation efficiency is similar across galaxies with stellar masses of ~10^7-10^10 Msun at z>5. The age spread of cluster stars is typically a few Myrs and increases with cluster mass. The metallicity dispersion of cluster members is ~0.08 dex in [Z/H] and does not depend on cluster mass significantly. Our findings support the scenario that present-day old globular clusters (GCs) were formed during relatively normal star formation in high-redshift galaxies. Simulations with a stricter/looser star formation model form a factor of a few more/fewer bound clusters per stellar mass formed, while the shape of the mass function is unchanged. Simulations with a lower local star formation efficiency form more stars in bound clusters. The simulated clusters are larger than observed GCs due to finite resolution. Our simulations are among the first cosmological simulations that form bound clusters self-consistently in a wide range of high-redshift galaxies.

Tables

Table 1. Table 1: Initial conditions of simulations studied in this paper.

Simulation	$M_{\mathrm{halo}}$ ($z=5$)	$M_{\ast }$ ($z=5$)	$m_b$	${ϵ}_{\mathrm{gas}}$	${ϵ}_{\mathrm{star}}$	$m_{\mathrm{DM}}$	${ϵ}_{\mathrm{DM}}$	Notes
	[$M_{☉}$]	[$M_{☉}$]	[$M_{☉}$]	[pc]	[pc]	[$M_{☉}$]	[pc]
z5m10a_hr	$6.6\times {10}^9$	$1.5\times {10}^7$	119.3	0.14	0.7	6.5e2	10
z5m11c_hr	$7.6\times {10}^{10}$	$9.4\times {10}^8$	890.8	0.28	1.4	4.9e3	21	8 times better mass resolution than z5m11c
z5m11c	$7.8\times {10}^{10}$	$7.8\times {10}^8$	7126.5	0.42	2.1	3.9e4	42	8 times lower mass resolution than z5m11c_hr
z5m12b	$8.6\times {10}^{11}$	$2.2\times {10}^{10}$	7126.5	0.42	2.1	3.9e4	42

Table 2. Table 2: Simulation restarts. From each parent simulation, a starburst is re-started and run for the cosmic time and redshift labeled (see Fig. 2 ) with different star formation models. We consider the following criteria for star formation: (i) molecular ( mol ), (ii) self-gravitating ( sg ), (iii) density threshold ( den ), and (iv) converging flow ( cf ) (as defined in Sections 2.1 and 2.2 ). Each model name refers to a combination of star formation criteria with certain choice of local star formation efficiency ϵ ff subscript italic-ϵ ff \epsilon_{\rm ff} listed below (see Section 2.2 for details). We refer to these re-simulations by Parent simulation_Model name (e.g. z5m11c_hr_G18). More detailed comparison and discussion between different star formation models are presented in Section 5 .

Parent	Redshift	Cosmic time	$M_{\mathrm{halo}}$	Model	Star formation	${ϵ}_{\mathrm{ff}}$	$M_{\ast }$ formed	$f_{\mathrm{bound}}$
simulation		[Gyr]	[$M_{☉}$]	name	criteria		[$M_{☉}$]
z5m10a_hr	6.605–6.143	0.820–0.901	$5.8\times {10}^9$	G18	mol+sg+cf	1	$1.57\times {10}^7$	0.26
z5m11c_hr	5.562–5.024	1.023–1.163	$7.6\times {10}^{10}$	FIRE	mol+sg+den	1	$2.56\times {10}^8$	0.28
				no $n_{\mathrm{th}}$	mol+sg	1	$3.50\times {10}^8$	0.17
				G18	mol+sg+cf	1	$4.44\times {10}^8$	0.26
				G18_e50	mol+sg+cf	0.5	$3.00\times {10}^8$	0.39
z5m11c	5.562–5.024	1.023–1.163	$7.6\times {10}^{10}$	G18	mol+sg+cf	1	$1.55\times {10}^8$	0.12
z5m12b	6.333–6.080	0.866–0.913	$4.5\times {10}^{11}$	G18	mol+sg+cf	1	$3.31\times {10}^9$	0.33

Equations

$$\alpha \equiv \frac{\Vert \nabla \otimes \mathbf{v}{\Vert }_i^2+(c_{\mathrm{s},i}/h_i{)}^2}{8\pi G{\rho }_i}<1,$$

Full text

††pagerange: Self-consistent proto-globular cluster formation in cosmological simulations of high-redshift galaxies–References

1 Introduction

Globular clusters (GCs) are spherical, tightly bound collections of stars that are typically 5–13 Gyr old, massive ($10^{4}$–$10^{6}\,M_{\sun}$), metal-poor ($\rm-2.5\lesssim[Fe/H]\lesssim 0.3$), and compact (with half-light radii of 0.5–10 pc) (see e.g. Harris, 1991; Brodie & Strader, 2006; Bastian & Lardo, 2018, for a series of reviews). GCs are found in almost all galaxies above $M_{\ast}\sim 10^{9}\,M_{\sun}$ in the local Universe.

The origin of GCs remains a longstanding puzzle in the context of star formation and galaxy formation. A number of theories have been proposed to explain their formation. For example, Peebles & Dicke (1968) suggested that GCs were the first bound structures formed in the universe when the Jean mass following recombination is $10^{5}$–$10^{6}\,M_{\sun}$. Similarly, Fall & Rees (1985) suggested that GCs were formed via thermal instabilities of metal-poor gas in galactic halos where the Jeans mass is $10^{6}\,M_{\sun}$ (see also Kang et al., 1990; Shapiro et al., 1992). Peebles (1984) has suggested that GCs were formed in dark matter (DM) halos at high redshift that merged into more massive halos at later times. Recent cosmological simulations have shown that GC-like objects can form in $10^{6}$–$10^{8}\,M_{\sun}$ halos in the very early universe (e.g. Trenti et al., 2015; Kimm et al., 2016; Ricotti et al., 2016). Naoz & Narayan (2014) suggested GCs can form without a DM halo in the early Universe due to the relative stream velocity between DM and baryons induced at the epoch of recombination (see also e.g. Chiou et al., 2019). Some GCs may also be the nuclei or nuclear clusters of former dwarf galaxies that were tidally disrupted (e.g. Mackey & van den Bergh, 2005).

Since the 1980s, a population of stellar systems called young massive clusters (YMCs) have been found in some extreme environments in the local Universe, including galaxy mergers, starburst galaxies, and galactic nuclei (see Portegies Zwart et al., 2010, for a recent review). The YMCs are younger than 1 Gyr, with masses of $10^{4}$–$10^{8}\,M_{\sun}$ and half-light radii only a few pc (e.g. Larsen, 2004). The YMCs broadly follow a mass function of ${\rm d}N/{\rm d}M\sim M^{-2}$ with a truncation at the high-mass end (e.g. Zhang & Fall, 1999; Gieles et al., 2006; Larsen, 2009). Given the similarities between GCs and YMCs, many authors have suggested that present-day GCs were in fact YMCs formed at early times (see Longmore et al., 2014; Kruijssen, 2014, for recent reviews).

YMCs are preferentially formed in highly gas-rich, turbulent conditions, with gas surface densities of $\Sigma_{\rm gas}\sim 10^{2}$–$10^{4}\,M_{\sun}\,{\rm pc}^{-2}$ and turbulent velocities of $\sigma_{v}\sim 10$–$100\,{\rm km}\,{\rm s}^{-1}$, translating to pressures of $P/k\sim 10^{6}$–$10^{8}$ K ${\rm cm}^{-3}$ (see e.g. Longmore et al., 2014).111Throughout this paper, when we refer to a ‘high-pressure’ region where bound clusters form, by ‘pressure’ we mean the gravitational force per unit area $P_{\rm grav}\sim G\Sigma_{\rm gas}^{2}$. This can be balanced by turbulent ram pressure $P_{\rm turb}\sim\rho\sigma_{v}^{2}$, which is comparable to $P_{\rm grav}$ (following Elmegreen & Efremov, 1997). We do not mean thermal pressure as the gas is usually cold. Such pressures are almost three orders of magnitude higher than those in our Milky Way (MW) and typical disk galaxies in the local Universe (see e.g. Kruijssen & Longmore, 2013), but are commonly seen in high-redshift galaxies (e.g. Kruijssen, 2014, and references therein). Both theories and numerical simulations showed that in such high-pressure regions, the integrated star formation efficiency over a cloud (i.e. the fraction of the cloud mass that turns into stars) can reach above 50% and hence a large fraction of stars from in bound clusters (e.g. Fall et al., 2010; Skinner & Ostriker, 2015; Kim et al., 2016; Grudić et al., 2018a; Tsang & Milosavljević, 2018; Li et al., 2019, for recent studies).222For example, Grudić et al. (2018a) found the integrated star formation efficiency of a cloud is set by its surface density, $\epsilon_{\rm int}\sim(1+\Sigma_{\rm crit}/\Sigma_{\rm gas})^{-1}$. This can be derived assuming star formation in the cloud stops when feedback from newly formed stars becomes larger than the self-gravity of the cloud. They found $\Sigma_{\rm crit}\sim 3000\,M_{\sun}\,{\rm pc}^{-2}$ using high-resolution cloud-scale simulations. Most stars form in one free-fall time of the cloud, so $\epsilon_{\rm int}$ is also the star formation efficiency per cloud free-fall time. We will refer back to these results later in this paper. This further supports the idea that GCs form in regular star formation in high-redshift galaxies where bound clusters can form efficiently (see e.g. Elmegreen & Efremov, 1997; Elmegreen et al., 2012; Kruijssen, 2012).

Resolving the formation of YMCs or proto-GCs in cosmological simulations of galaxy formation is a challenging problem. This is due to the multi-scale nature of this task: a wide dynamic range is required to simultaneously capture cosmic structure formation, hierarchical fragmentation of the interstellar medium (ISM), and star formation and feedback physics in giant molecular clouds (GMCs). State-of-the-art cosmological simulations specifically designed to study GC formation and evolution down to $z=0$ cannot yet resolve substructures in the GMCs. A commonly adopted approach is to spawn tracer particles to represent bound clusters. When and where cluster particles form are determined empirically according to local gas properties (see e.g. Kravtsov & Gnedin 2005; Li et al. 2017; the E-MOSAICS simulations, Pfeffer et al. 2018; Kruijssen et al. 2019). Such models may include free parameters for physics on under-resolved scales (e.g. Li et al., 2018).

A different family of models for GC formation and evolution across cosmic time use semi-analytic approach built upon DM halo merger trees (e.g. Beasley et al., 2002; Muratov & Gnedin, 2010; Li & Gnedin, 2014; Kruijssen, 2015; Choksi et al., 2018; El-Badry et al., 2019). These models usually determine GC formation rates by halo properties and/or galaxy-scale gas conditions, such as merger rates (e.g. Li & Gnedin, 2014), halo growth rates (e.g. Choksi et al., 2018), and average gas surface densities (e.g. El-Badry et al., 2019). These models are complementary to cosmological simulations.

Thanks to the dramatically increased speed of supercomputers over the past few years, cosmological simulations now can achieve unprecedented resolution. In the meanwhile, models used in state-of-the-art simulations have also been substantially improved, making it possible to better capture the multi-phase ISM, star formation in GMCs, disruption of GMCs by stellar feedback, and the launch and propagation of galactic winds, at least in high-resolution cosmological zoom-in simulations (e.g. Ceverino et al., 2014; Hopkins et al., 2014, 2018c; Wang et al., 2015; Agertz & Kravtsov, 2016). Recently, Kim et al. (2018) reported a proto-GC formed during a gas-rich merger of two low-mass $z\sim 7$ galaxies in a cosmological zoom-in simulation from the Feedback in Realistic Environments (FIRE) project.333https://fire.northwestern.edu Mandelker et al. (2018) found that GC-like objects formed via fragmentations of cold streams from the cosmic web in a cosmological zoom-in simulation of a $z\sim 6$ galaxy.

In this paper, we study a suite of high-resolution cosmological zoom-in simulations of galaxies at $z\geq 5$, where we use physically motivated models for the multi-phase ISM, star formation, and stellar feedback. The simulations are only run to $z=5$, so that we can use sufficiently high mass resolution at a reasonable computational cost. The goal of this paper is to show that YMCs/proto-GCs form in our simulations over a wide range of galaxy masses and identify the key processes resulting in cluster formation. We will show that newly formed clusters broadly follow a power-law mass function of ${\rm d}N/{\rm d}M\sim M^{-2}$ as arises naturally from hierarchical structure formation (e.g. Elmegreen & Falgarone, 1996; Fujii & Portegies Zwart, 2013; Guszejnov et al., 2018). Our simulations are among the first to form bound clusters in ‘regular’ star formation in gas-rich high-redshift galaxies self-consistently.

Being able to explicitly resolve the formation of YMCs/proto-GCs in cosmological simulations is important for the following reasons. First of all, it complements high-resolution cloud-scale simulations of cluster formation (e.g. Mac Low & Klessen, 2004; McKee & Ostriker, 2007; Krumholz, 2014, and references above) by bridging small- and large-scale simulations together. Secondly, it is also complementary to previous studies on cluster formation in analytic models (e.g. Elmegreen & Efremov, 1997; Kruijssen, 2012; Grudić et al., 2018b) and idealized simulations of isolated galaxies and gas-rich mergers (e.g. Li et al., 2004; Kruijssen et al., 2011, 2012; Hopkins et al., 2012; Hopkins et al., 2013a; Maji et al., 2017; Lahén et al., 2019; Moreno et al., 2019) by adding a cosmological environment. Finally, doing so will provide insights and calibrate empirical models implemented in other simulations and semi-analytic models.

Understanding cluster formation in $z\geq 5$ galaxies is also important for understanding the galaxy populations at these redshifts and the reionization history. It has been suggested that the progenitors of present-day GCs may have contributed significantly to the cosmic star formation rate densities at $z\geq 5$ and thus the ionizing photon budgets available for cosmic reionization (e.g. Boylan-Kolchin, 2018; Zick et al., 2018). Recently, Bouwens et al. (2017b) reported that dwarf galaxies at $z\sim 6$–8 found in the Hubble Frontier Fields (HFFs) have very small sizes, with individual sources reaching as small as $\sim 10$–30 pc, similar to star cluster complexes and/or super star clusters seen in the local Universe (see also e.g. Bouwens et al., 2017a; Kawamata et al., 2018). Vanzella et al. (2017a, b, 2019) also reported a superdense star-forming region at $z=6.143$ with an effective radius less than 13 pc and a star formation rate density of $\sim 1000\,M_{\sun}\,{\rm yr}^{-1}\,{\rm kpc}^{-2}$ behind one of the HFF clusters. It is quite possible that these small objects in the early Universe are newly formed star clusters/YMCs/proto-GCs that may eventually become present-day GCs. These clusters might bias the interpretation of the faint-end galaxy ultraviolet luminosity functions (UVLFs) at $z\geq 5$, as they can be preferentially selected in surface-brightness-limited surveys (see e.g. Bouwens et al., 2017b; Ma et al., 2018a).

This paper is organized as follows. In Section 2, we introduce our simulation sample and describe the star formation and feedback physics adopted in these simulations. In Section 3, we present case studies to show the key processes involved in cluster formation. We also present the mass functions for newly formed clusters in these simulations. In Section 4, we show the properties of these clusters. In Section 5, we show how our results are affected by the details in our star formation model. We discuss the implications of our results in Section 6 and conclude in Section 7.

We adopt a standard flat $\Lambda$CDM cosmology with Planck 2015 cosmological parameters $H_{0}=68{\rm\,km\,s^{-1}\,Mpc^{-1}}$, $\Omega_{\Lambda}=0.69$, $\Omega_{m}=1-\Omega_{\Lambda}=0.31$, $\Omega_{b}=0.048$, $\sigma_{8}=0.82$, and $n=0.97$ (Planck Collaboration et al., 2016). We use a Kroupa (2002) initial mass function (IMF) from 0.1–$100\,M_{\sun}$, with IMF slopes of $-1.30$ from 0.1–$0.5\,M_{\sun}$ and $-2.35$ from 0.5–$100\,M_{\sun}$.

2 Methods

2.1 Cosmological zoom-in simulations

Our study is based on a sample of cosmological zoom-in simulations from the FIRE project, which we summarize in Table 1. The initial conditions for halos z5m10a, z5m11c, and z5m12b are first introduced in Ma et al. (2018b). These three halos are selected at $z=5$ with approximate halo mass $10^{10}$, $10^{11}$, and $10^{12}\,M_{\sun}$, respectively, at this redshift. We have run halos z5m10a and z5m11c with 8 times better mass resolution than those used in our previous work, which we refer to as z5m10a_hr and z5m11c_hr. We include halo z5m11c run at two different resolution for comparison. Halos z5m11c and z5m12b are identical to those used in Ma et al. (2019). In Table 1, we list the mass resolution and force softening lengths for baryonic and dark matter (DM) particles adopted in these simulations. All the multi-scale cosmological zoom-in initial conditions are generated with the code music (Hahn & Abel, 2011).

The simulations are run using an identical version of the code gizmo444http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html (Hopkins, 2015) in the mesh-less finite-mass (MFM) mode. We use the FIRE-2 models for the multi-phase ISM, star formation, and stellar feedback. We briefly review and highlight the key ingredients that are important to this study below, but refer to Hopkins et al. (2018b, c) for more details of the numerical implementation and tests. Gas follows an ionized+atomic+molecular cooling curve between 10–$10^{10}$ K, including metallicity-dependent fine-structure and molecular cooling at low temperatures and metal-line cooling at high temperatures. The ionization states and cooling rates for H and He are computed following Katz et al. (1996) and cooling rates for heavy elements are calculated from a compilation of cloudy runs (Ferland et al., 2013), applying a uniform, redshift-dependent ionizing background from Faucher-Giguère et al. (2009) and heating from local sources. Self-shielding is accounted for with a local Sobolev/Jeans-length approximation.

Star formation is allowed only if the following criteria are met:

(i) Molecular (mol). The self-shielded molecular fraction of each gas particle ($f_{\rm mol}$) is estimated following Krumholz & Gnedin (2011). Star formation only takes place in molecular gas.

(ii) Self-gravitating (sg). This requires that the gravitational potential energy dominates over kinetic plus thermal energy at the resolution scale, as described by the virial parameter

[TABLE]

where $\otimes$ is the outer product, $c_{\rm s}$ is the sound speed, $h$ is the resolution scale, and the subscript $i$ implies that the quantities are evaluated for individual gas particles (Hopkins et al., 2013b).

(iii) Density threshold (den). The number density of hydrogen exceeds a threshold of $n_{\rm H}>n_{\rm th}=1000\,{\rm cm}^{-3}$.

If all criteria above are met, a gas particle will turn into a star particle at a rate $\dot{\rho}_{\ast}=\epsilon_{\rm ff}\,f_{\rm mol}\,\rho/t_{\rm ff}$, where $\epsilon_{\rm ff}$ is the local star formation efficiency and $t_{\rm ff}$ is the free-fall time of the particle. We adopt $\epsilon_{\rm ff}=1$ by default. We note that this efficiency reflects the rate at which locally self-gravitating clumps fragment, while the realized star formation efficiency of GMCs is regulated by feedback at $\sim 1$–10% per cloud free-fall time for typical MW GMC conditions (see e.g. Hopkins et al., 2018c, and reference therein). Only in extremely high-surface-density regions ($\gg 10^{3}\,M_{\sun}\,{\rm pc}^{-2}$) where feedback is no longer sufficiently strong to overcome gravity can the cloud-scale star formation efficiency reach near unity (e.g. Grudić et al., 2018a). We will vary $\epsilon_{\rm ff}$ and consider criteria for star formation different from the FIRE-2 model in Section 2.2.

Every star particle is treated as a single stellar population with known age, mass, and metallicity (inherited from its parent gas particle). All feedback quantities are calculated directly from standard stellar population synthesis models starburst99 (Leitherer et al., 1999) assuming a Kroupa (2002) IMF. The simulations account for the following feedback mechanisms: (1) photoionization and photoelectric heating, (2) radiation pressure for UV/optical single scattering and multiple scattering of re-radiated IR photons (the latter is usually subdominant in the simulations), and (3) energy, momentum, mass, and metal injection from discrete supernovae and continuous stellar winds (from both OB and AGB stars). More details on the numerics of the radiative feedback (1 and 2) and mechanical feedback (3) are presented in Hopkins et al. (2018a, b). We include metal yields from Type-II SNe, Type-Ia SNe, and AGB winds (see Hopkins et al., 2018c, Appendix A).

The simulations adopt the sub-resolution turbulent metal diffusion and mixing model described in Su et al. (2017) and Escala et al. (2018). We do not model the chemistry appropriate for primordial gas nor the formation and evolution of Pop iii stars, but assume a metallicity floor at $Z=10^{-4}\,Z_{\sun}$.

2.2 Simulation restarts

Fig. 2 shows the stellar surface density maps for the four simulated galaxies in Table 1 at the final redshift of these simulations ($z=5$). Every galaxy contains a number of small, cluster-like objects which survive at least to $z=5$. The central surface density of the clusters in galaxy z5m10a_hr is $\sim 100$–$300\,M_{\sun}\,{\rm pc}^{-2}$, while those in galaxy z5m12b typically have central surface densities over $10^{4}\,M_{\sun}\,{\rm pc}^{-2}$. Galaxy z5m11c_hr shows more smaller clusters than z5m11c owing to its higher mass resolution. Fig. 2 shows the star formation history of the four galaxies. All galaxies exhibit bursty star formation to some extent on time-scales of $\sim 100\,{\rm Myr}$ (see also Ma et al., 2018b, fig. 5). We note that when and where a gas particle turns into a star particle and a SN occurs are sampled stochastically from the SFR and SN rates in our simulations, so running the same simulation twice does not lead to identical results (Hopkins et al., 2018c). The difference between z5m11c and z5m11c_hr in their morphologies and star formation histories is likely due to stochastic effects. They should be regarded as two realizations of the same halo.

Most of the stars and hence clusters in the galaxies are formed during these bursts of star formation. The time-scale of cluster formation is of the same order as the free-fall time of its parent cloud, which is at most a few Myrs for typical densities of GCs (e.g. Kim et al., 2018; Grudić et al., 2018a). The default time interval between consecutive snapshots in our simulations ($\sim 20\,{\rm Myr}$) is too large for studying the formation of GCs, so for each simulation, we restart the run from a snapshot prior to a starburst with more frequent outputs ($\sim 0.5$ Myr between snapshots) until the end of the burst. The grey shaded regions in Fig. 2 illustrate the starbursts we re-simulate for each galaxy. The starting and final cosmic time and redshift for each burst are listed in Table 2. Galaxies z5m11c and z5m11c_hr both exhibit a starburst slightly before $z=5$, so we re-run them for the same redshift interval. We point out that the starburst we select for z5m12a_hr is triggered by a major merger, while those in other simulations are not associated with mergers.

We consider alternative star formation models and parameters when re-running the starbursts to understand the robustness of cluster formation in our simulations. In addition to the star formation criteria (i)–(iii) above, we also consider a fourth criterion, requiring a locally convergent flow:

(iv) Converging flow (cf). Star formation is only allowed in converging flows where $\nabla\cdot{\bf v}<0$.

Each star formation model combines some of the above criteria and certain choice of $\epsilon_{\rm ff}$. In Table 2, we list the star formation models we explore in this paper. The standard FIRE-2 star formation model consists of criteria mol+sg+den with local star formation efficiency $\epsilon_{\rm ff}=1$. The ‘no $n_{\rm th}$’ model only includes criteria mol+sg, i.e. no density threshold for star formation compared to the FIRE-2 model. The ‘G18’ model is adopted in Grudić et al. (2018a), which consists of criteria mol+sg+cf. The ‘G18_e50’ model is identical to ‘G18’ model, except that we adopt a lower local star formation efficiency of $\epsilon_{\rm ff}=0.5$. In Sections 3 and 4, we choose model ‘G18’ as our default model when studying GC formation during these starbursts. Grudić et al. (2019) showed that this model reproduces the observed fraction of stars forming in bound clusters in a suite of high-resolution simulations of individual clouds similar to MW and M51 GMCs. This model also gives good convergence with mass resolution. In Section 5, we will study how different star formation models affect our results.555The virial parameter $\alpha$ defined in Equation 1 usually exhibits some time variability due to random fluctuations in local gas motion. To avoid spurious star formation, Grudić et al. (2018a) also required $\alpha$ to change smoothly with time. This makes a small difference in the cluster mass function shown in Section 10, within stochastic effects (less than a factor of 2). We find that the converging-flow criterion plays a more dominant and systematic effect.

2.3 Cluster Identification

When the re-simulations end, we run the cluster finding algorithm666The code is publicly available at https://github.com/omgspace/Phinder. developed by Grudić et al. (2018b) on all star particles to identify star clusters in the simulated galaxies. The algorithm is summarized as follows. First, we calculate the gravitational potential generated by all star particles. Next, for each particle, we search its 32 nearest neighbors for the particle at the lowest gravitational potential. We move to this particle and do the same search until a local potential minimum is identified. Therefore, every particle is associated with a local potential minimum and all the particles are grouped into particle associations. Lastly, for each group, we re-compute the gravitational potential generated only by the group members and remove those unbound to this group to get a list of bound clusters and their member particles. We also require all bound clusters to have at least 32 particles. The same force softening lengths as in the simulations are used when we calculate the gravitational potential.

In our analysis below, we regard bound clusters with half-mass radius smaller than $R_{\rm half}=100\,{\rm pc}$ as proto-GC candidates. This radius cut is to exclude the galaxies themselves and multiplicities in our analysis and only affects a small number of objects. We only focus on clusters that are newly formed during the starbursts to get a clean comparison between different star formation models (Section 5). As our re-simulations only last for $<100\,{\rm Myr}$, we ignore cluster destruction in our analysis. We expect this only has a small effect on low-mass and under-resolved clusters that may be disrupted faster than this time-scale. All cluster properties are calculated at the end of the re-simulations. We list the total stellar mass formed in each restart and the fraction of stars ending up in bound clusters in Table 2. We refer to bound clusters, YMCs, and proto-GCs interchangeably in the text below.

3 The formation of GC candidates

3.1 GC formation in high-pressure clouds

In this section, we investigate in what conditions and by what physical processes proto-GCs can form. We do so by studying examples of proto-GC candidates formed in the simulation z5m11c_hr_G18, where we re-simulate the starburst in galaxy z5m11c_hr using the star formation model in Grudić et al. (2018b). Approximately $1/4$777For reference, only 4% of the stars formed prior to the burst still remain in bound clusters by $z=5$, probably due to a combinations of reasons below. First, a lower fraction of stars were formed in bound clusters at early time, as much fewer massive clusters can form in low-mass galaxies (cf. z5m10_hr in Fig. 10). Second, clusters that are marginally resolved by $\lesssim 100$ particles will be numerically dissolved over 100 Myr (cf. discussion in Section 6.3). Lastly, cluster disruption is likely efficient because of strong tidal shocks in the progenitor of this galaxy. The oldest cluster found at $z=5$ was formed at $z\sim 10$ (over 600 Myr ago) with metallicity $\rm[Z/H]\sim-2$. of the stars formed in the burst belong to bound clusters with at least 32 particles ($10^{4.5}\,M_{\sun}$) by the end of the re-simulation. Proto-GCs in other simulations are formed in similar ways.

We begin with an overview of the galactic-scale ISM properties during the starburst. In Fig. 3, we show projected gas images at six epochs from different stages of the starburst. Each image represents a physical length of 16 kpc along each dimension. Panel A shows a rapid gas infall from nearly all directions prior to the starburst, building up a large amount of gas in the central region of the galaxy (panel B). At epoch C, a patch of the ISM right next to the center has collapsed to sufficiently high densities to trigger star formation. The galaxy maintains an extremely high gas fraction (close to 80% within the central 5 kpc) for the first $50\,{\rm Myr}$ of the starburst before the SFR decreases dramatically. The gas is highly turbulent with typical velocity dispersion $\sigma_{v}\sim 100\,{\rm km}\,{\rm s}^{-1}$, corresponding to Mach number ${\cal M}\sim 10$–30 for gas of $10^{3}$–$10^{4}$ K. The average gas surface density is over $\Sigma_{\rm gas}\sim 10^{2}\,M_{\sun}\,{\rm pc}^{-2}$, with certain regions in the galaxy reaching $10^{3}$–$10^{4}\,M_{\sun}\,{\rm pc}^{-2}$.

Panel D illustrates the process where feedback from stars that formed earlier starts to launch galactic winds, which sweep up the gas nearby and create shell-like structures in the ISM. Remarkably, the winds compress the gas to high pressure and (surface) density where star formation takes place rapidly, as shown in panels D and E. This is like the ‘positive feedback’ driven by active galactic nucleus (AGN) outflows (e.g. Silk, 2013; Bieri et al., 2015, though the engine here is massive stars rather than AGNs) and the ‘triggered star formation’ model888The ‘triggered star formation’ model is usually discussed in the context of massive star formation within GMCs, but it has also been proposed for cluster formation (e.g. Elmegreen, 2002). (e.g. Elmegreen & Lada, 1977; Palous et al., 1994; Elmegreen et al., 2002; Whitworth & Francis, 2002; Tan, 2000, 2005; Walch, 2014). Finally, feedback from the starburst eventually blows out almost all the gas from the central region and thus star formation is temporarily suppressed in the galaxy (panel F).

Now we study how proto-GCs form in our simulations. Every particle in the simulation has a unique particle ID. By design, each star particle inherits the same ID from its parent gas particle. Thus we are able to trace the particles belonging to a proto-GC identified at the end of the simulation back to the time before they formed. In Fig. 4, we illustrate the formation process of a cluster with mass $M_{\rm cl}=1.4\times 10^{7}\,M_{\sun}$ identified at $z=5$. The formation time ($t=0$) is defined as the epoch when half of the cluster’s member particles have formed (i.e. cosmic time $t=1.121$ Gyr, redshift $z=5.172$ for this cluster). Each column shows a slab ($2\,{\rm kpc}\times 2\,{\rm kpc}\times 500\,{\rm pc}$, the third dimension refers to the direction perpendicular to the page) in the ISM at the given time centered on the median coordinate of all cluster members (and/or their parent gas particles). The color (grey) maps in the first/third (second/fourth) rows represent the projected gas surface density ($\Sigma=\int\rho{\rm d}l$). These images show the gas distribution in the ISM and highlight the dense structures in this region. We checked our results are not sensitive to the thickness of the slab as long as it is larger than the typical size of these dense structures (a few 10 pc to 100 pc). The vector fields in the first/third rows show the direction of gas motion with respect to the center of mass of all cluster members (and/or their progenitors). The blue points in the second/fourth rows show all star particles in the slab, regardless of their formation time. The red points highlight the cluster members (after they turn into stars).

This cluster is formed in a high-pressure cloud compressed by feedback from recent star formation nearby, which is an important channel for forming bound clusters in our simulations. This can be seen from panel D of Fig. 3, where the white dashed box shows a $2\,{\rm kpc}\times 2\,{\rm kpc}$ region centered on this cluster at its formation time. Fig. 4 presents the formation process of this cluster in more detail. Panels A and B show a gas flow starting to sweep through this region from the upper right to the bottom left 5–10 Myr before the cluster formed. In the central region of panel B, we find that the flow has entrained a large amount of dense gas, which will merge shortly with the cloud complex at the center of the frame, as shown by the velocity field. As the gas flow keeps sweeping across and compressing the ISM, a dense cloud has formed by $t=-2\,{\rm Myr}$ at the center of panel C, where stars that eventually end up in the cluster start to form at the density peaks. During the time when the bulk of cluster stars form (panels C–E), the cloud maintains a super high surface density exceeding $10^{4}\,M_{\sun}\,{\rm pc}^{-2}$. At such high surface density, young stars newly born in the cloud cannot provide sufficient pressure support from feedback to overcome the self-gravity of the cloud, so the cloud converts almost all of its mass into stars in a cloud free-fall time. In the meanwhile, gas from the left side of the frame keeps feeding the cloud, so the cloud can survive for a few Myrs. At $t=5\,{\rm Myr}$ (panel E), feedback from the cluster eventually blows out the surrounding gas and terminates cluster formation. These gas flows will likewise drive cluster formation at nearby locations in the same way as described above.

Fig. 5 shows the process by which another cluster formed at redshift $z=5.282$ (cosmic time $t=1.092$ Gyr). The cluster has a final mass $M_{\rm cl}=1.6\times 10^{6}\,M_{\sun}$ by $z=5$. The white box in panel C of Fig. 3 shows where it forms in the galaxy. Again, each panel represents a slab of $2\,{\rm kpc}\times 2\,{\rm kpc}\times 500\,{\rm pc}$ in the ISM. The symbols and colors are identical to those in Fig. 4. We find that at $t=-10\,{\rm Myr}$ (panel A), two cloud complexes (left and right of the image center) are moving toward each other and about to collide. A third complex above the center is also falling to the central region. In panel B, the two clouds at the center have already collided and become a larger complex in which star formation is triggered (see the blue points in the bottom panel). In the meanwhile, the cloud above center keeps moving toward the center and eventually merges with the central cloud complex in panel C. In panel D, the complex breaks into two pieces, as feedback acting at the lower left is pushing them apart (see the velocity fields). The piece at the center of frame reaches a high pressure due to a combination of feedback-driven pressure from the bottom left and cloud collision from the top. The cloud maintains a surface density $\sim 10^{4}\,M_{\sun}\,{\rm pc}^{-2}$ during the time the majority of cluster stars form (panels D and E). Finally, at $t=6\,{\rm Myr}$, the cluster blows out the surrounding gas. High-quality images and animations for Figs. 3–5 are available at this URL999http://www.tapir.caltech.edu/~xchma/HiZFIRE/globular/.

Our simulations suggest a scenario where bound clusters form efficiently in gas-rich, turbulent galaxies with strong star formation activity. The large amount of gas needed can be supplied by rapid gas accretion as commonly happened in the high-redshift universe and/or gas-rich mergers (e.g. Li et al., 2004; Maji et al., 2017). The turbulence is driven by rapid gas inflows (converting gravitational energy to kinetic energy) and feedback processes. In regions where bound cluster form, external pressure provided by feedback-driven winds and cloud/complex/stream collision must be sufficiently high to compress the gas into high-pressure clouds with surface densities $\gg 10^{3}\,M_{\sun}\,{\rm pc}^{-2}$ (see also Elmegreen, 2002). Note that this is much larger than the average gas surface density in the galaxy, indicating that compression by external pressure is important to achieve such conditions. Once the clouds become self-gravitating, they will turn nearly all of their mass into stars in a cloud free-fall time and stars formed in these clouds are gravitationally bound at birth.

In Fig. 6, we further demonstrate that bound clusters preferentially form in high-pressure, high-density regions. For every star particle formed in simulation z5m11c_hr_G18, we find the density of its parent gas particle in the final snapshot prior to its birth (the time separation between two successive snapshots is 0.5 Myr). The solid curve in Fig. 6 shows the probability density function for the ‘birth density’ of all stars formed during the starburst.101010It is highly non-trivial to associate individual star particles with a well-defined ‘birth cloud’ in our simulations, so here we use a local quantity, the birth density of each star particle, rather than cloud-scale properties (which are also non-trivial to define) for diagnostics. We expect the density of star-forming gas in a cloud to correlate with cloud pressure in our simulations, hence stars forming at high densities also means they form in high-pressure (high-surface-density) regions (see also Kruijssen, 2012). The dashed and dotted curves break the distribution function into stars belonging to a bound cluster and stars not bound to any cluster by the end of the simulation. Stars that end up in bound clusters tend to form in one-order-of-magnitude denser gas than other stars. We checked that the birth density distribution of cluster stars does not strongly depend on cluster mass, at least in the range of $M_{\rm cl}\sim 10^{5}$–$10^{7}\,M_{\sun}$. For each cluster, we compute the median birth density of all member stars as a substitute for the average density of its birth cloud. In Fig. 6, the grey histogram shows the distribution of this density for all 277 clusters formed in the re-simulation. This is consistent with the density distribution of all cluster-forming gas (the dashed line). Stars formed at much lower densities (lower-pressure regions) may be unbound in the first place owing to lower cloud-scale star formation efficiencies ($\sim 1$–10% as oppose to near unity in high-pressure clouds) or form marginally bound objects that will be shortly disrupted by feedback. Our findings are in line with previous analytic models that suggest bound clusters form in the most high-pressure, high-density regions in a turbulent ISM (e.g. Elmegreen & Efremov, 1997; Elmegreen, 2002; Elmegreen et al., 2012; Kruijssen, 2012).

The clusters in Figs. 4 and 5 do not morphologically resemble GCs in the Universe when they just formed. In Fig. 7, we show the morphological evolution of these two clusters after their formation. The top (bottom) row shows the cluster in Fig. 4 (5) at $t=2$, 8, and 14 Myr. Each panel shows the projected image of all star particles (blue points) in a $(500\,{\rm pc})^{3}$ box centered on the cluster. The cluster members are shown by the red points. Both clusters are irregular at $t=2\,{\rm Myr}$, but they quickly become spherical over the next 10 Myr. We note that this morphological relaxation process happens faster in our simulations than in reality because a cluster is only resolved by a moderate number of particles (e.g. the cluster of mass $M_{\rm cl}\sim 1.5\times 10^{7}\,M_{\sun}$ in Fig. 4 contains $\sim 2\times 10^{4}$ particles). The half-mass radii of the two clusters are $R_{\rm half}=36.9$ and 23.2 pc, respectively, by $t=14\,{\rm Myr}$. This suggests that clusters formed in high-pressure clouds tend to be gravitationally bound at birth.

3.1.1 Hierarchical cluster formation

In the examples shown in Figs. 4 and 5, almost all stars in the cluster form collectively in a high-pressure gas cloud. Nonetheless, a cloud can form a star cluster complex that may break into several bound clusters at a later time. On the other hand, multiple clusters formed in the same cloud complex can merge together to become a more massive cluster. This is similar to the ‘top-down-and-bottom-up’ hierarchical cluster formation picture discussed in Grudić et al. (2018b, and references therein). We also identify such processes in our cosmological simulations, which we will explicitly show here.

Fig. 9 presents an example of multiple clusters forming ‘top-down’ in one cloud. The left panel is 2 kpc along each side, showing the projected gas image of a 500 pc-thick slab in the galaxy, similar to the top panels in Figs. 4 and 5. The middle panel shows the image zoomed in on the central $500\,{\rm pc}\times 500\,{\rm pc}$. The grey scale shows the gas surface density and the blue points show all the star particles in the $(500\,{\rm pc})^{3}$ box. The central cloud complex is compressed by a gas flow in a similar way as in Fig. 4. The complex forms a YMC where we identify six components labeled by C1–C6 in the middle panel, which become six bound clusters by $z=5$. We define $t=0$ as the formation time of cluster C1. At this time, all six components are close to each other, within a 150 pc-radius sphere from C1. The right panel shows the images of the six clusters 20 Myr later. Each panel is $500\,{\rm pc}\times 500\,{\rm pc}$ centered on the cluster of interest. The blue points show all star particles in the $(500\,{\rm pc})^{3}$ box, with cluster stars highlighted by the red points. Although clusters C1–C4 are still in close proximity to each other (see e.g. the bottom-left panel), the distances between them increase relative to 20 Myr ago, indicating that they are falling apart. Clusters C5 and C6 are even further away from C1. The masses of the six clusters range in $M_{\rm cl}=1.5\times 10^{5}$–$3\times 10^{6}\,M_{\sun}$ (with decreasing mass from C1 to C6).

Fig. 9 provides an example of ‘bottom-up’ cluster formation, where a bound cluster of mass $M_{\rm cl}=6\times 10^{6}\,M_{\sun}$ identified at $z=5$ is formed via multiple cluster mergers. The left panel shows the gas image of a $2\,{\rm kpc}\times 2\,{\rm kpc}\times 500\,{\rm pc}$ slab in the ISM. The middle panel shows the image zoomed in on the central $500\,{\rm pc}\times 500\,{\rm pc}$. The blue points show all star particles in the $(500\,{\rm pc})^{3}$ box, with cluster stars highlighted by red points. The right panel shows the images of this cluster at six subsequent epochs. Each panel is 500 pc on each side. There are at least three clusters seen at $t=9$ and 16 Myr that merge together and form a more massive cluster by $t=36\,{\rm Myr}$. Note that the three progenitor clusters are formed in the same cloud complex compressed by a gas flow, so the final cluster is formed ‘top-down’ first and then ‘bottom-up’.

Although we identify both ‘top-down’ and ‘bottom-up’ cluster formation in our simulations, it is worth noting that these processes are not physically different. Whether a complex forms clusters ‘top-down’ or ‘bottom-up’ is set by the velocity structure of individual cluster-forming clumps in the complex, which determines whether the clusters will fly apart or become bound and merge together. We also stress that ‘top-down’ or ‘bottom-up’ process does not necessarily happen every time a bound cluster forms.

3.2 Mass function of newly formed clusters

In Fig. 10, we present the mass function of bound clusters formed in the starbursts re-simulated with our default ‘G18’ star formation model. For every galaxy, we show the number of clusters per logarithmic mass, ${\rm d}N/{\rm d}\log M_{\rm cl}$, divided by the total stellar mass formed during the starburst. To guide the eye, we also show the power-law mass function, ${\rm d}N/{\rm d}M_{\rm cl}\propto M_{\rm cl}^{-2}$, with the black dashed line.111111We stress that here we study the mass function of bound clusters newly formed within a short time period. We do not attempt to connect this to the mass function of present-day GCs in this paper. Doing so requires evolving the clusters over cosmic time (including internal dynamics), which is not possible at our current resolution but worth future studies (see Section 6.3 for more discussion). This is consistent with the observed YMC mass functions (e.g. Zhang & Fall, 1999; Portegies Zwart et al., 2010, and reference therein) and in agreement with cosmological simulations (e.g. Kravtsov & Gnedin, 2005; Li et al., 2017). A power-law slope of $-2$ is generically arises from scale-free structure formation (e.g. Guszejnov et al., 2018). We only consider clusters with at least 32 particles. The mass function is shown every 0.5 dex in logarithmic mass. We also fit our data by ${\rm d}N/{\rm d}M_{\rm cl}\propto M_{\rm cl}^{-\alpha}$ and label the actual slopes $\alpha$ in Fig. 10.

In all four galaxies, the newly formed bound clusters follow a mass function broadly consistent with ${\rm d}N/{\rm d}M_{\rm cl}\propto M_{\rm cl}^{-2}$. Moreover, we find that all galaxies have similar cluster formation efficiencies, which we refer to the number of bound clusters with a given mass formed per star formation in the galaxy over the same period. This is different from the fraction of stars formed in bound clusters as in Table 2, because the latter does not fully converge with resolution. For example, forming a cluster of mass $M_{\rm cl}\sim 10^{6\pm 0.25}\,M_{\sun}$ requires a total stellar mass $\sim 2\times 10^{7}\,M_{\sun}$ formed in the galaxy. The stellar mass required to form a cluster of mass $M_{\rm cl}$ increases linearly with $M_{\rm cl}$ following the power-law mass function with a slope of $-2$. The variation between different galaxies is within a factor of 2–3. This suggests that our results on the cluster mass function and formation efficiency are not sensitive to the resolution of these simulations, at least when they are run with the same adaptive, resolution-free star formation model. In Section 5, we will study how the details of the star formation prescriptions affect our results.

All the four galaxies are highly gas-rich and turbulent during the starbursts we study here. The bound clusters are formed in similar ways as we present in Section 3.1. We caution that the tentative universal cluster formation efficiency found above may only apply to starburst galaxies with similar gas fraction and the extent of turbulent support to those in our simulations. Galaxies in the early universe preferentially meet such conditions, while low-mass galaxies at later times are likely less efficient in forming bound clusters (e.g. Elmegreen, 2018, and reference therein). We defer a detailed investigation on cluster formation efficiency over a broad range of galaxy mass, redshift, and ISM properties to a future study.

4 Properties of bound clusters

In Section 3.1, we show how bound clusters form in high-pressure, high-surface-density clouds in gas-rich, turbulent ISM. We present the properties of these clusters in this section. We only show results from simulations using the ‘G18’ star formation model. All cluster properties are measured at the end of the re-simulations.

4.1 Age and metallicity distribution

Fig. 11 shows the stacked age distribution of cluster stars in bound clusters formed in simulation z5m11c_hr_G18 for several cluster mass intervals as labeled. The spread of stellar ages increases with cluster mass from $\sim 1\,{\rm Myr}$ for the central 68% stars around the median stellar age in clusters below $M_{\rm cl}\sim 10^{5.5}\,M_{\sun}$ to roughly 3 Myr in clusters more massive than $M_{\rm cl}\sim 10^{6.5}\,M_{\sun}$. This is broadly consistent with the fact that the spread of ages in nearby YMCs is found to be a few Myrs (e.g. Blum et al., 2001; Mac Low & Klessen, 2004; Hollyhead et al., 2015) and with the picture where stars form rapidly at near unity efficiency within a free-fall time of the high-pressure clouds (e.g. Grudić et al., 2018a, and references in Section 1). Note that Fig. 11 essentially shows the SFR of the cluster. The increasing SFR in the early stage of cluster formation means the cluster grows its mass super-linearly, consistent with previous work (e.g. Murray & Chang, 2015; Murray et al., 2017; Li et al., 2018).

Fig. 12 shows the stacked metallicity distribution of cluster stars in simulation z5m11c_hr_G18 for the same cluster mass bins. The stacked distribution shows a $1\sigma$ dispersion of stellar metallicity in $\rm[Z/H]$ of $\sim 0.08\,{\rm dex}$ and does not depend on cluster mass significantly. For individual clusters, the $1\sigma$ dispersion in $\rm[Z/H]$ ranges in 0.03–0.12 dex. This is because most stars in the cluster are born in the same parent cloud with strong turbulence, where metal mixing is expected to be efficient. Our results suggest that cluster stars tend to have near uniform abundance at birth. The observed abundance variations in old GCs may be resulted from stellar evolution effects or subsequent generation of star formation, which still remains an open question (e.g. Bastian & Lardo, 2018, for a recent review). The metallicity of star clusters broadly traces the gas-phase metallicity when and where they formed. We find that the mean metallicity of all bound clusters formed during the starburst in z5m11c_hr_G18 has a median $\rm[Z/H]\sim-0.90$ and $1\sigma$ dispersion 0.12 dex, following the gas-phase mass–metallicity relation at $z=5$ at stellar mass $M_{\ast}\sim 10^{9}\,M_{\sun}$ (e.g. Ma et al., 2016a).

We have also checked z5m10a_hr_G18 and z5m12b_G18. We find that in all simulations, the spread in stellar ages increases with cluster mass and is of order a few Myrs. We caution that the trend where the age spread increases with cluster mass is valid only for clusters formed in the same galaxy. We do not compare the age spread for clusters in different simulations, because they are run at different mass resolution. Moreover, we find the metallicity spread is independent of cluster mass in all simulations and the $1\sigma$ dispersion in $\rm[Z/H]$ is $\sim 0.08\,{\rm dex}$ in all galaxies.

4.2 The size–mass relation

In Fig. 13, we present the size–mass relation for all bound clusters formed during the starburst we re-simulate for each galaxy with the G18 star formation model, where $R_{\rm half}$ is the half-mass radius. The black dotted lines represent constant-density lines on this diagram, $\rho=3M_{\rm cl}/8\pi R_{\rm half}^{3}$, from $0.1\,m_{\rm H}\,{\rm cm}^{-3}$ on the top left (where $m_{\rm H}$ is the mass of a hydrogen atom) to $10^{5}\,m_{\rm H}\,{\rm cm}^{-3}$ on the bottom right. We only show clusters with more than 32 star particles. Most clusters have a half-mass radius $\sim 6$–40 pc. This is systematically larger than present-day YMCs/GCs and those formed in high-resolution cloud-scale simulations (e.g. Grudić et al. 2019). We suspect this is mostly because we have finite mass resolution and force softening lengths in our simulations such that the cuspy density profile in the inner region of the clusters are not well resolved.

We note that the density of clusters differs from simulation to simulation, depending on the galaxy mass and mass resolution. For example, clusters in z5m12b have systematically higher densities than those in other simulations, while clusters in z5m11c have the lowest densities on average. To understand this further, in Fig. 14, we show the distribution function of density for all star-forming gas particles (weighted by SFR) in these galaxies.121212Here we take all gas particles within the galaxy from all snapshots during the starbursts we re-simulate and show their SFR-weighted density distribution. This is technically different from Fig. 6, where we take all star particles formed during the starburst and show the non-weighted density distribution of their progenitor gas particles right before their formation. However, given that star particles are created stochastically from the SFRs of gas particles, the results between Figs. 6 and 14 are conceptually identical. We reiterate that cluster stars tend to form at the high-density end of this distribution (Fig. 6). We find that the trend in the average cluster density among these simulations broadly follow that in the density of star-forming gas. This can be qualitatively understood as cluster density broadly traces the average density of its birth cloud131313The stellar density of the cluster is a factor of a few lower than the average gas density of its birth cloud, because not all gas in the cloud turns into stars and the stars expand in radius when the remaining gas is expelled by feedback (see e.g. panel F in Figs. 4–5)., which correlates with the density of star-forming gas (where the cloud fragments). At the same resolution (e.g. z5m11c vs. z5m12b), more massive galaxies form stars in denser gas and hence clusters are more compact. This is because z5m12b has a higher average ISM density than the other galaxies. Moreover, better resolution tends to restrict star formation to higher gas densities and the clusters are smaller in sizes (z5m11c vs. z5m11c_hr). Since we can resolve fragmentation until the local turbulent Jeans mass $M_{J}\sim\sigma_{v}^{3}/(G^{3/2}\rho^{1/2})$ becomes comparable to the mass of a particle $m_{b}$. Combined with $\sigma_{v}\sim h^{1/2}$ following energy cascade in supersonic turbulence and $m_{b}\sim\rho\,h^{3}$ ($h$ is the resolution scale), this leads to $\rho\sim m_{b}^{-1/2}$. This is why stars tend to form at a factor of $\sim 2\sqrt{2}\approx 2.83$ higher densities in z5m11c_hr than in z5m11c following a factor of 8 difference in mass resolution. We will further discuss the effects of resolution in Section 5.

At the low-mass end, where the clusters are only marginally resolved by $\sim 100$ particles or less, Fig. 13 shows a large scatter in $R_{\rm half}$. This is likely because of artificial N-body relaxation owing to a small number of particles and/or these clusters being disrupted too easily (also artificial due to under-resolved densities). In short, we caution that cluster sizes in our simulations depend on the mass resolution, especially for clusters than only contain a small number particles, whose sizes are heavily affected by numerical effects.

5 Effects of star formation model

In this section, we investigate how the star formation model used in our simulations affects bound cluster formation. We note again the star formation criteria studied here but refer to Section 2 for details: (i) molecular (mol), (ii) self-gravitating (sg), based on the local virial parameter, (iii) a density threshold of $n_{\rm th}=10^{3}\,{\rm cm}^{-3}$ (den), and (iv) converging flow (cf). If some or all of these criteria are met, a gas particle is eligible to convert to a star particle at a rate $\dot{\rho}_{\ast}=\epsilon_{\rm ff}\,f_{\rm mol}\,\rho/t_{\rm ff}$, where $\epsilon_{\rm ff}$ is the local star formation efficiency of the star-forming particle. Our default choice is $\epsilon_{\rm ff}=1$.

We compare four different star formation prescriptions as follows: (1) the FIRE model, which includes criteria mol+sg+den, (2) the ‘no $n_{\rm th}$’ model, which only includes criteria mol+sg, (3) the ‘G18’ model, which consists of criteria mol+sg+cf, and (4) the ‘G18_e50’ model, which is the same as the ‘G18’ model, except that we adopt $\epsilon_{\rm ff}=0.5$ (50%). Models (1)–(3) adopt the default $\epsilon_{\rm ff}=1$. Model (2, ‘no $n_{\rm th}$’) tends to adaptively pick up overdense structures we are able to resolve for star formation. Model (1, ‘FIRE’) adds a density threshold to model (2), but the effects of the threshold are not resolution independent (see Hopkins et al., 2013b, and discussion below). In practice, a higher $n_{\rm th}$ is usually applied at ultra-high resolution (e.g. Wheeler et al., 2018, $m_{b}\sim 30\,M_{\sun}$). Model (3, ‘G18’) is also stricter than model (2), but in contrast to model (1), it is also adaptive to resolution. Model (4, ‘G18_e50’) forms stars at a factor of 2 lower rate at a given density compared to model (3), which means star formation is delayed. We choose model (3) as our default model in Sections 3 and 4 as it reproduces the fraction of stars formed in bound clusters as observed in a suite of high-resolution star-forming cloud simulations similar to MW and M51 GMCs. It also converges fast with resolution (see Grudić et al. 2019, in preparation).

The first diagnostic we analyze here is the probability density function of gas density for star-forming gas141414Here we compare simulations of the same starburst in galaxy z5m11c_hr re-run from identical initial condition, so the difference is entirely due to the star formation model, not resolution or the mean ISM density (see discussion above in Section 4.2). (weighted by the SFR of each gas particle). We re-run the starburst in galaxy z5m11c_hr using models (1)–(4) and present the results in Fig. 15. As we have discussed above, model (2) is the least strict, where stars may form in gas below $100\,{\rm cm}^{-3}$ and the peak density for star formation is at $\sim 10^{3}\,{\rm cm}^{-3}$. Model (3) is stricter than model (2), which pushes star formation to a factor of 3 higher densities. In model (4), stars form at lower rates while gas keeps collapsing and thus star formation is further restricted to higher densities.151515We have also tested $\epsilon_{\rm ff}=0.1$, where star formation occurs at even higher densities. The simulation then becomes too expensive to finish, since there are a large number of gas particles at $n\gg 10^{5}\,{\rm cm}^{3}$ and the timestep is very small for these particles. Model (1) is also stricter than model (2), since gas below $10^{3}\,{\rm cm}^{-3}$ cannot form stars until it collapses to reach the density threshold. Compared to model (3), model (1) makes it harder for star formation in low-density gas and easier at higher densities.

In Section 3.1, we find that bound clusters preferentially form in high-density, high-pressure gas. Stars formed in low-density gas are likely born unbound in the beginning or in clusters that will be disrupted by feedback shortly. Intuitively, we would expect cluster formation is more (less) efficient when the star formation model is stricter (looser). This is confirmed by the normalized cluster mass functions as shown in Fig. 16. Model (2) produces the least clusters at all masses (per stellar mass formed). Models (1) and (3) produce more clusters than model (2) at all masses. Model (4) produces the most clusters at nearly all masses. Compared to model (3), model (1) has more clusters below $\sim 10^{5}\,M_{\sun}$, but less clusters at higher masses. All of these models produce a cluster mass function consistent with ${\rm d}N/{\rm d}M_{\rm cl}\propto M_{\rm cl}^{-2}$. We also note again that stochastic effects in the simulations, such as when and where a star particle can form and a SN can occur (stochastically sampled from the rates), can also generate random variations in the (normalized) number of clusters formed in each mass bin. Nevertheless, we find these two findings are robust to our star formation model: (a) bound clusters form efficiently in high-density, high-pressure clouds in a gas-rich, turbulent galaxy and (b) the newly formed clusters broadly follow a power-law mass function with a slope of $-2$. The fraction of stars formed in bound clusters do, however, depend on the details in the star formation criteria. Models (2), (3), and (4) produce 17%, 26%, and 39% of the stars in bound clusters during the starburst in galaxy z5m11c_hr (see Table 2). As expected, this fraction increases when the star formation criteria become stricter, as star formation is restricted to higher density gas.

Fig. 16 shows that a lower $\epsilon_{\rm ff}$ actually increases the cluster formation efficiency. This might seem counterintuitive at first glance. We emphasize that $\epsilon_{\rm ff}$ here is the local star formation efficiency in the densest gas following the fragmentation of star-forming clouds down to the smallest resolvable scales in the simulations. Using a lower $\epsilon_{\rm ff}$ delays star formation in a cloud, so the cloud can collapse further to higher densities. As the surface densities of star-forming clumps become higher, it effectively enhances the cloud-scale star formation efficiencies, resulting in a larger fraction of stars formed in bound clusters. We point out that this trend is opposite to that in Li et al. (2018), where the authors find that lowering $\epsilon_{\rm ff}$ reduces the cloud-scale star formation efficiencies and thereby the bound fraction. Admittedly, the simulations in Li et al. (2018) use very different star formation and feedback models from ours and the physical meanings of $\epsilon_{\rm ff}$ are likely different between the two studies. Tracking down the exact cause of this discrepancy is beyond the scope of this paper, but it is worth further investigations.

Models (2)–(4) are all adaptive to resolution, while model (1) contains a density criterion which is not necessarily independent of resolution. In Fig. 14, we can see that using the same adaptive star formation model (G18), the density distribution of star-forming gas depends both on galaxy mass (probably the average ISM density if more physical) and mass resolution. Galaxies z5m11c and z5m12b are run at the same resolution ($m_{b}\sim 7000\,M_{\sun}$), but star formation happens in denser gas in z5m12b, likely due to a higher mean ISM density in more massive systems. Galaxies z5m11c and z5m11c_hr have similar masses (they should be treated as different realizations of the same initial condition) but are run using different resolution ($m_{b}\sim 900$ and $\sim 7000\,M_{\sun}$). Star formation in z5m11c_hr tends to happen in higher-density gas, because we can resolve turbulence on smaller scales than in z5m11c. Therefore, adding a density threshold may have a simulation-by-simulation effect on star/cluster formation. We do not use model (1) as our default in this paper.

6 Discussion

We have shown that our simulations are able to form bound clusters self-consistently in $z\geq 5$ galaxies over a broad range of masses. We show that in a highly gas-rich, turbulent ISM, bound clusters preferentially form in high-pressure regions compressed by feedback-driven winds and collisions between gas clouds/streams. The newly formed clusters broadly follow a power-law mass function with a slope of $-2$ (equal mass per decade in cluster mass). Our results are in line with previous analytic models and numerical experiments, but obtained in cosmological simulations self-consistently for the first time. The clusters formed in our simulations are small in size ($\sim 6$–40 pc in $R_{\rm half}$, although not as small as real GCs), with cluster stars typically spanning a few Myrs in stellar age and $<0.2$ dex in metallicity ($\sim 0.08$ dex dispersion in $\rm[Z/H]$), in broad agreement with observed cluster properties. Our findings support the idea that present-day GCs were star clusters/YMCs formed in high-redshift galaxies. We restrict the scope of this paper primarily to how bound clusters form, but briefly discuss some other related topics below.

6.1 Cluster formation across cosmic time

We have focused on a sample of cosmological zoom-in simulations run to $z=5$. However, the key processes we identify for cluster formation are not exclusive to the early universe. As long as the galaxy is gas-rich and turbulent, there is always a large probability to form high-pressure regions due to gas stream collision and/or compression by feedback-driven winds. Such conditions are ubiquitous in high-redshift galaxies, but it may also happen at later times, for example, in gas-rich mergers161616Mergers are also common at high redshifts (e.g. Kim et al. 2018 studied a proto-GC formation in a galaxy merger at $z\sim 7$), but they are not required to bring in a large amount of gas and stir turbulence. Normal gas infall is sufficient to create a gas-rich, turbulent environment. and starburst galaxies (e.g. Portegies Zwart et al., 2010; Oklopčić et al., 2017). In the left panel of Fig. 17, we show the projected stellar image of simulation m12i from Wetzel et al. (2016) (also see e.g. Hopkins et al., 2018c) at $z=1$, when the galaxy has just undergone a starburst. The halo mass and stellar mass at this epoch are $M_{\rm halo}\sim 10^{12}\,M_{\sun}$ and $M_{\ast}\sim 2\times 10^{10}\,M_{\sun}$, respectively, comparable to galaxy z5m12b at $z=5$. This galaxy becomes a MW-mass disk galaxy by $z=0$. This simulation uses the FIRE-2 model for star formation (i.e. model 1 in Section 5) and has the same mass resolution as z5m12b. We identify a number of bound clusters formed in this galaxy during the starburst at $z\sim 1$ in similar ways as we present in Section 3.1.

On the other hand, bound clusters are indeed more difficult to form at lower redshifts. In the right panel of Fig. 17, we present the projected stellar image of another galaxy, m11q, at $z\sim 1$, from the FIRE simulation suite (see table 1 in Hopkins et al., 2018c). Its halo mass at this epoch is $M_{\rm halo}\sim 10^{11}\,M_{\sun}$, similar to that of z5m11c_hr at $z=5$, but it has a factor of 2 lower stellar mass. It is also run with the FIRE star formation model and has the same mass resolution as z5m11c_hr ($m_{b}\sim 900\,M_{\sun}$). This galaxy has also just undergone a starburst at this epoch. However, we do not find any bound cluster with at least 32 particles formed during the starburst in this galaxy. The peak SFR during the burst is $\sim 1\,M_{\sun}\,{\rm yr}^{-1}$, almost a factor of 10 lower than that in z5m11c_hr (cf. Fig. 2). Galaxy m11q also has a much lower gas fraction and gas surface density than z5m11c_hr at the peak of the starburst.

These results suggest that even at $z\sim 1$, cluster formation in low-mass galaxies is less efficient than in the early universe. Below $z\sim 1$, it is expected that the gas fraction in galaxies decreases with time and the ISM becomes less turbulent (see e.g. Faucher-Giguère et al., 2013; Hayward & Hopkins, 2017). At the same time, MW-mass galaxies can maintain a stable gas disk (e.g. Ma et al., 2017a, b) in which GMCs are formed mainly via gravitational instability (Toomre, 1964). This is a different mode of star formation from that presented in Section 3.1. All of these effects make bound clusters more difficult to form at later times, consistent with the fact that YMCs tend to form in extreme environments at the present day (e.g. Portegies Zwart et al., 2010). This is in line with the scenario where present-day GCs preferentially formed at high redshifts, so they tend to be old and metal-poor.

6.2 How does feedback affect cluster formation?

We do not explicitly study the effects of feedback strength (per stellar mass) on cluster formation in this paper, but briefly discuss them here. It has been suggested that the cluster formation efficiency will decrease when feedback becomes stronger (e.g. Li et al., 2018). This is plausible since stronger feedback may reduce the star formation efficiency in individual clouds by blowing away more gas and unbinding the stellar association formed in the cloud.

A higher feedback strength might thus limit bound clusters to form only in the highest-density clouds. On the other hand, it can drive turbulence and generate fast gas flows more efficiently. Such processes may expedite forming high-pressure clouds (via compression) and actually enhance cluster formation in the galaxy. The two effects outlined above are competing with each other and it is not obvious what the net effect is. In our simulations, all feedback quantities (e.g. luminosities, SNe and mass loss rates, etc.) are IMF-averaged and directly computed from stellar population synthesis models. Artificially boosting or decreasing feedback strengths will result in violations of certain observational constraints such as the Kennicutt–Schmidt relation, at least in disk galaxies at $z\sim 0$ (see fig. 35 in Hopkins et al., 2018c). We do not expect the input feedback strength from stars to change by a large factor, but there are certain effects that may enhance the feedback strength per stellar mass, including a top-heavy IMF (e.g. Marks et al., 2012), stellar rotation and binarity (e.g. Ma et al., 2016b), and radiation pressure due to Lyman-$\alpha$ resonance scattering (e.g. Kimm et al., 2018). The impact of these on cluster formation is beyond the scope of the current paper but worth studying in a future work.

6.3 Future directions

Our high-resolution cosmological zoom-in simulations are able to form bound clusters self-consistently and systematically. However, studying the dynamic evolution of these clusters on cosmological time-scales over 10 Gyr is a more challenging problem. First of all, at finite mass resolution, clusters in the simulations are resolved by a much smaller number of particles, where $N$-body relaxation happens faster than it should. A cluster may lose a large fraction of its mass or even dissolve on time-scales comparable to 10 relaxation times, which is of order 50 Myr for a cluster with 100 particles and density $10^{2}$–$10^{3}\,m_{\rm H}\,{\rm cm}^{-3}$. Secondly, using a finite force softening length, we cannot resolve the central density profile of these clusters, so they can be tidally disrupted more easily in the simulations (see also e.g. van den Bosch et al., 2018). Moreover, as all particles have nearly equal mass in the simulations, some internal dynamic processes due to unequal mass in $N$-body interactions, such as mass segregation, cannot be captured. To reliably trace the dynamic evolution, disruption, and survival of bound clusters over cosmic time in cosmological simulations, it might be necessary to use a hybrid approach where the formation of clusters is explicitly resolved but the dynamics of these clusters are followed by tracer particles (e.g. Li et al., 2017). Encouraged by previous studies using this method, we expect some clusters formed in the inner region of high-redshift galaxies to be kicked to large radii and survive to $z=0$ (e.g. Kruijssen et al., 2019; Li & Gnedin, 2019).

We find that the cluster formation efficiency, or the fraction of stars formed in bound clusters of a given mass, in these simulations depends on the star formation criteria. Our default choice (the G18 model) is motivated by high-resolution cloud-scale simulations, but it is not guaranteed that this model is still the ‘best’ application for cosmological simulations. Note that the models studied in Section 5 do not produce a big difference in global galaxy properties, such as stellar mass, SFR, and star formation history, since they are mainly regulated by feedback (e.g. Hopkins et al., 2013b). The formation of bound clusters, however, may serve as a test bed for star formation models in cosmological simulations (see also e.g. Li et al., 2018).

Last but not least, the results from our simulations are likewise useful for inspiring and calibrating sub-grid cluster formation models for low-resolution cosmological simulations and semi-analytic models (e.g. Choksi et al., 2018; El-Badry et al., 2019). This is worth studying in the future.

7 Conclusion

In this paper, we present the self-consistent formation of bound star clusters in a sample of high-resolution cosmological zoom-in simulations of $z\geq 5$ galaxies. For each galaxy, we focus on a starburst which we re-simulate with different star formation models. Our default simulations adopt the FIRE-2 stellar feedback model in Hopkins et al. (2018c) and a star formation prescription motivated by high-resolution cloud-scale simulations from Grudić et al. (2018a). We identify the key processes that result in the formation of bound star clusters and study the mass function and basic properties of the bound clusters formed during the starbursts. By comparing simulations of the same starburst with different star formation models, we are able to test the robustness of these results to our star formation model. Our main findings are the following:

(i) We find that bound star clusters preferentially form in high-pressure clouds with gas surface density above $10^{4}\,M_{\sun}\,{\rm pc}^{-2}$, where the cloud-scale star formation efficiencies are nearly unity and the clusters tend to be gravitationally bound at birth (Section 3.1, Figs. 4–7). These high-pressure clouds are formed due to compression by feedback-driven winds and collisions of smaller clouds/gas streams in highly gas-rich, turbulent environments (Figs. 3–5).

(ii) In some cases, clusters may form hierarchically. Multiple clusters may form in one cloud complex and then fall apart (Fig. 9, ‘top-down’). On the other hand, several small clusters born in the same cloud complex may merge to become a more massive cluster (Fig. 9, ‘bottom-up’).

(iii) The newly formed bound clusters broadly follow a power-law mass function of ${\rm d}N/{\rm d}M_{\ast}\propto M_{\rm cl}^{-2}$. In our default simulations, the efficiency of forming a cluster of a given mass is comparable in galaxies in a broad range of mass: one cluster of mass $M_{\rm cl}$ ($\pm 0.25$ dex) forms out of every $20\,M_{\rm cl}$ of stellar mass formed in the galaxy (Section 3.2, Fig. 10). About 25% of the stars form in bound clusters (Table 2). These results only refer to clusters formed during the starburst we re-simulate for each galaxy.

(iv) The age spread of stars in a cluster is typically a few Myrs and increases with cluster mass (Section 4.1, Fig. 11). The metallicity dispersion of cluster stars is $\sim 0.08$ dex and does not depend strongly on cluster mass (Fig. 12). The mean metallicity of the clusters broadly traces the gas-phase metallicity of the galaxy when and where the clusters form.

(v) The bound clusters in our simulations typically have half-mass radii in $\sim 6$–40 pc (Section 4.2, Fig. 13), which are considerably larger than real GCs. This means our simulations cannot reliably trace the dynamical evolution of clusters across cosmic time to study present-day GC populations.

(vi) The cluster formation efficiency relies on the star formation criteria used in the simulations. In general, simulations with a stricter (looser) star formation model form more (less) bound clusters per stellar mass formed, though we only find a factor of a few variation for the star formation criteria we consider. A lower local star formation efficiency results in a larger fraction of stars formed in bound clusters (Section 5, Fig. 16). However, the mass function of newly formed clusters follows a power-law function of slope $-2$ regardless of star formation model.

(vii) Bound clusters can form at lower redshifts, but they are more difficult to form efficiently due to lower gas fraction, turbulent support, and surface density, especially in low-mass systems (Section 6.1, Fig. 17). This is in line with the scenario that present-day GCs formed at high redshifts, so they tend to be old and metal-poor.

(viii) The fraction of stars formed in bound clusters and their sizes in our simulations are not fully converged. However, the fact that bound clusters preferentially form in high-pressure regions, the power-law mass function with a slope of $-2$ for newly formed clusters, and the small age and metallicity spread of cluster stars are all robust.

Our findings support the picture that present-day GCs were in fact regular YMCs formed in high-redshift galaxies during normal star formation activity.

Acknowledgement

We acknowledge our referee, Oleg Gnedin, for a constructive report that helps us correct a few errors and inaccurate statements and improve the clarity of this paper. We also thank many colleagues, especially Paul Shapiro and Hui Li, for sending us helpful comments and discussions after this paper appeared on arXiv. The simulations used in this paper were run on XSEDE computational resources (allocations TG-AST120025, TG-AST130039, TG-AST140023, TG-AST140064, and TG-AST190028). This work was supported in part by a Simons Investigator Award from the Simons Foundation (EQ) and by NSF grant AST-1715070. EQ thanks the Princeton Astrophysical Sciences department and the theoretical astrophysics group and Moore Distinguished Scholar program at Caltech for their hospitality and support. MG and PFH was supported by an Alfred P. Sloan Research Fellowship, NASA ATP Grant NNX14AH35G, and NSF Collaborative Research Grant #1411920 and CAREER grant #1455342. CAFG was supported by NSF through grants AST-1517491, AST-1715216, and CAREER award AST-1652522, by NASA through grant 17-ATP17-0067, by STScI through grant HST-GO-14681.011, and by a Cottrell Scholar Award from the Research Corporation for Science Advancement. MBK acknowledges support from NSF grant AST-1517226 and CAREER grant AST-1752913 and from NASA grants NNX17AG29G and HST-AR-14282, HST-AR-14554, HST-AR-15006, and HST-GO-14191 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS5-26555. AW was supported by NASA, through ATP grant 80NSSC18K1097 and HST grants GO-14734 and AR-15057 from STScI. DK was supported by NSF grant AST-1715101 and the Cottrell Scholar Award from the Research Corporation for Science Advancement.