Globular Cluster Formation from Colliding Substructure

The origin of globular clusters (GCs) remains an unsolved problem in star and galaxy formation studies. With masses in the range 10⁴–10⁶ $\,{M}_{\odot }$, half-light radii of a few parsecs, a bimodal metallicity distribution, internal (star-to-star) variation in their light-element abundances, and typical ages >10 Gyr, these remarkably compact stellar systems are a common feature of galaxies in the local universe (see, e.g., the reviews by Harris 2001; West et al. 2004; Brodie & Strader 2006; Gratton et al. 2012; Kruijssen 2014; Forbes et al. 2018b). Young massive stellar clusters, with masses and densities comparable to those of GCs, are observed in galaxy mergers throughout the local universe (e.g., Whitmore et al. 1999), suggesting that the progenitors of metal-rich globular clusters could still be forming today in unusually high pressure environments. The mass of the entire GC system of a galaxy and its radial extent correlate tightly with the dark matter (DM) mass of the host halo (see, e.g., Harris et al. 2015, 2017; Forbes et al. 2016; Hudson & Robison 2018, and references therein), suggestive of a picture in which star formation in dense proto-GCs was relatively immune to the feedback mechanisms that hamper most star formation in the field.

The Milky Way (MW) contains about 160 known GCs. The blue (metal-poor) halo population is more extended, with a spatial distribution that falls off as ${r}^{-3.6\pm 0.2}$ at Galactocentric distances $\gtrsim 3.5\,\mathrm{kpc}$ (Bica et al. 2006). Only six known GCs inhabit the Galaxy at distances ≳90 kpc (Harris 1996; Laevens et al. 2014). The red (metal-rich) globular clusters are more spatially concentrated than the blue clusters: they lie within the solar circle and form a flattened, rotating population. Blue GCs outnumber red globular clusters by a ratio of 3:1. Despite their different location and kinematics, the red and blue clusters appear to have similar internal properties, such as masses, sizes, and ages. The age spread among the bulk of globular clusters is about 2 Gyr, but there are a few outliers (e.g., Leaman et al. 2013). The youngest GC, Whiting 1, is metal-rich and has an estimated age of 5.7 ± 0.3 Gyr (Valcheva et al. 2015), similar to those of Pal 1 and Terzan 7. The newly discovered metal-poor GC in Crater (Belokurov et al. 2014; Laevens et al. 2014), located at a distance of $145\pm 0.3\,\mathrm{kpc}$, is well described by a simple stellar population with an age of 7.5 ± 0.4 Gyr (Weisz et al. 2016). While most Galactic GCs formed before the peak of the cosmic star formation rate density at redshift 2 (Madau & Dickinson 2014), a handful were clearly assembled as late as z ≲ 1.

It has been an enduring challenge to explain such features, and despite a wide variety of proposed phenomenological models and a flurry of work on the present-day properties and evolution of GCs in galaxy halos, a self-consistent scenario of GC formation is yet to be constructed. It has long been suggested, for example, that a significant fraction of MW's GCs formed ex situ during the very early stages of galaxy formation, typically in dwarf-sized systems that were subsequently accreted into the Galactic potential well (Searle & Zinn 1978; Peebles 1984; Rosenblatt et al. 1988; Cen 2001; Bromm & Clarke 2002; Diemand et al. 2005; Moore et al. 2006; Boley et al. 2009; Griffen et al. 2010; Trenti et al. 2015; Kimm et al. 2016; Ricotti et al. 2016; Boylan-Kolchin 2017; Creasey et al. 2019, but see Chiou et al. 2019). By often relying on ad hoc prescriptions for the formation of globular clusters inside their own DM halos, these "pregalactic" models have been somewhat successful in reproducing some puzzling properties of the GC populations. The observed spatial coincidence between GCs and multiple tidal debris streams in the outer halo of M31 does indeed suggest a direct connection between some GCs and dwarf galaxy remnants (Mackey et al. 2010). In situ mechanisms have focused instead on the formation of GCs inside the main progenitor of their present-day host, perhaps from cooling-induced fragmentation of low-metallicity infalling gas (Fall & Rees 1985), during the merger with a gas-rich massive companion (Ashman & Zepf 1992; Muratov & Gnedin 2010; Choksi et al. 2018), or within the dense cores of giant molecular clouds in early galactic disks (Kravtsov & Gnedin 2005; Kruijssen 2015; Pfeffer et al. 2018).

Most clusters in the MW contain no measurable amounts of DM. The faint tidal tails observed around some GCs (Grillmair et al. 1995; Odenkirchen et al. 2003) provide strong constraints on their mass-to-light ratios and indicate that their mass distribution does not extend beyond the optical radius (Moore 1996). The absence of DM could be the result of stripping in the strong tidal field of the Galaxy (Mashchenko & Sills 2005; Saitoh et al. 2006), and does not preclude by itself the formation of GCs at the center of DM substructure. Ex situ scenarios in which GCs form at the center of their DM halos do predict, however, detectable amounts of DM in the outskirts of GCs that have always resided within the weak tidal field of the outer MW halo. This is not borne out by the data: recent observations of stellar kinematics in NGC 2419, located at 90 kpc from the center of our Galaxy, and in MGC1, located at 200 kpc from the center of M31, show that these GCs cannot be deeply embedded within dark halos having a virial mass greater than 10⁶ $\,{M}_{\odot }$ (Baumgardt et al. 2009; Conroy et al. 2011; Ibata et al. 2013).

In this paper we propose and investigate a possible scenario for the formation of GCs within the framework of ΛCDM hierarchical galaxy assembly. Many GC formation models leave largely unanswered the fundamental question of how, at early cosmological times, self-bound and extremely dense aggregates of ∼10⁵ stars, largely mono-metallic in iron-peak elements, were able to form and survive without retaining much of the DM of their former host. GCs are the result of a special mode of star formation that requires extremely high gas pressures, p ∼ 10⁷ k_B cm⁻³ K, some 3 dex higher than typical interstellar values. It is these pressures that hinder the dispersal of star-forming material, protect dense proto-GCs from the feedback processes that regulate star formation in the field, and produce high star formation efficiencies (e.g., Elmegreen & Efremov 1997). They are the kind of pressure that would result, e.g., if typical atomic clouds in the interstellar medium (ISM) with densities of a few atoms cm⁻³ were to collide at the Galactic orbital speeds of 200 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$, which is precisely the outcome of the mechanism proposed here: early high-speed collisions between subhalos infalling onto massive galaxy hosts, leading to the formation of dark-matter-less GCs.

Numerical simulations in ΛCDM have shown that massive galaxy halos are populated by a rich spectrum of substructures that collapsed at early times, the fossil remnants of a hierarchical merging process that is never complete (Klypin et al. 1999; Moore et al. 1999; Diemand et al. 2007a, 2008; Springel et al. 2008). About one thousand subhalos with preinfall masses m_peak ≳ 10⁸ $\,{M}_{\odot }$ are predicted to have been accreted by our Galaxy over cosmic time (see Han et al. 2018, and references therein).⁹ While cold gas condensation and normal, Population II star formation can take place at high redshift in these "atomic-cooling" subhalos because of hydrogen atomic radiative processes, their overall efficiency of converting baryons into stars must remain low in order to avoid overproducing the observed abundance of dwarf galaxy satellites of the MW (Madau et al. 2008; Boylan-Kolchin et al. 2014). Mechanisms such as early reionization of the intergalactic medium, supernova feedback, and H₂-regulation have all been invoked in order to reduce the star formation efficiencies of low-mass field halos and even prevent the smallest ones from forming stars altogether, and "dark dwarfs" with long gas depletion timescales are often produced in cosmological simulations (e.g., Okamoto et al. 2008; Kuhlen et al. 2013; Shen et al. 2014; Benítez-Llambay et al. 2015; Sawala et al. 2015; Ricotti et al. 2016). Reionization heating is expected to photoevaporate low-density gas from these shallow potential wells, but to have little effect on relatively high-density gas in their cores that can self-shield from UV background radiation. The detailed mapping betwen the mass spectrum of DM substructure and the population of Galactic satellites remains, however, a matter of debate, and so is the smallest mass halo that is capable of hosting an observable baryonic counterpart. A recent probabilistic comparison (abundance matching) between the luminosity function of MW's dwarfs and N-body zoom-in simulations sets a 1σ upper bound to the preinfall peak subhalo mass of Segue I—the faintest MW satellite—of m_peak < 2.4 × 10⁸ $\,{M}_{\odot }$ (Jethwa et al. 2018). Similarly, a fresh reassessment of the "missing satellites problem," cast in terms of number counts, suggests that luminous satellites of the MW inhabit subhalos with infall masses as small as 10⁷–10⁸ $\,{M}_{\odot }$ (Kim et al. 2018).

