Formation of globular cluster candidates in merging proto-galaxies at high redshift: a view from the FIRE cosmological simulations

Abstract

Using a state-of-the-art cosmological simulation of merging proto-galaxies at high redshift from the FIRE project, with explicit treatments of star formation and stellar feedback in the interstellar medium, we investigate the formation of star clusters and examine one of the formation hypotheses of present-day metal-poor globular clusters. We find that frequent mergers in high-redshift proto-galaxies could provide a fertile environment to produce long-lasting bound star clusters. The violent merger event disturbs the gravitational potential and pushes a large gas mass of ≳ 10^5–6 M_⊙ collectively to high density, at which point it rapidly turns into stars before stellar feedback can stop star formation. The high dynamic range of the reported simulation is critical in realizing such dense star-forming clouds with a small dynamical time-scale, t_ff ≲ 3 Myr, shorter than most stellar feedback time-scales. Our simulation then allows us to trace how clusters could become virialized and tightly bound to survive for up to ∼420 Myr till the end of the simulation. Because the cluster's tightly bound core was formed in one short burst, and the nearby older stars originally grouped with the cluster tend to be preferentially removed, at the end of the simulation the cluster has a small age spread.

1 INTRODUCTION

By the interaction of gravity and pressure, gas becomes unstable and collapses to turn into stars. Observationally, we know that most stars form in ‘clustered fashion’, with some stars forming en masse in dense molecular clouds of ≳ 10⁴ M_⊙ and others in relatively loose associations (for reviews, see Portegies Zwart, McMillan & Gieles 2010 and Bressert et al. 2010). The fraction of stars that form in bound clusters, often referred to as Γ, is affected by the host galaxy environment, with galaxies with higher star formation activity such as starbursts having higher Γs than local spiral and dwarf galaxies (e.g. Kruijssen 2012; Adamo & Bastian 2015). Simulations of an isolated star clusters in an idealized setup have a long history and can now resolve the formation of individual stars (for recent reviews, see Kruijssen 2013 and Dale 2015). However, they cannot be used to study the formation of a dense environment in a galactic context, in which those star clusters spawn. It is also difficult to utilize this type of calculations in studying how bound versus unbound star clusters form, and how they evolve and survive in a galactic environment. A high-resolution, galaxy-wide simulation is required to address these questions.

However, star clusters have rarely been simulated in a galactic context, with numerical accuracy high enough to resolve their dynamical evolution. Note that in most previous galaxy simulations the mass of each resolution element approximately equates to that of a star cluster. Only after we reach a mass resolution of ≲ 10³ M_⊙ can we begin to reliably resolve the kinematics amongst the member star particles of a massive ∼10^5 − 6 M_⊙ cluster in a galaxy-scale simulation.^1,2 It is also required that realistic physics models be in place which describe star formation and stellar feedback at such high resolution. Most previous galaxy-scale simulations lack such sophisticated physics models.

Recently, however, this is becoming increasingly possible with modern cosmological simulations. In this paper, we use simulations from the FIRE project (Feedback In Realistic Environments; Hopkins et al. 2014)³ to study how star clusters form and evolve in high-redshift merging proto-galaxies. This has been proposed (Whitmore & Schweizer 1995; Schweizer et al. 1996; Schweizer 2006) as a formation channel for the present-day old, metal-poor ‘blue’ globular clusters (GCs) which tend to be distributed throughout the galactic halo (as opposed to metal-rich ‘red’ GCs which are found mostly in galactic bulges). For reviews of their properties and various hypotheses of their origin, see Peebles & Dicke (1968), Harris & Racine (1979), Fall & Rees (1985), Harris (2001), and Brodie & Strader (2006). Our investigation is complementary to other recent studies which have explored alternative GC formation scenarios: for example, in binary merging galaxies at low redshift (e.g. Li, Mac Low & Klessen 2004; Bournaud, Duc & Emsellem 2008; Kruijssen et al. 2011, 2012; Renaud, Bournaud & Duc 2015; Maji et al. 2017), or in primordial mini-haloes at redshift z ≳ 10 (e.g. Kimm et al. 2016; Ricotti, Parry & Gnedin 2016).

Specifically, we carry out an investigation using a cosmological ‘zoom-in’ simulation from Ma et al. (2015) at z ∼ 6, with ∼ pc-scale resolution, run with the FIRE physics that includes star formation only in dense, self-gravitating gas and stellar feedback from supernovae (SNe), stellar winds, photoionization and photoelectric heating, and radiation pressure. Our setup allows us to explicitly follow the formation and evolution of at least the most massive star clusters. We show that, in the simulation, most stars form in unbound associations but some form in resolved, bound clusters. Frequent mergers in high-redshift proto-galaxies provide a fertile environment to produce the latter population by pushing large gas masses (≳ 10^5 − 6 M_⊙) collectively to high density, at which point it turns into stars before stellar feedback can disrupt the clouds. We explore subsequent dynamical evolution, virialization, and the age and metallicity spreads of the resulting clusters.

The remainder of this article is organized as follows. In Section 2 we detail the simulation code and methods. Section 3 presents results of the simulation focusing on the formation and evolution of bound clusters. Finally in Section 4 we summarize our findings and conclusions.

