Globular Cluster Formation from Colliding Substructure

TL;DR

This paper proposes a novel model where globular clusters form from high-speed collisions of infalling subhalos during galaxy assembly, leading to observable properties consistent with real globular cluster populations.

Contribution

It introduces a new formation mechanism for globular clusters triggered by colliding subhalos, supported by simulations and theoretical analysis.

Findings

Produces a linear GC-halo mass relation matching observations

Generates old globular populations with realistic age distributions

Results in a bimodal metallicity distribution similar to observed galaxies

Abstract

We investigate a scenario where the formation of Globular Clusters (GCs) is triggered by high-speed collisions between infalling atomic-cooling subhalos during the assembly of the main galaxy host, a special dynamical mode of star formation that operates at high gas pressures and is intimately tied to LCDM hierarchical galaxy assembly. The proposed mechanism would give origin to "naked" globulars, as colliding dark matter subhalos and their stars will simply pass through one another while the warm gas within them clashes at highly supersonic speed and decouples from the collisionless component, in a process reminiscent of the Bullet galaxy cluster. We find that the resulting shock-compressed layer cools on a timescale that is typically shorter than the crossing time, first by atomic line emission and then via fine-structure metal-line emission, and is subject to gravitational instability and fragmentation. Through a combination of kinetic theory approximation and high-resolution $N$-body simulations, we show that this model may produce: (a) a GC number-halo mass relation that is linear down to dwarf galaxy scales and agrees with the trend observed over five orders of magnitude in galaxy mass; (b) a population of old globulars with a median age of 12 Gyr and an age spread similar to that observed; (c) a spatial distribution that is biased relative to the overall mass profile of the host; and (d) a bimodal metallicity distribution with a spread similar to that observed in massive galaxies.

Equations

$$\frac{d{N}_{\mathrm{sub}}}{d\xi }={\xi }^{-1}\left(a_1{\xi }^{-{\alpha }_1}+a_2{\xi }^{-{\alpha }_2}\right)\exp (-b{\xi }^{\beta }),$$

$$n_{\mathrm{sub}}(\xi ,x)\propto \frac{1}{cx(1+cx{)}^2},$$

$$4\pi R_{\mathrm{vir}}^3{n}_{\mathrm{sub}}(>\xi ,x)=\frac{c^3}{gcx(1+cx{)}^2}{\int }_{\xi }^1\frac{d{N}_{\mathrm{sub}}}{d{\xi }^{\prime }}d{\xi }^{\prime },$$

$$f(v_{\mathrm{rel}})=\frac{1}{2\sqrt{\pi }{\sigma }_v^3}{v}_{\mathrm{rel}}^2\exp \left(-\frac{v_{\mathrm{rel}}^2}{4{\sigma }_v^2}\right),$$

$$\begin{array}{rl}{\sigma }_v^2(x)=& \left(\frac{G{M}_{\mathrm{vir}}}{R_{\mathrm{vir}}}\right)g(1+cx{)}^{2}x\\ & \times {\int }_x^{\infty }\left[\frac{\ln (1+cs)}{s^3(1+cs{)}^2}-\frac{c}{s^2(1+cs{)}^3}\right]ds.\end{array}$$

$$\begin{array}{rl}{N}_{\mathrm{coll}}=& \Delta t4\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}$$

$$N_{\mathrm{coll}}\simeq 59\left(\frac{\Delta t}{\mathrm{Gyr}}\right).$$

$$N_{\mathrm{coll}}\propto R_{\mathrm{vir}}^{-3}{M}_{\mathrm{vir}}^{1.9}(M_{\mathrm{vir}}/R_{\mathrm{vir}}{)}^{1/2}{c}^4.$$

$$N_{\mathrm{in}-\mathrm{s}}^{\mathrm{GC}}=A(M_{\mathrm{vir}}/10^{12}{\mathrm{M}}_{\odot }{)}^{0.87},$$

$$N_{\mathrm{ex}-\mathrm{s}}^{\mathrm{GC}}=N_{\mathrm{in}-\mathrm{s}}^{\mathrm{GC}}{\int }_{{\xi }_{\min }}^1{\xi }^{0.87}\frac{d{N}_{\mathrm{sub}}}{d\xi }d\xi ,$$