The central idea of this work is to suggest that high-speed collisions between infalling substructures during the assembly of the main galaxy host may result in the formation of "naked" GCs. Like in the Bullet galaxy gluster (Clowe et al. 2006), colliding DM subhalos and their stars will simply pass through one another, largely continuing on their original orbital path (close encounters are more effective at disrupting the colliding satellites at low speeds). The atomic gas within them, however, will collide at highly supersonic speed and decouple from the collisionless component. The collision will compress and shock-heat the gas to characteristic temperatures ≳10⁵ K, creating the very high pressure region that is conducive to GC formation. We will show that the low-metallicity gas will cool via atomic line emission on a timescale that is much shorter than the crossing time. Under the appropriate conditions, shock dissipation of the relative kinetic energy will lead to the coalescence of the colliding gas clouds and the formation of a Jeans unstable slab, rather than cloud disruption. A high star formation efficiency may result from the short dynamical timescale and high binding energy of the splash remnant.

2.1. Kinetic Theory Approach

An estimate of the number of close collisions can be made within the framework of kinetic theory (e.g., Makino & Hut 1997). Let us start by selecting subhalos by their preinfall mass, i.e., by using the "unevolved" substructure mass function—the distribution of peak bound masses m_peak of subhalos accreted by the host at all previous redshifts (e.g., van den Bosch et al. 2005; Jiang & van den Bosch 2016). In ΛCDM cosmological simulations, this universal redshift-independent function is well fitted by a double Schechter function of the form (Han et al. 2018)

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{\mathrm{sub}}}{d\xi }={\xi }^{-1}\left({a}_{1}{\xi }^{-{\alpha }_{1}}+{a}_{2}{\xi }^{-{\alpha }_{2}}\right)\exp (-b{\xi }^{\beta }),\end{eqnarray} \tag{ 1 }$$

where ξ = m_peak/M_vir is the ratio of the peak mass of a subhalo to the host halo mass, and (dN_sub/dξ)dξ is the number of subhalos with masses between ξ and ξ + dξ. When adopting for the host mass the virial definition corresponding to the overdensity of the spherical collapse model, the best-fit parameters are a₁ = 0.11, a₂ = 0.32, α₁ = 0.95, α₂ = 0.08, b = 8.9, and β = 1.9 (Han et al. 2018). In the case of small subhalos, the unevolved spatial distribution (the spatial distribution at fixed infall mass, see Han et al. 2016) follows the spherical Navarro et al. (1997); (NFW) density profile of the host,

$$\begin{eqnarray}&&{n}_{\mathrm{sub}}(\xi ,x)\propto \displaystyle \frac{1}{{cx}{\left(1+{cx}\right)}^{2}},\end{eqnarray} \tag{ 2 }$$

where x ≡ r/R_vir, R_vir is the virial radius of the host, and c is its "concentration parameter."¹⁰ The subhalo unevolved mass function and spatial distributions are approximately separable (i.e., the spatial distributions of different peak mass subhalos are similar except for a change in amplitude), and the number density of small subhalos with peak masses >ξ M_vir, n_sub(>ξ, x), can be written as

$$\begin{eqnarray}&&4\pi {R}_{\mathrm{vir}}^{3}{n}_{\mathrm{sub}}(\gt \xi ,x)=\displaystyle \frac{{c}^{3}}{{gcx}{\left(1+{cx}\right)}^{2}}{\int }_{\xi }^{1}\displaystyle \frac{{{dN}}_{\mathrm{sub}}}{d\xi ^{\prime} }d\xi ^{\prime} ,\end{eqnarray} \tag{ 3 }$$

where $g\equiv 1/[\mathrm{ln}(1+c)-c/(1+c)]$.

Let us further make the simplifying assumption that the velocity distribution of subhalos is a local Maxwellian with one-dimensional velocity dispersion σ_v. The distribution of relative encounter speeds is then

$$\begin{eqnarray}&&f({v}_{\mathrm{rel}})=\displaystyle \frac{1}{2\sqrt{\pi }{\sigma }_{v}^{3}}{v}_{\mathrm{rel}}^{2}\exp \left(-\displaystyle \frac{{v}_{\mathrm{rel}}^{2}}{4{\sigma }_{v}^{2}}\right),\end{eqnarray} \tag{ 4 }$$

and the mean encounter speed is ${\bar{v}}_{\mathrm{rel}}=4{\sigma }_{v}/\sqrt{\pi }$. For isotropic orbits in an NFW potential, the velocity dispersion obtained by solving the Jeans equation is (Łokas & Mamon 2001)

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{v}^{2}(x) & = & \left(\displaystyle \frac{{{GM}}_{\mathrm{vir}}}{{R}_{\mathrm{vir}}}\right)g{\left(1+{cx}\right)}^{2}x\\ & & \times \,{\int }_{x}^{\infty }\left[\displaystyle \frac{\mathrm{ln}(1+{cs})}{{s}^{3}{\left(1+{cs}\right)}^{2}}-\displaystyle \frac{c}{{s}^{2}{\left(1+{cs}\right)}^{3}}\right]{ds}.\end{array}\end{eqnarray} \tag{ 5 }$$

Under the assumption that collisions occur on random orbits, we can finally estimate the mean number of encounters between subhalo pairs of peak masses (in units of M_vir) $\gt {\xi }_{1}$ and >ξ₂ in a time interval Δt as

$$\begin{eqnarray}\begin{array}{rcl}{N}_{\mathrm{coll}} & = & {\rm{\Delta }}t\,4\pi {R}_{\mathrm{vir}}^{3}{\int }_{0}^{1}{\bar{v}}_{\mathrm{rel}}(x){x}^{2}{dx}\\ & & \times {\int }_{{\xi }_{1}}^{1}{\int }_{{\xi }_{2}}^{1}{n}_{\mathrm{sub}}(\xi ,x){n}_{\mathrm{sub}}(\xi ^{\prime} ,x)(\pi {b}_{\max }^{2})d\xi d\xi ^{\prime} ,\end{array}\end{eqnarray} \tag{ 6 }$$

where b_max is the maximum impact parameter for a close collision.

Consider now, for illustrative purposes, a MW-sized host halo of mass M_vir = 10¹² $\,{M}_{\odot }$ and concentration c = 10 (e.g., Duffy et al. 2008). Its isotropic velocity dispersion profile in Equation (5) reaches a maximum value of 102 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$ at x = 0.08, corresponding to a mean collision speed of 230 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$. These relative orbital velocities are much greater than the internal velocities of subhalos, and few such encounters will actually result in mergers. We can therefore neglect gravitational focusing and, in order to ensure that the region of overlap between the gaseous cores of the interacting clumps is extensive, define penetrating collisions as those with b_max = 0.1(r_vir,1 + r_vir,2), where r_vir,1 and r_vir,2 are the preinfall virial radii of the two clumps. The mean number of impacts within R_vir between subhalos having peak masses m_peak > 10⁸ $\,{M}_{\odot }$ is then

$$\begin{eqnarray}&&{N}_{\mathrm{coll}}\simeq 59\left(\displaystyle \frac{{\rm{\Delta }}t}{\mathrm{Gyr}}\right).\end{eqnarray} \tag{ 7 }$$

This is clearly an interesting figure—comparable after ${\rm{\Delta }}t\sim$ a few Gyr to the number of GCs observed today in the halo of the MW—if many such encounters were to result in the formation of one or more globular clusters. Because of the steep substructure mass function, most collisions involve subhalos at the small-mass end of the distribution: the mean frequency of close encounters between atomic-cooling subhalos with m_peak > 10⁸ $\,{M}_{\odot }$ and satellites that are ten times more massive is 3.6 times smaller than Equation (7). There are only N_coll ∼ 3(Δt/Gyr) close collisions between massive satellites with m_peak > 10⁹ $\,{M}_{\odot }$, but these may give origin to multiple GCs per encounter. The same calculation predicts about 1 collision Gyr⁻¹ between atomic-cooling subhalos in a dwarf galaxy-sized host with M_vir = 10¹⁰ $\,{M}_{\odot }$.

In reality, of course, many of the simplified assumptions behind our integration of Equation (6) do not hold in practice: (1) subhalos are often accreted at similar times and locations as members of groups along filaments, and this causes an enhancement of encounters at small angles (Benson 2005; Zentner et al. 2005; González & Padilla 2016); (2) after accretion, subhalos are susceptible to dynamical friction and tidal stripping, and their mass and spatial distributions evolve away from their form at infall (Diemand et al. 2007b; Springel et al. 2008); (3) the abundance of subhalos above a given peak mass evolves with time, and so does the background gravitational potential in which they move; (4) subhalos are predicted to have a non-Maxwellian orbital velocity function, with centrally rising velocity anisotropy (Sawala et al. 2017); and (5) ram pressure stripping will remove the outer gaseous component of subhalos on orbits with small pericenter radii, a process that is most effective close to pericenter passage (Mayer et al. 2006; Grcevich & Putman 2009). In Section 4 we shall tackle some of these issues by using the densely time-sampled snapshots of the Via Lactea N-body simulation to track collisions between infalling substructure in a massive MW-sized host that grows and evolves in a fully cosmological context.