2 METHODOLOGY

The simulation studied in this work is one of a suite of high-resolution cosmological simulations at high redshift from the FIRE project, and is presented and described in detail in Ma et al. (2015, z5m10h run therein). We briefly review its important features for completeness. The simulation was run using the Gizmo code (Hopkins 2015)⁴ which solves gravity using a tree-particle mesh (TreePM) method with fully adaptive gravitational softenings (scaled to the inter-particle separation), and solves the hydrodynamics using the Lagrangian ‘pressure-energy’ (P-sph) formulation of smoothed particle hydrodynamics with various improvements to alleviate known issues with fluid mixing and shock-capturing in older SPH formulations (see Hopkins 2013; Hopkins et al. 2014). The simulation follows a high-resolution Lagrangian region around a ‘target halo’ of virial mass ∼10¹⁰ M_⊙ (stellar mass of the ‘target galaxy’ ∼4.7 × 10⁷ M_⊙) at z ∼ 6, with fixed mass resolution 800 h⁻¹M_⊙ and minimum force resolution 1.4 h⁻¹ pc (proper) for gas (and fixed to 1.4 h⁻¹ pc for star particles). This resolution is high enough to resolve the dynamics between member particles of a relatively massive star cluster studied in this paper (but see also footnote 1). At z = 5 of the original Ma et al. (2015) simulation, we identify a group of long-lasting bound clusters, one of which we hereafter call the cluster ‘A’ (see Figs 1–3). We specifically re-ran the interval from z = 7 to z = 5 of the original simulation to generate more snapshots at finer intervals during ∼80 Myr before and after the cluster ‘A’ formed, in order to study when, where, and how the cluster forms.⁵

Baryonic physics is treated using the FIRE-1 model, described in detail in Hopkins et al. (2014). Briefly, radiative heating and cooling in 10 − 10¹⁰ K take into account molecular, atomic, ionized, and metal-line processes (with 11 independently tracked species), and include photoheating by local sources and redshift-dependent ultraviolet background (Faucher-Giguère et al. 2009), and self-shielding. Star formation occurs via spawning of star particles only from gas which meets a series of criteria: it must be locally self-gravitating and Jeans-unstable according to a sink particle criterion (including both thermal and turbulent support; see Hopkins, Narayanan & Murray 2013), molecular (following Krumholz & Gnedin 2011), and denser than a star formation threshold n_th = 500 cm⁻³.⁶ A new star particle inherits the mass, metallicity, gravitational softening length, and particle ID number from its progenitor gas particle. Once a star particle forms, the simulation explicitly tracks feedback from (1) local and long-range radiation pressure (including single-scattering and multiple-scattering of re-radiated infrared photons), (2) energy, momentum, mass and metal injection from SNe (Types Ia and II) and stellar mass-loss (OB and AGB-star winds), and (3) photoionization and photoelectric heating. The rates for each channel are calculated as a function of the stars’ age and metallicity using a stellar population synthesis model Starburst99 (Leitherer et al. 1999) assuming the Kroupa (2002) initial mass function (IMF).

For more details on how each item above is implemented, we refer the interested readers to section 2 of Hopkins, Quataert & Murray (2012), and section 3 and appendix A of Hopkins et al. (2014). A series of FIRE simulations using cosmological and isolated initial conditions have reported reasonable star formation histories, stellar mass--halo mass relation, Kennicutt–Schmidt relation (Hopkins et al. 2014), mass--metallicity relation (Ma et al. 2016), multiphase interstellar medium (ISM; Hopkins, Quataert & Murray 2011), galactic outflows (Muratov et al. 2015), dense neutral hydrogen content of galactic haloes (Faucher-Giguère et al. 2015, 2016), galaxy structures and metallicity gradients (Ma et al. 2017), among others. These studies validate the ‘realistic’ baryonic physics we adopt in the present study.

3 RESULTS

We start by giving an overview of the simulation result and the two distinct populations of star clusters. Then we discuss the formation, evolution, and composition of long-lived bound star clusters.

3.1 Overview of simulation results: from the host galaxy to the target cluster

Fig. 1 shows the z = 5 snapshot of the target galaxy and the star cluster of interest studied in the present paper (cluster ‘A’ defined in Section 2). Images are rendered as in the way described in Hopkins et al. (2005, 2014). In the left panel, molecular cloud complexes in cold filaments in the upper side of the panel indicate fresh flows of dense gas that has not turned into stars yet. The right panel shows multiple star clusters of masses ranging from 10⁵ to 10⁷ M_⊙ scattered around the galaxy. Each of these clusters is dynamically resolved with 10² − 10³ star particles within 10 − 10² pc half-mass radii (Fig. 4). Among them, marked with a white circle of radius 300 pc is cluster ‘A’ that forms at z = 6.92 and survives ∼420 Myr afterwards.

Fig. 2 reveals the movement of cluster ‘A’ throughout its host galaxy.^7,8 The second column captures the moment right before most star particles in the cluster ‘A’ form (red, young star particles inside the black circle) during a major proto-galaxy merger. After the cluster forms, it moves to an orbit with a relatively large radius, piercing through the galactic nucleus at z = 6.30 (fifth column), but later reaching far out into the extended galactic halo at z = 5.15 (sixth column).