$$N_{\mathrm{GC}}\propto M_{\mathrm{vir}}^{0.98}.$$

$$m_{\mathrm{peak}}>1.6\times 10^8{\mathrm{M}}_{\odot }[\mu (1+z)/5{]}^{-3/2}{\mathcal{B}}^{-1/2}.$$

$$r_{\mathrm{vir}}\simeq 4.4\mathrm{kpc}{\left(\frac{m_{\mathrm{peak}}}{10^{8.5}{\mathrm{M}}_{\odot }}\right)}^{1/3}\left(\frac{5}{1+z}\right)$$

$$\begin{array}{rl}2m_c=& \frac{8\pi }{3}{\rho }_1{r}_c^3\simeq 10^7{\mathrm{M}}_{\odot }{n}_1{\left(\frac{r_c}{0.1r_{\mathrm{vir}}}\right)}^3\times \\ & m_{8.5}{\left(\frac{5}{1+z}\right)}^3,\end{array}$$

$$T_2=\frac{3\mu m_p}{16k_B}{V}_s^2\simeq 1.3\times 10^5\mathrm{K}{V}_{100}^2,$$

$$t_{\mathrm{cool}}=\frac{3k_B{T}_2{n}_{\mathrm{tot}}}{2n_{\mathrm{H},2}^2\Lambda }\simeq 0.03\mathrm{Myr}\left(\frac{V_{100}^{3.6}}{n_1}\right).$$

$$t_{\mathrm{cross}}\equiv \frac{2r_c}{V_s}\simeq 7.8\mathrm{Myr}(V_{100}^{-1}{r}_{0.4}),$$

$$\mathcal{M}\equiv \frac{V_s}{c_s}\simeq 13.6(T_{3.9}^{-1/2}{V}_{100}),$$

$$n_{\mathrm{H},2}={\mathcal{M}}^2{n}_{\mathrm{H},1}\simeq 140\mathrm{c}{\mathrm{m}}^{-3}(n_1{T}_{3.9}^{-1}{V}_{100}^2),$$

$$\begin{array}{rl}{M}_{\mathrm{BE}}=& 1.182\frac{c_s^3}{\sqrt{G^3{\mathcal{M}}^2{\rho }_1}}\\ \simeq & 7.7\times 10^5{\mathrm{M}}_{\odot }(n_1^{-1/2}{T}_{3.9}^2{V}_{100}^{-1}).\end{array}$$

$$n_{\mathrm{H},2}^2{\Lambda }_{\mathrm{O}}=n_{\mathrm{H},2}^2{\mathcal{A}}_{\mathrm{O}}{k}_{\mathrm{OI}}(k_B{T}_{\mathrm{OI}}),$$

$$t_{\mathrm{cool}}^{\mathrm{O}}=\frac{5k_{B}T{n}_{\mathrm{tot}}}{2n_{\mathrm{H},2}^2{\Lambda }_{\mathrm{O}}}\simeq 0.4\mathrm{Myr}\left(\frac{T_{3.9}^{1.59}}{V_{100}^2{n}_1{Z}_{-1.3}}\right),$$

$${\lambda }_G=\frac{2c_s^2}{G\Sigma },$$

$$t_G=\frac{c_s}{\pi G\Sigma }.$$

$$t_G={\left(\frac{c_s}{2\pi G{\rho }_1{V}_s}\right)}^{1/2}\simeq 3.8\mathrm{Myr}{\left(\frac{\mathcal{C}}{n_1{V}_{100}}\right)}^{1/2},$$

$${\Sigma }_G={\left(\frac{2{\rho }_1{V}_s{c}_s}{\pi G}\right)}^{1/2}\simeq 19{\mathrm{M}}_{\odot }\mathrm{p}{\mathrm{c}}^{-2}(n_1{V}_{100}\mathcal{C}{)}^{1/2},$$

$$c_s\le \frac{8\pi G{\rho }_1{r}_c^2}{V_s}\simeq 4.3\mathrm{km}{\mathrm{s}}^{-1}(V_{100}^{-1}{n}_1{r}_{0.4}^2)$$