2.2. Scaling with Host Mass and Ex situ Globulars

One of the most intriguing aspects of GC phenomenology is their relationship to DM halos: the total mass of the GC population of a galaxy and the number of GCs correlate linearly with the DM mass of the host halo. This trend is valid over five orders of magnitude in galaxy mass (see, e.g., Harris et al. 2015; Forbes et al. 2016; Burkert & Forbes 2019, and references therein) and contrasts markedly with the nonlinear relation between total stellar mass and DM-halo mass (e.g., Behroozi et al. 2013). The observed correlation may point to a fundamental GC system–DM connection rooted in the cluster formation physics, or simply be the inevitable consequence of hierarchical assembly and the central limit theorem (Boylan-Kolchin 2017; Burkert & Forbes 2019; El-Badry et al. 2019).

It is of obvious interest at this stage to use kinetic theory and derive a scaling relation between the number of collisions and the virial mass of the host following Equation (6). The collision frequency between subhalos more massive than a given m_peak ≪ M_vir at fixed concentration parameter is proportional to ${n}_{\mathrm{sub}}^{2}\propto {R}_{\mathrm{vir}}^{-6}\,{M}_{\mathrm{vir}}^{1.9}$ (from Equations (1) and (3)) times ${\bar{v}}_{\mathrm{rel}}\,\propto {\left({M}_{\mathrm{vir}}/{R}_{\mathrm{vir}}\right)}^{1/2}$ (see Equation (5)). Including a dependence on halo concentration and assuming that the collisions of interest continue for a timescale Δt that is independent of host mass, one infers from Equation (6):

$$\begin{eqnarray}&&{N}_{\mathrm{coll}}\propto {R}_{\mathrm{vir}}^{-3}\,{M}_{\mathrm{vir}}^{1.9}{\left({M}_{\mathrm{vir}}/{R}_{\mathrm{vir}}\right)}^{1/2}{c}^{4}.\end{eqnarray} \tag{ 8 }$$

Noting that ${R}_{\mathrm{vir}}\propto {M}_{\mathrm{vir}}^{1/3}$ and using the halo concentration–mass relation of Duffy et al. (2008), $c\propto {M}_{\mathrm{vir}}^{-0.09}$, we obtain the following scaling ${N}_{\mathrm{coll}}\propto {M}_{\mathrm{vir}}^{0.87}$. Under the assumption that GCs populate present-day DM halos in direct proportion to ${N}_{\mathrm{coll}}$, let us write

$$\begin{eqnarray}&&{N}_{\mathrm{in}-{\rm{s}}}^{\mathrm{GC}}=A{\left({M}_{\mathrm{vir}}/{10}^{12}{M}_{\odot }\right)}^{0.87},\end{eqnarray} \tag{ 9 }$$

where A is an arbitrary normalization factor. In this scenario, ${N}_{\mathrm{in}-{\rm{s}}}^{\mathrm{GC}}$ tracks the population of GCs formed in situ within a massive galaxy host. The same substructure collision mechanism, at work in satellites prior to infall, would also produce a significant population of globular cluster that formed ex situ and were subsequently accreted. The ${N}_{\mathrm{in}-{\rm{s}}}^{\mathrm{GC}}-{M}_{\mathrm{vir}}$ relation given above and the substructure mass function in Equation (1) allow a straightforward calculation of the mean number of accreted GCs as

$$\begin{eqnarray}&&{N}_{\mathrm{ex}-{\rm{s}}}^{\mathrm{GC}}={N}_{\mathrm{in}-{\rm{s}}}^{\mathrm{GC}}\,{\int }_{{\xi }_{\min }}^{1}{\xi }^{0.87}\displaystyle \frac{{{dN}}_{\mathrm{sub}}}{d\xi }d\xi ,\end{eqnarray} \tag{ 10 }$$

where ξ_min ≡ m_min/M_vir is the critical subhalo mass below which GCs cannot form. Obviously, collisions between m_peak > 10⁸ $\,{M}_{\odot }$ clumps as a mechanism for cluster formation require at least a few (say ∼4) atomic-cooling subhalos in a given host; this trivial condition yields, after integrating Equation (1), m_min ≃ 10^9.6 $\,{M}_{\odot }$. This critical mass could be larger—thereby decreasing the fraction of accreted globular clusters—if the relative orbital velocity v_rel required to trigger the formation of a GC by impact had to exceed ∼50 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$. Note, however, that lower velocity encounters, as those expected in a less massive host, may still produce the high pressures that are conducive to GC formation if the colliding gaseous cores were typically denser, e.g., at very high redshifts.

Assuming here that all halos below m_min = 10^9.6 $\,{M}_{\odot }$ do not contain a GC, Equation (10) gives ${N}_{\mathrm{ex}-{\rm{s}}}^{\mathrm{GC}}/{N}_{\mathrm{in}-{\rm{s}}}^{\mathrm{GC}}=0.35,0.73,1.17$ for M_vir = 10¹¹, 10¹², 10¹³ $\,{M}_{\odot }$, i.e., most GC formation occurs in situ for the lowest-mass host halos and most GCs in massive elliptical systems are actually accreted.¹¹ With the addition of the accreted globular clusters, the scaling of the total number of GCs, ${N}_{\mathrm{GC}}={N}_{\mathrm{ex}-{\rm{s}}}^{\mathrm{GC}}+{N}_{\mathrm{in}-{\rm{s}}}^{\mathrm{GC}}$, with host virial mass becomes very close to linear,

$$\begin{eqnarray}&&{N}_{\mathrm{GC}}\propto {M}_{\mathrm{vir}}^{0.98}.\end{eqnarray} \tag{ 11 }$$

Therefore, in a collision-driven scenario, a constant GC number-to-halo mass ratio is the result of encounter probability calculations (Equation (8)) and hierarchical clustering (Equation (10)): GCs are the result of a distinctive mode of star formation that tracks the total mass of their host galaxy rather than its stellar mass.

Figure 1 shows the observed GC number–halo mass relation from a recent compilation by Burkert & Forbes (2019); (see also Harris et al. 2017; Forbes et al. 2018a), compared with our theoretical prediction for an assumed normalization A = 95 in Equation (9). The model appears to reproduce the observed trend reasonably well, with the expected average number of GC systems dropping below a few in dwarf galaxies with M_vir ≲ 10¹⁰ $\,{M}_{\odot }$. With the adopted parameters, an MW-mass system with M_vir = 10¹² $\,{M}_{\odot }$ today would host a grand total of 165 massive GCs, of which 70 formed ex situ. Estimates of the number of accreted GCs within the MW today range from 30 to 90 (Forbes & Bridges 2010; Leaman et al. 2013), so it is conceivable that a substantial fraction of globular cluster in massive hosts may have an external origin. The significant scatter in the observed relation on dwarf galaxy scales may reflect variations in merger histories, uncertainties in dark halo mass determinations based on kinematical tracers, or environmental effects such as tidal stripping (Burkert & Forbes 2019).

Zoom In Zoom Out Reset image size

While it is reassuring that a scenario in which GCs are the result of collisions between subhalos may be able to accomodate a uniform GC production rate per unit host halo mass as implied by the data, it may be premature to read too much into this comparison given the above-mentioned limitations of kinetic theory. In particular, our modeling so far has not provided any information on the age and age spread of GCs, on the epoch at which the GC-to-halo mass relation may actually be established, on the normalization of such relation, and on the baryonic physics leading to the formation of massive globular clusters. In the following sections, we shall discuss the conditions under which high-speed impacts may lead to cloud coalescence, gravitational instability, and the formation of "naked," bound GCs, and use numerical simulations to argue that a collision-based framework may fulfill several key observational constraints on GCs that have emerged over the last two decades.

Observations of interacting galaxies like the Antennae and Mice pairs show widespread stellar cluster formation generated by the collision, and high-resolution simulations of encounters between disk galaxies with realistic ISM reveal the presence of shock-induced gas compression and star formation at the collision interface (e.g., Saitoh et al. 2009). There is a vast literature on supersonic cloud–cloud collisions as a possible triggering mechanism for star formation in the ISM (e.g., Stone 1970; Smith 1980; Gilden 1984; Nagasawa & Miyama 1987; Habe & Ohta 1992; Anathpindika 2009; Arreaga-García et al. 2014; Balfour et al. 2017), and cloud–cloud collisions in the young Galaxy have already been invoked in the context of GC formation (e.g., Murray & Lin 1989; Kang et al. 1990; Lin & Murray 1991; Kumai et al. 1993; Hartwick 2009). Below, we show that, when a similar collision occurs between subhalos, the postshock gas remnant will cool down on a timescale that is much shorter than the compression time, and may, under the right conditions, become gravitationally unstable.