To closely inspect the formation site of the cluster ’A’, in Fig. 3 we zoom in on the 2 kpc × 2 kpc region around the cluster ’A’. Before the cluster ’A’ forms, we see a dense gas clump forming in a galactic merger event (first and second columns; Section 3.4). This efficiently forms new star particles in a short time, which then violently disperse the gas via stellar feedback, but only after a large number of star particles have formed (third and fourth columns). Later, we see the cluster pass through the centre of the host galaxy, triggering a mass loss of the cluster (fifth column; Section 3.3).

3.2 Two populations of star clusters

We now examine the population of star clusters formed in our simulation between z = 7 and z = 5. To identify star clusters in the simulated galaxy, we utilize the six-dimensional phase-space halo finder Rockstar (Behroozi, Wechsler & Wu 2013)⁹ modified to run on star particles instead of dark matter particles.¹⁰ To ensure that the identified star clusters are dynamically resolved and are not numerical artefacts, we only consider clusters with more than 100 member particles (or >8.74 × 10⁴ M_⊙). We then directly extract the stellar mass, M_cl, and spherical half-mass radius, R_h, of each cluster.¹¹

In Fig. 4 we plot the M_cl–R_h relation of all star clusters identified in 76 snapshots finely spaced between z = 7 and z = 5. We plot only those that survived at least one output interval (e.g. Δ t_out ∼ 1.6 Myr at 6.5 < z < 7; see footnote 5). In other words, we plot clusters that were found in two or more snapshots to eliminate spurious transient associations. Each data point is then coloured by the fraction of member particle masses that are bound to the cluster, i.e. bound mass fraction f_bd derived from the standard Rockstar output (but see footnote 11). Additionally, four star clusters of mass above 10^5.5 M_⊙ that survived more than 300 Myr are annotated with blue circles. These long-lasting, bound star clusters are identified at z = 5 by considering the mass M_cl, bound mass fraction f_bd, and lineage information (see Section 3.3). These clusters are grouped with few gas particles and almost no dark matter particles, thus not associated with dark matter over-densities or galactic subhaloes. Among them, the cluster ’A’, studied extensively here, is marked with a thicker blue circle at z = 5. The force resolution of the simulation is well below the y-axis range (see Section 2), indicating that clusters shown here are not affected by force softening, although they may of course be affected by the finite particle number (i.e. mass resolution).

In this figure various features are prominent. First, a large fraction of the clusters congregate around M_cl ∼ 10^{5.0 − 5.5} M_⊙ and R_h ∼ 10^{1.5 − 2.5} pc, giving mean stellar densities, defined as |$\rho _{\rm cl} \equiv 3M_{\rm cl}/(8\pi R_{\rm h}^3)$|⁠, of ∼10⁻² M_⊙ pc⁻³. The majority of these clusters are loose ’unbound associations’ with low f_bd. These associations, according to our lineage tree analysis (Section 3.3), often dissolve into an extended stellar disc or bulge in 10–30 Myr by rapidly losing their unbound member particles. This behaviour is expected since most of the unbound associations inherit the properties of normal molecular clouds that are only marginally gravitationally bound after turning just a few per cent of the cloud mass into stars. Indeed, their observational counterparts such as OB associations tend to drift apart in a ≳ 10 Myr time-scale. Readers should note that the density of a cluster/association (or its gas progenitor) can be much less than the star formation threshold (500 cm⁻³ or 12 M_⊙ pc⁻³), because stars form in a small self-gravitating sub-clumps above the threshold density, distributed widely over a lower-mean-density giant molecular cloud (GMC) complex.

A much smaller fraction of the stellar mass is represented by self-gravitating bound clusters that often survive for a longer time. For example, the long-lasting bound clusters identified at z = 5 (blue circles in Fig. 4 and their ancestors) have radii as small as R_h ∼ 10 pc with mean stellar densities as large as ∼10² M_⊙ pc⁻³. Their mean densities are often comparable to or larger than the star formation threshold, and their typical bound mass fractions are higher than 80 per cent (red face colour). In addition, they tend to maintain their locations on the M_cl − R_h plane for ≳ 100 Myr (as indicated by groups of red points crowded around the four blue circles). This second group of star clusters could be categorized as ’bound clusters’, which could potentially become candidates for present-day metal-poor ’blue’ GCs. We note that these numerically formed clusters are similar in their sizes and masses to the objects observed at z ≳ 6 by Vanzella et al. (2017) and conjectured as proto-GCs.

To illustrate the two populations of star clusters, Fig. 5 shows the normalized probability distribution function (PDF) of mean cluster stellar densities. We see a clear bimodal distribution with the majority of mostly unbound associations forming one peak at ∼10⁻² M_⊙ pc⁻³, while bound clusters forming another peak at ∼30 M_⊙ pc⁻³. The corresponding surface density for this second bimodal group is ∼10³ M_⊙ pc⁻², as can be seen in the surface density PDF of Fig. 6. It is about an order-of-magnitude larger than that of a typical GMC (∼10² M_⊙ pc⁻²), implying a rather unusual dynamical process by which these bound clusters may have formed and evolved. The bound cluster's surface density also corresponds to the surface density scale at which star formation over Myr time-scales may become highly efficient despite strong stellar feedback (see Section 3.4 and Grudić et al. 2016). In fact, most of these dense objects can be traced through multiple snapshots and are ancestors of just a few bound clusters (see Section 3.3). Note that the second bimodal peak around ∼30 M_⊙ pc⁻³ in Fig. 5 corresponds well with the maximum gas density when the cluster ’A’ formed (red dotted line at z = 6.92). This high gas density was realized when two proto-galaxies merge as can be seen in the second column of Fig. 3, but typically not before or after the merger (black dotted line in Fig. 5 at z = 6.99; representative gas density PDF for a ’normal’ galaxy). We will come back to this discussion in Section 3.4.