$$\frac{{\lambda }_G}{2}=\pi c_s{t}_G\simeq 107\mathrm{pc}\left(\frac{r_{0.4}^3{n}_1}{V_{100}^2}\right)$$

$$M_G=\pi {\Sigma }_G{\left(\frac{{\lambda }_G}{2}\right)}^2\simeq 10^{6.15}{\mathrm{M}}_{\odot }\left(\frac{r_{0.4}^7{n}_1^3}{V_{100}^4}\right).$$

$$d=\frac{4r_c{T}_3}{{\mathcal{M}}^{2}T}\simeq 3\mathrm{pc}(r_{0.4}^5{n}_1^2{V}_{100}^{-4}),$$

Full text

1. Introduction

The origin of globular clusters (GCs) remains an unsolved problem in star and galaxy formation studies. With masses in the range $10^{4}$–$10^{6}\,\,{\rm M_{\odot}}$, half-light radii of a few pc, 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 globulars 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; Forbes et al., 2016; Harris et al., 2017; Hudson & Robison, 2018, and references therein), a trend the former 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 \mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 318\relax}}3.5 kpc (Bica et al., 2006). Only six known GCs inhabit the Galaxy at distances \mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 318\relax}}90 kpc (Harris, 1996; Laevens et al., 2014). The red (metal-rich) globulars are more spatially concentrated than the blue clusters: they lie within the solar circle and form a flattened, rotating population. Blue GCs outnumber red globulars 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 globulars 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\pm 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$ kpc, is well-described by a simple stellar population with an age of $7.5\pm 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\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 316\relax}}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 globulars 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^{6}\,\,{\rm 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 $\Lambda$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 $\sim 10^{5}$ 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\sim 10^{7}\,k_{B}$ cm${}^{-3}\,$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*-3* were to collide at the Galactic orbital speeds of $200\,\,{\rm km\,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.

2. Subhalo-Subhalo Collisions

Numerical simulations in $\Lambda$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 (Moore et al., 1999; Klypin et al., 1999; Diemand et al., 2007a, 2008; Springel et al., 2008). About one thousand subhalos with pre-infall masses m_{\rm peak}\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 318\relax}}10^{8}\,\,{\rm M_{\odot}} are predicted to have been accreted by our Galaxy over cosmic time (see Han et al., 2018, and references therein).111To be more precise, $m_{\rm peak}$ is the maximum mass attained by a subhalo over its entire merger history. This mass is greater than the subhalo mass at infall, as most subhalos start being stripped at \mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 316\relax}}2 host virial radii, regardless of host mass (Behroozi et al., 2014). Infalling halos generally become accreted shortly after reaching peak mass.  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 H2-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 photo-evaporate 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\sigma$ upper bound to the pre-infall peak subhalo mass of Segue I – the faintest MW satellite – of $m_{\rm peak}<2.4\times 10^{8}\,\,{\rm 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^{7}-10^{8}\,\,{\rm 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 \mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 318\relax}}10^{5}\,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 pre-infall mass, i.e. by using the “unevolved” substructure mass function – the distribution of peak bound masses $m_{\rm 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 $\Lambda$CDM cosmological simulations, this universal redshift-independent function is well fitted by a double Schechter function of the form (Han et al., 2018)

[TABLE]

where $\xi=m_{\rm peak}/M_{\rm vir}$ is the ratio of the peak mass of a subhalo to the host halo mass, and $(dN_{\rm sub}/d\xi)d\xi$ is the number of subhalos with masses between $\xi$ and $\xi+d\xi$. When adopting for the host mass the virial definition corresponding to the overdensity of the spherical collapse model, the best-fit parameters are $a_{1}=0.11,a_{2}=0.32,\alpha_{1}=0.95,\alpha_{2}=0.08,b=8.9,$ and $\beta=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,

[TABLE]

where $x\equiv r/R_{\rm vir}$, $R_{\rm vir}$ is the virial radius of the host, and $c$ is its “concentration parameter”.222In reality, the dynamics of subhalos differs somewhat from that of DM particles because of dynamical friction, but this is a small effect for the majority of subhalos that are located in the outer regions of the main host and at the low-mass end (Han et al., 2016).  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 $>\xi M_{\rm vir}$, $n_{\rm sub}(>\xi,x)$, can be written as

[TABLE]

where $g\equiv 1/[\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 $\sigma_{v}$. The distribution of relative encounter speeds is then

[TABLE]

and the mean encounter speed is $\bar{v}_{\rm 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)

[TABLE]

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_{\rm vir}$) $>{\xi}_{1}$ and $>{\xi}_{2}$ in a time interval $\Delta t$ as

[TABLE]

where $b_{\rm max}$ is the maximum impact parameter for a close collision.

Consider now, for illustrative purposes, a MW-sized host halo of mass $M_{\rm vir}=10^{12}\,\,{\rm 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\,\,{\rm km\,s^{-1}}$ at $x=0.08$, corresponding to a mean collision speed of $230\,\,{\rm km\,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_{\rm max}=0.1(r_{\rm vir,1}+r_{\rm vir,2})$, where $r_{{\rm vir},1}$ and $r_{\rm vir,2}$ are the pre-infall virial radii of the two clumps. The mean number of impacts within $R_{\rm vir}$ between subhalos having peak masses $m_{\rm peak}>10^{8}\,\,{\rm M_{\odot}}$ is then

[TABLE]

This is clearly an interesting figure – comparable after $\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 globulars. 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_{\rm peak}>10^{8}\,\,{\rm M_{\odot}}$ and satellites that are ten times more massive is 3.6 times smaller than Equation (7). There are only $N_{\rm coll}\sim 3\,(\Delta t/{\rm Gyr})$ close collisions between massive satellites with $m_{\rm peak}>10^{9}\,\,{\rm M_{\odot}}$, but these may give origin to multiple GCs per encounter. The same calculation predicts about 1 collision Gyr*-1* between atomic-cooling subhalos in a dwarf galaxy-sized host with $M_{\rm vir}=10^{10}\,\,{\rm 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 § 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 non-linear 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; El-Badry et al., 2019; Burkert & Forbes, 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_{\rm peak}\ll M_{\rm vir}$ at fixed concentration parameter is proportional to $n_{\rm sub}^{2}\propto R_{\rm vir}^{-6}\,M_{\rm vir}^{1.9}$ (from Eqs. 1 and 3) times $\bar{v}_{\rm rel}\propto(M_{\rm vir}/R_{\rm vir})^{1/2}$ (see Eq. 5). Including a dependence on halo concentration and assuming that the collisions of interest continue for a timescale $\Delta t$ that is independent of host mass, one infers from Equation (6):

[TABLE]

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

[TABLE]

where $A$ is an arbitrary normalization factor. In this scenario, $N^{\rm GC}_{\rm in-s}$ 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 globulars that formed ex-situ and were subsequently accreted. The $N^{\rm GC}_{\rm in-s}-M_{\rm vir}$ relation given above and the substructure mass function in Equation (1) allow a straightforward calculation of the mean number of accreted GCs as

[TABLE]

where $\xi_{\rm min}\equiv m_{\rm min}/M_{\rm vir}$ is the critical subhalo mass below which GCs cannot form. Obviously, collisions between $m_{\rm peak}>10^{8}\,\,{\rm M_{\odot}}$ clumps as a mechanism for cluster formation require at least a few (say $\sim 4$) atomic-cooling subhalos in a given host; this trivial condition yields, after integrating Equation (1), $m_{\rm min}\simeq 10^{9.6}\,\,{\rm M_{\odot}}$. This critical mass could be larger – thereby decreasing the fraction of accreted globulars – if the relative orbital velocity $v_{\rm rel}$ required to trigger the formation of a GC by impact had to exceed $\sim 50\,\,{\rm km\,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_{\rm min}=10^{9.6}\,\,{\rm M_{\odot}}$ do not contain a GC, Equation (10) gives $N^{\rm GC}_{\rm ex-s}/N^{\rm GC}_{\rm in-s}=0.35,0.73,1.17$ for $M_{\rm vir}=10^{11},10^{12},10^{13}\,\,{\rm 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.333Similar trends are seen in recent phenomenological models of GC formation based on DM merger trees (e.g. Boylan-Kolchin, 2017).  With the addition of the accreted globulars, the scaling of the total number of GCs, $N_{\rm GC}=N^{\rm GC}_{\rm ex-s}+N^{\rm GC}_{\rm in-s}$, with host virial mass becomes very close to linear,

[TABLE]

Therefore, in a collision-driven scenario, a constant GC number-to-halo mass ratio is the result of encounter probability calculations (Eq. 8) and hierarchical clustering (Eq. 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 Forbes et al. 2018a; Harris et al. 2017), 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_{\rm vir}\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 316\relax}}10^{10}\,\,{\rm M_{\odot}}. With the adopted parameters, a MW-mass system with $M_{\rm vir}=10^{12}\,\,{\rm 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 globulars 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).

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 globulars. 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.

3. GC Formation in High-Speed Impacts

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 post-shock 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_{\rm vir}>8000$ K,

[TABLE]

Here, $\mu$ is the mean molecular weight per particle ($\mu=1.23$ for neutral primordial gas), ${\cal B}\equiv\Delta_{c}/(\Omega_{m}^{z}\,18\pi^{2})$, $\Delta_{c}=18\pi^{2}+82y-39y^{2}$ is the redshift-dependent density contrast at virialization (Bryan & Norman, 1998), $y\equiv\Omega_{m}^{z}-1$, $\Omega_{m}^{z}=\Omega_{m}(1+z)^{3}/[\Omega_{m}(1+z)^{3}+\Omega_{\Lambda}]$, and $z$ is the infall redshift. In evaluating these expressions, we have assumed a Planck Collaboration et al. (2018) flat cosmology with parameters $(\Omega_{m},\Omega_{\Lambda},h)=(0.315,0.685,0.674)$, and ${\cal 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 generation 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}\,\,{\rm cm^{-3}}\equiv n_{1}\,\,{\rm cm^{-3}}$ within radii $r_{c}\sim 0.1\,r_{\rm vir}$ (roughly independent of halo mass) (Ricotti, 2009; Prieto et al., 2013; Visbal et al., 2014a), and follows an $r^{-2}$ profile in the outer regions. Here,

[TABLE]

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 et al., 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_{\rm H,1}=0.75n_{c}$, mass density $\rho_{1}=n_{c}m_{p}$, temperature $T_{1}=8000$ K, and radius $r_{c}=0.1\,r_{\rm vir}$. These are the density and thermal pressure conditions of the warm diffuse component of the interstellar medium. 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,

[TABLE]

where $m_{8.5}\equiv m_{\rm peak}/10^{8.5}\,\,{\rm M_{\odot}}$, is about 10 per cent 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}\simeq 2v_{\rm rel}/3$ relative to the unpertubed medium, heating the gas to

[TABLE]

where $V_{100}\equiv V_{s}/100\,\,{\rm km\,s^{-1}}$ and we have set the molecular weight to $\mu=0.59$. The density is enhanced by a factor 4, the post-shock pressure is $p_{2}=(3/4)\rho_{1}V_{s}^{2}$, and the gas cools via hydrogen, helium, and metal line emission on a (isobaric) timescale

[TABLE]

Here, $n_{\rm tot}$ is the total particle number density ($n_{\rm tot}/n_{\rm H,2}=9/4$ for completely ionized gas of primordial composition), $n_{\rm H,2}=4\,n_{\rm H,1}$ is the post-shock mean hydrogen density, $\Lambda$ is the cooling function (e.g., Wang et al., 2014), and we have assumed collisional ionization equilibrium and ignored gas clumping (see Fig. 2). Since $t_{\rm cool}$ is typically much shorter than the crossing time,

[TABLE]

where $r_{0.4}\equiv r_{c}/0.4\,{\rm kpc}$, the internal energy increase caused by the shock is radiated away and the impact is effectively isothermal ($T_{2}=T_{1}=T$), characterized by a sound speed $c_{s}=\sqrt{k_{B}T/\mu m_{p}}$ and a Mach number

[TABLE]

where $T_{3.9}\equiv T/8000\,{\rm K}$ and we have assumed $\mu=1.23$. The post-shock hydrogen density is

[TABLE]

and the shock velocity relative to the unpertubed medium is now $V_{s}\simeq v_{\rm rel}/2$. The collision is also approximately one-dimensional, since the lateral rarefaction timescale is $\sim{\cal M}\,t_{\rm cross}$. For a strictly isothermal collision, when the compression phase is completed after a time $t_{\rm 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_{\rm cross}/{\cal M}$. If the collision is to result in gravitational instability and collapse, there must be density perturbations that can grow on a timescale \mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 316\relax}}t_{\rm cross}. We shall see below, however, that gas further downstream will actually cool below its pre-shock temperature via fine-structure metal lines, i.e. the collision is “more dissipative than isothermal” (e.g. Whitworth & Clarke, 1997).

Ly$\alpha$ cooling becomes inefficient below 8000 K, and the fast shock will dissociate most pre-existing 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 in the main host – will act to suppress H2 formation in the post-shock gas, and the temperature will then plateau around 8000 K (e.g., Kang et al., 1990). The Bonnor-Ebert mass in such warm, dense, and pressurized medium is

[TABLE]

a factor of $1/{\cal M}$ smaller than its value in the pre-shock gas (where $M_{\rm BE}>m_{c}$, i.e. the pre-collision clumps are gravitionally stable), and comparable to the observed mass of GCs.

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

[TABLE]

where ${\cal A}_{\rm O}$ is the abundance of O relative to hydrogen, ${\cal 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,

[TABLE]

is still shorter than the crossing time, but it is considerably longer than $t_{\rm 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_{\rm cool}^{\rm O}$, the flow will cool below its pre-shock 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 non-equilibrium chemistry package krome (Grassi et al., 2014) to construct more sophisticated models of the cooling and thermal properties of the post-shock gas. The thermal model includes non-equilibrium cooling for H and He at all temperatures and for metals below $10^{4}$ 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 post-shock temperature for a case with $V_{s}=100\,\,{\rm km\,s^{-1}}$ and $n_{\rm H,1}=0.75\,\,{\rm 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_{3}<T$ and hydrogen densities $n_{\rm H,3}=n_{\rm H,2}T/T_{3}$ (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.

The one-dimensional compression will create a cold layer of thickness-to-diameter ratio $\sim 2n_{\rm H,1}/n_{\rm H,3}<2/{\cal 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 $\Sigma$ 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

[TABLE]

and

[TABLE]

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 $\Sigma(t)=2\rho_{1}V_{s}t$, 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 $\Sigma_{G}\equiv\Sigma(t_{G})=2\rho_{1}V_{s}t_{G}$ (Kang et al., 1990). Using Equation (24) one derives then

[TABLE]

and

[TABLE]

where ${\cal C}\equiv c_{s}/(1\,\,{\rm km\,s^{-1}})$. Non-linear fragmentation occurs before the start of the rarefaction era when $t_{G}\leq t_{\rm cross}-t_{\rm cool}^{\rm O}\simeq t_{\rm cross}$. This inequality is fulfilled when

[TABLE]

(corresponding with the choosen scalings to $T_{3}\leq 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

[TABLE]

and masses

[TABLE]

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,444Note that some scenarios for the formation of multiple stellar populations within GCs require globulars to have been initially much more massive than they are today, (e.g. D’Ercole et al., 2008).  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^{-2}\,Z_{\odot}$ (see Fig. 3). This may provide a plausible explanation for the threshold metallicity, [Fe/H]$=-2.4$ dex, 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,

[TABLE]

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

[TABLE]

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 post-shock 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 non linear 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 $\Lambda$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

[TABLE]

and $t_{\rm cool}/t_{\rm cross}<1$ for $V_{s}<343\,\,{\rm km\,s^{-1}}~{}(n_{1}\,r_{0.4})^{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 sub-clusters instead. When the clouds have pre-collision 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.

4. N-body Simulation

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_{\rm vir,1}+r_{\rm 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 dark matter-only simulation that follows the formation of a MW-sized halo in a $\Lambda$CDM cosmology. The high-resolution region was sampled with 234 million particles of mass $2.1\times 10^{4}\,\,{\rm M_{\odot}}$ and evolved with a force resolution of 90 pc starting from redshift 50 (Diemand et al., 2007a, b). 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 ($\Delta 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 viceversa.555The phase-space group-finder 6DFOF finds peaks in phase-space density and links the particles within those peaks into groups. These groups correspond to the cores of halos and subhalos, whose total extent stretches well beyond these cores and is found by 6DFOF in a second step (Diemand et al., 2006).

VLI initial conditions were generated with the original version of the GRAFIC2 package (Bertschinger, 2001), which incorrectly used the baryonic instead of the dark matter 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_{\rm peak}$, before infall.

We have identified all collisions between atomic cooling systems ($T_{\rm vir}>8000\,$K) that, at the instant of closest approach, are unbound and interpenetrating, i.e. involve subhalos with relative velocities $v_{\rm rel}>3\,V_{\rm max,1}$ and whose centers are separated by a distance $b<b_{\rm max}=0.1\,(r_{\rm vir,1}+r_{\rm vir,2})$. Here $V_{\rm max,1}$ is the maximum circular velocity of the larger (“1”) of the two subhalos at impact, $3\,V_{\rm max}$ is the escape velocity from the center of a spherically-averaged NFW density profile, and the choosen maximum impact parameter $b_{\rm 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_{\rm vir}=288\,$kpc (note that VLI virial radius drops below 40 kpc at z\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 318\relax}}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 $\sim$ 150 globulars 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_{\rm vir}<8000$ K subhalos.666By way of illustration, removing the constraints on the virial temperature of interacting subhalos and counting all multiple encounters would result in nearly 1000 unbound close collisions.  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 globulars becomes blurred.

The median masses at infall of the lighter and heavier member of the colliding pairs are $10^{8.3}\,\,{\rm M_{\odot}}$ and $10^{9.7}\,\,{\rm M_{\odot}}$, respectively. The median relative velocity is $160\,\,{\rm km\,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_{\rm rel}\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 318\relax}}350\,\,{\rm km\,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\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 316\relax}}b\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 316\relax}}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.

5. Summary and Implications for GC Formation Models

The presence of substructure within dark matter 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 dark-matter-less globulars, 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 a 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 $\sim 10^{5}$ 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\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 318\relax}}10^{-2.5}\,Z_{\odot}, 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${}^{2}\times$cross section$\times$velocity factor – and of hierarchical assembly – globulars being brought in by the accretion of smaller satellites. With the addition of ex-situ globulars, the scaling of the total number of GCs with host virial mass is very close to linear, $N_{\rm GC}\propto M_{\rm 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 $\sim 10^{9.6}\,\,{\rm M_{\odot}}$.

(iii) The Ages of GCs: Our model differs for 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\,\,{\rm km\,s^{-1}}$, a situation that may not be conducive to the cooling and fragmention of the splash remnant (see § 3.1). The others may give origin to a population of young metal-poor globulars like Crater (Weisz et al., 2016). The dot-dashed histogram on the same panel shows the relative frequency histogram of the redshift of formation for 55 Milky Way GCs with Hubble Space Telescope photometry (VandenBerg et al., 2013), augmented by the age determinations for the young globulars 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]$\simeq-1.55$, and a high-metallicity tail at $\simeq-0.55$ (Harris et al., 2016). Only 30% of MW globulars 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_{*}/\,{\rm 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}\,\,{\rm M_{\odot}}$ (see Fig. 5).

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

[TABLE]

The zero-point was chosen to give $12+\log({\rm O/H})=8.18$ at $\log(m_{*}/\,{\rm 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 globulars, 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 globulars, 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 $\log(Z/Z_{\odot})\geq-1$: the red globulars 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 towards 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 towards even lower metallicities. We notice here two effects that have been ignored and may skew the predicted distribution in the opposite direction, towards 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 globulars (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 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 Fig. 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\sim 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 $\mu$. The two orbits have the same specific angular momentum $h$ and the subhalos collide at the semiparameter location of the ellipses, $p=h^{2}/\mu$. 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 $(\mu e_{1}/h,\mu/h)$ and $(\mu e_{2}/h,\mu/h)$, respectively, where $e_{1}$ and $e_{2}$ are the orbital eccentricities, with $e_{1}<e_{2}$. 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 $\vec{v}_{f}$ at $p$ of the combined remnant as

[TABLE]

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

[TABLE]

The new specific orbital energy is

[TABLE]

and the new eccentricity is

[TABLE]

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

[TABLE]

and the change in orbital energy is again

[TABLE]

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 globulars 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 3-D 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.).