We shall focus on collisions between subhalos where both members of the pair are above the atomic-cooling threshold at infall, i.e., have a mass corresponding to a virial temperature of T_vir > 8000 K,

$$\begin{eqnarray}&&{m}_{\mathrm{peak}}\gt 1.6\times {10}^{8}\,\,{M}_{\odot }{[\mu (1+z)/5]}^{-3/2}\,{{ \mathcal B }}^{-1/2}.\end{eqnarray} \tag{ 12 }$$

Here, μ is the mean molecular weight per particle (μ = 1.23 for neutral primordial gas), ${ \mathcal B }\,\equiv \,{{\rm{\Delta }}}_{c}/({{\rm{\Omega }}}_{m}^{z}\,18{\pi }^{2})$, Δ_c = 18π² + 82y − 39y² is the redshift-dependent density contrast at virialization (Bryan & Norman 1998), $y\equiv {{\rm{\Omega }}}_{m}^{z}-1$, ${{\rm{\Omega }}}_{m}^{z}\,={{\rm{\Omega }}}_{m}{\left(1+z\right)}^{3}/[{{\rm{\Omega }}}_{m}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}]$, and z is the infall redshift. In evaluating these expressions, we have assumed a Planck Collaboration et al. (2018) flat cosmology with parameters $({{\rm{\Omega }}}_{m},{{\rm{\Omega }}}_{{\rm{\Lambda }}},h)=(0.315,0.685,0.674)$, and ${ \mathcal B }=1$. Atomic-cooling subhalos are likely to be polluted by metals through inefficient star formation, which is key since GCs will inherit the gas-phase metallicity of the colliding subhalo population. According to cosmological radiation hydrodynamics simulations by Wise et al. (2014) that follow the buildup of dwarf galaxies from their early Population III progenitors (see also Ricotti et al. 2016), most atomic-cooling halos host metal-enriched stars by redshift 6, when infalling substructures begin colliding at high speed. The gas shock heated in a close encounter will then be pre-enriched by previous generations of massive star formation.

Gas distributed isothermally in atomic-cooling halos develops warm central cores of constant baryon density ${n}_{c}\,\approx 0.6\,{[(1+z)/5]}^{3}\,\,{\mathrm{cm}}^{-3}\equiv {n}_{1}\,\,{\mathrm{cm}}^{-3}$ within radii r_c ∼ 0.1 r_vir (roughly independent of halo mass); (Ricotti 2009; Prieto et al. 2013; Visbal et al. 2014a), and follows an r⁻² profile in the outer regions. Here,

$$\begin{eqnarray}&&{r}_{\mathrm{vir}}\simeq 4.4\,\mathrm{kpc}{\left(\displaystyle \frac{{m}_{\mathrm{peak}}}{{10}^{8.5}{M}_{\odot }}\right)}^{1/3}\left(\displaystyle \frac{5}{1+z}\right)\end{eqnarray} \tag{ 13 }$$

is the virial radius of the system. Prior to infall, the susceptibility of such systems to internal and external feedback mechanims may decrease their gas content to 5%–10% of their virial mass (Wise et al. 2014). After infall, their centrally concentrated gaseous cores are expected to survive disruption by ram pressure before the first pericenter passage (Mayer et al. 2006; Visbal & Haiman 2014b).

3.1. Shock Heating, Cooling, and Gravitational Instability

To highlight the dominant physical processes at work during a high-speed close interaction, let us consider the idealized situation of a head-on collision between two identical gaseous cores of characteristic hydrogen number density n_H,1 = 0.75n_c, mass density ρ₁ = n_cm_p, temperature T₁ = 8000 K, and radius r_c = 0.1 r_vir. These are the density and thermal pressure conditions of the warm diffuse component of the ISM. In the absence of molecular cooling, even when gas is heated by shocks and compression above the 8000 K limit, atomic radiative losses together with photoelectric heating from dust grains, cosmic ray heating, and soft X-ray heating will keep the gas temperature near this value (Wolfire et al. 1995). The total gas mass involved in the impact,

$$\begin{eqnarray}\begin{array}{rcl}2{m}_{c} & = & \displaystyle \frac{8\pi }{3}{\rho }_{1}{r}_{c}^{3}\simeq {10}^{7}\,{M}_{\odot }\,{n}_{1}{\left(\displaystyle \frac{{r}_{c}}{0.1{r}_{\mathrm{vir}}}\right)}^{3}\\ & & \times {m}_{8.5}{\left(\displaystyle \frac{5}{1+z}\right)}^{3},\end{array}\end{eqnarray} \tag{ 14 }$$

where m_8.5 ≡ m_peak/10^8.5 $\,{M}_{\odot }$, is about 10% of the mass corresponding to the universal baryon fraction. Following a highly supersonic interpenetrating impact, most of the gas decouples from the collisionless component. Two shock waves arise and propagate from the contact discontinuity with velocity (in the adiabatic phase) V_s ≃ 2v_rel/3 relative to the unpertubed medium, heating the gas to

$$\begin{eqnarray}&&{T}_{2}=\displaystyle \frac{3\mu {m}_{p}}{16{k}_{B}}{V}_{s}^{2}\simeq 1.3\times {10}^{5}\,{\rm{K}}\,{V}_{100}^{2},\end{eqnarray} \tag{ 15 }$$

where ${V}_{100}\equiv {V}_{s}/100\,\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and we have set the molecular weight to μ = 0.59. The density is enhanced by a factor of 4, the postshock pressure is p₂ = (3/4)ρ₁V_s², and the gas cools via hydrogen, helium, and metal-line emission on a (isobaric) timescale

$$\begin{eqnarray}&&{t}_{\mathrm{cool}}=\displaystyle \frac{3{k}_{B}{T}_{2}{n}_{\mathrm{tot}}}{2{n}_{{\rm{H}},2}^{2}{\rm{\Lambda }}}\simeq 0.03\,\mathrm{Myr}\,\left(\displaystyle \frac{{V}_{100}^{3.6}}{{n}_{1}}\right).\end{eqnarray} \tag{ 16 }$$

Here, n_tot is the total particle number density (n_tot/n_H,2 = 9/4 for completely ionized gas of primordial composition), n_H,2 = 4 n_H,1 is the postshock mean hydrogen density, Λ is the cooling function (e.g., Wang et al. 2014), and we have assumed collisional ionization equilibrium and ignored gas clumping (see Figure 2). Since t_cool is typically much shorter than the crossing time,

$$\begin{eqnarray}&&{t}_{\mathrm{cross}}\,\equiv \,\displaystyle \frac{2{r}_{c}}{{V}_{s}}\simeq 7.8\,\mathrm{Myr}\,({V}_{100}^{-1}\,{r}_{0.4}),\end{eqnarray} \tag{ 17 }$$

where r_0.4 ≡ r_c/0.4 kpc, the internal energy increase caused by the shock is radiated away and the impact is effectively isothermal (T₂ = T₁ = T), characterized by a sound speed ${c}_{s}=\sqrt{{k}_{B}T/\mu {m}_{p}}$ and a Mach number

$$\begin{eqnarray}&&{ \mathcal M }\,\equiv \,\displaystyle \frac{{V}_{s}}{{c}_{s}}\simeq 13.6\,({T}_{3.9}^{-1/2}\,{V}_{100}),\end{eqnarray} \tag{ 18 }$$

where T_3.9 ≡ T/8000 K and we have assumed μ = 1.23. The postshock hydrogen density is

$$\begin{eqnarray}&&{n}_{{\rm{H}},2}={{ \mathcal M }}^{2}\,{n}_{{\rm{H}},1}\simeq 140\,{\mathrm{cm}}^{-3}\,({n}_{1}\,{T}_{3.9}^{-1}\,{V}_{100}^{2}),\end{eqnarray} \tag{ 19 }$$

and the shock velocity relative to the unpertubed medium is now V_s ≃ v_rel/2. The collision is also approximately one-dimensional, since the lateral rarefaction timescale is $\sim { \mathcal M }\,{t}_{\mathrm{cross}}$. For a strictly isothermal collision, when the compression phase is completed after a time t_cross, a rarefaction wave will propagate through the shock-compressed gas slab at the local sound speed, and the medium will rapidly expand into the surroundings on a timescale $\sim {t}_{\mathrm{cross}}/{ \mathcal M }$. If the collision is to result in gravitational instability and collapse, there must be density perturbations that can grow on a timescale  ≲ t_cross. We shall see below, however, that gas further downstream will actually cool below its preshock temperature via fine-structure metal lines, i.e., the collision is "more dissipative than isothermal" (e.g., Whitworth & Clarke 1997).