Figs 5 and 6 indicate that these bound clusters are in a separate group of a distinctive evolutionary process, not simply a tail of a single, unimodal distribution. Of course, because we adopted a finite threshold density for star formation, a patch of gas could have collapsed further – if we had infinite resolution – before forming stars. So the densities of the bound clusters seen here are probably lower limits. This only strengthens our inference of a bimodal population.

3.3 Evolution of a long-lasting bound star cluster

3.3.1 Method and overview

We now focus on how the long-lasting bound star clusters evolve in time. Before moving to discuss the simulation result, a remark on numerical resolution and the ’intra-cluster’ dynamics would be timely. The long-lasting bound clusters such as cluster ’A’ are resolved in mass (with ≳ 10³ particles), and the force resolution is also well below their half-mass radii. Therefore, the kinematics amongst the cluster member star particles is numerically resolved. However, this does not mean that our simulation reproduces the phase-space distribution of the actual stars in observed clusters. Partly due to the mass resolution – each particle represents a stellar population of ∼10³ M_⊙ – we simply cannot resolve a variety of internal, ’intra-cluster’ evolution processes such as stellar two-body relaxation, mass segregation, evaporation, stellar mass loss, or binary interactions (for a comprehensive review, see e.g. Benacquista 2006; Portegies Zwart et al. 2010). Thus, it would be prudent to assume that the reported simulation at best marginally resolves the detailed intra-cluster evolution. With this caveat in mind, in all subsequent analyses, we focus only on globally averaged characteristics of a star cluster, but not on e.g. radial profiles of stellar density or metallicity inside a cluster, which are more severely affected by the intra-cluster evolution.¹²

In order to trace star clusters evolving in time, we combine the descendant information in standard Rockstar outputs, and supplementary information such as member particle IDs of star clusters selected at z = 5. The selected cluster's ’lineage tree’ – or the main branch of a merger tree – is built by linking the ’main ancestors’ determined by Rockstar. (When a cluster P bequeaths the most particles to a cluster D in the next snapshot, P is defined as D's parent, D as P's descendant. D may have many parents, the most massive of which is defined as the main ancestor.)¹³ We then draw the lineage tree of the cluster ’A’ in Fig. 7. The ancestor of the cluster ’A’ starts out in the upper right corner, then loses its mass at early times to dramatically evolve to a compact – i.e. smaller R_h and higher mean density – and tightly bound – i.e. higher bound mass fraction – cluster at z = 5. The following section describes this evolution in detail.

3.3.2 Three phases of cluster evolution

For a more quantitative evaluation of the simulation results, evolution of the cluster's various characteristics is shown in thick black solid lines in Figs 8–11. For example, Fig. 8 depicts the cluster's mass M_cl, bound mass fraction f_bd, and virial ratio E_kin/|E_pot|. As in other figures, 0 Myr corresponds to the moment the cluster ’A’ forms. The grey shaded region represents the ∼80 Myr period during which the simulation produced outputs at finer intervals (Section 2).

From Figs 7 to 11, three phases of evolution are noticeable:

(1) Phase 1: formation and initial evolution of the cluster. In the first ∼100 Myr of its life (blue to first red circle in Fig. 7), the cluster's mass and radius frequently change (Figs 8 and 9), partly because of its proximity to its host galaxy and neighbouring clusters. In other words, according to our lineage tree analysis in the first ∼100 Myr, the cluster's mass and size change because it loses its loosely bound member particles, or it captures neighbouring particles as it passes by them (but negligible in situ star formation; see Section 3.5). The cluster gradually becomes tightly bound with f_bd increasing from <10 per cent to ∼80 per cent, while its virial ratio asymptotes to a virial equilibrium (Fig. 8).¹⁴ However, as the bottom panel of Fig. 9 demonstrates, the cluster's core is already very dense at its formation. The evolution in mass, size, and f_bd is primarily due to changes among loosely bound old stars initially associated with the newly formed cluster.¹⁵

