Globular Cluster Formation from Colliding Substructure
Piero Madau, Alessandro Lupi, Juerg Diemand, Andreas Burkert, Douglas N. C. Lin
aa r X i v : . [ a s t r o - ph . GA ] J a n Preprint typeset using L A TEX style emulateapj v. 12/16/11
GLOBULAR CLUSTER FORMATION FROM COLLIDING SUBSTRUCTURE
Piero Madau , Alessandro Lupi , J¨urg Diemand , Andreas Burkert , and Douglas N. C. Lin Department of Astronomy & Astrophysics, University of California, 1156 High Street, Santa Cruz, CA 95064, USA Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy Cisco Systems (Switzerland) GmbH, Richtistrasse 7, 8304 Wallisellen, Z¨urich, Switzerland University Observatory Munich (USM), Scheinerstrasse 1, D-81679 Munich, Germany Max-Planck-Institut f¨ur Extraterrestrische Physik (MPE), Giessenbachstr. 1, D-85748 Garching, Germany
ABSTRACTWe 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 galaxyhost, a special dynamical mode of star formation that operates at high gas pressures and is intimatelytied to ΛCDM 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 whilethe warm gas within them clashes at highly supersonic speed and decouples from the collisionlesscomponent, 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 atomicline emission and then via fine-structure metal-line emission, and is subject to gravitational instabilityand 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 islinear down to dwarf galaxy scales and agrees with the trend observed over five orders of magnitudein galaxy mass; (b) a population of old globulars with a median age of 12 Gyr and an age spreadsimilar to that observed; (c) a spatial distribution that is biased relative to the overall mass profile ofthe host; and (d) a bimodal metallicity distribution with a spread similar to that observed in massivegalaxies. Subject headings: cosmology: theory — dark matter — galaxies: formation — globular clusters:general — methods: numerical INTRODUCTION
The origin of globular clusters (GCs) remains an un-solved problem in star and galaxy formation studies.With masses in the range 10 –10 M ⊙ , half-light radiiof a few pc, a bimodal metallicity distribution, inter-nal (star-to-star) variation in their light-element abun-dances, and typical ages >
10 Gyr, these remarkably com-pact stellar systems are a common feature of galaxiesin the local Universe (see, e.g., the reviews by Harris2001; West et al. 2004; Brodie & Strader 2006; Grat-ton et al. 2012; Kruijssen 2014; Forbes et al. 2018b).Young massive stellar clusters, with masses and densi-ties comparable to those of GCs, are observed in galaxymergers throughout the local Universe (e.g., Whitmoreet al. 1999), suggesting that the progenitors of metal-rich globulars could still be forming today in unusuallyhigh pressure environments. The mass of the entire GCsystem of a galaxy and its radial extent correlate tightlywith 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), atrend the former suggestive of a picture in which star for-mation in dense proto-GCs was relatively immune to thefeedback mechanisms that hamper most star formationin 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 − . ± . atGalactocentric distances ∼ > . ∼ >
90 kpc (Harris 1996; Laevens et al. 2014). The red (metal-rich) globulars are more spatially concentratedthan the blue clusters: they lie within the solar circleand form a flattened, rotating population. Blue GCsoutnumber red globulars by a ratio of 3:1. Despite theirdifferent location and kinematics, the red and blue clus-ters appear to have similar internal properties, such asmasses, sizes, and ages. The age spread among the bulkof globulars is about 2 Gyr, but there are a few outliers(e.g., Leaman et al. 2013). The youngest GC, Whiting1, is metal-rich and has an estimated age of 5 . ± . ± . . ± . z ∼ < ex-situ dur-ing the very early stages of galaxy formation, typicallyin dwarf-sized systems that were subsequently accretedinto the Galactic potential well (Searle & Zinn 1978; Pee-bles 1984; Rosenblatt et al. 1988; Cen 2001; Bromm &Clarke 2002; Diemand et al. 2005; Moore et al. 2006; Bo- Formation of GCs in Subhalo-Subhalo Collisionsley et al. 2009; Griffen et al. 2010; Trenti et al. 2015;Kimm et al. 2016; Ricotti et al. 2016; Boylan-Kolchin2017; Creasey et al. 2019, but see Chiou et al. 2019).By often relying on ad-hoc prescriptions for the forma-tion of globulars inside their own DM halos, these “pre-galactic” models have been somewhat successful in repro-ducing some puzzling properties of the GC populations.The observed spatial coincidence between GCs and mul-tiple tidal debris streams in the outer halo of M31 doesindeed suggest a direct connection between some GCsand dwarf galaxy remnants (Mackey et al. 2010). In-situ mechanisms have focused instead on the formationof GCs inside the main progenitor of their present-dayhost, perhaps from cooling-induced fragmentation of low-metallicity infalling gas (Fall & Rees 1985), during themerger 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 inearly galactic disks (Kravtsov & Gnedin 2005; Kruijssen2015; Pfeffer et al. 2018).Most clusters in the MW contain no measurableamounts of DM. The faint tidal tails observed aroundsome GCs (Grillmair et al. 1995; Odenkirchen et al. 2003)provide strong constraints on their mass-to-light ratiosand indicate that their mass distribution does not ex-tend beyond the optical radius (Moore 1996). The ab-sence of DM could be the result of stripping in the strongtidal field of the Galaxy (Mashchenko & Sills 2005; Saitohet al. 2006), and does not preclude by itself the forma-tion of GCs at the center of DM substructure.
Ex-situ scenarios in which GCs form at the center of their DMhalos do predict, however, detectable amounts of DM inthe outskirts of GCs that have always resided within theweak tidal field of the outer MW halo. This is not borneout by the data: recent observations of stellar kinemat-ics in NGC 2419, located at 90 kpc from the center ofour Galaxy, and in MGC1, located at 200 kpc from thecenter of M31, show that these GCs cannot be deeplyembedded within dark halos having a virial mass greaterthan 10 M ⊙ (Baumgardt et al. 2009; Conroy et al. 2011;Ibata et al. 2013).In this paper we propose and investigate a possiblescenario for the formation of GCs within the frameworkof ΛCDM hierarchical galaxy assembly. Many GC for-mation models leave largely unanswered the fundamen-tal question of how, at early cosmological times, self-bound and extremely dense aggregates of ∼ stars,largely mono-metallic in iron-peak elements, were ableto form and survive without retaining much of the DMof their former host. GCs are the result of a specialmode of star formation that requires extremely high gaspressures, p ∼ k B cm − K, some 3 dex higher thantypical interstellar values. It is these pressures that hin-der the dispersal of star-forming material, protect denseproto-GCs from the feedback processes that regulate starformation in the field, and produce high star formationefficiencies (e.g., Elmegreen & Efremov 1997). They arethe kind of pressure that would result, e.g., if typicalatomic clouds in the interstellar medium (ISM) with den-sities of a few atoms cm − were to collide at the Galacticorbital speeds of 200 km s − , which is precisely the out-come of the mechanism proposed here: early high-speedcollisions between subhalos infalling onto massive galaxy hosts, leading to the formation of dark-matter-less GCs . SUBHALO-SUBHALO COLLISIONS
Numerical simulations in ΛCDM have shown that mas-sive galaxy halos are populated by a rich spectrum ofsubstructures that collapsed at early times, the fossilremnants of a hierarchical merging process that is nevercomplete (Moore et al. 1999; Klypin et al. 1999; Diemandet al. 2007a, 2008; Springel et al. 2008). About one thou-sand subhalos with pre-infall masses m peak ∼ > M ⊙ arepredicted to have been accreted by our Galaxy over cos-mic time (see Han et al. 2018, and references therein). While cold gas condensation and normal, Population IIstar formation can take place at high redshift in these“atomic-cooling” subhalos because of hydrogen atomicradiative processes, their overall efficiency of convertingbaryons into stars must remain low in order to avoidoverproducing the observed abundance of dwarf galaxysatellites of the MW (Madau et al. 2008; Boylan-Kolchinet al. 2014). Mechanisms such as early reionization ofthe intergalactic medium, supernova feedback, and H -regulation have all been invoked in order to reduce thestar formation efficiencies of low-mass field halos andeven prevent the smallest ones from forming stars al-together, and “dark dwarfs” with long gas depletiontimescales 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 tophoto-evaporate low-density gas from these shallow po-tential wells, but to have little effect on relatively high-density gas in their cores that can self-shield from UVbackground radiation. The detailed mapping betwen themass spectrum of DM substructure and the population ofGalactic satellites remains, however, a matter of debate,and so is the smallest mass halo that is capable of host-ing an observable baryonic counterpart. A recent prob-abilistic comparison (abundance matching) between theluminosity function of MW’s dwarfs and N -body zoom-insimulations sets a 1 σ upper bound to the pre-infall peaksubhalo mass of Segue I – the faintest MW satellite – of m peak < . × M ⊙ (Jethwa et al. 2018). Similarly,a fresh reassessment of the “missing satellites problem”,cast in terms of number counts, suggests that luminoussatellites of the MW inhabit subhalos with infall massesas small as 10 − M ⊙ (Kim et al. 2018).The central idea of this work is to suggest that high-speed collisions between infalling substructures duringthe assembly of the main galaxy host may result in theformation of “naked” GCs. Like in the Bullet galaxygluster (Clowe et al. 2006), colliding DM subhalos andtheir stars will simply pass through one another, largelycontinuing on their original orbital path (close encoun-ters are more effective at disrupting the colliding satel-lites at low speeds). The atomic gas within them, how-ever, will collide at highly supersonic speed and decou-ple from the collisionless component. The collision will To be more precise, m peak is the maximum mass attained by asubhalo over its entire merger history. This mass is greater than thesubhalo mass at infall, as most subhalos start being stripped at ∼ < adau et al. 3compress and shock-heat the gas to characteristic tem-peratures ∼ > K, creating the very high pressure regionthat is conducive to GC formation. We will show that thelow-metallicity gas will cool via atomic line emission ona timescale that is much shorter than the crossing time.Under the appropriate conditions, shock dissipation ofthe relative kinetic energy will lead to the coalescenceof the colliding gas clouds and the formation of a Jeansunstable slab, rather than cloud disruption. A high starformation efficiency may result from the short dynamicaltimescale and high binding energy of the splash remnant.
Kinetic Theory Approach
An estimate of the number of close collisions canbe made within the framework of kinetic theory (e.g.,Makino & Hut 1997). Let us start by selecting subhalosby their pre-infall mass, i.e. by using the “unevolved”substructure mass function – the distribution of peakbound masses m peak of subhalos accreted by the host atall previous redshifts (e.g., van den Bosch et al. 2005;Jiang & van den Bosch 2016). In ΛCDM cosmologicalsimulations, this universal redshift-independent functionis well fitted by a double Schechter function of the form(Han et al. 2018) dN sub dξ = ξ − (cid:0) a ξ − α + a ξ − α (cid:1) exp( − bξ β ) , (1)where ξ = m peak /M vir is the ratio of the peak mass of asubhalo to the host halo mass, and ( dN sub /dξ ) dξ is thenumber of subhalos with masses between ξ and ξ + dξ .When adopting for the host mass the virial definitioncorresponding to the overdensity of the spherical col-lapse model, the best-fit parameters are a = 0 . , a =0 . , α = 0 . , α = 0 . , b = 8 . , and β = 1 . n sub ( ξ, x ) ∝ cx (1 + cx ) , (2)where x ≡ r/R vir , R vir is the virial radius of the host,and c is its “concentration parameter”. The subhalounevolved mass function and spatial distributions areapproximately separable (i.e., the spatial distributionsof different peak mass subhalos are similar except for achange in amplitude), and the number density of smallsubhalos with peak masses > ξM vir , n sub ( > ξ, x ), can bewritten as4 πR n sub ( > ξ, x ) = c gcx (1 + cx ) Z ξ dN sub dξ ′ dξ ′ , (3)where g ≡ / [ln(1 + c ) − c/ (1 + c )].Let us further make the simplifying assumption thatthe velocity distribution of subhalos is a local Maxwellian In reality, the dynamics of subhalos differs somewhat fromthat of DM particles because of dynamical friction, but this is asmall effect for the majority of subhalos that are located in theouter regions of the main host and at the low-mass end (Han et al.2016). with one-dimensional velocity dispersion σ v . The distri-bution of relative encounter speeds is then f ( v rel ) = 12 √ πσ v v exp (cid:18) − v σ v (cid:19) , (4)and the mean encounter speed is ¯ v rel = 4 σ v / √ π . Forisotropic orbits in an NFW potential, the velocity dis-persion obtained by solving the Jeans equation is ( Lokas& Mamon 2001) σ v ( x ) = (cid:18) GM vir R vir (cid:19) g (1 + cx ) x × Z ∞ x (cid:20) ln(1 + cs ) s (1 + cs ) − cs (1 + cs ) (cid:21) ds. (5)Under the assumption that collisions occur on randomorbits, we can finally estimate the mean number of en-counters between subhalo pairs of peak masses (in unitsof M vir ) > ξ and > ξ in a time interval ∆ t as N coll =∆ t πR Z ¯ v rel ( x ) x dx × Z ξ Z ξ n sub ( ξ, x ) n sub ( ξ ′ , x )( πb ) dξdξ ′ , (6)where b max is the maximum impact parameter for a closecollision.Consider now, for illustrative purposes, a MW-sizedhost halo of mass M vir = 10 M ⊙ and concentration c = 10 (e.g., Duffy et al. 2008). Its isotropic velocitydispersion profile in Equation (5) reaches a maximumvalue of 102 km s − at x = 0 .
08, corresponding to a meancollision speed of 230 km s − . These relative orbital ve-locities are much greater than the internal velocities ofsubhalos, and few such encounters will actually resultin mergers. We can therefore neglect gravitational fo-cusing and, in order to ensure that the region of over-lap between the gaseous cores of the interacting clumpsis extensive, define penetrating collisions as those with b max = 0 . r vir , + r vir , ), where r vir , and r vir , are thepre-infall virial radii of the two clumps. The mean num-ber of impacts within R vir between subhalos having peakmasses m peak > M ⊙ is then N coll ≃ (cid:18) ∆ t Gyr (cid:19) . (7)This is clearly an interesting figure – comparable after∆ t ∼ a few Gyr to the number of GCs observed today inthe halo of the MW – if many such encounters were toresult in the formation of one or more globulars. Becauseof the steep substructure mass function, most collisionsinvolve subhalos at the small-mass end of the distribu-tion: the mean frequency of close encounters betweenatomic-cooling subhalos with m peak > M ⊙ and satel-lites that are ten times more massive is 3.6 times smallerthan Equation (7). There are only N coll ∼ t/ Gyr)close collisions between massive satellites with m peak > M ⊙ , but these may give origin to multiple GCs perencounter. The same calculation predicts about 1 colli-sion Gyr − between atomic-cooling subhalos in a dwarfgalaxy-sized host with M vir = 10 M ⊙ . Formation of GCs in Subhalo-Subhalo Collisions Fig. 1.—
The GC system number-host halo mass relation. Thedata points are observational estimates from a recent compilationby Burkert & Forbes (2019). The lines show the close-to-linear rela-tion predicted by a kinetic theory model of cluster formation follow-ing subhalo-subhalo collisions (the normalization is arbitrary, seetext for details). Solid curve: Mean number of GCs formed eitherin the main progenitor of a halo of mass M vir ( in-situ globulars)or in a companion system and later accreted ( ex-situ globulars).Dashed curve: Ex-situ globulars only. Dot-dashed curve:
In-situ globulars only.
In reality, of course, many of the simplified assump-tions behind our integration of Equation (6) do not holdin practice: 1) subhalos are often accreted at similartimes and locations as members of groups along fila-ments, and this causes an enhancement of encounters atsmall angles (Benson 2005; Zentner et al. 2005; Gonz´alez& Padilla 2016); 2) after accretion, subhalos are suscep-tible to dynamical friction and tidal stripping, and theirmass and spatial distributions evolve away from theirform at infall (Diemand et al. 2007b; Springel et al. 2008);3) the abundance of subhalos above a given peak massevolves with time, and so does the background gravita-tional potential in which they move; 4) subhalos are pre-dicted to have a non-Maxwellian orbital velocity func-tion, with centrally rising velocity anisotropy (Sawalaet al. 2017); and 5) ram-pressure stripping will removethe outer gaseous component of subhalos on orbits withsmall pericenter radii, a process that is most effectiveclose to pericenter passage (Mayer et al. 2006; Grcevich& Putman 2009). In § Via Lactea N -body simulation to track collisions be-tween infalling substructure in a massive MW-sized hostthat grows and evolves in a fully cosmological context. Scaling with Host Mass and Ex-situ Globulars
One of the most intriguing aspects of GC phenomenol-ogy is their relationship to DM halos: the total mass ofthe GC population of a galaxy and the number of GCscorrelate linearly with the DM mass of the host halo.This trend is valid over five orders of magnitude in galaxymass (see, e.g., Harris et al. 2015; Forbes et al. 2016;Burkert & Forbes 2019, and references therein), and con-trasts markedly with the non-linear relation between to- tal stellar mass and DM-halo mass (e.g., Behroozi et al.2013). The observed correlation may point to a funda-mental GC system-DM connection rooted in the clus-ter formation physics, or simply be the inevitable conse-quence of hierarchical assembly and the central limit the-orem (Boylan-Kolchin 2017; El-Badry et al. 2019; Burk-ert & Forbes 2019).It is of obvious interest at this stage to use kinetictheory and derive a scaling relation between the num-ber of collisions and the virial mass of the host followingEquation (6). The collision frequency between subhalosmore massive than a given m peak ≪ M vir at fixed con-centration parameter is proportional to n ∝ R − M . (from Eqs. 1 and 3) times ¯ v rel ∝ ( M vir /R vir ) / (see Eq.5). Including a dependence on halo concentration andassuming that the collisions of interest continue for atimescale ∆ t that is independent of host mass, one infersfrom Equation (6): N coll ∝ R − M . ( M vir /R vir ) / c . (8)Noting that R vir ∝ M / and using the haloconcentration-mass relation of Duffy et al. (2008), c ∝ M − . , we obtain the following scaling N coll ∝ M . .Under the assumption that GCs populate present-dayDM halos in direct proportion to N coll , let us write N GCin − s = A ( M vir / M ⊙ ) . , (9)where A is an arbitrary normalization factor. In thisscenario, N GCin − s tracks the population of GCs formed in-situ within a massive galaxy host. The same substructurecollision mechanism, at work in satellites prior to infall,would also produce a significant population of globularsthat formed ex-situ and were subsequently accreted. The N GCin − s − M vir relation given above and the substructuremass function in Equation (1) allow a straightforwardcalculation of the mean number of accreted GCs as N GCex − s = N GCin − s Z ξ min ξ . dN sub dξ dξ, (10)where ξ min ≡ m min /M vir is the critical subhalo mass be-low which GCs cannot form. Obviously, collisions be-tween m peak > M ⊙ clumps as a mechanism for clus-ter formation require at least a few (say ∼
4) atomic-cooling subhalos in a given host; this trivial conditionyields, after integrating Equation (1), m min ≃ . M ⊙ .This critical mass could be larger – thereby decreasingthe fraction of accreted globulars – if the relative orbitalvelocity v rel required to trigger the formation of a GCby impact had to exceed ∼
50 km s − . Note, however,that lower velocity encounters, as those expected in a lessmassive host, may still produce the high pressures thatare conducive to GC formation if the colliding gaseouscores were typically denser, e.g. at very high redshifts.Assuming here that all halos below m min = 10 . M ⊙ do not contain a GC, Equation (10) gives N GCex − s /N GCin − s =0 . , . , .
17 for M vir = 10 , , M ⊙ , i.e. mostGC formation occurs in-situ for the lowest-mass host ha-los and most GCs in massive elliptical systems are actu-ally accreted. With the addition of the accreted glob-ulars, the scaling of the total number of GCs, N GC = Similar trends are seen in recent phenomenological models adau et al. 5 N GCex − s + N GCin − s , with host virial mass becomes very closeto linear, N GC ∝ M . . (11) Therefore, in a collision-driven scenario, a constant GCnumber-to-halo mass ratio is the result of encounter prob-ability calculations (Eq. 8) and hierarchical clustering(Eq. 10): GCs are the result of a distinctive mode ofstar formation that tracks the total mass of their hostgalaxy rather than its stellar mass.
Figure 1 shows the observed GC number-halo mass re-lation 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 assumednormalization A = 95 in Equation (9). The model ap-pears to reproduce the observed trend reasonably well,with the expected average number of GC systems drop-ping below a few in dwarf galaxies with M vir ∼ < M ⊙ .With the adopted parameters, a MW-mass system with M vir = 10 M ⊙ today would host a grand total of 165massive GCs, of which 70 formed ex-situ . Estimates ofthe number of accreted GCs within the MW today rangefrom 30 to 90 (Forbes & Bridges 2010; Leaman et al.2013), so it is conceivable that a substantial fraction ofglobulars in massive hosts may have an external origin.The significant scatter in the observed relation on dwarfgalaxy scales may reflect variations in merger histories,uncertainties in dark halo mass determinations based onkinematical tracers, or environmental effects such as tidalstripping (Burkert & Forbes 2019).While it is reassuring that a scenario in which GCs arethe result of collisions between subhalos may be able toaccomodate a uniform GC production rate per unit hosthalo mass as implied by the data, it may be prematureto 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 onthe age and age-spread of GCs, on the epoch at which theGC-to-halo mass relation may actually be established, onthe normalization of such relation, and on the baryonicphysics leading to the formation of massive globulars.In the following sections, we shall discuss the conditionsunder which high-speed impacts may lead to cloud co-alescence, gravitational instability, and the formation of“naked”, bound GCs, and use numerical simulations toargue that a collision-based framework may fulfill severalkey observational constraints on GCs that have emergedover the last two decades. GC FORMATION IN HIGH-SPEED IMPACTS
Observations of interacting galaxies like the Antennaeand Mice pairs show widespread stellar cluster formationgenerated by the collision, and high-resolution simula-tions of encounters between disk galaxies with realisticISM reveal the presence of shock-induced gas compres-sion and star formation at the collision interface (e.g.,Saitoh et al. 2009). There is a vast literature on super-sonic cloud-cloud collisions as a possible triggering mech-anism 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 of GC formation based on DM merger trees (e.g. Boylan-Kolchin2017). et al. 2014; Balfour et al. 2017), and cloud-cloud colli-sions in the young Galaxy have already been invoked inthe 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 col-lision occurs between subhalos, the post-shock gas rem-nant will cool down on a timescale that is much shorterthan the compression time, and may, under the rightconditions, become gravitationally unstable.We shall focus on collisions between subhalos whereboth members of the pair are above the atomic-coolingthreshold at infall, i.e. have a mass corresponding to avirial temperature of T vir > m peak > . × M ⊙ [ µ (1 + z ) / − / B − / . (12)Here, µ is the mean molecular weight per particle ( µ =1 .
23 for neutral primordial gas),
B ≡ ∆ c / (Ω zm π ),∆ c = 18 π + 82 y − y is the redshift-dependent densitycontrast at virialization (Bryan & Norman 1998), y ≡ Ω zm −
1, Ω zm = Ω m (1+ z ) / [Ω m (1+ z ) +Ω Λ ], and z is theinfall redshift. In evaluating these expressions, we haveassumed a Planck Collaboration et al. (2018) flat cosmol-ogy with parameters (Ω m , Ω Λ , h ) = (0 . , . , . B = 1. Atomic cooling subhalos are likely to be pol-luted by metals through inefficient star formation, whichis key since GCs will inherit the gas-phase metallicityof the colliding subhalo population. According to cos-mological radiation hydrodynamics simulations by Wiseet al. (2014) that follow the buildup of dwarf galaxiesfrom their early Population III progenitors (see also Ri-cotti et al. 2016), most atomic-cooling halos host metal-enriched stars by redshift 6, when infalling substructuresbegin colliding at high speed. The gas shock heated ina close encounter will then be pre-enriched by previousgeneration of massive star formation.Gas distributed isothermally in atomic-cooling halosdevelops warm central cores of constant baryon density n c ≈ . z ) / cm − ≡ n cm − within radii r c ∼ . 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, r vir ≃ . (cid:18) m peak . M ⊙ (cid:19) / (cid:18)
51 + z (cid:19) (13)is the virial radius of the system. Prior to infall, the sus-ceptibility of such systems to internal and external feed-back mechanims may decrease their gas content to 5-10%of their virial mass (Wise et al. 2014). After infall, theircentrally concentrated gaseous cores are expected to sur-vive disruption by ram pressure before the first pericenterpassage (Mayer et al. 2006; Visbal et al. 2014b). Shock Heating, Cooling, and GravitationalInstability
To highlight the dominant physical processes at workduring a high-speed close interaction, let us consider theidealized situation of a head-on collision between twoidentical gaseous cores of characteristic hydrogen numberdensity n H , = 0 . n c , mass density ρ = n c m p , temper-ature T = 8000 K, and radius r c = 0 . r vir . These arethe density and thermal pressure conditions of the warmdiffuse component of the interstellar medium. In the ab-sence of molecular cooling, even when gas is heated by Formation of GCs in Subhalo-Subhalo Collisionsshocks and compression above the 8000 K limit, atomicradiative losses together with photoelectric heating fromdust grains, cosmic ray heating, and soft X-ray heatingwill keep the gas temperature near this value (Wolfireet al. 1995). The total gas mass involved in the impact,2 m c = 8 π ρ r c ≃ M ⊙ n (cid:18) r c . r vir (cid:19) × m . (cid:18)
51 + z (cid:19) , (14)where m . ≡ m peak / . M ⊙ , is about 10 per cent ofthe mass corresponding to the universal baryon fraction.Following a highly supersonic interpenetrating impact,most of the gas decouples from the collisionless compo-nent. Two shock waves arise and propagate from the con-tact discontinuity with velocity (in the adiabatic phase) V s ≃ v rel / T = 3 µm p k B V s ≃ . × K V , (15)where V ≡ V s /
100 km s − and we have set the molec-ular weight to µ = 0 .
59. The density is enhanced bya factor 4, the post-shock pressure is p = (3 / ρ V s ,and the gas cools via hydrogen, helium, and metal lineemission on a (isobaric) timescale t cool = 3 k B T n tot n , Λ ≃ .
03 Myr (cid:18) V . n (cid:19) . (16)Here, n tot is the total particle number density( n tot /n H , = 9 / n H , = 4 n H , is the post-shock meanhydrogen density, Λ is the cooling function (e.g., Wanget al. 2014), and we have assumed collisional ionizationequilibrium and ignored gas clumping (see Fig. 2). Since t cool is typically much shorter than the crossing time, t cross ≡ r c V s ≃ . V − r . ) , (17)where r . ≡ r c / . T = T = T ), characterized bya sound speed c s = p k B T /µm p and a Mach number M ≡ V s c s ≃ . T − / . V ) , (18)where T . ≡ T / µ = 1 . n H , = M n H , ≃
140 cm − ( n T − . V ) , (19)and the shock velocity relative to the unpertubedmedium is now V s ≃ v rel /
2. The collision is also ap-proximately one-dimensional, since the lateral rarefac-tion timescale is ∼ M t cross . For a strictly isothermalcollision, when the compression phase is completed aftera time t cross , a rarefaction wave will propagate throughthe shock-compressed gas slab at the local sound speed,and the medium will rapidly expand into the surround-ings on a timescale ∼ t cross / M . If the collision is to re-sult in gravitational instability and collapse, there must Fig. 2.—
Radiative cooling timescale t cool times the pre-shockhydrogen density n H , , as a function of shock speed for atomic gaswith metallicity Z = 0 . Z ⊙ (top curve) and Z = 0 . Z ⊙ (bottomcurve). The cooling function and electron fraction in collisionalionization equilibrium have been taken from Wang et al. (2014).Strong shock jump conditions have been assumed. The power-lawapproximation used in Equation (16) for low metallicity gas, isplotted with the dashed line in the range 40 < V s <
300 km s − . be density perturbations that can grow on a timescale ∼ < t cross . We shall see below, however, that gas furtherdownstream will actually cool below its pre-shock tem-perature via fine-structure metal lines, i.e. the collisionis “more dissipative than isothermal” (e.g. Whitworth &Clarke 1997).Ly α cooling becomes inefficient below 8000 K, and thefast shock will dissociate most pre-existing molecules. Inthe absence of an ambient UV radiation field, enoughmolecular hydrogen may reform behind the shock via theH − channel or on the surface of dust grains to rapidlycool the compressed gas to 100 K. It seems likely, how-ever, that photodissociation by local stellar sources ofLyman-Werner photons – either within the colliding sub-halos or in in the main host – will act to suppress H for-mation in the post-shock gas, and the temperature willthen plateau around 8000 K (e.g., Kang et al. 1990). TheBonnor-Ebert mass in such warm, dense, and pressurizedmedium is M BE =1 . c s p G M ρ ≃ . × M ⊙ ( n − / T . V − ) . (20)a factor of 1 / M smaller than its value in the pre-shockgas (where M BE > m c , i.e. the pre-collision clumps aregravitionally stable), and comparable to the observedmass of GCs.GCs are generally metal poor. While the lowest metal-licity globular listed in the 2010 edition of the Harris(1996) catalog has [Fe/H]= − . − . Fig. 3.—
Post-shock temperature vs. time (in Myr, mea-sured since the fluid element was shocked) for a case with V s =100 km s − and n H , = 0 .
75 cm − . Strong shock jump conditionsand isobaric cooling have been assumed, and the non-equilibriumchemistry network has been taken from the krome package (Grassiet al. 2014). The gas is also assumed to be purely atomic, dust-free, and optically thin. The thermal model includes atomic cool-ing, metal cooling following Bovino et al. (2016), and pdV work.Four curves are shown for different metallicities: Z = 0 (solid line), Z = 10 − . Z ⊙ (dot-dashed line), Z = 10 − Z ⊙ (long-dashed line),and Z = 10 − Z ⊙ (short-dashed line). of C II (158 µ m) and O I (63 µ m). The rate coeffi-cient for the collisional excitation of O I by H atoms is k OI ≃ . × − cm s − T . . e − T OI /T where T OI = 228K is the energy of the excited O I level over k B (Draine2011). At densities far below the critical density of theline, this process removes energy at the rate n , Λ O = n , A O k OI ( k B T OI ) , (21)where A O is the abundance of O relative to hydrogen, A O = 4 . × − Z − . (Anders & Grevesse 1989), Z − . is the gas metallicity in units of 10 − . solar, and alloxygen is assumed to be neutral and in the gas phase.The isobaric cooling timescale, t Ocool = 5 k B T n tot n , Λ O ≃ . (cid:18) T . . V n Z − . (cid:19) , (22)is still shorter than the crossing time, but it is consider-ably longer than t cool in Equation (16) as a result of themismatch in timescales between fine-structure metal-linecooling below 8000K and atomic cooling at higher tem-peratures. After a time ∼ t Ocool , the flow will cool belowits pre-shock temperature and contract further, generat-ing a growing, cold and very dense layer at the center ofthe shock-bounded slab.The arguments above are based on simple timescaleestimates, and to check the robustness of our results wehave used the non-equilibrium chemistry package krome (Grassi et al. 2014) to construct more sophisticated mod-els of the cooling and thermal properties of the post-shock gas. The thermal model includes non-equilibriumcooling for H and He at all temperatures and for met-als below 10 K Bovino et al. (2016), and compressionalheating. Strong shock jump conditions and isobaric cool- ing were adopted for gas at different metallicities (rang-ing from 0 to 0.1 solar), and the cooling function wascomputed assuming dust-free and optically thin condi-tions, and neglecting molecular-phase processes. Figure3 shows the time evolution of the post-shock temperaturefor a case with V s = 100 km s − and n H , = 0 .
75 cm − .The figure is consistent with the qualitative estimatesgiven above: low-density gas, shock-heated to high tem-peratures, cools fast via atomic hydrogen and heliumline emission to about 6500 K, and then levels off for atimescale that is longer at lower metallicities. Eventually,metal cooling via fine-structure lines of C and O takesover. The end result is a runaway cooling phenomenonthat drives the shocked medium to temperatures T < T and hydrogen densities n H , = n H , T /T (assuming iso-baric cooling at constant mean molecular weight) thatare, respectively, well below and above those of the col-liding gas clouds. Note that these calculations do notexplicitly include internal (within the colliding substruc-tures) and external (within the host galaxy) sources ofUV and X-ray photoheating.The one-dimensional compression will create a coldlayer of thickness-to-diameter ratio ∼ n H , /n H , < / M ∼ .
01 ( T . V − ) that depends on the entity ofradiative losses. In this case, gravitational accelerationsare most important in the dynamics of flows transverseto the collision axis, and it is more appropriate to dis-cuss the criterion for the gravitational instability of athin, 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 ofconstant surface density Σ and isothermal sound speed c s , the dispersion relation for modes transverse to thecollision axis gives the wavelength and timescale of thefastest growing mode as λ G = 2 c s G Σ , (23)and t G = c s πG Σ . (24)Perturbations smaller than λ G are gravitationally sta-ble, while those with longer wavelengths are unstablebut grow more slowly. The surface density of the coldlayer between the two shocks increases with time asΣ( t ) = 2 ρ V s t , where t is measured from the onset ofmetal cooling below 8000 K. The shocked gas is first li-able to gravitational fragmentation at time t = t G , whenΣ G ≡ Σ( t G ) = 2 ρ V s t G (Kang et al. 1990). Using Equa-tion (24) one derives then t G = (cid:18) c s πGρ V s (cid:19) / ≃ . (cid:18) C n V (cid:19) / , (25)andΣ G = (cid:18) ρ V s c s πG (cid:19) / ≃
19 M ⊙ pc − ( n V C ) / , (26)where C ≡ c s / (1 km s − ). Non-linear fragmentation oc-curs before the start of the rarefaction era when t G ≤ t cross − t Ocool ≃ t cross . This inequality is fulfilled when c s ≤ π Gρ r c V s ≃ . − ( V − n r . ) (27) Formation of GCs in Subhalo-Subhalo Collisions(corresponding with the choosen scalings to T ≤ λ G πc s t G ≃
107 pc (cid:18) r . n V (cid:19) (28)and masses M G = π Σ G (cid:18) λ G (cid:19) ≃ . M ⊙ (cid:18) r . n V (cid:19) . (29)While the latter is again comparable to the observedcharacteristic mass of GCs after accounting for mass lossby stellar evolution and tidal disruption after birth, itis also somewhat ill-defined – scaling strongly with cloudsize, gas density, and shock velocity. Neverthless, it is en-couraging that collisions between atomic cooling, metal-poor subhalos may offer a mechanism for imprinting thesignature of the GC mass scale on the collapsing shell.It should be noted that the conditions for gravitationalinstability given in Equations (27-29) can only be met aslong as the gas metallicity Z is not much below 10 − Z ⊙ (see Fig. 3). This may provide a plausible explanationfor the threshold metallicity, [Fe/H]= − . d along the lineconnecting the cloud centers, d = 4 r c T M T ≃ r . n V − ) , (30)which is comparable to the typical half-light radius ofGCs. The free-fall time for a uniform, pressure-freesphere of such gas, t ff = r π Gρ = 1 M s πT Gρ T ≃ . V − n / r . ) , (31)is shorter than the several-Myr evolutionary timescalefor massive stars.The above analysis is highly idealized, and is onlymeant to provide an idea of the general conditions underwhich the splash remnant may become Jeans unstable. Itis clear from the previous discussion that one expects sig-nificant spatial structures in the post-shock region bothin temperature and density as a function of distance fromthe shock fronts (e.g., Smith 1980; Kang et al. 1990; Ku-mai et al. 1993), and the thin shell approximation maybe inadequate in the case of very inhomogeneous flows(Yamada & Nishi 1998). A shocked slab is known tobe susceptible to a number of hydrodynamical instabil-ities like the non linear thin shell instability (Vishniac1994), which may compete with the gravitational insta-bility and produce substructures on the scale of the slabthickness. The thermodynamic treatment of the prob-lem has been simplified by neglecting molecular cooling, Note that some scenarios for the formation of multiple stellarpopulations within GCs require globulars to have been initiallymuch more massive than they are today, (e.g. D’Ercole et al. 2008). heating from external X-ray radiation and cosmic rays,dust processes, and gas clumping. The simple modelof head-on collisions between identical, uniform-densityclumps is obviously unrealistic, and should be extendedto off-center impacts between subhalos with internal sub-structure.Nevertheless, our calculations may elucidate the condi-tions for a special, dynamical mode of star formation fol-lowing substructure collisions, a mode that is intimatelytied to ΛCDM hierarchical galaxy assembly. It seemsplausible that GC formation by impact requires the rela-tive velocities of the colliding subhalos to be in a specificrange. If the collision velocity is too low, shocks may notbe able to produce the extremely high pressure environ-ments that are a prerequisite to the formation of denseand tightly bound clusters. Conversely, if the velocity istoo high and the shock too violent, the interacting cloudswill expand and disperse before significant radiative cool-ing can occur. Equations (16) and (17) give t cool t cross ≃ . × − (cid:18) V . n r . (cid:19) , (32)and t cool /t cross < V s <
343 km s − ( n r . ) . .This upper velocity limit is only weakly dependent on theproperties (gas density and size) of the original collidingclumps.Simulations of interstellar cloud-cloud collisions mayalso 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 atmoderately supersonic speeds, star formation operatesin a global hub-and-spoke mode that produces a cen-tral monolithic stellar cluster. At higher collision veloc-ities, a spider’s-web mode operates and delivers a loosedistribution of independent, small sub-clusters instead.When the clouds have pre-collision substructure, how-ever, the collision velocity becomes less critical (Balfouret al. 2017). In these numerical experiments the gas isevolved with a barotropic equation of state. Ultimately,the detailed fate of subhalo-subhalo high-speed collisionsshould be addressed with the help of cosmological hydro-dynamic simulations that include all the relevant heatingand cooling processes. N-BODY SIMULATION
High-speed close interactions between satellites orbit-ing within a parent halo have been advocated as a majormechanism for the morphological evolution of galaxiesin clusters (Moore et al. 1996; Baushev 2018), and col-lisions between protogalaxies have been proposed as anew pathway to form supermassive black holes at veryhigh redshifts (Inayoshi et al. 2015). Cosmological N -body simulations by Tormen et al. (1998) showed thatfast satellite-satellite encounters with impact parameter b < ( r vir , + r vir , ) are fairly common and can lead to sig-nificant mass loss and disruption. In this section we makeuse of the Via Lactea I (VLI) N -body simulation to ar-gue that rarer, nearly-central collisions between atomic-cooling subhalos still occur frequently enough at highredshifts to represent a plausible pathway to the forma-tion of GCs.VLI is a dark matter-only simulation that follows theformation of a MW-sized halo in a ΛCDM cosmology.adau et al. 9 Fig. 4.—
Frequency distribution of redshifts (top left panel), peak masses (top right panel), relative velocities (middle left panel), impactparameters (middle right panel), impact angles (bottom left panel), and first pericenter distances (bottom right panel) for all unboundinterpenetrating collisions between atomic cooling subhalos in the
Via Lactea simulation (see text for details). The top right panel shows thepeak mass frequency histograms for both the lighter (dashed line) and heavier (solid line) members of all colliding pairs, while the bottomright panel depics the first pericenter distance for all colliding subhalos that reach pericenter before disruption. The dashed histogram inthe left top panel delineates the redshifts of formation for a sample of Milky Way GCs. The age determinations by VandenBerg et al.(2013), Sarajedini et al. (2007), Dotter et al. (2010), Valcheva et al. (2015), and Weisz et al. (2016) have been converted to redshift usinga Planck Collaboration et al. (2018) cosmology. The distribution is poorly known at z ∼ > . × M ⊙ and evolved with a forceresolution of 90 pc starting from redshift 50 (Diemandet al. 2007a,b). We stored and analyzed 200 outputsfrom redshift 16 to z = 0. The simulation has sufficientmass resolution to follow metal and atomic-cooling sub-halos through many orbits and severe mass stripping,and sufficient output time resolution (∆ t = 68 . VLI initial conditions were generated with the originalversion of the GRAFIC2 package (Bertschinger 2001),which incorrectly used the baryonic instead of the darkmatter power spectrum for the refinement levels, leadingto reduced small-scale power. Compared to
Via LacteaII (Diemand et al. 2008), the abundance of subhalos inVLI is suppressed by a factor of 1.7, and one shouldtherefore view the derived collision frequency strictly as alower limit. This is even more so as even state-of-the-artcosmological simulations still suffer from significant over-merging, and many subhalos will be artificially disruptedbefore a collision by numerical resolution effects (van denBosch et al. 2018). The distance of closest approach be-tween subhalos was found by linearly interpolating dis-tances amid snapshots. Once a subhalo is accreted by themain host, its diffuse outer layers are rapidly stripped offby tidal forces, with tidal mass losses being more signifi-cant for more massive subhalos (Diemand et al. 2007b).The accurate tracking of the accretion history of sub-structure allows us to define the epoch when the satellitemass reached a maximum value, denoted here as m peak ,before infall.We have identified all collisions between atomic cool-ing systems ( T vir > v rel > V max , andwhose centers are separated by a distance b < b max =0 . r vir , + r vir , ). Here V max , is the maximum circu-lar velocity of the larger (“1”) of the two subhalos atimpact, 3 V max is the escape velocity from the centerof a spherically-averaged NFW density profile, and thechoosen maximum impact parameter b max correspondsto the sum of the gas core radii of the pair. The impactparameter condition ensures that the region of overlapbetween clumps colliding off-center is extensive and sois the resulting gas splash, and excludes events wheresmall, dense clumps may plow through the tenuous haloof larger subhalos without much resistance, i.e. withouttheir gas becoming dislodged. In order to minimize the The 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, whosetotal extent stretches well beyond these cores and is found by6DFOF in a second step (Diemand et al. 2006). effect of ram-pressure stripping on the gaseous contentof interacting substructure, the collisions of interest hereare those that occur before the first pericenter passage ofeach subhalo. Clumps with correlated infall histories of-ten undergo multiple collisions: for simplicity, we assumehere that the smaller member of a colliding pair depletesits gas reservoir without replenishment during the firstinterpenetrating encounter, and do not count any closeinteractions it may be involved in afterwards.Figure 4 shows the frequency distributions of redshifts,masses, relative velocities, impact parameters, angles ofimpact, and first pericenter distances for all encountersthat satisfy the above criteria. We tally 133 collisionevents within today’s VLI virial volume, correspondingto R vir = 288 kpc (note that VLI virial radius drops be-low 40 kpc at z ∼ > ∼
150 globulars observed today in the halo of the MW.If, on the other hand, many more low-mass GCs formedinitially than is currently observed, and were selectivedestroyed by dynamical processes acting over a Hubbletime (e.g., Gnedin & Ostriker 1997; Vesperini 1998), thenour model would require either the formation of multipleGCs per impact , or a higher collision frequency, perhapsinvolving the more numerous population of T vir < Collisions occur at early times in the transi-tion region between the main host and the field: we countonly 10 events within a Galactocentric distance of 50 kpc.For the most part, these collisions involve subhalos thatare already dynamically associated before accretion intothe main host, i.e. they are either part of the same in-falling halo or two separate clumps descending along afilament and organized into small groups with correlatedtrajectories (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 heaviermember of the colliding pairs are 10 . M ⊙ and 10 . M ⊙ ,respectively. The median relative velocity is 160 km s − ,and there is a weak negative correlation (with corre-lation coefficient r = − .
3) between relative velocityand redshift: the few extreme-velocity impacts with v rel ∼ >
350 km s − all occur at z < .
5, when the depthof the gravitational potential of the main host is larger.Collisions have typical impact parameters in the range0 . ∼ < b ∼ < . . SUMMARY AND IMPLICATIONS FOR GC FORMATIONMODELS
The presence of substructure within dark matter ha-los is a unique signature of a Universe where systems By way of illustration, removing the constraints on the virialtemperature of interacting subhalos and counting all multiple en-counters would result in nearly 1000 unbound close collisions. adau et al. 11grow hierarchically through the accretion of smaller-massunits. We have investigated a scenario where the forma-tion of GCs is triggered by high-speed collisions betweeninfalling, atomic-cooling subhalos during the early assem-bly of the main galaxy host. This is a special, dynam-ical mode of star formation that operates at extremelyhigh gas pressures, is relatively immune from the feed-back processes that regulate star formation in the field,and tracks the total mass of the main host rather than itsstellar mass. The proposed mechanism would give ori-gin to dark-matter-less globulars, as colliding DM subha-los and their stars will simply pass through one anotherwhile the warm gas within them shocks and decouples.The well-known Bullet galaxy cluster provides a strik-ing illustration of the process, albeit on different scales.In a MW-sized host, where the relative orbital veloci-ties between subhalos substantially exceed their internalvelocities, most close encounters are unbound collisionsbetween satellites that are accreted at similar early timesand locations, and occur in the outer regions of the mainhost.Below, we summarize our main results and discusssome implications for GC formation models. (i) Fragmentation of Shocked Gas and the Masses ofGCs:
Shock heating and cooling, the encounter geom-etry, and the complexities of multiphase gaseous innerhalos are all key factors in determining the outcome ofa subhalo-subhalo impact. We have shown that, underidealized conditions, the low metallicity warm gas in thecores of interacting subhalos will be shock-heated to char-acteristic temperatures ∼ K, and will cool rapidlyfirst via atomic line emission and, further downstream,via fine-structure lines of C II and O I. Because the col-lision is more dissipative than isothermal, the resultingshock-compressed slab is liable to gravitational instabil-ity. An idealized analysis in the thin shell approximationyields the conditions under which the imprint of the GCmass scale could be present on the cooling, collapsingshell. The requirements for gravitational fragmentationcan only be satisfied for gas metallicities Z ∼ > − . Z ⊙ ,thus offering a natural interpretation for the observedminimum metallicity of GCs. We caution, however, thatthese findings are based on a simplified thermodynamictreatment of the problem that neglects molecular cool-ing, heating from external X-ray radiation and cosmicrays, dust processes, and gas clumping. (ii) The GC system-DM Connection: An analysisof the scaling behavior of the encounter frequencywithin the kinetic theory approximation points to a GCnumber-halo mass relation that is the result of bothencounter probability calculations – the subhalo colli-sion rate per unit volume being the usual density × crosssection × velocity factor – and of hierarchical assembly –globulars being brought in by the accretion of smallersatellites. With the addition of ex-situ globulars, thescaling of the total number of GCs with host virial mass isvery close to linear, N GC ∝ M . , in agreement with thetrend observed over five orders of magnitude in galaxymass. This uniform GC production rate per unit hosthalo mass is predicted to break down on dwarf galaxyscales, perhaps below a critical mass of ∼ . M ⊙ . (iii) The Ages of GCs: Our model differs for much pre-vious work as it does not assume an arbitrary value forthe redshift when metal-poor GC formation is shut-off.The details of the redshift distribution in top left panelof Figure 4 reflect the mass assembly history of the sim-ulated MW-sized host system, but it is again noteworthythat a scenario in which GCs are the result of collidingsubstructures would produce a population of old clusterswith typical ages >
10 Gyr, a median age of 12 Gyr (cor-responding in the adopted cosmology to a median red-shift of 3.5), and an age spread that is similar to the oneobserved. In contrast to many pregalactic scenarios (e.g.,Katz & Ricotti 2013; Kimm et al. 2016; Boylan-Kolchin2017), in our model GCs have extended formation his-tories and typically form after the epoch of reionization:only about 38% of all close encounters occur at redshiftsgreater than 4. Seven collision events in our sample takeplace at z <
2. Four of these “late” impacts have relativevelocities in excess of 350 km s − , a situation that maynot be conducive to the cooling and fragmention of thesplash remnant (see § HubbleSpace Telescope photometry (VandenBerg et al. 2013),augmented by the age determinations for the young glob-ulars Crater, Pal 1, Terzan 7, and Whiting 1 (Sarajediniet al. 2007; Dotter et al. 2010; Valcheva et al. 2015; Weiszet al. 2016). The distribution is poorly known at z > (iv) The Metallicity of GCs:
The GC population inthe MW is observed to be clearly bimodal, with a low-metallicity component peaking at [Fe/H] ≃ − .
55, anda high-metallicity tail at ≃ − .
55 (Harris et al. 2016).Only 30% of MW globulars have [Fe/H] > −
1, but manymassive galaxies possess strongly bimodal GC systems,with nearly equal numbers of metal-rich and metal-poorclusters (Peng et al. 2006). Since stars within most GCsdo not show an internal spread in iron-peak elements,there must exist a mechanism that chemically homoge-nizes the gas within a protocluster before the onset ofstar formation. The dominant mode of chemical mix-ing is thought to be turbulent diffusion (Murray & Lin1990), which has been shown to produce a stellar abun-dance 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 inter-acting subhalo pair that triggers their formation, and auseful perspective can be obtained by assigning a stel-lar mass to each subhalo at infall following the medianstellar-to-halo mass relation from Behroozi et al. (2013).This redshift-dependent prescription, extrapolated to thesmall scales of interest here, leads to a broad range ofstellar masses, 10 . < m ∗ / M ⊙ < . for our sample ofcolliding substructures, a distribution that extends overnearly 6 decades in mass with a median value equal to10 . M ⊙ (see Fig. 5).A general tendency of decreasing metallicity towardslower stellar masses is commonly accepted, but the exactform of the stellar mass ( m ∗ ) vs. gas-phase metallicity( Z ) relation (hereafter MZR) and its evolution with red-2 Formation of GCs in Subhalo-Subhalo Collisions Fig. 5.—
Stellar mass (left panel) and gas-phase mean metallicity distribution (right panel) in our sample of colliding subhalo pairs (seetext for details). The dashed histogram shows the observed distribution of stellar metallicities in Galactic globular clusters. shifts are currently poorly known as a result of the pres-ence of strong systematic uncertainties affecting metal-licity diagnostics. A few studies, mostly at z ∼
0, havetried to extend the MZR to the low-mass dwarf galaxyregime, deriving power-law relations ( Z ∝ m α ∗ ) withslopes α = 0 . ± .
03 (Berg et al. 2012), ≃ . . ± . . ± .
08 (Blanc et al. 2019).Here, we adopt the intermediate-range value of α = 0 . Y -interceptequal to 5.50, i.e.log( Z/Z ⊙ ) = 12 + log(O / H) − .
69= 5 .
50 + 0 .
37 log( m ∗ / M ⊙ ) − . . (33)The zero-point was chosen to give 12 + log(O / H) = 8 . m ∗ / M ⊙ ) = 7 .
25, in agreement with the DEEP2low-mass MZR (Zahid et al. 2012; Blanc et al. 2019). Lit-tle is known about the redshift evolution of the MZR ondwarf-galaxy scales; in order to minimize the number offree parameters and for ease of interpretation, we assumeno evolution with time in the following (see also Hidalgo2017).The expected GC metallicity, computed by averag-ing the gas metallicities of each colliding pair, is shownin 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 spreadin metallicities that is similar to that observed, witha distribution that is strongly bimodal. In contrast tosome previous work, we do not explicitly assume sepa-rate pathways for the formation of blue and red globulars,and simply predict the metallicity of each GC based onthe enrichment level of the interacting subhalo pair thattriggered its formation. About 30% of the colliding pairshave metallity log(
Z/Z ⊙ ) ≥ −
1: the red globulars arethe result of impacts that involve at least one massive,more chemically evolved satellite, and the metallicity bi-modality reflects a bimodality in stellar and peak halomasses for the colliding subhalos. The age spread is sim-ilar in both the blue and red populations, but the red and blue peaks are shifted towards lower values compared tothe MW GC data.Given the intrinsic uncertainties in the stellar-to-halomass relation, the MZR, and their evolution with redshifton small mass-galaxy scales, it seems again ill-advised todraw definite conclusions from this comparison. Evolu-tionary corrections on the MZR will shift the predictedvalues towards even lower metallicities. We notice heretwo effects that have been ignored and may skew thepredicted distribution in the opposite direction, towardshigher metallicities: 1) collisions involving massive, en-riched subhalos that, albeit rarer, could produce severalGCs per event. The impact of this (unknown) multi-plicity factor has been neglected in Figure 5; and 2)high-speed bound collisions between subhalos falling inon radial orbits and the central most massive progenitor,which may result in a more centrally concentrated sub-population of metal-rich globulars (see also Griffen et al.2010). We count only 17 such collisions, compared to the133 subhalo-subhalo impacts that satisfy our criteria. (v) The Spatial Distribution of GCs:
A common fea-ture of many globular formation models is the reliance onsome ad hoc assumptions to identify the sites where GCsform. The median Galactocentric distance of all knownMW (blue+red) GCs is 5 kpc (Harris 1996), with themetal-poor, [Fe/H] < − z ∼
10 (Moore et al. 2006; Katz & Ricotti 2014).Our collision-driven scenario may offer a new mechanismadau et al. 13for biasing the spatial distribution of GCs relative to theoverall mass profile. This is because, in an inelastic colli-sion, the splash remnant will lose orbital energy and falldeeper into the Galactic potential rather than sharingthe orbits of the progenitor subhalos. It is interesting tobriefly examine here the impact of kinetic energy dissipa-tion on orbital parameters. Consider, for simplicity, twosubhalos moving on coplanar, coaxial, prograde ellipticalorbits in a Keplerian potential of gravitational param-eter µ . The two orbits have the same specific angularmomentum h and the subhalos collide at the semiparam-eter location of the ellipses, p = h /µ . Let us decomposethe velocity vectors along the outward radial direction(ˆ r ) and the prograde azimuthal direction (ˆ θ ) at this po-sition. The radial and azimuthal velocity components ofthe two subhalos just prior to impact are ( µe /h, µ/h )and ( µe /h, µ/h ), respectively, where e and e are theorbital eccentricities, with e < e . In the case of a per-fectly inelastic encounter between two bodies of equalmass m c , conservation of linear momentum determinesthe postcollision instantaneous orbital velocity vector ~v f at p of the combined remnant as ~v f = 12 ( ~v + ~v ) = µ h [( e + e )ˆ r + 2ˆ θ ] . (34)In the absence of external torques the total angular mo-mentum of the system remains unchanged, but the dis-sipation of kinetic energy leads to a net loss of orbitalenergy ∆ E = − m c µ h ( e − e ) ≤ . (35) The new specific orbital energy is ǫ = 12 v f − µ h = µ h (cid:20) ( e + e ) − (cid:21) , (36)and the new eccentricity is e = q ǫ f h /µ = e + e . (37)Another simple situation to analyze is the case of twoclumps moving on coplanar, coaxial prograde orbits withdifferent angular momenta but same pericenter distance r p . Conservation of angular momentum in this caseyields for the new eccentricity of the collision remnant √ e = √ e + √ e , (38)and the change in orbital energy is again∆ E = − m c µ r p ( e + e − e ) ≤ . (39)The above examples emphasize the fact that the remnantwill initially be moving on a less energetic orbit with aneccentricity that is intermediate between those of the col-liding pair. One would expect globulars to progressivelylose memory of their initial infall direction as they or-bit in a host halo that is clumpy and triaxial. Detailedcalculations of the expected radial profile and 3-D dis-tribution of GCs are beyond the scope of this paper andare postponed to future work.We thank Z. Haiman and M. Krumholz for very help-ful comments and suggestions on various aspects of thismanuscript. Support for this work was provided byNASA through a contract to the WFIRST-EXPO Sci-ence Investigation Team (15-WFIRST15-0004), admin-istered by the GSFC (P.M.). REFERENCESAnathpindika, S. 2009, A&A, 504, 437Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53,197Andrews, B. H., & Martini, P. 2013, ApJ, 765, 140Angulo, R. E., Lacey, C. G., Baugh, C. M., & Frenk, C. S. 2009,MNRAS, 399, 983Arreaga-Garc´ıa, G., Klapp, J., & Morales, J. S. 2014,International Journal of Astronomy and Astrophysics, 4, 192Ashman, K. M., & Zepf, S. E. 1992, ApJ, 384, 50Balfour, S. K., Whitworth, A. P., & Hubber, D. A. 2017,MNRAS, 465, 3483Balfour, S. K., Whitworth, A. P., Hubber, D. A., & Jaffa, S. E.2015, MNRAS, 453, 2471Baumgardt, H., Cˆot´e, P., Hilker, M., et al. 2009, MNRAS, 396,2051Baushev, A. N. 2018, NewA, 60, 69Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57Behroozi, P. S., Wechsler, R. H., Lu, Y., et al. 2014, ApJ, 787, 156Belokurov, V., Irwin, M. J., Koposov, S. E., et al. 2014, MNRAS,441, 2124Ben´ıtez-Llambay, A., Navarro, J. F., Abadi, M. G., et al. 2015,MNRAS, 450, 4207Benson, A. J. 2005, MNRAS, 358, 551Berg, D. A., Skillman, E. D., Marble, A. R., et al. 2012, ApJ, 754,98 Bertschinger, E. 2001, ApJS, 137, 1Bica, E., Bonatto, C., Barbuy, B., & Ortolani, S. 2006, A&A,450, 105Blanc, G. A., Lu, Y., Benson, A., Katsianis, A., & Barraza, M.2019, arXiv e-prints, arXiv:1904.02721Boley, A. C., Lake, G., Read, J., & Teyssier, R. 2009, ApJ, 706,L192Bovino, S., Grassi, T., Capelo, P. R., Schleicher, D. R. G., &Banerjee, R. 2016, A&A, 590, A15Boylan-Kolchin, M. 2017, MNRAS, 472, 3120Boylan-Kolchin, M., Bullock, J. S., & Garrison-Kimmel, S. 2014,MNRAS, 443, L44Brodie, J. P., & Strader, J. 2006, ARA&A, 44, 193Bromm, V., & Clarke, C. J. 2002, ApJ, 566, L1Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80Burkert, A., & Forbes, D. 2019, arXiv e-prints, arXiv:1901.00900Cen, R. 2001, ApJ, 560, 592Chiou, Y. S., Naoz, S., Burkhart, B., Marinacci, F., &Vogelsberger, M. 2019, arXiv e-prints, arXiv:1904.08941Choksi, N., Gnedin, O. Y., & Li, H. 2018, MNRAS, 480, 2343Clowe, D., Bradaˇc, M., Gonzalez, A. H., et al. 2006, ApJ, 648,L109Conroy, C., Loeb, A., & Spergel, D. N. 2011, ApJ, 741, 72Creasey, P., Sales, L. V., Peng, E. W., & Sameie, O. 2019,MNRAS, 482, 219
D’Ercole, A., Vesperini, E., D’Antona, F., McMillan, S. L. W., &Recchi, S. 2008, MNRAS, 391, 825Diemand, J., Kuhlen, M., & Madau, P. 2006, ApJ, 649, 1—. 2007a, ApJ, 657, 262—. 2007b, ApJ, 667, 859Diemand, J., Kuhlen, M., Madau, P., et al. 2008, Nature, 454, 735Diemand, J., Madau, P., & Moore, B. 2005, MNRAS, 364, 367Dotter, A., Sarajedini, A., Anderson, J., et al. 2010, ApJ, 708, 698Draine, B. T. 2011, Physics of the Interstellar and IntergalacticMedium (Princeton University Press)Duffy, A. R., Schaye, J., Kay, S. T., & Dalla Vecchia, C. 2008,MNRAS, 390, L64El-Badry, K., Quataert, E., Weisz, D. R., Choksi, N., &Boylan-Kolchin, M. 2019, MNRAS, 482, 4528Elmegreen, B. G., & Efremov, Y. N. 1997, ApJ, 480, 235Elmegreen, B. G., & Elmegreen, D. M. 1978, ApJ, 220, 1051Fall, S. M., & Rees, M. J. 1985, ApJ, 298, 18Feng, Y., & Krumholz, M. R. 2014, Nature, 513, 523Forbes, D. A., Alabi, A., Romanowsky, A. J., et al. 2016,MNRAS, 458, L44Forbes, D. A., & Bridges, T. 2010, MNRAS, 404, 1203Forbes, D. A., Read, J. I., Gieles, M., & Collins, M. L. M. 2018a,MNRAS, 481, 5592Forbes, D. A., Bastian, N., Gieles, M., et al. 2018b, Proceedingsof the Royal Society of London Series A, 474, 20170616Gilden, D. L. 1984, ApJ, 279, 335Gnedin, O. Y., & Ostriker, J. P. 1997, ApJ, 474, 223Gonz´alez, R. E., & Padilla, N. D. 2016, ApJ, 829, 58Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, MNRAS,439, 2386Gratton, R. G., Carretta, E., & Bragaglia, A. 2012, A&A Rev.,20, 50Grcevich, J., & Putman, M. E. 2009, ApJ, 696, 385Griffen, B. F., Drinkwater, M. J., Thomas, P. A., Helly, J. C., &Pimbblet, K. A. 2010, MNRAS, 405, 375Grillmair, C. J., Freeman, K. C., Irwin, M., & Quinn, P. J. 1995,AJ, 109, 2553Habe, A., & Ohta, K. 1992, PASJ, 44, 203Han, J., Cole, S., Frenk, C. S., Benitez-Llambay, A., & Helly, J.2018, MNRAS, 474, 604Han, J., Cole, S., Frenk, C. S., & Jing, Y. 2016, MNRAS, 457,1208Harris, W. E. 1996, AJ, 112, 1487Harris, W. E. 2001, in Saas-Fee Advanced Course 28: StarClusters, ed. L. Labhardt & B. Binggeli, 223Harris, W. E., Blakeslee, J. P., & Harris, G. L. H. 2017, ApJ, 836,67Harris, W. E., Blakeslee, J. P., Whitmore, B. C., et al. 2016, ApJ,817, 58Harris, W. E., Harris, G. L., & Hudson, M. J. 2015, ApJ, 806, 36Hartwick, F. D. A. 2009, ApJ, 691, 1248Hidalgo, S. L. 2017, A&A, 606, A115Hudson, M. J., & Robison, B. 2018, MNRAS, 477, 3869Ibata, R., Nipoti, C., Sollima, A., et al. 2013, MNRAS, 428, 3648Inayoshi, K., Visbal, E., & Kashiyama, K. 2015, MNRAS, 453,1692Jethwa, P., Erkal, D., & Belokurov, V. 2018, MNRAS, 473, 2060Jiang, F., & van den Bosch, F. C. 2016, MNRAS, 458, 2848Jimmy, Tran, K.-V., Saintonge, A., et al. 2015, ApJ, 812, 98Kang, H., Shapiro, P. R., Fall, S. M., & Rees, M. J. 1990, ApJ,363, 488Katz, H., & Ricotti, M. 2013, MNRAS, 432, 3250—. 2014, MNRAS, 444, 2377Kim, S. Y., Peter, A. H. G., & Hargis, J. R. 2018, PhysicalReview Letters, 121, 211302Kimm, T., Cen, R., Rosdahl, J., & Yi, S. K. 2016, ApJ, 823, 52Klypin, A., Kravtsov, A. V., Valenzuela, O., & Prada, F. 1999,ApJ, 522, 82Kravtsov, A. V., & Gnedin, O. Y. 2005, ApJ, 623, 650Kruijssen, J. M. D. 2014, Classical and Quantum Gravity, 31,244006—. 2015, MNRAS, 454, 1658Kuhlen, M., Madau, P., & Krumholz, M. R. 2013, ApJ, 776, 34Kumai, Y., Basu, B., & Fujimoto, M. 1993, ApJ, 404, 144Laevens, B. P. M., Martin, N. F., Sesar, B., et al. 2014, ApJ, 786,L3Larson, R. B. 1985, MNRAS, 214, 379 Leaman, R., VandenBerg, D. A., & Mendel, J. T. 2013, MNRAS,436, 122Li, Y.-S., & Helmi, A. 2008, MNRAS, 385, 1365Lin, D. N. C., & Murray, S. D. 1991, in Astronomical Society ofthe Pacific Conference Series, Vol. 13, The Formation andEvolution of Star Clusters, ed. K. Janes, 55–72 Lokas, E. L., & Mamon, G. A. 2001, MNRAS, 321, 155Ma, X., Kasen, D., Hopkins, P. F., et al. 2015, MNRAS, 453, 960Mackey, A. D., Huxor, A. P., Ferguson, A. M. N., et al. 2010,ApJ, 717, L11Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415Madau, P., Kuhlen, M., Diemand, J., et al. 2008, ApJ, 689, L41Makino, J., & Hut, P. 1997, ApJ, 481, 83Mashchenko, S., & Sills, A. 2005, ApJ, 619, 258Mayer, L., Mastropietro, C., Wadsley, J., Stadel, J., & Moore, B.2006, MNRAS, 369, 1021Moore, B. 1996, ApJ, 461, L13Moore, B., Diemand, J., Madau, P., Zemp, M., & Stadel, J. 2006,MNRAS, 368, 563Moore, B., Ghigna, S., Governato, F., et al. 1999, ApJ, 524, L19Moore, B., Katz, N., Lake, G., Dressler, A., & Oemler, A. 1996,Nature, 379, 613Muratov, A. L., & Gnedin, O. Y. 2010, ApJ, 718, 1266Murray, S. D., & Lin, D. N. C. 1989, ApJ, 339, 933—. 1990, ApJ, 357, 105Nagasawa, M., & Miyama, S. M. 1987, Progress of TheoreticalPhysics, 78, 1250Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490,493Odenkirchen, M., Grebel, E. K., Dehnen, W., et al. 2003, AJ,126, 2385Okamoto, T., Gao, L., & Theuns, T. 2008, MNRAS, 390, 920Peebles, P. J. E. 1984, ApJ, 277, 470Peng, E. W., Jord´an, A., Cˆot´e, P., et al. 2006, ApJ, 639, 95Pfeffer, J., Kruijssen, J. M. D., Crain, R. A., & Bastian, N. 2018,MNRAS, 475, 4309Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2018,arXiv e-prints, arXiv:1807.06209Prieto, J., Jimenez, R., & Haiman, Z. 2013, MNRAS, 436, 2301Ricotti, M. 2009, MNRAS, 392, L45Ricotti, M., Parry, O. H., & Gnedin, N. Y. 2016, ApJ, 831, 204Rosenblatt, E. I., Faber, S. M., & Blumenthal, G. R. 1988, ApJ,330, 191Saitoh, T. R., Daisaka, H., Kokubo, E., et al. 2009, PASJ, 61, 481Saitoh, T. R., Koda, J., Okamoto, T., Wada, K., & Habe, A.2006, ApJ, 640, 22Sarajedini, A., Bedin, L. R., Chaboyer, B., et al. 2007, AJ, 133,1658Sawala, T., Pihajoki, P., Johansson, P. H., et al. 2017, MNRAS,467, 4383Sawala, T., Frenk, C. S., Fattahi, A., et al. 2015, MNRAS, 448,2941Searle, L., & Zinn, R. 1978, ApJ, 225, 357Shen, S., Madau, P., Conroy, C., Governato, F., & Mayer, L.2014, ApJ, 792, 99Smith, J. 1980, ApJ, 238, 842Springel, V., Wang, J., Vogelsberger, M., et al. 2008, MNRAS,391, 1685Stone, M. E. 1970, ApJ, 159, 277Tormen, G., Diaferio, A., & Syer, D. 1998, MNRAS, 299, 728Trenti, M., Padoan, P., & Jimenez, R. 2015, ApJ, 808, L35Valcheva, A. T., Ovcharov, E. P., Lalova, A. D., et al. 2015,MNRAS, 446, 730van den Bosch, F. C., Ogiya, G., Hahn, O., & Burkert, A. 2018,MNRAS, 474, 3043van den Bosch, F. C., Tormen, G., & Giocoli, C. 2005, MNRAS,359, 1029VandenBerg, D. A., Brogaard, K., Leaman, R., & Casagrande, L.2013, ApJ, 775, 134Vesperini, E. 1998, MNRAS, 299, 1019Visbal, E., Haiman, Z., & Bryan, G. L. 2014a, MNRAS, 442, L100—. 2014b, MNRAS, 445, 1056Vishniac, E. T. 1983, ApJ, 274, 152—. 1994, ApJ, 428, 186Wang, Y., Ferland, G. J., Lykins, M. L., et al. 2014, MNRAS,440, 3100 adau et al. 15adau et al. 15