Zoom In Zoom Out Reset image size

Lyα cooling becomes inefficient below 8000 K, and the fast shock will dissociate most preexisting molecules. In the absence of an ambient UV radiation field, enough molecular hydrogen may reform behind the shock via the H⁻ channel or on the surface of dust grains to rapidly cool the compressed gas to 100 K. It seems likely, however, that photodissociation by local stellar sources of Lyman–Werner photons—either within the colliding subhalos or in the main host—will act to suppress H₂ formation in the postshock gas, and the temperature will then plateau around 8000 K (e.g., Kang et al. 1990). The Bonnor–Ebert mass in such a warm, dense, and pressurized medium is

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{BE}} & = & 1.182\,\displaystyle \frac{{c}_{s}^{3}}{\sqrt{{G}^{3}\,{{ \mathcal M }}^{2}\,{\rho }_{1}}}\\ & \simeq & 7.7\times {10}^{5}\,\,{M}_{\odot }({n}_{1}^{-1/2}\,{T}_{3.9}^{2}\,{V}_{100}^{-1}),\end{array}\end{eqnarray} \tag{ 20 }$$

a factor of $1/{ \mathcal M }$ smaller than its value in the preshock gas (where M_BE > m_c, i.e., the precollision clumps are gravitionally stable), and comparable to the observed mass of GCs.

GCs are generally metal-poor. While the lowest metallicity globular cluster listed in the 2010 edition of the Harris (1996) catalog has [Fe/H] = −2.4, the median of the Galactic GC population lies at [Fe/H] = −1.3. Even at the low metallicities of blue globular clusters, however, the postshock dense medium will eventually cool below 8000 K via the collisional excitation of the fine-structure lines of C ii (158 μm) and O i (63 μm). The rate coefficient for the collisional excitation of O i by H atoms is ${k}_{\mathrm{OI}}\simeq 1.3\times {10}^{-9}\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}\,{T}_{3.9}^{0.41}{e}^{-{T}_{\mathrm{OI}}/T}$ where T_OI = 228 K is the energy of the excited O i level over k_B (Draine 2011). At densities far below the critical density of the line, this process removes energy at the rate

$$\begin{eqnarray}&&{n}_{{\rm{H}},2}^{2}{{\rm{\Lambda }}}_{{\rm{O}}}={n}_{{\rm{H}},2}^{2}\,{{ \mathcal A }}_{{\rm{O}}}\,{k}_{\mathrm{OI}}({k}_{B}{T}_{\mathrm{OI}}),\end{eqnarray} \tag{ 21 }$$

where ${{ \mathcal A }}_{{\rm{O}}}$ is the abundance of O relative to hydrogen, ${{ \mathcal A }}_{{\rm{O}}}=4.3\times {10}^{-5}\,{Z}_{-1.3}$ (Anders & Grevesse 1989), Z_−1.3 is the gas metallicity in units of 10^−1.3 solar, and all oxygen is assumed to be neutral and in the gas phase. The isobaric cooling timescale,

$$\begin{eqnarray}&&{t}_{\mathrm{cool}}^{{\rm{O}}}=\displaystyle \frac{5{k}_{B}{{Tn}}_{\mathrm{tot}}}{2{n}_{{\rm{H}},2}^{2}{{\rm{\Lambda }}}_{{\rm{O}}}}\simeq 0.4\,\mathrm{Myr}\,\left(\displaystyle \frac{{T}_{3.9}^{1.59}}{{V}_{100}^{2}{n}_{1}{Z}_{-1.3}}\right),\end{eqnarray} \tag{ 22 }$$

is still shorter than the crossing time, but it is considerably longer than t_cool in Equation (16) as a result of the mismatch in timescales between fine-structure metal-line cooling below 8000K and atomic cooling at higher temperatures. After a time $\sim {t}_{\mathrm{cool}}^{{\rm{O}}}$, the flow will cool below its preshock temperature and contract further, generating a growing, cold, and very dense layer at the center of the shock-bounded slab.

The arguments above are based on simple timescale estimates, and to check the robustness of our results we have used the nonequilibrium chemistry package krome (Grassi et al. 2014) to construct more sophisticated models of the cooling and thermal properties of the postshock gas. The thermal model includes nonequilibrium cooling for H and He at all temperatures and for metals below 10⁴ K (Bovino et al. 2016), and compressional heating. Strong shock jump conditions and isobaric cooling were adopted for gas at different metallicities (ranging from 0 to 0.1 solar), and the cooling function was computed assuming dust-free and optically thin conditions, and neglecting molecular-phase processes. Figure 3 shows the time evolution of the postshock temperature for a case with V_s = 100 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and n_H,1 = 0.75 $\,{\mathrm{cm}}^{-3}$. The figure is consistent with the qualitative estimates given above: low-density gas, shock-heated to high temperatures, cools fast via atomic hydrogen and helium line emission to about 6500 K, and then levels off for a timescale that is longer at lower metallicities. Eventually, metal cooling via fine-structure lines of C and O takes over. The end result is a runaway cooling phenomenon that drives the shocked medium to temperatures T₃ < T and hydrogen densities n_H,3 = n_H,2T/T₃ (assuming isobaric cooling at constant mean molecular weight) that are, respectively, well below and above those of the colliding gas clouds. Note that these calculations do not explicitly include internal (within the colliding substructures) and external (within the host galaxy) sources of UV and X-ray photoheating.

Zoom In Zoom Out Reset image size

The one-dimensional compression will create a cold layer of thickness-to-diameter ratio $\sim 2{n}_{{\rm{H}},1}/{n}_{{\rm{H}},3}\lt 2/{{ \mathcal M }}^{2}\,\sim 0.01\,({T}_{3.9}\,{V}_{100}^{-2})$ that depends on the entity of radiative losses. In this case, gravitational accelerations are most important in the dynamics of flows transverse to the collision axis, and it is more appropriate to discuss the criterion for the gravitational instability of a thin, ram pressure–bounded slab of shocked gas (e.g., Stone 1970; Elmegreen & Elmegreen 1978; Vishniac 1983; Gilden 1984; Larson 1985). For an infinitely thin sheet of constant surface density Σ and isothermal sound speed c_s, the dispersion relation for modes transverse to the collision axis gives the wavelength and timescale of the fastest growing mode as

$$\begin{eqnarray}&&{\lambda }_{G}=\displaystyle \frac{2{c}_{s}^{2}}{G{\rm{\Sigma }}},\end{eqnarray} \tag{ 23 }$$

and

$$\begin{eqnarray}&&{t}_{G}=\displaystyle \frac{{c}_{s}}{\pi G{\rm{\Sigma }}}.\end{eqnarray} \tag{ 24 }$$

Perturbations smaller than ${\lambda }_{G}$ are gravitationally stable, while those with longer wavelengths are unstable but grow more slowly. The surface density of the cold layer between the two shocks increases with time as Σ(t) = 2ρ₁V_st, where t is measured from the onset of metal cooling below 8000 K. The shocked gas is first liable to gravitational fragmentation at time t = t_G, when Σ_G ≡ Σ(t_G) = 2ρ₁V_st_G (Kang et al. 1990). Using Equation (24) one derives then

$$\begin{eqnarray}&&{t}_{G}={\left(\displaystyle \frac{{c}_{s}}{2\pi G{\rho }_{1}{V}_{s}}\right)}^{1/2}\simeq 3.8\,\mathrm{Myr}\,{\left(\displaystyle \frac{{ \mathcal C }}{{n}_{1}{V}_{100}}\right)}^{1/2},\end{eqnarray} \tag{ 25 }$$

and

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{G}={\left(\displaystyle \frac{2{\rho }_{1}{V}_{s}{c}_{s}}{\pi G}\right)}^{1/2}\simeq 19\,{M}_{\odot }{\mathrm{pc}}^{-2}\,{\left({n}_{1}{V}_{100}{ \mathcal C }\right)}^{1/2},\end{eqnarray} \tag{ 26 }$$

where ${ \mathcal C }\equiv {c}_{s}/(1\,\,\mathrm{km}\,{{\rm{s}}}^{-1})$. Nonlinear fragmentation occurs before the start of the rarefaction era when ${t}_{G}\leqslant {t}_{\mathrm{cross}}\,-{t}_{\mathrm{cool}}^{{\rm{O}}}\simeq {t}_{\mathrm{cross}}$. This inequality is fulfilled when

$$\begin{eqnarray}&&{c}_{s}\leqslant \displaystyle \frac{8\pi \,G{\rho }_{1}{r}_{c}^{2}}{{V}_{s}}\simeq 4.3\,\,\mathrm{km}\,{{\rm{s}}}^{-1}\,({V}_{100}^{-1}\,{n}_{1}\,{r}_{0.4}^{2})\end{eqnarray} \tag{ 27 }$$