(2) Phase 2: mass loss of the cluster. At ∼100 Myr after the cluster's formation (between z = 6.4 and z = 6.3; two red circles in Fig. 7), the cluster's mass decreases by >70 per cent (Fig. 8). This event refers to the period when the cluster pierces through the host galaxy's gravitational centre (fifth columns in Figs 2 and 3), a pericentric approach to the Galactic centre since the cluster's formation with a very high relative velocity (Fig. 10). Strong tidal shocking at the pericentre and tidal interaction with the host galaxy remove a large fraction of the cluster's mass in a short time, leaving the most tightly bound core of the clusters, and bringing its bound mass fraction from ∼80 per cent to >90 per cent (Fig. 8). It showcases that the energy gain by a tidal shock can alter the cluster's binding energy (e.g. Kruijssen et al. 2011). This process preferentially strips the old star particles which existed before the formation of the cluster itself, but were associated with the cluster by Rockstar when the cluster formed. The process thus leaves the core of tightly bound young stars with narrower distribution in particle ages (Sections 3.4 and 3.5). We speculate that a tidal shock-induced process like this may be essential to transform a population of star clusters to leave only tightly bound clusters (e.g. Kruijssen 2015). It is worth noting that the cluster's evolution is driven not only by dynamic relaxation effects of the cluster itself, but also by the embedding galactic contexts such as the cluster passing through strong tidal fields and experiencing mass losses.

(3) Phase 3: settling of the cluster. After the mass loss event, the cluster eventually settles in at M_cl ∼ 10^6.0 M_⊙ and R_h ∼ 10^1.4 pc. The cluster maintains these properties for the next ∼320 Myr. It also preserves most other characteristics such as the bound mass fraction, virial ratio, density, metallicity, and velocity dispersion (Figs 8, 9 and 11). One of the reasons that the cluster could maintain its characteristics is because the cluster has now moved to an orbit with a relatively large radius after a fast velocity increase around the major mass loss event in Phase 2 (Fig. 10). Our observation is broadly consistent with Kruijssen (2015) who argued that external perturbation events such as galaxy mergers could cause clusters to migrate into the halo, and thus limit the duration of the shock-induced disruption phase. The cluster may survive for an extended period of time since it is now less likely to undergo tidal shocks that could disrupt it. It is also interesting that the cluster is ejected from the main concentration of stars into the halo, where Milky Way's metal-poor ’blue’ GCs are indeed found.¹⁶

3.4 Formation of a long-lasting bound star cluster

3.4.1 Method and overview

In this section, we go back to the moment when these long-lasting bound star clusters form, and examine the formation conditions of the bound clusters. In particle-based hydrodynamics codes like Gizmo a newly born star particle inherits the unique particle ID number from its progenitor gas particle. Using these unique IDs, we can identify the ’gas progenitor’ of a cluster in an earlier snapshot before the cluster formed. The cluster's ’formation time’ is defined as the moment when 50 per cent of the gas progenitor member particles have turned into star particles. Using such information, in Figs 9 to 11 we present the evolution of the cluster ’A’ and its gas progenitor. In addition to the gas progenitor of the cluster ’A’, progenitors of three other long-lasting bound clusters identified at z = 5 are shown (blue dashed lines in the bottom panels; recall that these three clusters are more massive than 10^5.5 M_⊙ and survived > 300 Myr at z = 5; marked by thin blue circles in Fig. 4). Also shown are the progenitors of nine unbound associations at z = 5 – that are more massive than 1.5 × 10⁵ M_⊙ at z = 5, but fail to meet the criteria of a long-lasting bound cluster (thin red lines).

At early times, more than 100 Myr before the formation, the gas progenitor particles of the cluster ’A’ are spread over ≳ kpc, effectively a random sample of the entire galaxy (Fig. 9). About 30 Myr before the formation time, the gas progenitor's mean density reaches ∼10 cm⁻³ (or ∼10⁻¹ M_⊙ pc⁻³), a typical value of a dense ISM in this relatively dense high-redshift galaxy. Then, for the next 30 Myr, a galaxy merger event disturbs the gravitational potential around the gas progenitor sitting between the two merging proto-galaxies (first and second columns in Figs 2 and 3). This event forces the gas cloud into strong compressive shocks (Fig. 10), forming a self-gravitating cloud with its velocity dispersion decreased to ≲ 20 km s⁻¹ (Fig. 11). The mean density of the gas progenitor accordingly increases quickly to ≳ 10³ cm⁻³ (or ≳ 10 M_⊙ pc⁻³; see Fig. 9, or the red dotted line of Fig. 5). Consequently, the gas cloud collapses and turns a significant portion of its mass (≳ 10^5 − 6 M_⊙, comparable to that of a star cluster) simultaneously into stars in the gas free-fall time-scale, t_ff ≲ 3 Myr (see also Section 3.5). Due to its high density and compactness, the newly formed cluster remains gravitationally bound. Note that this formation scenario applies also to other long-lasting bound star clusters.

In order to illustrate the nature of gas compression triggering star cluster formation, in Fig. 12 we plot the radial profile of the mass-weighted gas infall velocity when the cluster ’A’ starts to form. The plot is centred on the gas progenitor of the cluster – the location of densest gas within 100 pc from the centre of the gas progenitor – at 10.5 Myr before the cluster's formation time (when the first stars in the cluster begin to form). The negative velocity values indicate the inward motion of the gas. The merger event disturbs the gravitational potential around the progenitor and compresses the gas cloud to high density where it is strongly self-gravitating. Therefore, the actual gas infall speed is of the order of the purely gravitational infall speed at each radius R, |$v_{\rm ff} (R) = \sqrt{8GM_{\rm tot,{\, }<}(R)/(\pi ^2R)}$|⁠, where M_{tot,  <}(R) is the total mass within R. The gas infall speed does occasionally become larger than the purely gravitational collapse speed. This is indicative of a massive gas cloud with little support from rotation or shear, which collapses under self-gravity until it undergoes efficient star formation at high density.

