https://orcid.org/0000-0003-4394-6085
ABSTRACT
We perform cosmological hydrodynamical simulations to study the formation of proto-globular cluster candidates in progenitors of present-day dwarf galaxies |$(M_{\rm vir} \approx 10^{10}\, {\rm M}_\odot$| at z = 0) as part of the ‘Feedback in Realistic Environment’ (FIRE) project. Compact (r1/2 < 30 pc), relatively massive (0.5 × 105 ≲ M⋆/M⊙ ≲ 5 × 105), self-bound stellar clusters form at 11 ≳ z ≳ 5 in progenitors with |$M_{\rm vir} \approx 10^9\, \mathrm{M}_{\odot }$|. Cluster formation is triggered when at least |$10^7\, \mathrm{M}_{\odot }$| of dense, turbulent gas reaches |$\Sigma _{\rm gas} \approx 10^4\, {\rm M}_\odot \, {\rm pc}^{-2}$| as a result of the compressive effects of supernova feedback or from cloud–cloud collisions. The clusters can survive for |$2-3\, {\rm Gyr}$|; absent numerical effects, they could possibly survive substantially longer, perhaps to z = 0. The longest lived clusters are those that form at significant distance – several hundreds of pc – from their host galaxy. We therefore predict that globular clusters forming in progenitors of present-day dwarf galaxies will be offset from any pre-existing stars within their host dark matter haloes as opposed to deeply embedded within a well-defined galaxy. Properties of the nascent clusters are consistent with observations of some of the faintest and most compact high-redshift sources in Hubble Space Telescope lensing fields and are at the edge of what will be detectable as point sources in deep imaging of non-lensed fields with JWST. By contrast, the star clusters’ host galaxies will remain undetectable.
1 INTRODUCTION
It has been over 200 years since William Herschel declared that globular clusters (GCs) ‘are generally but little known and are undoubtedly the most interesting objects in the heavens’ (Herschel 1814), and while GCs have been the subject of intense and detailed study, many aspects of their formation and evolution remain but little known. While GCs are ubiquitous in massive (|$L \gtrsim 0.1\, L^{\star }$|) galaxies at z = 0, exactly how and when they typically form are topics of considerable debate. Some observations of both individual GCs and GC systems around galaxies are naturally reproduced if metal-poor GC formation is connected to specific conditions present only in the high-redshift Universe (Peebles & Dicke 1968; Moore et al. 2006). On the other hand, the existence of metal-rich GCs and observations of dense, massive star clusters forming in extreme settings such as mergers in the low-redshift Universe point to a connection between GC formation and the high pressure, high surface density tail of the distribution of star-forming gas in galactic disks (Ashman & Zepf 1992; Elmegreen & Efremov 1997).
Broadly speaking, these can be thought of as pre-galactic and galactic models, respectively. Although both classes of model are capable of reproducing many bulk properties of the GC population at z = 0, they differ in the typical epoch of cluster formation – in or near the epoch of reionization, z ∼ 6–10, for pre-galactic models, and near, but prior to, the peak of the cosmic star formation history (SFH) at z ∼ 2–3 for the galactic models. As a result, the abundance and properties of GCs forming in the reionization era have the potential to definitively discriminate between formation models.
It is therefore momentous that we appear to be on the cusp of directly observing the formation of GCs in the high-redshift Universe. Recent results indicate that many high-z ‘galaxies’ have properties similar to nascent star clusters (or star cluster complexes; Ishigaki et al. 2018; Bouwens et al. 2021, 2022), and individual systems magnified by gravitational lensing have revealed clear candidates for GC-like objects in formation at z ∼ 3 − 6 (Johnson et al. 2017; Vanzella et al. 2017, 2019, 2021a). JWSTwill likely unveil huge numbers of GCs in formation (Carlberg 2002; Boylan-Kolchin 2017, 2018; Renzini 2017; Pozzetti, Maraston & Renzini 2019), thereby providing a wealth of information about the ‘how, when, and where’ of GC formation.
Recent progress in modelling the formation of GCs has also been substantial, with a variety of approaches providing frameworks for understanding the formation of GCs both within their host galaxies and in the broader context of galaxy formation theory. Ongoing work on this front includes (1) simulations – at very high resolution but without a full cosmological context – of cluster formation within galaxies or molecular cloud complexes (He, Ricotti & Geen 2019; Li et al. 2019, 2022; Lahén et al. 2020; Lee, Shin & Kim 2021; Hislop et al. 2022; Lahén, Naab & Kauffmann 2022); (2) numerical or seminumerical models of GC formation (which track GCs in cosmological context without directly resolving their formation; e.g. Katz & Ricotti 2014; Ricotti, Parry & Gnedin 2016; Li et al. 2017; Renaud, Agertz & Gieles 2017; Pfeffer et al. 2018; Creasey et al. 2019; El-Badry et al. 2019; Carlberg 2020; Halbesma et al. 2020; Phipps et al. 2020; Reina-Campos et al. 2022); (3) simulations linking GC formation to specific conditions in the high-redshift Universe (e.g. Mandelker et al. 2018; Madau et al. 2020; Lake et al. 2021); and (4) empirical models connecting GCs to high-redshift haloes (e.g. Trenti, Padoan & Jimenez 2015; Boylan-Kolchin 2017; Valenzuela et al. 2021).
Given the inexorable increase in computing power, it is also now possible to directly resolve GC formation in cosmological simulations of galaxy formation (Kimm et al. 2016; Kim et al. 2018; Ma et al. 2020). Most of these focus on galaxies that are fairly massive relative to a typical galaxy [i.e. M⋆(z)], in large part because the ubiquity of GCs in |$L \gtrsim 0.1\, L^{\star }$| galaxies at z = 0 guarantees the conditions for GC formation are universally met at some time in such objects. The presence of GCs in many nearby dwarf (|$M_{\star } \lesssim 3\times 10^9\, {\rm M}_{\odot }$|) galaxies may be an important clue to the origin of GCs more generally.
GCs are present even in galaxies as faint as |$M_{\star } \approx 10^5\, {\rm M}_{\odot }$| (Eridanus II; Bechtol et al. 2015; Koposov et al. 2015; Crnojević et al. 2016; Simon et al. 2021), a regime in which standard models of galaxy formation predict that the majority of star formation should have occurred by the end of the reionization era, z ∼ 6 (Bullock, Kravtsov & Weinberg 2000; Benson et al. 2002; Somerville 2002; Ricotti & Gnedin 2005). The fraction of stars contained in GCs in low-mass (|$M_\star \lesssim 10^{7}\, \mathrm{M}_{\odot }$|) galaxies is |$1\!-\!10~{{\ \rm per\ cent}}$| (Georgiev et al. 2010; Hudson, Harris & Harris 2014; Larsen 2017), which is much higher than in more massive systems and indicates that cluster-related star formation plays an important role in the growth of these systems; this is especially true at early times, as the clusters can contain |${\sim }25~{{\ \rm per\ cent}}$| of the metal-poor ([Fe/H] < −2) stars in dwarf galaxies (Larsen et al. 2014). However, the fraction of galaxies hosting at least one GC drops off strongly towards low stellar masses (Georgiev et al. 2010; Burkert & Forbes 2020; Eadie, Harris & Springford 2022), meaning that conditions for GC formation are met in only a subset of dwarf galaxies; an understanding of which dwarf galaxy progenitors achieve the necessary conditions for GC formation may prove essential for understanding GC formation more broadly.
In this paper, we perform a series of cosmological zoom simulations to study the formation of bound stellar clusters in the progenitors of present-day dwarf galaxies (|$M_{\rm halo}(z=0)\approx 10^{10}\, {\rm M}_{\odot }$|). We refer to the star clusters of interest as proto-globular cluster candidates (GCCs). In these proof-of-principle simulations, we focus on the conditions that lead to the formation of dense star clusters at high redshift in such galaxies and explore when this cluster formation occurs, the masses and lifetimes of these clusters, and the properties of the clusters relative to their host galaxies (e.g. the mass of the cluster relative to the host galaxy and the location of formation of the clusters with respect to the size of the galaxies). This paper is organized as follows: in Section 2, we discuss our simulation set-up. Section 3 describes the formation and properties of the clusters and connections to their larger scale environments. Section 4 discusses implications for the detectability of clusters in situ in the high-redshift Universe and for GC formation models as well as the sensitivity of our results to variations in the treatment of galaxy formation physics and numerics in the simulations. In Section 5, we present our conclusions. Unless otherwise noted, all lengths quoted in this paper are physical, not comoving.
2 SIMULATIONS
Our simulation suite is part of the ‘Feedback In Realistic Environment’ project (FIRE, Hopkins et al. 2014, 2018b).1 We select seven realizations of haloes with virial masses of |${\sim }10^{10}\, \text{M}_\odot$| at z = 0 – hosts of present-day dwarf galaxies – from Fitts et al. (2017) and re-simulate them with an updated version of the gizmo2 code (Hopkins et al. 2014) and FIRE-2 galaxy formation prescriptions (Hopkins et al. 2018b). These simulated haloes comprise the six galaxies with the highest z = 0 stellar mass in the Fitts et al. (2017) suite (m10h, m10i, m10j, m10k, m10l, m10m) along with the one of the lowest M⋆(z = 0) galaxies (m10b). The gravity solver in gizmo is a descendant of gadget3 (first described in Springel et al. 2008) and the hydrodynamical equations are treated via the mesh-free finite mass (MFM) Lagrangian Godunov method, which provides adaptive spatial resolution while maintaining conservation of mass, energy, and momentum. We adopt a flat dark energy + cold dark matter (ΛCDM) cosmology with h = 0.71, Ωm = 0.266 = 1 − ΩΛ, and Ωb = 0.0449, all consistent with 7-yr data from WMAP (Komatsu et al. 2011). These differ slightly from the latest constraints based on full-mission Planck observations (Planck Collaboration VI 2018), but these differences are unimportant for the topics considered here.
The initial conditions (ICs)3 are generated using the zoom technique (Katz & White 1993; Oñorbe et al. 2014) at z = 125, embedded within periodic cosmological boxes of |$L=25\, {\rm Mpc}/h$|, and are computed via the code music (Hahn & Abel 2011). The ICs have particle masses of |$m_{\rm gas}=500\, {\rm M}_\odot$| and |$m_{\rm dm}=2500\, {\rm M}_\odot$|, and we perform the simulations with the Plummer-equivalent adaptive gravitational softening lengths of ϵgas = 2 pc, ϵ⋆ = 3.7 pc, and ϵdm = 35 pc (physical). For each simulation, we output snapshots every 10 Myr – comparable to the expected formation period of massive star clusters (Bastian, Hollyhead & Cabrera-Ziri 2014; Hollyhead et al. 2015) – over the period z = 15–4 (yielding 134 snapshots in this redshift range) to enable detailed study of the formation of any bound stellar clusters within this redshift range. For z < 4, we output snapshots every 250 Myr, resulting in 180 snapshots in total. We identify and characterize haloes and subhaloes using a modified version of the code rockstar (Behroozi, Wechsler & Wu 2013; Wetzel & Garrison-Kimmel 2020a, b).
Our simulations explicitly resolve the multiphase interstellar medium (ISM), with heating and cooling of gas modelled in the temperature range of T = 10–1010 K. Star formation happens in self-shielding, self-gravitating, Jeans-unstable gas clouds above a density threshold of |$n_{\rm thresh}=1000\, {\rm cm}^{-3}$|. The star formation efficiency per free-fall time of such gas is 100 per cent, though, as described in Hopkins et al. (2018b), this does not imply a large global star formation efficiency. Our stellar feedback model includes SNe Ia and II, multiwavelength photoheating, cosmic ray heating, stellar winds, and radiation pressure, all adopted from starburst99 stellar evolutionary model (Leitherer et al. 1999) assuming a Kroupa initial mass function (Kroupa 2002). The ultraviolet background radiation model is an updated version of the original model presented in Faucher-Giguère et al. (2009)4 and completes the reionization around z = 6. Relative to earlier FIRE simulations, including the Fitts et al. (2017) suite, our fiducial FIRE-2 model includes metal-diffusion physics, a slightly updated ultraviolet background radiation model, and a correction to the cosmic ray heating source term that avoids spurious heating at very early times in the simulations. We also adopt an updated stellar mass-loss algorithm (as we discuss further in Section 4.2) that differs from the default FIRE-2 implementation outlined in Hopkins et al. (2018b) but follows the default treatment in FIRE-3 (Hopkins et al. 2023). In a second set of simulations, we use the default FIRE-2 algorithm for the stellar mass-loss processes (as described in Hopkins et al. 2018b and in Section 4.2) or vary the density threshold from our standard choice of |$n_{\rm thresh}=1000\, {\rm cm}^{-3}$|, keeping all other fiducial FIRE-2 physics unchanged, to study the impact on star formation activity and any potential difference in efficiency of galaxies in forming GCCs.
We identify star clusters using the Phinder algorithm first described in Grudić et al. (2018b). phinder searches for local minima of the stellar gravitational potential and identifies bound particles within each group. We only keep clusters with at least 32 bound members, resulting in a minimum cluster mass of |$M_{\rm cl}\approx 1.6\times 10^4\, {\rm M}_\odot$|, in our cluster catalogues. However, even this conservative choice typically results in bound objects with very large half-mass radii R1/2 > 100 pc. These objects get destroyed very quickly in the tidal forces of the ISM, meaning they are not good candidates for long-lived and self-bound star clusters. In our analysis, we only consider bound clusters with initial stellar mass of |$M_{\rm cl}=100\, m_\star \approx 5\times 10^4\, {\rm M}_\odot$| and 3D half-mass radii of r1/2 < 50 pc. This choice results in the selection of long-lived, bound stellar clusters, as we show below. Given that our chosen force softenings for baryonic particles are comparable to half-light radii of the most compact GCs, we do not expect to resolve the true internal structure of our bound clusters in this proof-of-concept work. In a future paper, we will present more detailed results on GCC formation using updated FIRE3 physics and smaller force softenings; our preliminary analysis shows that these GCCs have sizes that are significantly smaller (∼5–7 pc) than those presented here. Since our simulations employ softened rather than direct gravitational force calculations, the internal dynamics of the clusters could not be accurately tracked over cosmological times even had they formed with the smaller sizes of present-day GCs. The metallicity of each cluster is computed from the metallicity of its star particles, each of which is modelled as a single stellar population with known age, mass, and metallicity that is inherited from its parent gas particle.
Our analysis results in several GCCs across our simulation suite. Most form within the most massive halo in the simulation volume; the exceptions are the m10k cluster and the second cluster in the m10i realization. Though many lower mass cluster candidates are identified using this procedure, we focus only on clusters with |$M_{\star }\gt 5\times 10^{4}\, \mathrm{M}_{\odot }$| (>100 star particles at birth) for the remainder of this paper. Of the seven haloes from the Fitts et al. (2017) suite that we have resimulated, six contain such a star cluster, and each of these haloes forms only a single cluster with |$M_{\star } \gt 10^5\, \mathrm{M}_{\odot }$|. The lone exception, halo m10b, has a central galaxy with a much lower stellar mass at z = 0 based on the Fitts et al. (2017) suite, a factor of 15 less than the other haloes simulated here, providing a tentative indication of a connection between total stellar mass formed and the presence of massive star clusters at fixed z = 0 halo mass; we return to this point in Section 4. In the following sections, we discuss the formation of the clusters, their properties, and connections with their host galaxies and dark matter haloes.
3 FORMATION AND PROPERTIES OF STELLAR CLUSTERS
The basic properties of each cluster are listed in the first several columns of Table 1. The clusters form between redshift 4.5 and 11.0 and have stellar masses that range from 5 × 104 to |$5\times 10^5\, \mathrm{M}_{\odot }$| at formation. The 3D half-mass sizes of the clusters are all smaller than 30 pc; as noted in Section 2, the sizes we find here should be considered upper limits, as our adopted force softening for baryonic particles (|$2\!-\!4\, {\rm pc}$|) precludes the formation of significantly denser systems. The clusters typically form with very low metallicities (|$[{\rm Fe/H}] \approx -2.5~{\rm to}\, -3$|), with dispersion in the iron metallicity of 0.15–0.2 dex. These metallicities are comparable to or slightly lower than those of the most metal-poor clusters observed in low-mass galaxies (Beasley et al. 2019), which we discuss further in Section 4.3; the spread in metallicity is somewhat larger than observed in MW GCs (0.045 dex; Carretta et al. 2009; Bailin 2019). The iron spread in the simulated clusters is inherited from the progenitor gas clouds as opposed to resulting from self-enrichment during the star formation process. We define the formation time tform of the cluster as the time when the first 10 members of the cluster form. Our results are insensitive to this precise definition, as the full duration of star formation is less than 10 Myr in all of the clusters studied here.
Properties of the stellar clusters and their host galaxies, all computed at the epoch of cluster formation. Columns specify: (1) zf: redshift the cluster has formed; (2) r1/2: 3D stellar half-mass radius of the cluster at the formation time; (3) Mcl: mass of each cluster at the formation epoch; (4) d: the distance from its host; (5) t50: the cosmological time since the formation epoch for each cluster that has lost 50 per cent of its mass; (6) Δtform: formation time window in which all the cluster members have formed; (7) [Fe/H]: the iron abundance of the cluster; (8) σ[Fe/H]: the spread in the iron abundance; (9) Mvir: host DM halo virial mass; (10) rvir: host’s virial radius; (11) Vm: host’s maximum circular velocity; (12) M⋆: stellar mass (excluding the stellar mass contained in the nascent clusters) within |$0.1\, r_{\rm vir}$| of the host halo’s centre; (13) r1/2, h: host’s 3D stellar half-mass radius; (14) [Fe/H]h: iron abundance for all stars associated with the host. The m10i simulations forms two clusters with |$M_{\rm cl}\gt 10^5\, {\rm M}_\odot$| at different redshifts, with the later-forming cluster (denoted m10i† in the table) hosted in a separate halo that is not part of the merger tree of the main halo (but is still within the high-resolution region of the simulation). The clusters denoted with an asterisk, m10l and m10m, each form two coeval clusters that form separately out of a single massive GMC within the main galaxy and shortly thereafter merge together to form one massive, long-lived cluster. One halo, m10b, does not form any bound clusters in excess of |$M_\star =5\times 10^4\, \mathrm{M}_{\odot }$|.
| cluster(s) | host dark matter halo and galaxy (at zf) | |||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Halo | zf | r1/2 | |$M_{\star ,\rm cl}$| | d | t50 | Δtform | [Fe/H] | σ[Fe/H] | Mvir | rvir | Vm | M⋆ | r1/2, h | [Fe/H]h |
| (pc) | (M⊙) | (pc) | (Myr) | (Myr) | (M⊙) | (kpc) | (km s−1) | (M⊙) | (pc) | |||||
| m10b | – | – | – | – | – | – | – | – | – | – | – | – | – | – |
| m10i | 11.0 | 26 | 2.3 × 105 | 325 | 197 | 6.9 | −2.5 | 0.20 | 5.1 × 108 | 2.2 | 34 | 3.9 × 104 | 124 | −3.1 |
| m10i† | 4.9 | 27 | 1.6 × 105 | 310 | 270 | 9.6 | −2.7 | 0.17 | 1.3 × 109 | 5.8 | 33 | 1.6 × 105 | 175 | −3.0 |
| m10l* | 8.1 | 17 | 2.2 × 105 | 802 | 1800 | 7.4 | −2.9 | 0.19 | 7.4 × 108 | 3.2 | 34 | 2.1 × 105 | 171 | −2.9 |
| m10j | 9.0 | 18 | 1.3 × 105 | 1280 | 1620 | 9.2 | −2.3 | 0.14 | 4.4 × 108 | 2.4 | 29 | 4.9 × 104 | 85 | −3.4 |
| m10k | 4.5 | 19 | 3.6 × 105 | 60 | 940 | 7.2 | −2.8 | 0.18 | 1.2 × 109 | 6.4 | 29 | 1.6 × 105 | 169 | −2.8 |
| m10m* | 6.0 | 15 | 4.5 × 105 | 628 | 2500 | 7.3 | −2.8 | 0.16 | 1.6 × 109 | 5.4 | 39 | 1.6 × 105 | 320 | −3.0 |
| m10h | 5.7 | 6 | 5.3 × 104 | 614 | 94 | 4.2 | −2.7 | 0.17 | 1.1 × 109 | 5.0 | 33 | 1.1 × 105 | 370 | −2.7 |
| cluster(s) | host dark matter halo and galaxy (at zf) | |||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Halo | zf | r1/2 | |$M_{\star ,\rm cl}$| | d | t50 | Δtform | [Fe/H] | σ[Fe/H] | Mvir | rvir | Vm | M⋆ | r1/2, h | [Fe/H]h |
| (pc) | (M⊙) | (pc) | (Myr) | (Myr) | (M⊙) | (kpc) | (km s−1) | (M⊙) | (pc) | |||||
| m10b | – | – | – | – | – | – | – | – | – | – | – | – | – | – |
| m10i | 11.0 | 26 | 2.3 × 105 | 325 | 197 | 6.9 | −2.5 | 0.20 | 5.1 × 108 | 2.2 | 34 | 3.9 × 104 | 124 | −3.1 |
| m10i† | 4.9 | 27 | 1.6 × 105 | 310 | 270 | 9.6 | −2.7 | 0.17 | 1.3 × 109 | 5.8 | 33 | 1.6 × 105 | 175 | −3.0 |
| m10l* | 8.1 | 17 | 2.2 × 105 | 802 | 1800 | 7.4 | −2.9 | 0.19 | 7.4 × 108 | 3.2 | 34 | 2.1 × 105 | 171 | −2.9 |
| m10j | 9.0 | 18 | 1.3 × 105 | 1280 | 1620 | 9.2 | −2.3 | 0.14 | 4.4 × 108 | 2.4 | 29 | 4.9 × 104 | 85 | −3.4 |
| m10k | 4.5 | 19 | 3.6 × 105 | 60 | 940 | 7.2 | −2.8 | 0.18 | 1.2 × 109 | 6.4 | 29 | 1.6 × 105 | 169 | −2.8 |
| m10m* | 6.0 | 15 | 4.5 × 105 | 628 | 2500 | 7.3 | −2.8 | 0.16 | 1.6 × 109 | 5.4 | 39 | 1.6 × 105 | 320 | −3.0 |
| m10h | 5.7 | 6 | 5.3 × 104 | 614 | 94 | 4.2 | −2.7 | 0.17 | 1.1 × 109 | 5.0 | 33 | 1.1 × 105 | 370 | −2.7 |
Properties of the stellar clusters and their host galaxies, all computed at the epoch of cluster formation. Columns specify: (1) zf: redshift the cluster has formed; (2) r1/2: 3D stellar half-mass radius of the cluster at the formation time; (3) Mcl: mass of each cluster at the formation epoch; (4) d: the distance from its host; (5) t50: the cosmological time since the formation epoch for each cluster that has lost 50 per cent of its mass; (6) Δtform: formation time window in which all the cluster members have formed; (7) [Fe/H]: the iron abundance of the cluster; (8) σ[Fe/H]: the spread in the iron abundance; (9) Mvir: host DM halo virial mass; (10) rvir: host’s virial radius; (11) Vm: host’s maximum circular velocity; (12) M⋆: stellar mass (excluding the stellar mass contained in the nascent clusters) within |$0.1\, r_{\rm vir}$| of the host halo’s centre; (13) r1/2, h: host’s 3D stellar half-mass radius; (14) [Fe/H]h: iron abundance for all stars associated with the host. The m10i simulations forms two clusters with |$M_{\rm cl}\gt 10^5\, {\rm M}_\odot$| at different redshifts, with the later-forming cluster (denoted m10i† in the table) hosted in a separate halo that is not part of the merger tree of the main halo (but is still within the high-resolution region of the simulation). The clusters denoted with an asterisk, m10l and m10m, each form two coeval clusters that form separately out of a single massive GMC within the main galaxy and shortly thereafter merge together to form one massive, long-lived cluster. One halo, m10b, does not form any bound clusters in excess of |$M_\star =5\times 10^4\, \mathrm{M}_{\odot }$|.
| cluster(s) | host dark matter halo and galaxy (at zf) | |||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Halo | zf | r1/2 | |$M_{\star ,\rm cl}$| | d | t50 | Δtform | [Fe/H] | σ[Fe/H] | Mvir | rvir | Vm | M⋆ | r1/2, h | [Fe/H]h |
| (pc) | (M⊙) | (pc) | (Myr) | (Myr) | (M⊙) | (kpc) | (km s−1) | (M⊙) | (pc) | |||||
| m10b | – | – | – | – | – | – | – | – | – | – | – | – | – | – |
| m10i | 11.0 | 26 | 2.3 × 105 | 325 | 197 | 6.9 | −2.5 | 0.20 | 5.1 × 108 | 2.2 | 34 | 3.9 × 104 | 124 | −3.1 |
| m10i† | 4.9 | 27 | 1.6 × 105 | 310 | 270 | 9.6 | −2.7 | 0.17 | 1.3 × 109 | 5.8 | 33 | 1.6 × 105 | 175 | −3.0 |
| m10l* | 8.1 | 17 | 2.2 × 105 | 802 | 1800 | 7.4 | −2.9 | 0.19 | 7.4 × 108 | 3.2 | 34 | 2.1 × 105 | 171 | −2.9 |
| m10j | 9.0 | 18 | 1.3 × 105 | 1280 | 1620 | 9.2 | −2.3 | 0.14 | 4.4 × 108 | 2.4 | 29 | 4.9 × 104 | 85 | −3.4 |
| m10k | 4.5 | 19 | 3.6 × 105 | 60 | 940 | 7.2 | −2.8 | 0.18 | 1.2 × 109 | 6.4 | 29 | 1.6 × 105 | 169 | −2.8 |
| m10m* | 6.0 | 15 | 4.5 × 105 | 628 | 2500 | 7.3 | −2.8 | 0.16 | 1.6 × 109 | 5.4 | 39 | 1.6 × 105 | 320 | −3.0 |
| m10h | 5.7 | 6 | 5.3 × 104 | 614 | 94 | 4.2 | −2.7 | 0.17 | 1.1 × 109 | 5.0 | 33 | 1.1 × 105 | 370 | −2.7 |
| cluster(s) | host dark matter halo and galaxy (at zf) | |||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Halo | zf | r1/2 | |$M_{\star ,\rm cl}$| | d | t50 | Δtform | [Fe/H] | σ[Fe/H] | Mvir | rvir | Vm | M⋆ | r1/2, h | [Fe/H]h |
| (pc) | (M⊙) | (pc) | (Myr) | (Myr) | (M⊙) | (kpc) | (km s−1) | (M⊙) | (pc) | |||||
| m10b | – | – | – | – | – | – | – | – | – | – | – | – | – | – |
| m10i | 11.0 | 26 | 2.3 × 105 | 325 | 197 | 6.9 | −2.5 | 0.20 | 5.1 × 108 | 2.2 | 34 | 3.9 × 104 | 124 | −3.1 |
| m10i† | 4.9 | 27 | 1.6 × 105 | 310 | 270 | 9.6 | −2.7 | 0.17 | 1.3 × 109 | 5.8 | 33 | 1.6 × 105 | 175 | −3.0 |
| m10l* | 8.1 | 17 | 2.2 × 105 | 802 | 1800 | 7.4 | −2.9 | 0.19 | 7.4 × 108 | 3.2 | 34 | 2.1 × 105 | 171 | −2.9 |
| m10j | 9.0 | 18 | 1.3 × 105 | 1280 | 1620 | 9.2 | −2.3 | 0.14 | 4.4 × 108 | 2.4 | 29 | 4.9 × 104 | 85 | −3.4 |
| m10k | 4.5 | 19 | 3.6 × 105 | 60 | 940 | 7.2 | −2.8 | 0.18 | 1.2 × 109 | 6.4 | 29 | 1.6 × 105 | 169 | −2.8 |
| m10m* | 6.0 | 15 | 4.5 × 105 | 628 | 2500 | 7.3 | −2.8 | 0.16 | 1.6 × 109 | 5.4 | 39 | 1.6 × 105 | 320 | −3.0 |
| m10h | 5.7 | 6 | 5.3 × 104 | 614 | 94 | 4.2 | −2.7 | 0.17 | 1.1 × 109 | 5.0 | 33 | 1.1 × 105 | 370 | −2.7 |
3.1 Cluster formation mechanisms
Fig. 1 shows the evolution of the progenitor gas cloud in the m10l halo prior the formation of a cluster at z = 8.1. The top portion contains six panels showing the time evolution of the gas surface density Σgas in the cluster’s environment prior to and immediately following the cluster formation epoch (times relative to the formation of the cluster, Δt ≡ t − tform, are given in each panel). Each panel spans |$(4\, \text{kpc})^2$| with a depth of 1 kpc and is centred on the centre of mass of the gas progenitors that eventually give birth to the cluster members; Σgas is computed in cells with areas of |$(4\, {\rm pc})^2$| and depths of 1 kpc. The larger subplot on the right is a zoomed-in view of the gas cloud in the midst of cluster formation (with dimension of |$(1\, \text{kpc})^2 \times 0.5~{\rm kpc}$|). The central galaxy of the halo, with centre coincident with the centre of the halo, is indicated in each subpanel by a white circle, with radius equal to the galaxy’s stellar half-mass radius at that time. The newly formed cluster is shown as a gold circle in the final panel, with size equal to its half-mas radius; it is significantly offset from the galaxy at this epoch.
Top: Time evolution of the gas surface density of the progenitor gas cloud that gives birth to a bound stellar cluster in m10l galaxy at z = 8.1. Each smaller sub-panel is |$(4\, {\rm kpc})^2$| with a depth of |$1\, \text{kpc}$|, centred on the centre of mass of the gas progenitors that eventually give birth to the cluster members. The time relative to the cluster formation epoch is noted in each frame and the pre-existing galaxy is indicated with a white circle with size equal to its stellar half-mass radius. The larger sub-panel on the right-hand side, spanning |$(1\, {\rm kpc})^2$| with a depth of |$0.5\, \text{kpc}$|, zooms into the region where the cluster forms just prior to its birth. At this epoch, a substantial amount of gas, |$\approx 10^7\, \mathrm{M}_{\odot }$|, reaches a surface density of |$\Sigma _{\rm gas} \sim 10^4\, \mathrm{M}_{\odot }\, {\rm pc}^{-2}$|, a necessary condition for cluster formation, as a result of a cloud–cloud collision. The cluster itself is shown at |$\Delta t=13\, {\rm Myr}$| as a gold circle with size equal to the cluster half-mass radius (20 pc). Bottom: |$(2\, {\rm kpc})^2$| gas surface density maps, with depth of |$1\, {\rm kpc}$|, at the last three times shown in the upper panels. Pre-existing stars are shown in cyan, stars that are part of the bound cluster that forms are shown in red, and stars forming from Δt = −17 Myr to the snapshot in question that are not bound to the cluster are shown in orange. The cloud–cloud collision results in a very high surface density outside of the pre-existing galaxy (whose stellar half-mass radius is shown in black in each panel). This gas is further compressed by supernovae resulting from stars born in the colliding gas (orange in left panel), leading to the formation of the stellar cluster. Many stars that are not formally identified as cluster members are formed at the same time in the immediate vicinity of the cluster as well (right panel), and feedback from the cluster formation launches a compressive shock that causes star formation within the pre-existing galaxy.
The panels show two dense patches of gas – originating from a somewhat earlier merger event – approaching one another with halo-centric speeds comparable to the |$v \approx 30\, {\rm km\,s^{-1}}$| virial velocity of the halo (which is also comparable to the turbulent velocity dispersion of the gas) and colliding, leading to a region of very high pressure and density by |$\Delta t=-7 \, {\rm Myr}$|. These are the requisite conditions for cluster formation. It is evident that feedback has cleared gas from the cluster’s immediate vicinity by |$\Delta t=13\, {\rm Myr}$|, leading to a compressive superbubble that stimulates high gas densities hundreds of pc from the cluster.
In the bottom panels of Fig. 1, we show grey-scale gas surface density maps, each with dimension |$(2 \, {\rm kpc})^2$| with a depth of 1 kpc, at the three snapshots closest to the epoch of cluster formation (|$\Delta t=-7, \, 3$|, and 13 Myr). The stars formed before |$\Delta t=-7\, {\rm Myr}$| are shown in cyan, while the stars formed during the epoch of cluster formation are shown in red (for particles identified by Phinder as cluster members) and gold (for stars identified as non-cluster members). The half-light radius of the pre-existing galaxy is shown as a black circle in each panel. The cloud–cloud collision results in very high surface densities near, but outside of, the pre-existing galaxy at |$-7\, {\rm Myr}$|. Feedback from a handful of stars formed in this gas causes further compression, initiating the formation of two coeval sub-clusters in very close proximity and a smattering of other stars. By |$\Delta t=13\, {\rm Myr}$|, the cluster stars have formed, along with a significant population of nearby stars that are not formally associated with the cluster but none the less form in its immediate environment, at distances of |$\lt 200\, {\rm pc}$|. Feedback from these combined populations clears the gas from the cluster region and launches a bubble, which collides with the pre-existing galaxy and initiates further star formation there as well; we discuss this extra-cluster star formation further in Section 3.3. Finally, the two adjacent sub-clusters merge within ∼40 Myr of their formation, with each sub-cluster contributing roughly half of the final stellar mass in the merged cluster; the cluster in halo m10m follows a very similar formation channel.
Fig. 2 is analogous to Fig. 1 but shows the gas surface density evolution immediately preceding the formation of a cluster in the m10i halo. In this case, there is no cloud–cloud collision; rather, filamentary gas accretion leads to high surface density near the centre of the halo. This high surface density initially promotes clustered star formation but not the formation of a cluster (|$\Delta t=-17 \, {\rm to}\, -7$| Myr). However, feedback from the subsequent SN explosions resulting from that star formation drives a bubble that compresses the already dense ISM even further, resulting in the formation of a cluster. Within 10 Myr, feedback from the lives and deaths of massive stars in the cluster has expelled all of the gas and star formation ceases.
Analogous to Fig. 1 but for the bound stellar cluster formed in m10i galaxy at z = 10.8. Unlike the cluster in m10l, there is no cloud–cloud collision inducing the formation of the cluster. In this case, it is filamentary gas accretion, coupled with compression due to supernova feedback, that pushes gas to |$\Sigma _{\rm gas}=10^4\, {\rm M}_\odot \, {\rm pc}^{-2}$|. In this case, the cluster formation happens at such an early epoch that there is barely even a ‘host galaxy’ to speak of, and the stellar mass formed as a result of the cluster formation far exceeds the pre-existing stellar content of the halo.
Figs 1 and 2 make it clear that there are multiple pathways to get to high densities and pressures that are conducive to cluster formation. To better quantify what constitutes ‘high density’ and how this is linked to the formation of stellar clusters, we compute Σgas in cells of |$4\times 4\times 1000\, {\rm pc}$|, centred on centre-of-mass of gas particles that give birth to the cluster, and plot the cumulative distribution of this quantity in 10 Myr intervals over a 50 Myr window spanning cluster formation, −37 ≤ Δt/Myr ≤ 13, in Fig. 3. The colours of the lines indicate the time relative to tform via the colourbar at the right; the thick solid line corresponds to the snapshot closest to cluster formation, Δt = 3 Myr. In both cases highlighted here, the maximum surface density is initially a few hundred |$\mathrm{M}_{\odot }\, {\rm pc}^{-2}$| and quickly rises until |${\sim }0.01~{{\ \rm per\ cent}}$| of the cells exceed |$\Sigma _{\rm gas} \approx 10^{4}\, {\rm pc}^{-2}$|. The bottom panels show the related quantity Mgas(> Σgas) for the same cells and indicate that |${\sim } 10^{7}\, \mathrm{M}_{\odot }$| of gas is contained in the cells exceeding |$\Sigma _{\rm gas} \approx 10^{4}\, \mathrm{M}_{\odot }\, {\rm pc}^{-2}$|, meaning that the region containing this high density has a surface area of order (20 pc)2. At this point, the clusters form rapidly: in each case studied in this paper, all star formation in the cluster takes place in less than 10 Myr.
Upper panels: The cumulative distribution of gas cells (|$4\times 4\times 1000\, {\rm pc}$| each) exceeding a surface density Σgas in 10 Myr increments immediately prior to and during the cluster formation epoch in cluster m10l (left) and m10i (right); line colours indicate the time relative to formation as shown in the colour bars on the right, with the thick solid line corresponding to Δt = 3 Myr, which is the epoch closest to tform. Lower panels: The cumulative distribution of gas mass in the same cells exceeding Σgas for the same snapshots. For halo m10l, we also show both quantities at Δt = 23 Myr (grey dotted line) and 33 Myr (black dotted line). Cluster formation occurs when |$\gtrsim 10^7\, \mathrm{M}_{\odot }$| of gas exceeds |$10^4\, {\rm M}_\odot \, {\rm pc}^{-2}$|, a result that holds true for all clusters in our simulation suite. The gas surface density drops precipitously after the cluster forms in all cases. Generally, the drop occurs very shortly after cluster formation (|$\Delta t = 10\, {\rm Myr}$|). In m10l, an even steeper drop occurs but is delayed by 20 Myr because of the additional star formation triggered by feedback from the cluster’s formation (see Fig. 1). However, in all cases, the cluster formation site itself is cleared of gas within 10 Myr of cluster formation.
The resulting feedback immediately removes all gas from newly formed cluster. Accompanying the gas removal from the cluster is a precipitous drop in both f(> Σgas) and M(> Σgas) at high gas surface densities (|$\Sigma _{\rm gas} \gtrsim 10^{3}\, \mathrm{M}_{\odot }\, {\rm pc}^{-2}$|) in m10i. Halo m10l does not see an immediate drop (and in fact sees a temporary increase) because of the star formation that is triggered within the galaxy by the feedback that clears the gas from the nascent cluster (see Fig. 1). However, m10l sees the same drop starting ∼20 Myr after cluster formation (grey dotted curves in Fig. 3), and after another 10 Myr (black dotted curves), feedback has removed most of the gas from the region considered here. We conclude that stellar feedback from cluster formation is very effective at clearing gas from the cluster formation site on the cluster formation time-scale of ≲10 Myr, in agreement with observations (Bastian et al. 2014; Hollyhead et al. 2015; Krumholz, McKee & Bland-Hawthorn 2019).
In the two cases highlighted here, and indeed in all the identified clusters, cluster formation occurs once giant molecular clouds (GMCs) reach surface densities exceeding |$\Sigma _{\rm thresh}=10^{4}\, \mathrm{M}_{\odot }\, {\rm pc}^{-2}$| (see also Murray, Quataert & Thompson 2010; Colín, Vázquez-Semadeni & Gómez 2013; Geen, Soler & Hennebelle 2017; Elmegreen 2018; Grudić et al. 2018a; Kim et al. 2018; Krumholz et al. 2019). This threshold surface density, which is a factor of ∼100 higher than is typical of GMCs in the Milky Way (Bolatto et al. 2008), is equivalent to a gravitational force per unit area (i.e. pressure) of |$P_{\rm thresh} \equiv G\, \Sigma _{\rm thresh}^2=2 \times 10^9 \, k_{\rm B}\, K \, {\rm cm}^{-3}$|. It is also what would be achieved by gas patches with density |$n\approx 10^{3}\, {\rm cm}^{-3}$| colliding at a relative velocity of |$v_{\rm rel} \approx 120\, {\rm km\,s^{-1}}$| (assuming |$P=\rho \, v_{\rm rel}^2$|), which is relevant for m10l (see Fig. 1): gas clumps on ballistic trajectories in a |$V_{\rm vir} \approx 30\, {\rm km\, s}^{-1}$| halo colliding near the halo’s centre will achieve a gravitational pressure that is close to Pthresh, and any additional force from SN feedback near the collision will likely be sufficient to achieve Pthresh. We do not find sufficiently high gas masses in halo m10b above this threshold, which likely explains its lack of GCCs at masses we can resolve.
The gas surface density can occasionally exceed Σthresh at times that are not linked to the formation of the massive and bound stellar clusters studied in this work. In all of those times, we find lower mass (|$M_\star \lt 5\times 10^4\, \mathrm{M}_{\odot }$|) clusters forming, and we also find that the amount of mass above Σthresh is commensurately lower: the mass of the cluster(s) that form is roughly a constant fraction of the high surface density gas, |$M_{\star } \approx 0.01\, M_{\rm gas}(\gt \Sigma _{\rm thresh})$|. We defer a more thorough investigation of the connection between gas surface density and cluster formation efficiency [and star formation rate (SFR)] to a future paper.
3.2 Cluster evolution after formation
While the formation of star clusters is an important topic in and of itself, the subsequent survival of those clusters is a crucial consideration when attempting to connect compact, bound stellar systems at high redshifts to GCs at z = 0. We study the evolution of the star clusters formed at high redshift in our simulations by tracking clusters across simulation snapshots. To do so, we run Phinder separately at each snapshot and then use the IDs of the each cluster at birth (the highest redshift snapshot at which it is identified) to find its descendant at each later time. The bound clusters formed in the simulations lose mass with time owing to both tidal interactions with their environment and numerical effects. We quantify cluster lifetimes via t50, the time relative to its birth at which a cluster has lost 50 per cent of its original stellar mass; the clusters typically disrupt completely soon after t50. We find a range of t50 values, from ∼100 Myr to 2.5 Gyr (see Table 1). The longest lived clusters are those that form at the greatest distances from the centre of their host haloes, indicating the importance of tides – coupled with numerical effects – in disrupting the clusters in our simulations. Notably, only one cluster forms within its host galaxy’s half-mass radius and that some form at |$\gt \, 4\, r_{\rm 1/2,h}$|, effectively completely outside of the galaxy.
In Fig. 4, we show the time evolution of the stellar mass profile of the (merged) m10l cluster over |$2.5\, {\rm Gyr}$| after its formation. For the first ∼1.3 Gyr, corresponding to over 100 crossing times, the cluster’s inner mass profile is remarkably stable, while secular mass-loss outside of r1/2 reduces the cluster mass by |${\sim 30~{{\ \rm per\ cent}}}$|. Over this period, the cluster completes approximately 13 regular orbits around its host galaxy, with stable apocentres of ≈1.1 kpc and pericentres of ≈200 pc. Subsequently, the mass-loss accelerates and the cluster eventually dissolves after |${\sim } 2.5\, {\rm Gyr}$|. The half-mass radius of the cluster remains nearly constant over the cluster’s entire evolution. Although our simulations do not have the ability to accurately track the orbits of stellar particles over a Hubble time – which would require a direct N-body integrator, not the softened force algorithm employed by gizmo – and the force softening adopted here precludes the formation of clusters with ∼pc-scale half-mass radii, these results indicate that the clusters forming in these simulations are long-lived and are plausible progenitors of present-day GCs that are observed in some dwarf galaxies.
The time evolution of the stellar mass profile of the GCC that forms in the m10l halo. The inner profile is very stable over the |${\sim }1.5\, {\rm Gyr}$|, a period over which the cluster completes nearly 15 orbits about its host galaxy. A combination of numerical and tidal effects lead to its dissolution after nearly 2.5 Gyr.
3.3 Connections between stellar clusters and their host haloes
The last six columns of Table 1 contain information about the dark matter halo within which each cluster forms and the properties of the central galaxy of each halo at the time of cluster formation. As expected for haloes that have |$M_{\rm vir}(z=0) =10^{10}\, \mathrm{M}_{\odot }$|, the progenitor haloes at high redshifts – when the clusters form – are approximately an order of magnitude lower in mass, |$M{\sim } (0.5\!-\!1.5)\times 10^{9}\, \mathrm{M}_{\odot }$|, and the maximum circular velocity of the host at the time of cluster formation is |$30\!-\!40~ {\rm km\,s^{-1}}$|. These numbers indicate that these GCCs are forming at high redshift in haloes that are a factor of ∼10 more massive than the atomic cooling limit of |$T_{\rm vir}=10^4\, {\rm K} \leftrightarrow V_{\rm halo} = 17\, {\rm km\,s^{-1}}$| (or |$M_{\rm vir} \approx 10^8\, \mathrm{M}_{\odot }$| at z ∼ 8).
At the time the clusters form, the haloes typically have very low stellar content: the most massive ‘galaxy’ – defined as the stellar mass within |$0.1\, r_{\rm vir}$| prior to the formation of the cluster – is |$10^5\, \mathrm{M}_{\odot }$|, and in some cases, only |$10^{4}\, \mathrm{M}_{\odot }$| of stars exist near the centre of the halo prior to the time of cluster formation. In this sense, the formation of the clusters is indeed pre-galactic, as there is essentially no galaxy present prior to the cluster formation episode. The pre-existing stellar populations tend to be fairly extended (|$r_{1/2,\star } \approx 150\!-\!400\, {\rm pc}$|) and metal-poor ([Fe/H] ∼ −3).
Fig. 5 shows the SFRs (top) and the archaeological SFHs5 (bottom) of the haloes in our fiducial suite; in both cases, we consider all stars within rvir and use time bins of 10 Myr. We separate the SFR in each cluster (shown in red) from the rest of the stars in the halo (shown in black). In the SFH plots, we include both stars in clusters and all other stars, but we mark the cluster formation epoch (red vertical bands) and the clusters’ contribution to the total SFHs (blue horizontal bands). The star formation in these simulated haloes is highly episodic, as has been noted before for simulations both at this mass scale (Stinson et al. 2007; Domínguez et al. 2015; Fitts et al. 2017; Sparre et al. 2017; Emami et al. 2019): many have prolonged periods of true quiescence punctuated by brief periods of relatively vigorous star formation (|$\dot{M}_{\star } \approx 0.01\!-\!0.1\, \mathrm{M}_{\odot }\, {\rm yr}^{-1}$| for tens of Myr; see also fig. 1 of Boylan-Kolchin et al. 2015).
Top: SFR within the virial radius of the main halo in each m10 realization with our fiducial FIRE2 physics, computed for 10 Myr time intervals. The formation of the stars formed in the massive bound clusters within each primary galaxy is plotted in red, while the black bars show all other star formation in the main halo. Bottom: the stellar mass assembly of the simulated galaxies over the cosmological evolution. The red vertical bands specify the cluster formation epochs for each realization, and the blue horizontal bands show the contribution of the cluster’s mass to the archaeological SFHs of each primary galaxy. The star formation in these haloes is very bursty and episodic, and the clusters often form at times of (relatively) intense but short-lived star formation. The star formation event leading to the formation of the cluster is often comprises most of the stellar content in the halo at the formation epoch, and in some cases, the cluster itself contains the majority of the stars within the halo at the time of its formation. We do not include the cluster that forms in the m10i box at z = 4.9, as it forms outside of the virial radius of the most massive galaxy.
The clusters contribute substantially to the total stellar mass within the virial radius at the time of cluster formation (see also Table 1). In fact, for all but the least massive cluster we consider here (m10h), the clusters outweigh the pre-existing stellar mass within the virial radius at the time of cluster formation; the earliest forming cluster (in m10i) is 10 times more massive than the pre-existing stars in that halo. Even in m10h, the cluster has a mass that is equal to 50 per cent of the total stellar mass that was present in the halo prior to cluster formation. These results demonstrate that high-redshift star clusters in progenitors of present-day dwarf galaxies may contain a substantial – or even dominant – fraction of the total stellar mass in the halo at the time the clusters form. As seen in Fig. 5, the clusters often form at times of star formation bursts in the dwarfs, which is not surprising given the importance of SN feedback in producing conditions conducive to star cluster formation (as discussed in Section 3.1). As a result, stars formed within ∼50 Myr of the cluster formation epoch often account for over 10 per cent (and, in the case of m10l, almost 70 per cent) of all stars residing within the halo’s virial radius at z = 3. The formation of both of the clusters themselves and of stars associate with the clusters’ birth clouds are defining events for the SFHs of these simulated dwarf galaxies. We return to observational implications of high fraction of a halo’s total stellar mass at the epoch of cluster formation attributable to the cluster in Section 4.1.
The formation sites of the GCs are also of substantial interest. Many ‘pre-galactic’ models for GCs – or standard interpretations of these models – assume that GCs form at the centres of dark matter mini-haloes, leading to the prediction that GCs should be immersed in a dense dark matter cocoon that should be detectable in the kinematics at the outskirts of GCs (Peebles 1984; Heggie & Hut 1996; Mashchenko & Sills 2005; Conroy, Loeb & Spergel 2011; Ibata et al. 2013; Peñarrubia et al. 2017; Boldrini, Mohayaee & Silk 2020). However, dense molecular clouds do not necessarily form directly at the centre of dark matter haloes, and so it is not at all obvious that even pre-galactic GC models necessarily produce GCs that form directly at the minimum of the dark matter gravitational potential and bear the indelible imprints of dark matter hosts. It is natural to assume that GCs form near the centres of their host galaxies, as the galaxies by definition trace the bulk of the star formation in the halo, but we have shown the GCCs simulated here can contain comparable stellar content to the galaxy in the early epochs corresponding to the GC formation times for these dwarf galaxies (which also form in atomic cooling haloes, not mini-haloes).
In fact, as Table 1 indicates, many of the simulated clusters form at substantial distances (|$3-5\, r_{1/2,\rm gal}$| or more) from their ‘host galaxies’ as defined by the pre-existing stars. While often there are connections between the galaxies and clusters, as evidenced by periods of increased star formation activity in the galaxy at the same time as cluster formation even for clusters that form well outside their host’s stellar half-mass radius, the clusters are generally not deeply embedded within the nascent galaxies, even though the haloes studied here all have cusped density profiles at their centres. This is likely due to a combination of the high redshift of formation, the (relatively) low masses of dwarf galaxy dark matter haloes, and conditions necessary for cluster formation: dense gaseous regions in the turbulent ISM of dwarf galaxy progenitors can span several hundred pc (see Figs 1 and 2), significantly exceeding the typical size of any pre-existing ‘galaxy’. As a result, stochastic events that trigger cluster formation – cloud collisions or compressive supernova shocks – often occur at locations that do not coincide with the bulk of the existing stars.
The large sizes (relative to the entire dark matter halo) of the GMC complexes that result in star clusters also results in clusters that do not typically form within the regions of the highest dark matter density, offering a natural explanation of the lack of observed dark matter around GCs even in scenarios where they form in low-mass haloes at high redshifts. In all cases studied here, the clusters’ mass is heavily dominated by stars, with dark matter contributing at most 15 per cent of the mass within r1/2. This dark matter contribution – equal to at most 10 dark matter particles – is consistent with the expected dark matter contribution of the main halo at the clusters’ formation location as opposed to these clusters forming at the centres of their own dark matter haloes.
4 DISCUSSION
4.1 Observational consequences
Given the recent advances in both theoretical and observational understanding of in situ star cluster formation in the distant Universe, along with the promise of forthcoming JWST observations, the potential observational implications of the results presented in Section 3 are worth exploring. Most of the clusters in our current sample – four of the primary six clusters – attain SFRs of |$0.02\!-\!0.03~\mathrm{M}_{\odot }\, {\rm yr}^{-1}$| (corresponding to |$(2\!-\!3)\times 10^{5}\, \mathrm{M}_{\odot }$| of stars formed in a 10 Myr window) at z ∼ 5–11, which would result in an absolute luminosity of MUV ≈ −14 mag (m ≈ 33 at z ≈ 8) using the stellar population modelling results of Boylan-Kolchin (2017) in the limit of no dust attenuation (which is likely appropriate for these very metal-poor, low-mass galaxies at high redshift). Given the sizes of these simulated systems, |$r_{1/2} \approx 20\, {\rm pc}$|, the corresponding surface brightness would be roughly |$\mu =22\!-\!23\, {\rm mag\, arcsec^2}$|, with an SFR surface density of |$\Sigma _{\rm SFR} \approx 10^{2}\, \mathrm{M}_{\odot }\, {\rm yr^{-1}\, kpc^{-2}}$| for |$\approx 10\, {\rm Myr}$|. These properties are in line with both predictions from GC formation models (Ricotti et al. 2016; Renzini 2017; Boylan-Kolchin 2018) and observations of possible star clusters in formation at high redshift (Vanzella et al. 2017, 2019, 2021).
The sizes of the clusters at formation – r1/2 ≈ 20 pc – are consistent with the most compact sources with comparable luminosities (MUV ≈ −15) seen in lensed Hubble Space Telescope images (Vanzella et al. 2019; Kikuchihara et al. 2020; Bouwens et al. 2021). As Figs 1, 2, and 5 demonstrate, there is typically substantial additional star formation that is coeval with cluster formation; accounting for this extra-cluster star formation, the total SFRs can be |$\approx 0.1\, \mathrm{M}_{\odot }\, {\rm yr}^{-1}$| over an area of several hundred pc. This would result in an absolute luminosity of MUV ≈ −15.5 mag (m ≈ 31.5 at z ≈ 8) for the full host system, which is at the edge of detectability in deep JWST fields. In lensing fields, the GCCs would likely stand out from the rest of the coeval star formation because of their much higher surface densities (Zick, Weisz & Boylan-Kolchin 2018).
Although all of the star clusters discussed here live for tens to hundreds of dynamical times, none survive to the present day. This is likely a reflection of the artificially large size of the clusters (owing to our adopted force softening) and the use of a treepm code, rather than a code that directly implements gravitational interactions, for our simulations. In the future, we plan to use hybrid schemes to more faithfully track the evolution of clusters formed at high redshift in progenitors of present-day dwarf galaxies to see if the clusters are indeed analogues of the old, metal-poor clusters observed in some low-mass dwarfs.
The results presented here are promising, if non-definitive, in this regard. Many properties of the clusters studied here are consistent with those expected of ancient GCs in Local Group dwarfs at earlier epochs, closer to their formation. In particular, the clusters here (1) contain a substantial fraction of the stars relative to the host galaxy; (2) have similar chemistry to the oldest stars in their host galaxies; and (3) form at locations that would likely allow them to survive to the present day. The conditions that lead to the formation of bound stellar clusters at high redshifts in our simulated dwarf galaxies therefore represent plausible formation scenarios for the growing population of observed ancient GCs in nearby dwarf galaxies. These clusters also form with the efficiency required for reproducing the observed relationship between GC system mass (MGCs) and dark matter halo mass (Mhalo) in the local Universe of MGCs/Mhalo ≈ 4 × 10−5 (Blakeslee, Tonry & Metzger 1997; Spitler & Forbes 2009; Georgiev et al. 2010; Hudson et al. 2014; Harris, Harris & Hudson 2015; Burkert & Forbes 2020): as shown in Boylan-Kolchin (2017), this is naturally achieved if a |$M_{\star }=2\times 10^5\, \mathrm{M}_{\odot }$| GC forms in a |$M_{\rm halo}\approx 10^{9}\, \mathrm{M}_{\odot }$| dark matter halo at z ∼ 8, which is remarkably similar to what we find in this work.
4.2 Sensitivity to simulation choices
The GCCs described in this paper form in very high surface density gas, often following compression caused by strong stellar feedback. It is therefore important to examine the effects of choices related to stellar feedback and star formation criteria in order to understand how robust our results are to details of the prescriptions we adopt. In this section, we explore the roles that work related to stellar mass-loss and the star formation density threshold play in driving our results.
As described and tested in detail in Hopkins et al. (2018a), we account for conversion of thermal energy into kinetic energy following supernova explosions (i.e. any ‘|$P\, \mathrm{ d}V$| work’ done) during the Sedov–Taylor phase of the expansion before the remnant reaches the resolution scale Δx at which the coupling to the nearest gas cells occurs, which is crucial for producing converged results. However, as discussed there and in Hopkins et al. (2018b), how to treat the ratio of thermal to kinetic energy injected for stellar mass-loss when the mass-loss is discretized into finite time-steps is more ambiguous if there is a continuously expanding bubble below the resolution scale. Moreover, it is clear that the Sedov–Taylor solution, which assumes a single instantaneous discrete energy injection event, is not the correct solution for a continuous wind. In our fiducial simulations therefore we ignore any unresolved |$P\, \mathrm{ d}V$| work from stellar mass-loss processes6, which is consistent with both observations and recent simulations (Harper-Clark & Murray 2009; Lancaster et al. 2021a, b).
However, it is interesting to ask what might happen if stellar mass-loss bubbles did undergo a prolonged energy-conserving phase during which substantial |$P\, \mathrm{ d}V$| work was done, converting almost all of the thermalized/shocked ejecta energy into kinetic energy (momentum) on large scales. We do this by treating each stellar mass-loss event (which injects some |$\Delta M \equiv \dot{M}_{\ast }\, \Delta t$|, with initial free-streaming kinetic luminosity/energy |$\Delta E \equiv \dot{E}\, \Delta t$|) as a ‘mini-supernova’ and applying the exact same treatment as we do for SNe following Hopkins et al. (2018a). The practical effect of this (given the various scaling for e.g. the cooling radii of SNe) is that most of total energy injection by stellar mass-loss is converted into momentum/kinetic energy on resolved scales (i.e. |$\sim 100~{{\ \rm per\ cent}}$| of the stellar mass-loss energy is converted into macroscopic momentum; see appendix D of Hopkins et al. 2018b), as compared to post-shock thermal energy that can be more efficiently radiated away. We emphasize that the particular functional form we adopt has no clear physical motivation, but it provides a useful comparison point for understand possible effects of unresolved stellar mass-loss physics on our results.
Surprisingly, the ‘stronger’ stellar mass-loss given by this assumption produces effectively weaker SN feedback (higher SFRs and stellar masses, by a factor of ∼2–3). This is related to the effects shown and discussed in Hopkins et al. (2020). If ‘early’ stellar feedback (processes that act before SNe explode in a young star-forming region, e.g. radiative heating and radiation pressure and stellar mass-loss) is much weaker, then that region collapses much further and produces many more stars in a denser configuration. When, approximately 40 Myr later, those stars begin to explode, the SNe (which carry energy that is an order-of-magnitude larger than the energy attributable to stellar mass-loss) are much more strongly clustered, making it easier for bubbles to overlap and driving much stronger outflows (as has also been seen in idealized experiments that vary the strength of SNe clustering, e.g. Martizzi, Faucher-Giguère & Quataert 2015; Walch & Naab 2015; Fielding, Quataert & Martizzi 2018). When early feedback is artificially made much stronger as in our experiment, clouds are disrupted earlier with much lower star formation efficiencies, producing much weaker SNe clustering and therefore weaker net large-scale SN feedback.
A comparison of SFHs between the ‘early’ stellar feedback model and our fiducial suite shows that the fiducial suite simulations have formed slightly fewer stars, in agreement with the argument outlined above. Furthermore, no bound stellar clusters above our chosen threshold cluster mass of M⋆ = 5 × 104 M⊙ form in the simulations with modified stellar mass-loss processes (stronger early feedback; less clustered SNe feedback). However, lower mass clusters |$(M_{\star } \lt 5\times 10^4\, \mathrm{M}_{\odot }$|) do form in the simulations with the alternate mass-loss treatment (which is the default treatment in FIRE-2 but not in FIRE-3). These low-mass clusters are always destroyed quickly (within 20 Myr). Kim et al. (2018) and Ma et al. (2020) both found that massive and long-lived clusters do form using this alternate mass-loss treatment in their simulations, which focus on significantly more massive systems (|$M_{\rm halo} \gtrsim 10^{10}\, \mathrm{M}_{\odot }$| at z ∼ 6) that have SFRs that are at least an order of magnitude higher than what we find here. This indicates that the larger scale impact of the treatment of |$P\, dV$| work from stellar mass-loss depends on the mass of the cluster: higher mass GMCs, which result in higher SFRs, are less dependent on the treatment of unresolved |$P\, dV$| work, likely because of their higher binding energies. One caveat to this conclusion is the relatively large force softening adopted in these simulations, which effectively limits the range of densities that star-forming gas can reach.
A related point is that cluster formation efficiency is also sensitive to the star formation criteria adopted in simulations (Grudić et al. 2018a, b; Ma et al. 2020). In our fiducial suite, we impose a star formation density threshold of |$n_{\rm H}\ge n_{\rm thresh}=1000\, {\rm cm}^{-3}$| (Hopkins et al. 2018b). We have also performed a simulation of m10m where we change only the star formation criterion (to |$n_{\rm thresh}=100\, {\rm cm}^{-3}$|), and leave everything else unmodified from our fiducial simulation set-up. In this simulation, stellar clusters form in a nearly identical manner as in the fiducial run. If we use this reduced density threshold, or alternately a flow convergence criteria, ∇ · v < 0, for the velocity field of the gas cells (Grudić et al. 2018a) instead of a density threshold criterion, along with the treatment of stellar mass-loss that is modified from our fiducial simulations, we find no bound clusters. These results strongly indicate that the treatment of unresolved |$P\, dV$| work from stellar mass-loss or other forms of unresolved early feedback is crucial in resolving the formation of star clusters in progenitors of z = 0 dwarf galaxies, while the choice of density threshold (if any) is not. Preliminary work using FIRE-3 galaxy formation physics (Hopkins et al. 2023), which follows the treatment of unresolved |$P\, \mathrm{ d}V$| work adopted here and does not adopt a density threshold for star formation – it requires star-forming gas to be self-gravitating, Jeans unstable, and within a converging flow – supports this conclusion (Sameie et al., in preparation).
4.3 Comparison with previous results
While our work covers a different mass regime from what has previously been studied in full cosmological simulations that resolve the formation of GCCs (e.g. Kimm et al. 2016; Kim et al. 2018; Ma et al. 2020), our results are broadly consistent with those from these previous numerical simulation and analytic arguments: GC candidates form preferentially at early cosmic times when the turbulent ISM of gas-rich galaxies can provide the requisite high pressures that are conducive to cluster formation.
Our results differ somewhat from those of semi-analytic models that aim to understand GC formation within the broader context of galaxy evolution across cosmic time. These models typically predict that GCs in z = 0 dwarf galaxies form at relatively late times (z ≲ 3), in part because the average metallicities of such galaxies are not predicted to reach the levels seen in most GCs until that point (Choksi, Gnedin & Li 2018; El-Badry et al. 2019; Kruijssen 2019; Reina-Campos et al. 2019). The clusters formed in our simulations form significantly earlier than this, at z ≳ 5. It is noteworthy that the GCCs found here have metallicities that are always at least as high as that of their host galaxies; in some cases, the GCCs are one full dex higher in metallicity. Allowing for GCs to have higher metallicities than the mean gas-phase metallicity of their progenitors could be an interesting and fruitful path forward for semi-analytic models of GC formation.
Although we are unable to follow the evolution of the GCCs to z = 0 and therefore cannot make definitive statements about whether these objects are truly GC analogues, many of their properties are grossly consistent with ancient GCs that are observed in some low-mass dwarf galaxies today. One clear difference from observations is the metallicity of the clusters in our simulations, which are slightly lower than metallicities of well-quantified GCs observed in dwarf galaxies (Beasley et al. 2019) and have larger iron metallicity spreads (0.15–0.2 dex; the metallicity distributions are reasonably well approximated by Gaussians with σ = 0.15 dex across the simulation suite) than are observed (0.045 dex; Carretta et al. 2009; Bailin 2019). This could be an indication that the clusters formed here would be disrupted and form the streams that are known to have lower metallicities (Martin et al. 2022), or it could mean that subtle aspects of the galaxy formation modelling in FIRE require refinements in this regime. Indeed, Wheeler et al. (2019) noted that low-mass dwarf galaxies from FIRE-2 lie slightly but systematically below the observed mass–metallicity relationship. An explicit treatment of Population III star formation and enrichment is likely to be important for understanding the properties of the lowest metallicity systems and their connections to GCs (see also Schauer et al. 2021).
More broadly, the existence of galaxies such as Eridanus II, which formed 80 per cent of its stars before z = 6 and hosts a GC with a stellar population indistinguishable from that of the galaxy (Simon et al. 2021), demonstrates that GC formation in very low-mass systems at high redshifts (z > 6) is possible and must be accounted for in models. The properties of the GGs and galaxies in our simulated dwarf galaxy haloes are in many ways similar to Eridanus II, with GC formation accompanying (or even preceding) the formation of the bulk of the stars in the galaxy, though our systems are somewhat higher in stellar mass. More detailed observational studies of GCs in dwarf galaxies, coupled with future simulations of such systems, hold the promise to reveal important aspects of star formation in metal-poor systems in the reionization era.
5 CONCLUSIONS
The recent discoveries of ancient GCs in low-mass (|$M_{\star }\sim 10^5\!-\!10^7\, \mathrm{M}_{\odot }$|) Local Group dwarf galaxies and of star cluster candidates in formation in the high-redshift Universe has reignited interest in a number of questions related to GCs in dwarf galaxies. We have used cosmological hydrodynamic simulations of seven haloes with virial masses of |$M_{\rm vir}(z=0)\sim 10^{10}\, \text{M}_\odot$| from the FIRE-2 project to investigate star cluster formation in the ancestors of present-day dwarf galaxies. We find that star cluster formation at high redshift (11 ≳ z ≳ 5) is indeed common in these systems, which is perhaps the most important high-level result from our study. In more detail, our principal conclusions include the following:
Relatively massive (|$M_{\star }\in [0.5\!-\!5]\times 10^5\, \mathrm{M}_{\odot }$|) and compact (|$6\, {\rm pc} \lesssim r_{1/2} \lesssim 30\, {\rm pc}$|) clusters form in haloes with virial masses of |$(0.5\!-\!2)\times 10^9\, \mathrm{M}_{\odot }$| – roughly a factor of 10 more massive than the atomic cooling threshold corresponding to |$T_{\rm vir}=10^4\, {\rm K}$| – in the redshift range 11 ≳ z ≳ 5. Of the seven systems studied here, five form one such cluster, while one forms two clusters and one forms no clusters.
The clusters form when |${\sim }10^7\, \mathrm{M}_{\odot }$| of dense and turbulent gas reaches a surface density of |$\Sigma _{\rm thresh} = 10^4\, \mathrm{M}_{\odot }\, {\rm pc}^{-2}$|. These conditions occur because of compressive shocks from nearby star formation, cloud–cloud collisions, or both.
Once the requisite conditions are achieved, cluster formation happens rapidly. Stellar feedback clears all gas from the cluster region by 10 Myr from the start of cluster star formation. In some cases, this feedback leads to shocks that trigger nearly coeval star formation hundreds of parsecs away from the cluster.
The clusters and star formation associated with the GMCs from which the clusters form exceed the pre-existing stellar content of the halo, meaning the clusters form before there is a well-defined host galaxy at these mass scales. The GCCs studied here therefore originate from a phase of galaxy formation that either predates or accompanies the star formation constituting the bulk of the host galaxy.
In several cases, the clusters constitute |${\sim } 10~{{\ \rm per\ cent}}$| of the total stellar mass in the host halo at z = 3. If the clusters were able to survive numerical heating and disruption, this would make them consistent with observations of GC-hosting dwarfs in the local Universe, where cluster stars typically constitute 1–10 per cent of the stellar mass of the galaxy (Georgiev et al. 2010; Hudson et al. 2014; Larsen 2017).
The clusters in our simulations form at higher redshift (11 ≳ z ≳ 5) than is predicted at this mass scale in many semi-analytic models of cluster formation (z ≲ 3; e.g. Choksi et al. 2018; El-Badry et al. 2019; Kruijssen 2019), in part because the models typically tie the mean gas-phase metallicity of a forming cluster to that of its host galaxy. By contrast, clusters in the simulations presented here form out of gas with higher-than-average metallicity (averaging over the gas in the halo at the cluster formation epoch) and often before a well-defined galaxy is even present.
The clusters live for tens to hundreds of dynamical times (0.2–2.5 Gyr), with clusters born far from the dynamical centre of the halo surviving the longest. The disruption of the clusters comes from a combination of physical and numerical effects; were it possible to accurately resolve the internal dynamics of the clusters, several might well survive to z = 0.
The clusters typically form outside of the half-light radius of any pre-existing galaxy (insofar as any such galaxy exists), well removed from the centre of the host halo. In some cases, the cluster formation sites are at |$0.25\!-\!0.5\, r_{\rm vir}$|. These offsets are natural since clusters are forming in regions of dense gas having turbulent velocity dispersion comparable to the halo virial velocity, which makes the cluster formation sites somewhat stochastic.
Given the formation sites, the clusters formed here are never deeply enshrouded in the centres of their own dark matter haloes, though they all form within the virial radius of a |$M_{\rm vir}\sim 10^9\, \mathrm{M}_{\odot }$| dark matter halo. The contribution of dark matter to the clusters’ mass profiles is minimal and consistent with the background halo density at the cluster formation location (typically hundreds of pc, or |${\sim }0.2\, r_{\rm vir}$|, from the halo centre).
Properties of these clusters are consistent with objects detected in HST observations at z ∼ 6–8 in lensing fields. They are also suggestive of similarities to clusters observed in Local Group dwarfs at z = 0.
The treatment of unresolved thermal energy deposition from stellar mass-loss is a primary numerical ambiguity affecting cluster formation physics in our suite.
While results from these ‘proof-of-principle’ simulations are both encouraging and intriguing, there are multiple avenues for improvement in the near term. Our choice of softening scale for baryonic particles, |$2\!-\!5\, {\rm pc}$|, leads to artificially large sizes for our clusters; running versions of these simulations with softenings that are roughly 10 times smaller should allow us to test whether the sizes of the clusters we form are realistic. A more challenging issue is faithfully tracing the internal dynamics of the clusters for cosmological times (∼10 or more Gyr). Hybrid numerical schemes that resolve cluster formation in cosmological simulations and then track cluster evolution with methods capable of tracking collisional stellar dynamics within a larger scale galactic environment (e.g. Rodriguez et al. 2023) are promising in this regard. Finally, the galaxy formation models employed here, which are part of the FIRE-2 suite, are subject to further refinement. Details of the treatment of various aspects of star formation will likely be important for accurately capturing the formation of bound, long-lived star clusters; updates incorporated into FIRE-3 (Hopkins et al. 2023) are a starting point in this direction. In the near future, however, it will be possible to run the kinds of simulations presented here with numerical parameters that allow us to form pc-scale clusters and to follow the evolution of the clusters and their host galaxies from birth in the high-redshift Universe to present day. Such simulations will be crucial for using JWST data to constrain models of the formation and evolution of GCs (Forbes et al. 2018; Renaud 2018; Adamo et al. 2020).
ACKNOWLEDGEMENTS
We thank the anonymous referee for a thorough and helpful report. MBK acknowledges support from NSF CAREER award AST-1752913, NSF grants AST-1910346 and AST-2108962, NASA grant NNX17AG29G, and HST-AR-15006, HST-AR-15809, HST-GO-15658, HST-GO-15901, HST-GO-15902, HST-AR-16159, and HST-GO-16226 from the Space Telescope Science Institute (STScI), which is operated by AURA, Inc., under NASA contract NAS5-26555. Support for PFH was provided by NSF Research Grants 1911233, 20009234, 2108318, NSF CAREER grant 1455342, and NASA grants 80NSSC18K0562 and HST-AR-15800. AW received support from: NSF grants CAREER 2045928 and 2107772; NASA ATP grant 80NSSC20K0513; HST grants AR-15809 and GO-15902 from STScI; a Scialog Award from the Heising-Simons Foundation; and a Hellman Fellowship. JSB was supported by NSF grant AST-1910346. EQ was supported in part by a Simons Investigator grant from the Simons Foundation and NSF AST grant 2107872. JS was supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2102729. DRW acknowledges support from HST-GO-15476, HST-GO-15901, HST-GO-15902, HST-AR-16159, and HST-GO-16226 from STScI.
This work used the Extreme Science and Engineering Discovery Environment (XSEDE; Towns et al. 2014), via allocation AST140080, and the Frontera computing project at the Texas Advanced Computing Center (via allocations AST21010 and AST20016), which are supported by National Science Foundation awards ACI-1548562 and OAC-1818253, respectively. The analysis in this paper is carried out by python packages numpy (Harris et al. 2020), matplotlib (Hunter 2007), scipy (Virtanen et al. 2020), and h5py (Collette 2013).
DATA AVAILABILITY
The data supporting the plots within this article are available on reasonable request to the corresponding author.
Footnotes
All ICs can be found at http://www.tapir.caltech.edu/~phopkins/publicICs.
Archaeological SFHs are computed by considering the formation times of all of the stars in the galaxy at z = 3 as opposed to summing the instantaneous star formation rate in any individual progenitor.