(corresponding with the choosen scalings to T₃ ≤ 2750 K). In the case of very fast collisions, the gas must cool to increasingly lower temperatures for gravitational instability to act on short enough timescales. As soon as condition (27) is satisfied, the unstable slab will break up into circular fragments of preferred radii

$$\begin{eqnarray}&&\displaystyle \frac{{\lambda }_{G}}{2}=\pi {c}_{s}{t}_{G}\simeq 107\,\mathrm{pc}\,\left(\displaystyle \frac{{r}_{0.4}^{3}{n}_{1}}{{V}_{100}^{2}}\right)\end{eqnarray} \tag{ 28 }$$

and masses

$$\begin{eqnarray}&&{M}_{G}=\pi {{\rm{\Sigma }}}_{G}{\left(\displaystyle \frac{{\lambda }_{G}}{2}\right)}^{2}\simeq {10}^{6.15}\,\,{M}_{\odot }\,\left(\displaystyle \frac{{r}_{0.4}^{7}{n}_{1}^{3}}{{V}_{100}^{4}}\right).\end{eqnarray} \tag{ 29 }$$

While the latter is again comparable to the observed characteristic mass of GCs after accounting for mass loss by stellar evolution and tidal disruption after birth,¹² it is also somewhat ill-defined—scaling strongly with cloud size, gas density, and shock velocity. Neverthless, it is encouraging that collisions between atomic-cooling, metal-poor subhalos may offer a mechanism for imprinting the signature of the GC mass scale on the collapsing shell.

It should be noted that the conditions for gravitational instability given in Equations (27)–(29) can only be met as long as the gas metallicity Z is not much below 10⁻² Z_⊙ (see Figure 3). This may provide a plausible explanation for the threshold metallicity, [Fe/H] = −2.4, below which GCs are not observed. It is also significant that, at the time when gravitational instability sets in, the gas will be crushed into a slab of thickness d along the line connecting the cloud centers,

$$\begin{eqnarray}&&d=\displaystyle \frac{4{r}_{c}{T}_{3}}{{{ \mathcal M }}^{2}T}\simeq 3\,\mathrm{pc}\,({r}_{0.4}^{5}\,{n}_{1}^{2}\,{V}_{100}^{-4}),\end{eqnarray} \tag{ 30 }$$

which is comparable to the typical half-light radius of GCs. The freefall time for a uniform, pressure-free sphere of such gas,

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{ff}} & = & \sqrt{\displaystyle \frac{3\pi }{32G{\rho }_{3}}}=\displaystyle \frac{1}{{ \mathcal M }}\,\sqrt{\displaystyle \frac{3\pi {T}_{3}}{32G{\rho }_{1}T}}\\ & \simeq & 2.2\,\mathrm{Myr}({V}_{100}^{-2}{n}_{1}^{1/2}\,{r}_{0.4}^{2}),\end{array}\end{eqnarray} \tag{ 31 }$$

is shorter than the several-Myr evolutionary timescale for massive stars.

The above analysis is highly idealized, and is only meant to provide an idea of the general conditions under which the splash remnant may become Jeans unstable. It is clear from the previous discussion that one expects significant spatial structures in the postshock region both in temperature and density as a function of distance from the shock fronts (e.g., Smith 1980; Kang et al. 1990; Kumai et al. 1993), and the thin shell approximation may be inadequate in the case of very inhomogeneous flows (Yamada & Nishi 1998). A shocked slab is known to be susceptible to a number of hydrodynamical instabilities like the nonlinear thin shell instability (Vishniac 1994), which may compete with the gravitational instability and produce substructures on the scale of the slab thickness. The thermodynamic treatment of the problem has been simplified by neglecting molecular cooling, heating from external X-ray radiation and cosmic rays, dust processes, and gas clumping. The simple model of head-on collisions between identical, uniform-density clumps is obviously unrealistic, and should be extended to off-center impacts between subhalos with internal substructure.

Nevertheless, our calculations may elucidate the conditions for a special, dynamical mode of star formation following substructure collisions, a mode that is intimately tied to ΛCDM hierarchical galaxy assembly. It seems plausible that GC formation by impact requires the relative velocities of the colliding subhalos to be in a specific range. If the collision velocity is too low, shocks may not be able to produce the extremely high pressure environments that are a prerequisite to the formation of dense and tightly bound clusters. Conversely, if the velocity is too high and the shock too violent, the interacting clouds will expand and disperse before significant radiative cooling can occur. Equations (16) and (17) give

$$\begin{eqnarray}&&\displaystyle \frac{{t}_{\mathrm{cool}}}{{t}_{\mathrm{cross}}}\simeq 3.5\times {10}^{-3}\left(\displaystyle \frac{{V}_{100}^{4.6}}{{n}_{1}\,{r}_{0.4}}\right),\end{eqnarray} \tag{ 32 }$$

and t_cool/t_cross < 1 for ${V}_{s}\lt 343\,\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\left({n}_{1}{r}_{0.4}\right)}^{0.217}$. This upper velocity limit is only weakly dependent on the properties (gas density and size) of the original colliding clumps.

Simulations of interstellar cloud–cloud collisions may also provide some insight on the fate of the shock-compressed layer. According to Balfour et al. (2015), when cold, uniform-density clouds collide head-on at moderately supersonic speeds, star formation operates in a global hub-and-spoke mode that produces a central monolithic stellar cluster. At higher collision velocities, a spider's-web mode operates and delivers a loose distribution of independent, small subclusters instead. When the clouds have precollision substructure, however, the collision velocity becomes less critical (Balfour et al. 2017). In these numerical experiments the gas is evolved with a barotropic equation of state. Ultimately, the detailed fate of subhalo–subhalo high-speed collisions should be addressed with the help of cosmological hydrodynamic simulations that include all the relevant heating and cooling processes.

High-speed close interactions between satellites orbiting within a parent halo have been advocated as a major mechanism for the morphological evolution of galaxies in clusters (Moore et al. 1996; Baushev 2018), and collisions between protogalaxies have been proposed as a new pathway to form supermassive black holes at very high redshifts (Inayoshi et al. 2015). Cosmological N-body simulations by Tormen et al. (1998) showed that fast satellite–satellite encounters with impact parameter b < (r_vir,1 + r_vir,2) are fairly common and can lead to significant mass loss and disruption. In this section we make use of the Via Lactea I (VLI) N-body simulation to argue that rarer, nearly central collisions between atomic-cooling subhalos still occur frequently enough at high redshifts to represent a plausible pathway to the formation of GCs.

VLI is a DM-only simulation that follows the formation of an MW-sized halo in a ΛCDM cosmology. The high-resolution region was sampled with 234 million particles of mass $2.1\times {10}^{4}\,\,{M}_{\odot }$ and evolved with a force resolution of 90 pc starting from redshift 50 (Diemand et al. 2007a, 2007b). We stored and analyzed 200 outputs from redshift 16 to z = 0. The simulation has sufficient mass resolution to follow metal and atomic-cooling subhalos through many orbits and severe mass stripping, and sufficient output time resolution (Δt = 68.5 Myr) to measure the orbital parameters of subhalos with good accuracy. Such output spacing provides many timesteps from infall to the first pericenter passage per subhalo. Substructure catalogs were constructed using the phase-space group-finder 6DFOF (Diemand et al. 2006), subhalos were linked across snapshots to follow their histories and trajectories backward and forward in time, and all substructure–substructure close encounters were recorded. When tracing halos backward in time, a subhalo was linked to its main progenitor only if the core of the latter contained at least 50% of the particles of the core of the former, and vice versa.¹³

VLI initial conditions were generated with the original version of the GRAFIC2 package (Bertschinger 2001), which incorrectly used the baryonic instead of the DM power spectrum for the refinement levels, leading to reduced small-scale power. Compared to Via Lactea II (Diemand et al. 2008), the abundance of subhalos in VLI is suppressed by a factor of 1.7, and one should therefore view the derived collision frequency strictly as a lower limit. This is even more so as even state-of-the-art cosmological simulations still suffer from significant overmerging, and many subhalos will be artificially disrupted before a collision by numerical resolution effects (van den Bosch et al. 2018). The distance of closest approach between subhalos was found by linearly interpolating distances amid snapshots. Once a subhalo is accreted by the main host, its diffuse outer layers are rapidly stripped off by tidal forces, with tidal mass losses being more significant for more massive subhalos (Diemand et al. 2007b). The accurate tracking of the accretion history of substructure allows us to define the epoch when the satellite mass reached a maximum value, denoted here as m_peak, before infall.