3.4.2 Merger-induced star cluster formation

Overall, our investigation paints a picture of bound star cluster formation when massive clouds reach very high densities – a condition which may be preferentially merger-induced (though in principle similar conditions could occur without mergers). Mergers can produce tidal shocks in regions of strong tidal acceleration (e.g. bridges, tidal tails, induced bars) with velocity jumps of the order of the circular velocity (Fig. 10), corresponding to Mach numbers of ∼50–100, which in turn produces extreme gas compression (e.g. Kim, Wise & Abel 2009). The resulting dense gas clumps efficiently cool and collapse. The free-fall time at the velocities shown in Fig. 12 reaches ∼0.5–1.5 Myr for the clump core with densities ∼10^3 − 4 cm⁻³. Therefore collapse occurs before the first SNe explode (which requires a time ∼3–10 Myr). Although radiative feedback and stellar winds are present in our simulation and act on shorter time-scales, the gas densities are sufficiently high such that they are unable to disrupt the cloud until a significant fraction – an order unity – of the gas has been converted into stellar mass (Grudić et al. 2016).¹⁷ The resulting, highly efficient star formation allows the cluster to remain bound. This is in line with previous analytic (Fall, Krumholz & Matzner 2010; Kruijssen 2012) and numerical (Grudić et al. 2016) models for the fraction of stars that form in bound clusters.

By disturbing and compressing the system, galaxy mergers provide unique opportunities to make the otherwise normal ISM shock to abnormally high densities. The violent environment is important to push a large gas mass of ≳ 10^5–6 M_⊙collectively to very high density in a short time, rather than letting it fragment into smaller dense pieces. Once at such high density, the cloud promptly turns a significant portion of itself into stars at a small dynamical time-scale, t_ff ≲ 3 Myr. Since proto-galaxies at high redshift frequently experience mergers, they could have provided a fertile environment to produce multiple long-lasting bound clusters – which we speculate to be candidates for present-day ’blue’ GCs.

We emphasize the critical difference between our merger-induced formation process of long-lasting bound star clusters and the normal turbulent fragmentation. All star particles require a density above n_th = 500 cm⁻³ to spawn in our simulation (Section 2), yet most of them end up in unbound associations. It is because, in a typical star-forming gas clump, only a small fraction of its mass in sub-clumps (with mass ≪10⁵ M_⊙) is sufficiently dense to form star particles. The newly born star particles then destroy their parent cloud complex. By contrast, in the cluster formation scenario identified in our simulation, a violent galaxy merger causes a large amount of self-gravitating gas mass – comparable to that of a star cluster – to simultaneously reach high density above the star formation threshold.

3.4.3 Caveats

It is worth discussing a few points about this process:

(1) Since the newly formed star cluster grouped by Rockstar includes not only the just born star particles but also nearby pre-existing ones, the newly formed cluster's mean stellar density could be smaller than that of the gas progenitor of the cluster ’A’ at z = 5. This results in the discontinuity in densities between the blue and black solid lines at 0 Myr in the middle/bottom panels of Fig. 9. In fact, at z = 6.90, just ∼3.2 Myr after its formation, the cluster ’A’ is composed of ∼80 per cent in pre-existing stars, and only ∼20 per cent in stars just born in the past 5 Myr (see Section 3.5 for more discussion). The pre-existing stars are typically in the outskirts of the cluster, and tend to be preferentially removed already by z = 5, leaving only the core of tightly bound young stars (but again, see footnote 15). But notice that if we consider the mean density of the core within 15 pc containing almost exclusively newly born star particles, the gas progenitor's density smoothly transitions to the newly formed cluster's core density at 0 Myr (thick black dashed line in the bottom panel of Fig. 9).

(2) Readers should note that the high-density cluster formation process described above is not a sufficient condition, but only a necessary condition for a long-lasting bound cluster. Among the nine unbound associations shown in Figs 9--11 (thin red lines), three of them broadly followed similar paths as that of the cluster ’A’ when formed. But eventually, they do not become long-lasting bound clusters at z = 5. Other processes after their formation – such as tidal disruption or gravitational capture – interrupt their lives as bound clusters at some points before z = 5. Obviously, however, we cannot exclude the possibility that they could have survived as bound clusters if better numerical resolution had been adopted.

It is also useful to comment on the numerical aspects of our simulation:

(1) As pointed out before (footnote 12), what may not be so robust in our calculation is the quantities like the cluster's final size and density, because the cluster's relaxation effects are only marginally resolved. In contrast, what is indeed physically robust is the fact that a large self-gravitating gas mass of ≳ 10^5 − 6 M_⊙ could reach very high density simultaneously in galaxy mergers, and should very efficiently form stars before any of the various stellar feedback channels intervenes. Based on our arguments above, capturing this in simulations requires the ability to resolve surface densities at least as high as ∼10³ M_⊙ pc⁻². In nature, gas clouds may reach even higher densities before they become self-gravitating and star-forming.