We have identified all collisions between atomic-cooling systems (T_vir > 8000 K) that, at the instant of closest approach, are unbound and interpenetrating, i.e., involve subhalos with relative velocities v_rel > 3 V_max,1 and whose centers are separated by a distance $b\lt {b}_{\max }=0.1({r}_{\mathrm{vir},1}+{r}_{\mathrm{vir},2})$. Here V_max,1 is the maximum circular velocity of the larger ("1") of the two subhalos at impact, 3 V_max is the escape velocity from the center of a spherically averaged NFW density profile, and the chosen maximum impact parameter b_max corresponds to the sum of the gas core radii of the pair. The impact parameter condition ensures that the region of overlap between clumps colliding off-center is extensive and so is the resulting gas splash, and excludes events where small, dense clumps may plow through the tenuous halo of larger subhalos without much resistance, i.e., without their gas becoming dislodged. In order to minimize the effect of ram pressure stripping on the gaseous content of interacting substructure, the collisions of interest here are those that occur before the first pericenter passage of each subhalo. Clumps with correlated infall histories often undergo multiple collisions: for simplicity, we assume here that the smaller member of a colliding pair depletes its gas reservoir without replenishment during the first interpenetrating encounter, and do not count any close interactions it may be involved in afterwards.

Figure 4 shows the frequency distributions of redshifts, masses, relative velocities, impact parameters, angles of impact, and first pericenter distances for all encounters that satisfy the above criteria. We tally 133 collision events within today's VLI virial volume, corresponding to R_vir = 288 kpc (note that VLI virial radius drops below 40 kpc at z ≳ 3). If most of these encounters were to result in the formation of one or more massive GCs, the predicted frequency would be comparable to the ∼150 globular clusters observed today in the halo of the MW. If, on the other hand, many more low-mass GCs formed initially than is currently observed, and were selective destroyed by dynamical processes acting over a Hubble time (e.g., Gnedin & Ostriker 1997; Vesperini 1998), then our model would require either the formation of multiple GCs per impact, or a higher collision frequency, perhaps involving the more numerous population of T_vir < 8000 K subhalos.¹⁴ Collisions occur at early times in the transition region between the main host and the field: we count only 10 events within a Galactocentric distance of 50 kpc. For the most part, these collisions involve subhalos that are already dynamically associated before accretion into the main host, i.e., they are either part of the same infalling halo or two separate clumps descending along a filament and organized into small groups with correlated trajectories (e.g., Li & Helmi 2008; Angulo et al. 2009). In this situation, the distinction between ex situ and in situ globular clusters becomes blurred.

Zoom In Zoom Out Reset image size

The median masses at infall of the lighter and heavier member of the colliding pairs are 10^8.3 and 10^9.7 $\,{M}_{\odot }$, respectively. The median relative velocity is 160 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$, and there is a weak negative correlation (with correlation coefficient r = −0.3) between relative velocity and redshift: the few extreme-velocity impacts with v_rel ≳ 350 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$ all occur at z < 3.5, when the depth of the gravitational potential of the main host is larger. Collisions have typical impact parameters in the range 0.3 ≲ b ≲ 1.7 kpc (median value 0.8 kpc), and the angles of impact between the initial velocity vectors of the two bodies is less than 50 degrees in about half of all encounters. Most interacting subhalos have highly radial orbits and plunge deep into their host halo (see also Wetzel 2011; González & Padilla 2016). In about 40% of all collisions, one or both of the interacting subhalos are tidally disrupted before falling within 50 kpc of the center of the main host.

The presence of substructure within DM halos is a unique signature of a universe where systems grow hierarchically through the accretion of smaller-mass units. We have investigated a scenario where the formation of GCs is triggered by high-speed collisions between infalling, atomic-cooling subhalos during the early assembly of the main galaxy host. This is a special, dynamical mode of star formation that operates at extremely high gas pressures, is relatively immune from the feedback processes that regulate star formation in the field, and tracks the total mass of the main host rather than its stellar mass. The proposed mechanism would give origin to DM-less globular clusters, as colliding DM subhalos and their stars will simply pass through one another while the warm gas within them shocks and decouples. The well-known Bullet galaxy cluster provides a striking illustration of the process, albeit on different scales. In an MW-sized host, where the relative orbital velocities between subhalos substantially exceed their internal velocities, most close encounters are unbound collisions between satellites that are accreted at similar early times and locations, and occur in the outer regions of the main host.

Below, we summarize our main results and discuss some implications for GC formation models.

(i) Fragmentation of shocked gas and the masses of GCs: shock heating and cooling, the encounter geometry, and the complexities of multiphase gaseous inner halos are all key factors in determining the outcome of a subhalo–subhalo impact. We have shown that, under idealized conditions, the low-metallicity warm gas in the cores of interacting subhalos will be shock heated to characteristic temperatures ∼10⁵ K, and will cool rapidly first via atomic line emission and, further downstream, via fine-structure lines of C ii and O i. Because the collision is more dissipative than isothermal, the resulting shock-compressed slab is liable to gravitational instability. An idealized analysis in the thin shell approximation yields the conditions under which the imprint of the GC mass scale could be present on the cooling, collapsing shell. The requirements for gravitational fragmentation can only be satisfied for gas metallicities Z ≳ 10^−2.5 Z_⊙, thus offering a natural interpretation for the observed minimum metallicity of GCs. We caution, however, that these findings are based on a simplified thermodynamic treatment of the problem that neglects molecular cooling, heating from external X-ray radiation and cosmic rays, dust processes, and gas clumping.

(ii) The GC system–DM connection: an analysis of the scaling behavior of the encounter frequency within the kinetic theory approximation points to a GC number–halo mass relation that is the result of both encounter probability calculations—the subhalo collision rate per unit volume being the usual density²×cross section × velocity factor—and of hierarchical assembly—globular clusters being brought in by the accretion of smaller satellites. With the addition of ex situ globular clusters, the scaling of the total number of GCs with host virial mass is very close to linear, ${N}_{\mathrm{GC}}\propto {M}_{\mathrm{vir}}^{0.98}$, in agreement with the trend observed over five orders of magnitude in galaxy mass. This uniform GC production rate per unit host halo mass is predicted to break down on dwarf galaxy scales, perhaps below a critical mass of ∼10^9.6 $\,{M}_{\odot }$.

(iii) The ages of GCs: our model differs from much previous work as it does not assume an arbitrary value for the redshift when metal-poor GC formation is shut off. The details of the redshift distribution in top left panel of Figure 4 reflect the mass assembly history of the simulated MW-sized host system, but it is again noteworthy that a scenario in which GCs are the result of colliding substructures would produce a population of old clusters with typical ages >10 Gyr, a median age of 12 Gyr (corresponding in the adopted cosmology to a median redshift of 3.5), and an age spread that is similar to the one observed. In contrast to many pregalactic scenarios (e.g., Katz & Ricotti 2013; Kimm et al. 2016; Boylan-Kolchin 2017), in our model GCs have extended formation histories and typically form after the epoch of reionization: only about 38% of all close encounters occur at redshifts greater than 4. Seven collision events in our sample take place at z < 2. Four of these "late" impacts have relative velocities in excess of 350 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$, a situation that may not be conducive to the cooling and fragmention of the splash remnant (see Section 3.1). The others may give origin to a population of young metal-poor globular clusters like Crater (Weisz et al. 2016). The dotted–dashed histogram on the same panel shows the relative frequency histogram of the redshift of formation for 55 MW GCs with Hubble Space Telescope photometry (VandenBerg et al. 2013), augmented by the age determinations for the young globular clusters Crater, Pal 1, Terzan 7, and Whiting 1 (Sarajedini et al. 2007; Dotter et al. 2010; Valcheva et al. 2015; Weisz et al. 2016). The distribution is poorly known at z > 4 as the typical uncertainty on the absolute age of GCs exceeds 1–1.5 Gyr.

(iv) The metallicity of GCs: the GC population in the MW is observed to be clearly bimodal, with a low-metallicity component peaking at [Fe/H] ≃ −1.55, and a high-metallicity tail at ≃ − 0.55 (Harris et al. 2016). Only 30% of MW globular cluster have [Fe/H] > −1, but many massive galaxies possess strongly bimodal GC systems, with nearly equal numbers of metal-rich and metal-poor clusters (Peng et al. 2006). Since stars within most GCs do not show an internal spread in iron-peak elements, there must exist a mechanism that chemically homogenizes the gas within a protocluster before the onset of star formation. The dominant mode of chemical mixing is thought to be turbulent diffusion (Murray & Lin 1990), which has been shown to produce a stellar abundance scatter that is much smaller than that of the star-forming gas (Feng & Krumholz 2014). In our model, GCs will inherit the gas-phase metallicity of the interacting subhalo pair that triggers their formation, and a useful perspective can be obtained by assigning a stellar mass to each subhalo at infall following the median stellar-to-halo mass relation from Behroozi et al. (2013). This redshift-dependent prescription, extrapolated to the small scales of interest here, leads to a broad range of stellar masses, 10^2.6 < m_*/$\,{M}_{\odot }$ < 10^8.6 for our sample of colliding substructures, a distribution that extends over nearly 6 decades in mass with a median value equal to 10^4.4 $\,{M}_{\odot }$ (see Figure 5).