(2) We argue that a numerical simulation with high resolution and a high dynamic range needs to be accompanied by appropriate subgrid physics to properly describe the clustered star formation scenario described here. For example, we find it harder to identify long-lasting bound clusters in the same runs with lower star formation threshold density (e.g. z5m10mr run with n_th = 50 cm⁻³ in Ma et al. 2015, see footnote 6) or without a self-gravity star formation criterion (e.g. z5m10e run in Ma et al. 2015), because in this cases stars artificially form at such low densities that clumps can never reach the critical mean densities described above (i.e. ≳ 10³ cm⁻³).¹⁸ This suggests that for the hydrodynamics code adopted in our study, a physically motivated star formation criterion is crucial to depict the Jeans gravitational collapse at the relevant mass and length resolution scales, and to numerically reproduce extraordinary conditions such as merger-induced star cluster formation.

3.5 Composition of a long-lasting bound star cluster at the end of the simulation (z = 5)

Finally we examine the cluster properties at z = 5 focusing on how the cluster formation event is imprinted in its composition. Fig. 13 shows the normalized PDF of the star cluster member particles’ creation times t_cr. A small dispersion in member particles’ ages for the cluster ’A’ (clustered at 0 Myr) reinforces the idea that the star formation episode lasted for a very short time. The two other small peaks located at −500 Myr and 100 Myr contribute negligibly to the total cluster mass. For comparison, we also plot the t_cr distribution for the cluster ’A’ at z = 6.90, just ∼3.2 Myr after its formation. Even though the cluster was initially grouped with nearby pre-existing stars of varying ages by Rockstar, they are typically in the outskirts of the cluster, and tend to be preferentially removed already by z = 5. This process left only the core of tightly bound young stars with very narrow distribution in particle ages. We also plot the PDF of one of the six unbound associations whose gas progenitors are shown in thin red lines in Figs 9–11, but do not resemble the evolution path of the cluster ’A’ progenitor. Indeed, its t_cr distribution reveals multiple generations of member star particles.

The width of the main spike in the cluster ’A’ members’ formation episode is of the order of the free-fall time (bottom panel of Fig. 13). Approximately 60 per cent of the entire member star particle population were born within an ∼3 Myr period in this main spike. That is, Δ t_cr ∼ 3 Myr ∼ t_ff, which is approximately consistent with the recent observations of young massive clusters (Bastian et al. 2013; Bastian, Hollyhead & Cabrera-Ziri 2014; Cabrera-Ziri et al. 2014).

In Fig. 14 we plot the normalized PDF of the star cluster member particles’ metallicities. This agrees reasonably well with internal abundance variations found in recent observations of metal-poor ’blue’ GCs (e.g. Gratton, Sneden & Carretta 2004; Da Costa et al. 2009; Marino et al. 2015; Roederer et al. 2016). It also shows that the star-forming gas was already enriched up to as high as 10⁻¹ Z_⊙, while approximately matching the observed mean abundance value. (One may also speculate that the cluster ’A’ exhibits a slightly narrower distribution than a representative unbound association does.) However, the large spread in the member particle metallicities for the unbound association (larger than the observed [Fe/H] spreads of ≲ 0.05 dex), or equivalently a large noise in the metallicity field, is to a great extent a numerical artefact because we do not include a sophisticated metal mixing scheme in the simulation – that is, presently metals are simply locked into individual particles once deposited by SNe. We plan to investigate if including a subgrid-scale turbulent metal diffusion scheme (e.g. Su et al. 2017) would reduce the artificially enhanced metallicity spread.

4 SUMMARY AND CONCLUSION

Using a state-of-the-art cosmological simulation of high-redshift merging proto-galaxies from the FIRE project, we have investigated the formation and evolution of star clusters, and in particular explored one formation hypothesis for present-day metal-poor ’blue’ GCs. In the simulation, two populations of clustered star formation emerge. Most stars form in unbound or loosely bound associations with mean stellar densities ∼10⁻² M_⊙ pc⁻³, or surface densities ∼1 M_⊙ pc⁻², corresponding to a typical GMC complex converting a few per cent of its mass into stars. This is the behaviour we have previously shown as ’typical’ in simulations with similar physics at low redshifts (Hopkins et al. 2012) and seen in observations (e.g. Evans et al. 2009). However, a few per cent of stars form in bound clusters with mean stellar densities ∼30 M_⊙ pc⁻³, or surface densities ∼10³ M_⊙ pc⁻² (Section 3.2). At these surface densities, Grudić et al. (2016) showed (in simulations of individual cluster-forming clouds, at sub-M_⊙ resolution but with the same physics) that stellar feedback begins to become unable to efficiently expel the cloud before/as it collapses into stars, leading to an order of unity fraction of gas turning into stars. Essentially, the self-gravity of the cloud becomes too large for small-scale feedback channels to overcome, particularly as the dynamical times for collapse become shorter than the massive star's evolution time-scale – and shorter than the time-scale for SNe explosion. A similar threshold was motivated by simple analytic arguments in Fall et al. (2010). Critically, the simulation here see both populations, with sufficient dynamic range to show that the same physics which produces low star formation rates in a galaxy-averaged sense (Hopkins et al. 2014; Orr et al. 2017) and most of the star formation in unbound associations can also produce a reasonable mass fraction in dense, bound, GC-like objects. It is also crucial that the simulation here includes stellar feedback not just from SNe but also from photoheating, radiation pressure, and stellar winds, since the cluster formation time-scales are shorter than ≲ 3 Myr, so SNe do not have the chance to explode before the clusters form (Section 2). In other words, simulations which do not include these other feedback mechanisms could therefore easily overestimate the formation efficiency of dense clusters.

We show that frequent mergers in high-redshift proto-galaxies can provide a fertile environment for the production of this GC-like population. The mergers disturb the gravitational potential and produce tidal shocks with Mach numbers ∼50–100, which rapidly and collectively push (pre-existing) large gas clouds of ≳ 10^5 − 6 M_⊙ to very high densities >10³ cm⁻³ (Section 3.4). Such clouds then cool and collapse to form stars extremely rapidly before stellar feedback can stop star formation. We typically see a dense ‘core’ in these systems, which is tightly bound, surrounded by a lower-density ‘envelope’ of loosely associated stars and nearby pre-existing stars. Those stars are quickly and preferentially stripped – not just by dynamical relaxation effects (which we do not necessarily resolve in the simulations) but also by the cluster passing through the tidal field of its host galaxy (Section 3.3). In the simulation we analysed, the tightly bound core is left behind and survives for ∼420 Myr, until the end of the simulation. This relic core has a very small age spread, and relatively small metallicity spread (Section 3.5).

Of course, the simulation here still suffers from various limitations. Our limited resolution, 800 h⁻¹M_⊙, means that we can only follow very massive clusters, and cannot follow detailed intra-cluster evolution (dynamical relaxation, mass segregation, etc; Section 3.3). Next-generation simulations in progress reach ∼30 M_⊙ resolution (Wheeler et al., in preparation), and in follow-up studies we can use these new runs to identify cluster-forming regions which are then refined to gain superior resolution. The limited evolution time of the present simulation – stopped at z = 5 – entails that we cannot yet predict the ultimate fate of the clusters (e.g. their locations in a z = 0 Milky Way-like galaxy, into which this progenitor halo should eventually evolve). Future simulations will also address this (Wetzel et al., in preparation), but also methods which replace clusters formed self-consistently with tracer particles which can be evolved for longer times (e.g. Renaud, Agertz & Gieles 2017; Li et al. 2017) could be useful. Improvements to the hydrodynamics and subgrid physics in our new FIRE-2 model (Hopkins et al. 2017) may also influence our predictions, although preliminary comparisons suggest that the hydrodynamic method (e.g. Gizmo’s mesh-free finite-mass mode versus the P-sph mode chosen here) does not significantly affect the conclusions presented here. Future work will include analysing different runs performed with different hydrodynamic solvers and star formation prescriptions to constrain their effects in the stellar distribution and cluster numbers. More accurate treatments of stellar yields (as opposed to the IMF-averaged yields adopted here) and subgrid metal diffusion in the dense ISM will be critical to make testable predictions for the internal abundance patterns within GCs, an obvious frontier for observations.

ACKNOWLEDGEMENTS

The authors thank Tom Abel, Nathan Bastian, Diederik Kruijssen, Joel Primack, and Eros Vanzella for useful discussions and insightful comments during the progress of this study. Ji-hoon Kim acknowledges support from NASA through an Einstein Postdoctoral Fellowship grant PF4-150147 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. He also acknowledges support from the Moore Center for Theoretical Cosmology and Physics at Caltech. He thanks for the support from the computational team at SLAC National Accelerator Laboratory, from Shawfeng Dong at the University of California Santa Cruz, and from Chris Mach at Caltech during the usage of the clusters for simulation analyses. The computing time used for the presented simulations was provided by Extreme Science and Engineering Discovery Environment (XSEDE) allocations TG-AST120025, TG-AST130039, TG-AST140023, and TG-AST140064. XSEDE is supported by the National Science Foundation (NSF) grant ACI-1053575. Philip Hopkins acknowledges support by the Gordon and Betty Moore Foundation through grant 776 to the Caltech Moore Center for Theoretical Cosmology and Physics, by the Alfred P. Sloan Foundation through Sloan Research Fellowship BR2014-022, by the NSF through grant AST-1411920 and CAREER grant 1455342, and by NASA through the ATP grant NNX14AH35G. The Flatiron Institute is supported by the Simons Foundation. Andrew Wetzel acknowledges support by a Caltech-Carnegie Fellowship, in part through the Moore Center for Theoretical Cosmology and Physics at Caltech, and by NASA through grant HST-GO-14734 from STScI. Claude-André Faucher-Giguère acknowledges support by the NSF through grants AST-1412836 and AST-1517491, by NASA through grant NNX15AB22G, and by STScI through grants HST-AR-14293.001-A and HST-GO-14268.022-A. Dušan Kereš acknowledges support by the NSF through grant AST-1412153, and by the Cottrell Scholar Award from the Research Corporation for Science Advancement. Shea Garrison-Kimmel acknowledges support from NASA through an Einstein Postdoctoral Fellowship grant PF5-160136. The publicly available yt code used in the analysis of this work is the product of collaborative efforts by many independent scientists from numerous institutions around the world. Their commitment to open science has helped make this work possible.

Footnotes

REFERENCES

Author notes