Zoom In Zoom Out Reset image size

A general tendency of decreasing metallicity toward lower stellar masses is commonly accepted, but the exact form of the stellar mass (m_*) versus gas-phase metallicity (Z) relation (hereafter MZR) and its evolution with redshift are currently poorly known as a result of the presence of strong systematic uncertainties affecting metallicity diagnostics. A few studies, mostly at z ∼ 0, have tried to extend the MZR to the low-mass dwarf galaxy regime, deriving power-law relations ($Z\propto {m}_{* }^{\alpha }$) with slopes α = 0.29 ± 0.03 (Berg et al. 2012), ≃0.5 (Andrews & Martini 2013), 0.44 ± 0.1 (Jimmy et al. 2015, high star formation rate bin), and 0.14 ± 0.08 (Blanc et al. 2019). Here, we adopt the intermediate-range value of α = 0.37 (Blanc et al. 2019; Ma et al. 2015) with Y-intercept equal to 5.50, i.e.,

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}(Z/{Z}_{\odot }) & = & 12+\mathrm{log}({\rm{O}}/{\rm{H}})-8.69\\ & = & 5.50+0.37\mathrm{log}({m}_{* }/\,{M}_{\odot })-8.69.\end{array}\end{eqnarray} \tag{ 33 }$$

The zero-point was chosen to give $12+\mathrm{log}({\rm{O}}/{\rm{H}})=8.18$ at $\mathrm{log}({m}_{* }/{M}_{\odot })=7.25$, in agreement with the DEEP2 low-mass MZR (Zahid et al. 2012; Blanc et al. 2019). Little is known about the redshift evolution of the MZR on dwarf galaxy scales; in order to minimize the number of free parameters and for ease of interpretation, we assume no evolution with time in the following (see also Hidalgo 2017).

The expected GC metallicity, computed by averaging the gas metallicities of each colliding pair, is shown in Figure 5, together with the observed distribution of [Fe/H] in Galactic globular clusters, both metal-poor and metal-rich (Harris 1996). Our simple scheme predicts a spread in metallicities that is similar to that observed, with a distribution that is strongly bimodal. In contrast to some previous work, we do not explicitly assume separate pathways for the formation of blue and red globular clusters, and simply predict the metallicity of each GC based on the enrichment level of the interacting subhalo pair that triggered its formation. About 30% of the colliding pairs have metallity $\mathrm{log}(Z/{Z}_{\odot })\geqslant -1$: the red globular clusters are the result of impacts that involve at least one massive, more chemically evolved satellite, and the metallicity bimodality reflects a bimodality in stellar and peak halo masses for the colliding subhalos. The age spread is similar in both the blue and red populations, but the red and blue peaks are shifted toward lower values compared to the MW GC data.

Given the intrinsic uncertainties in the stellar-to-halo mass relation, the MZR, and their evolution with redshift on small-mass galaxy scales, it seems again ill-advised to draw definite conclusions from this comparison. Evolutionary corrections on the MZR will shift the predicted values toward even lower metallicities. We notice here two effects that have been ignored and may skew the predicted distribution in the opposite direction, toward higher metallicities: (1) collisions involving massive, enriched subhalos that, albeit rarer, could produce several GCs per event. The impact of this (unknown) multiplicity factor has been neglected in Figure 5; and (2) high-speed bound collisions between subhalos falling in on radial orbits and the central most massive progenitor, which may result in a more centrally concentrated subpopulation of metal-rich globular clusters (see also Griffen et al. 2010). We count only 17 such collisions, compared to the 133 subhalo–subhalo impacts that satisfy our criteria.

(v) The spatial distribution of GCs: a common feature of many globular cluster formation models is the reliance on some ad hoc assumptions to identify the sites where GCs form. The median Galactocentric distance of all known MW (blue+red) GCs is 5 kpc (Harris 1996), with the metal-poor, [Fe/H] <−1, subpopulation being less spatially concentrated (median Galactocentric distance of 7.5 kpc). In our sample of colliding subhalos, we find that 80% of all clumps surviving complete disruption have first pericenter distances smaller than 20 kpc, with a median value of 12 kpc (see the bottom right panel in Figure 4). Models in which the bulk of the metal-poor GC subpopulation formed in satellite systems—many of which are now tidally disrupted—and were subsequently accreted onto the main galaxy tend to produce, however, clusters with a more extended spatial distribution than observed (Muratov & Gnedin 2010; Creasey et al. 2019), unless they are associated with rare, early progenitor halos at z ∼ 10 (Moore et al. 2006; Katz & Ricotti 2014). Our collision-driven scenario may offer a new mechanism for biasing the spatial distribution of GCs relative to the overall mass profile. This is because, in an inelastic collision, the splash remnant will lose orbital energy and fall deeper into the Galactic potential rather than sharing the orbits of the progenitor subhalos. It is interesting to briefly examine here the impact of kinetic energy dissipation on orbital parameters. Consider, for simplicity, two subhalos moving on coplanar, coaxial, prograde elliptical orbits in a Keplerian potential of gravitational parameter μ. The two orbits have the same specific angular momentum h and the subhalos collide at the semiparameter location of the ellipses, p = h²/μ. Let us decompose the velocity vectors along the outward radial direction ($\hat{r}$) and the prograde azimuthal direction ($\hat{\theta }$) at this position. The radial and azimuthal velocity components of the two subhalos just prior to impact are (μe₁/h, μ/h) and (μe₂/h, μ/h), respectively, where e₁ and e₂ are the orbital eccentricities, with e₁ < e₂. In the case of a perfectly inelastic encounter between two bodies of equal mass m_c, conservation of linear momentum determines the postcollision instantaneous orbital velocity vector ${{\boldsymbol{v}}}_{f}$ at p of the combined remnant as

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{f}=\displaystyle \frac{1}{2}({{\boldsymbol{v}}}_{1}+{{\boldsymbol{v}}}_{2})=\displaystyle \frac{\mu }{2h}\,[({e}_{1}+{e}_{2})\hat{r}+2\hat{\theta }].\end{eqnarray} \tag{ 34 }$$

In the absence of external torques the total angular momentum of the system remains unchanged, but the dissipation of kinetic energy leads to a net loss of orbital energy

$$\begin{eqnarray}&&{\rm{\Delta }}E=-\displaystyle \frac{{m}_{c}{\mu }^{2}}{4{h}^{2}}{\left({e}_{2}-{e}_{1}\right)}^{2}\leqslant 0.\end{eqnarray} \tag{ 35 }$$

The new specific orbital energy is

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{1}{2}{v}_{f}^{2}-\displaystyle \frac{{\mu }^{2}}{{h}^{2}}=\displaystyle \frac{{\mu }^{2}}{2{h}^{2}}\left[\displaystyle \frac{{\left({e}_{1}+{e}_{2}\right)}^{2}}{4}-1\right],\end{eqnarray} \tag{ 36 }$$

and the new eccentricity is

$$\begin{eqnarray}&&e=\sqrt{1+2{\epsilon }_{f}{h}^{2}/{\mu }^{2}}=\displaystyle \frac{{e}_{1}+{e}_{2}}{2}.\end{eqnarray} \tag{ 37 }$$

Another simple situation to analyze is the case of two clumps moving on coplanar, coaxial prograde orbits with different angular momenta but the same pericenter distance r_p. Conservation of angular momentum in this case yields for the new eccentricity of the collision remnant

$$\begin{eqnarray}&&\sqrt{1+e}=\displaystyle \frac{\sqrt{1+{e}_{1}}+\sqrt{1+{e}_{2}}}{2},\end{eqnarray} \tag{ 38 }$$

and the change in orbital energy is again

$$\begin{eqnarray}&&{\rm{\Delta }}E=-\displaystyle \frac{{m}_{c}\mu }{2{r}_{p}}({e}_{2}+{e}_{1}-2e)\leqslant 0.\end{eqnarray} \tag{ 39 }$$

The above examples emphasize the fact that the remnant will initially be moving on a less energetic orbit with an eccentricity that is intermediate between those of the colliding pair. One would expect globular clusters to progressively lose memory of their initial infall direction as they orbit in a host halo that is clumpy and triaxial. Detailed calculations of the expected radial profile and 3D distribution of GCs are beyond the scope of this paper and are postponed to future work.

We thank Z. Haiman and M. Krumholz for very helpful comments and suggestions on various aspects of this manuscript. Support for this work was provided by NASA through a contract to the WFIRST-EXPO Science Investigation Team (15-WFIRST15-0004), administered by the GSFC (P.M.).