Globular Cluster Formation Within The Aquarius Simulation
The Aquarius project is the first simulation that can resolve the full mass range of potential globular cluster formation sites. With a particle mass , Aquarius yields more than 100 million particles within the virial radius of the central halo which has a mass of , similar to that of the Milky Way. With this particle mass, dark matter concentrations (haloes) as small as 10 M will contain a minimum of 100 particles.
Here, we use this simulation to test a model of metal-poor globular cluster formation based on collapse physics. In our model, globular clusters form when the virial temperatures of haloes first exceed K as this is when electronic transitions allow the gas to cool efficiently. We calculate the ionising flux from the stars in these first clusters and stop the formation of new clusters when all the baryonic gas of the galaxy is ionised. This is achieved by adopting reasonable values for the star formation efficiencies and escape fraction of ionising photons which result in similar numbers and masses of clusters to those found in the Milky Way. The model is successful in that it predicts ages (peak age 13.3 Gyrs) and a spatial distribution of metal-poor globular clusters which are consistent with the observed populations in the Milky Way. The model also predicts that less than 5 of globular clusters within a radius of 100 kpc have a surviving dark matter halo, but the more distant clusters are all found in dark matter concentrations.
We then test a scenario of metal-rich cluster formation by examining mergers that trigger star formation within central gas disks. This results in younger ( 7–13.3 Gyrs), more centrally-located clusters (40 metal-rich GCs within 18 kpc from the centre of the host) which are consistent with the Galactic metal-rich population. We test an alternate model in which metal-rich globular clusters form in dwarf galaxies that become stripped as they merge with the main halo. This process is inconsistent with observed metal-rich globulars in the Milky Way because it predicts spatial distributions that are far too extended.
keywords:globular clusters: general – stars: formation – galaxy: formation – methods: numerical
Globular clusters (GCs) provide a remarkably rich source of information about galaxy formation. Unlike the diffuse stellar populations of galaxies, globular clusters mostly contain stellar populations with a narrow range of ages and are extremely homogenous, making them relatively simple to understand and model. The are extremely old, so they survive as a record of conditions and processes of the earliest stages of galaxy formation (West04). The properties of the Milky Way halo GC population led Searle78 to conclude that the halo of our Galaxy formed via the slow accretion of many small proto-galactic fragments, not via monolithic collapse as previously thought (Eggen62). In the 1990s, the (old) ages of globular clusters provided one of the motivations for considering cosmological models with non-zero cosmological constants (Ostriker95).
In the last two decades observations of extragalactic GC populations (notably with the Hubble Space Telescope) have revealed strong bimodality in the optical colours of GCs (e.g. Ashman92). The blue population is identified as metal-poor clusters and very old, whereas the red population is more metal-rich and not as old.
These observations have motivated a number of competing galaxy and GC formation scenarios (e.g. Forbes97; West04), which attempt to explain the bimodal colours, but no conclusive theory has emerged. However, the very existence of the bimodal colours is indicative of two epochs of star formation. In view of the large amount of data collected, there is a strong need for more detailed theoretical work that will provide specific predictions of where and when GCs formed (Ashman98; Brodie06).
In contrast with the extensive gains from observational studies of GCs, it has proven very difficult to predict the full range of observed GC properties in a self-consistent manner from theoretical models. The mass and spatial scales needed to study the physical conditions of GC formation are very difficult to simulate numerically in large models.
The problem of direct simulation was avoided in an early study by Beasley02 who used a semi-analytic model of galaxy formation in a cold dark matter (CDM) Universe. They assumed that globular clusters would form in each galaxy in numbers proportional to the numbers of stars forming (given by the semi-analytic model). They could then ascribe chemical properties to the GCs according to the same model. Their model successfully reproduce bimodal GC populations and other observations, but only if they invoke a truncation of metal-poor GC formation at redshift . This truncation was also investigated by Bekki05 in a collisionless dark matter simulation of the formation of a single galaxy (mass resolution ) with the simplifying assumption that GCs formed in all low-mass subhaloes forming before some truncation redshift. Bekki shows that the final () radial distribution of the objects is very sensitive to the truncation redshift which was set at to match the GC distribution of the Milky Way. In a more general study of structure formation Moore06 identified reionisation as the process responsible for the truncation, assuming that it took place by . They went on to suggest that the radial distribution of GCs and satellite galaxies could be used to constrain models of the reionisation process. These last two studies established the importance of reionisation and the truncation of globular cluster formation, but were not able to calculate the important parameters directly.
A more comprehensive approach to the problem was taken by Kravtsov05 who combined both gas and N-body codes to model a Milky Way sized galaxy, although the simulation only ran to a redshift of . Their model used the assumption that GCs formed in sufficiently massive giant molecular clouds—themselves in the disks of protogalaxies. This work was successful in matching several observed GC properties such as masses and sizes, but could not predict present-day positions. It was—like the earlier work—very reliant on assumptions about the conditions necessary for GC formation.
The N-body (dark matter) and semi-analytic approaches have recently been combined in a large scale cosmological simulation by Bekki08 to model the dynamical and chemical properties of GCs in a wide range of galaxies. The simulation covers a large volume, so the mass resolution is relatively low () compared to GC masses: they simulate GC formation by assuming every virialised dark matter halo of 10 or more particles will form a GC. They also rely on a truncation of metal-poor GC formation due to reionisation at (this value based on quasar data; see discussion below). From their model Bekki08 obtain old ages for both metal-rich GCs (peaking at z) and metal-poor GCs (peaking at z ) and also obtain bimodal metallicity distributions for about half the galaxies. They also find the galaxies without bimodal distributions tend to be small galaxies which lack metal-rich GCs. The model also produces more centrally concentrated distributions for the metal-rich GCs than for the metal-poor GCs although physical origin of this difference is not specifically discussed.
A common feature of the simulations described above is the reliance on ad hoc assumptions to identify the sites where GCs form. This was inevitable because the mass resolution was too low to resolve individual GCs. In this paper we present the first study of globular cluster formation based on a simulation in which globular cluster masses are well resolved.
We use the Aquarius suite of simulations (Springel08), the highest resolution simulations of Milky Way sized haloes done to date. It must be emphasised that this is the first work in which each GC formation site is directly resolved with a minimum of about 2000 particles. Although this first paper is based strictly on the dark matter components of the simulation, the exquisite resolution allows us to calculate the conditions for metal-poor GC formation directly. We also include a more qualitative model for the formation of metal-rich GCs which successfully predicts their centrally concentrated distributions.
Our model for metal-poor globular cluster formation is different from previous work in two major respects. First, we model the formation by directly identifying when and where the early haloes first reach a temperature of K, the threshold temperature above which rapid cooling takes place leading to collapse and star formation. Secondly, we do not assume an arbitrary value for the redshift when globular cluster formation is truncated. Instead, we directly estimate the number of ionising photons emitted by these early clusters and stop their formation when there are enough photons to ionise the entire galaxy.
We have used the simulation to test two models for the formation of the metal-rich GC population: (i) the stripping of GCs from disrupted satellite galaxies and (ii) the formation of star clusters in the gas disk of the forming galaxy, triggered by large merging events. We find that only the second model can produce the centralised distribution of metal-rich GCs as observed in the Milky Way.
The structure of this paper is as follows: In Section 2, we describe the Aquarius suite of cosmological simulations and our models for forming globular clusters in the simulations. In Sections 3 and 4, we present the results of our models in detail, notably comparing the spatial distributions with observations, first for the metal-poor GCs and then for the metal rich GCs. We then discuss our results in Section 5 and conclude by summarising up our results and future work in Section 6.
2 Models of Globular Cluster Formation
In this section we describe the methodology of our study: how we identify the globular cluster formation sites in the Aquarius simulation data. We start by summarising the relevant details of the simulations themselves. We then describe our models for metal-poor and metal-rich globular cluster formation in detail in Sections 2.2 and 2.3, respectively.
2.1 The Aquarius Simulations
The Aquarius Project actually consists of simulations of six different Milky Way-sized galaxies, one of which is analysed in this paper. Although a detailed description of the simulations can be found in Springel08; hereafter S08, we review the pertinent details here.
The starting cosmological parameters for Aquarius are the same as used in the Millennium Simulation (Springel05) project which are consistent with both the WMAP1 and WMAP5 results (Bennett03; Spergel07). Halo formation is tracked within a periodic box of side 100 Mpc in a cosmology with parameters = 0.25, = 0.75, = 0.9, = 1, and Hubble constant = 100 km s Mpc = 73 km s Mpc. We use a baryon fraction to convert from dark matter mass to baryonic mass.
Millennium simulation haloes of roughly Milky Way mass and without close neighbours at were selected for resimulation using particles in a box of dimension 10 Mpc. Identifying the Lagrangian region from when each halo formed, the mass distribution was rerun at a much higher spatial and mass resolution. Distant regions were sampled with more massive particles, but retained sufficient resolution to ensure an accurate representation of the tidal field at all times. For greater detail on the simulation, see S08.
A major feature of the simulations is the identification of substructure—the bound mass concentrations that will grow and merge over time to build structure. These are identified and measured in Aquarius using the same subfind (Springel01) algorithm used for the Millennium simulation. Although, subfind outputs are commonly called ‘subhaloes’, we use the words subhalo and halo interchangeably to represent potential GC structures. This is because a potential GC can be classed as either an object within much larger parent object (i.e. a subhalo) or an object outside the virial radius of any other halo (i.e. a halo). For future reference, wherever we state the word ‘halo’ and ‘subhalo’, we are referring to the same thing; the subfind outputs.
The way in which the haloes are linked between time steps works as follows. For each halo we take the most bound 10 of its particles and determine which halo they belong to at the next snapshot. The halo at the later time with the largest number of these particles is identified as the descendant. In most cases this is all that is required. However, there are occasional cases where subfind fails to identify a halo at one or more snapshots but picks it up again at a later time. For example, this can happen if a halo passes close to the centre of a larger halo. We attempt to deal with this by looking for descendants more than one snapshot later if a halo is not the most massive progenitor of the descendant identified using the procedure described above or if no descendant is identified. A halo which exists up to 3 snapshots later will be identified as the descendant if it contains more than half of the 10 most bound particles from the original halo and has no identified progenitors.
All these different data structures between time steps for each simulation are stored in a database very similar to that of the Millennium project111http://www.g-vo.org/Millennium. We also make use of the raw particle data in cases where we are unable to obtain the subfind halo information to track material from haloes that have merged. This is discussed in further detail in Section 2.2.2.
Each of the six Aquarius haloes was calculated using at least two different particle masses (‘resolutions’) to test for convergence. We selected our highest resolution halo (‘A’ halo) for our GC study in this paper and analysed it at two different resolutions so as to check our results for any dependence on simulation resolution. A summary of the two data sets is given in Table 1.
2.2 Formation of Metal-Poor Globular Clusters
2.2.1 Temperature Threshold
We use a relatively simple model to identify where the metal-poor GC-type objects would form within the Aquarius simulation based on the conditions necessary for the collapse of a proto-GC gas cloud and subsequent star formation. Recent simulations all rely on a similar model for GC formation (dating back to Peebles84), but the resolution of Aquarius allows us to measure the main parameter – temperature – directly.
The gas clouds cannot collapse without an efficient cooling mechanism: in the absence of significant amounts of heavy elements, the main cooling processes are the collisional excitation of hydrogen and helium, radiative recombination of hydrogen, and bremsstrahlung (e.g. Nishi02). The typical cooling function for primordial gas in the equilibrium state reveals an extremely rapid increase in cooling rate as the temperature rises through 10K. We therefore adopt 10K as a temperature threshold, above which gas clouds can efficiently cool and collapse to form GCs.
Given the very large mass of gas required to form a globular cluster (see Section 3.1), the early proto-cluster gas clouds must form in the potential wells of the dark matter subhaloes identified by subfind. By assuming the gas is in quasi-static equilibrium with the dark matter, we can use the virial theorem to relate the 1-D internal velocity dispersion measured of the dark matter subhaloes () to the (inferred) virial temperature of the gas, :
where is the mass of a hydrogen atom and is the mean molecular weight of the gas. We adopt molecular weight of , appropriate for a fully-ionised, primordial gas. Our 10 K temperature threshold therefore corresponds to a 1-D velocity dispersion of 11.9 kms.
We identify metal-poor GC formation sites by searching the entire merger tree of the simulation for any subhaloes that exceed the K threshold for the first time.
2.2.2 The Final Positions of Disrupted Subhaloes
Although we have excellent resolution to directly locate where candidates first form, it is no longer possible to locate these structures as distinct subhaloes at redshift zero if they have undergone merging events with either each other or the central host halo since it leads to their complete destruction. We address this problem by using the most bound particle of each subhalo as a marker to track its final position. This is equivalent to assuming that the collapsed baryonic gas forming a GC at the centre of each subhalo will follow a similar trajectory to the most-bound particle. We simply search the final simulation snapshot at to locate where each of these uniquely identified particles reside. This directly allows us to follow GC formation sites through to the present day.
As we show below, the majority of the haloes containing GC formation sites merge with the central halo by the present day, but some survive in separate haloes at a range of distances. We include both groups in our model (both merged and surviving haloes) but we restrict our discussion to the GC candidates that end up associated222By associated we mean either fully merged with the main halo or in surviving haloes that are within two times the half-mass radius (2 kpc) of the main halo at redshift zero. with the main Milky Way halo in the simulation at redshift zero. There are a few objects in the outer halo that do not satisfy this condition: we discuss these further in Section 5.
2.2.3 Reionisation within Aquarius from the first GCs
Previous simulations of GC formation have found it necessary to truncate the formation process after a certain redshift in order to avoid producing unrealistically large numbers of clusters (Kravtsov05, Bekki05; Bekki08). The truncation redshift has generally been set on the basis of external estimates of the redshift of reionisation. In this paper we take a different approach, based on an internal calculation of the ionisation from the first star clusters themselves.
Current estimates of when the Universe as a whole was reionised are based on absorption line studies of high-redshift quasars. Observations of Ly systems in high-redshift quasars indicate the inter-galactic medium (IGM) was fully ionised at (Gnedin02). The relative scarcity of quasars at redshifts greater than , points to another source of ionisation. Unless far more quasars are found, the photoionising contributions of high- massive stars seem to be the only plausible reason to account for the missing radiation. The key parameter of measuring their contribution observationally (from the quasar spectra) is the escape fraction of Ly photons from galaxies into the IGM, . Unfortunately this approach can not yet be applied: there is a wide range of estimates of (Madau96; Bianchi01; Ricotti02) and at best they are upper limits, as the absorption of ionising radiation from the molecular cloud and dust extinction are ignored.
In this work we estimate the ionising contribution of massive stars within the forming galaxy directly from the simulation. We assume that the local ionising radiation is dominated by flux from the first GCs to form. We calculate the number of ionising photons from each GC using a Salpeter initial mass function (IMF) based on Population II star formation models from the Starburst99 code by Leitherer99. We do not aim to ionise the entire IGM; instead we specifically ask if local ionising photons from the first GCs are sufficient to ionise the intra-galactic medium of the forming galaxy by some redshift, . Our full calculation in Section 3.1 shows that this does happen, and at earlier redshifts than given by the quasar estimates.
2.3 Formation of Metal-Rich Globular Clusters
It is more difficult to model the formation of the metal-rich GC population in our simulation. First, the defining chemical make-up of this population can only be modelled if we include baryonic gas processes. Secondly, although these clusters are slightly younger than the metal-poor GCs, they are actually the more centrally concentrated of the two populations in the Milky Way. This is counter to what happens in hierarchical merging processes where the first objects to form end up most centrally concentrated. We have considered two very simple models that can be tested qualitatively with our current simulation data.
2.3.1 Model 1: Tidal stripping of satellite dwarf galaxies
Based on previous models of tidal interactions removing the outer envelopes of dwarf elliptical (dE) galaxies (Bekki05), we investigated a stripping model of GC formation. The basic idea here is that the dwarf galaxies have formed their own (metal-rich) GCs that survive to join the main halo after the dwarf is disrupted. We tested this model by identifying as “stripped” any halo which merges with a more massive halo (typically, the mass increase in such a merger is a factor of 10 or more). The maximum mass of the progenitor halo, before merging, is taken to be proportional to the mass of the nuclear globular cluster that forms within it.
As we show below in Section 4, this results in far too extended a distribution. This is the problem mentioned earlier that, in hierarchical models, younger objects have less centrally-concentrated distributions than older ones.
2.3.2 Model 2: Major mergers with the central halo
The proposal that metal-rich GCs are formed in gas-rich mergers of interacting galaxies has been around for quite some time. Toomre72 were the first to investigate these merger events in detail but it was not until much later that the formation of GCs in this framework was discussed Schweizer87. Hubble Space Telescope observations of young massive star clusters forming in the merging Antennae galaxies (Whitmore95) provided strong motivation for this model (see also Holtzman92) which remains the subject of many studies (Zepf06).
We adopt the following simple model of metal-rich GC formation by gas-rich mergers. We search the merger tree for any halo (above some mass threshold) which merges with the central halo. We adopt the premise that during the merger event, stars will form via perturbations in and around the gas disk of the central halo at that particular redshift. Since the merging haloes are destroyed in this process, the location (radius with respect to the central halo) of where the stars form can be approximated by the radius of the gas disk. This we assume is a fraction of the half-mass radius of the central halo at the redshift of the merger event. Since we have access to the half-mass radius of the central halo at each time step, we essentially count the number of merger events at a given time step, and say that all of those in-falling haloes will create stars at some fixed fraction of the half-mass radius.
We adopt the major-merger scenario here because it is motivated not only by a rich range of observational results, but it is the only model we can find which can make metal-rich GCs with a more concentrated distribution than the metal-poor GCs, as is observed in the Milky Way.
3 Metal-Poor Globular Clusters
There are two things to note about this distribution. Firstly, there are far too many GC candidates for a Milky-Way sized galaxy, and secondly the distribution extends down to low redshift. We address these concerns by estimating, in the following section, the redshift at which the globular clusters that form can reionise the protogalaxy. Those haloes that form after reionisation are shown as hollow bars in the figure.
3.1 Reionisation Contributions
We define an ionising efficiency, , to determine the mass of gas ionised by each globular cluster in our model. is the mass of baryons ionised by the GC divided by the baryonic mass of the halo in which the GC first forms (assuming each of our GC formation sites forms only one GC). Each globular cluster can therefore ionise a region greater in baryonic mass than itself by a factor of
where is the star-formation efficiency (the fraction of the protocluster baryonic mass that forms stars), is the mean number of ionising photons per baryon locked up in stars, and is the fraction of those photons that escapes the cluster. In writing down this expression we assume that one ionising photon per baryon is sufficient to ionise the surrounding gas out to (where kpc is the dark matter half-mass radius of the main halo in both simulations). This radius contains the majority of the GCs and of the baryonic mass of the galaxy. This is a reasonable approximation given the uncertainties in the other factors.
For a self-gravitating gas cloud, star-formation efficiencies of order one third or more are required in order to form a bound cluster (Baumgardt07). However if, if not all the gas in a halo ends up in the proto-cluster the SFE could be considerably lower. Here we take because that gives masses for the GCs that seem to agree with observations - see Section 3.2.
The number of ionising photons emitted per baryon of stellar material depends strongly upon the stellar mass. We use the the Population II efficiency curve of Fig. 2 from Tumlinson04 (calculated using the STARBURST99 code of Leitherer99). Averaging over a Salpeter IMF (see Appendix) gives . This figure could be raised by moving towards a more top-heavy IMF. Note that the lifetime of a 10 star is around 40 million years, of order the dynamical time of the first globular cluster haloes that form and less than the dynamical time in later ones. It is a reasonable approximation, therefore, to assume instantaneous feedback.
As discussed in Section 2.2.3, the escape fraction of photons is the most uncertain factor. To obtain the correct number density of metal-poor GCs, we require , which is not unreasonable given the escape fraction in the early Universe would be considerably higher than we observe locally.
Putting all these factors together, we have . Figure 2 shows the cumulative mass of ionised gas for the AqA2 halo, starting at high redshift and moving towards the present. Note that we only include the contribution from GCs that end up within kpc of the central halo in the present day. The two curves correspond to the two different numerical resolutions and show good agreement – at the difference between the two curves represents less than one output time in the merger tree construction. The dashed line shows the total baryonic mass of the halo at the present day (i.e. the mass within a sphere centred on the most-bound particle and enclosing a mean density of 200 relative to the critical density).
The radiation produced by the metal-poor globular clusters is sufficient to ionise a mass equal to that of the present-day galaxy by a redshift of about 13 (Fig. 2). There is a lot of uncertainty in this estimate and varying by a factor of two could give an estimated ionisation redshift in the range 10 to 15. However, using = 210 (z = 13) gives a number of clusters that agrees well with observations and so we stick with that for the rest of our analysis.
Using this value for the reionisation cut in AqA2, we define a total sample of 173 metal-poor globular clusters, of which 125 lie within the 2 of the central halo. The majority of these (105) are no longer in separate dark matter haloes, but have merged with the central halo. The number of GCs we form is quite consistent with the number of metal-poor GCs observed in the Milky Way (103 with [Fe/H] ; Harris91).
In order to calculate the present day masses of our metal-poor GCs, we first need to know what fraction of the baryonic mass of a gaseous protocluster becomes locked up in stars. Baumgardt07 and Weidner07 found that in order to form a bound star cluster, a protocluster requires a star formation efficiency of approximately 0.3. However, this can be much lower if not all the gas within the halo cools to become part of the protocluster: for example, the cluster may form from a small fraction of the gas that has condensed at the centre of the halo. In this paper we adopt a value of 0.07 as this gives current-day GC masses that seem to agree with those in the Milky Way.
We then require an estimate of how much mass loss due to stellar evolution (winds, SNe) and dynamical evolution from tidal stripping and evaporation affects the cluster. Kruijssen08 showed that a typical globular cluster loses 70 of its initial stellar mass and that is the value that we adopt here. Overall, this means we have from an initial dark matter mass of M, a final present day stellar mass of .
Figure 3 compares both our candidates identified before to those of the Milky Way with [Fe/H] (assuming for 13 Gyr old stellar population (Maraston05). We obtain a range of masses between – . The masses are consistent with the mean of, but have a narrower spread than, the observed Milky Way GCs mass distribution. This may be due to observational errors in calculating the mass and/or scatter in the mass-loss from individual clusters.
With the reionisation cut, all the metal-poor GCs form in the redshift range 22-13, corresponding to look-back times of 13.3-13.5 Gyr. The precise age of the globular clusters presumably mimics the assembly time of the galactic halo, but the prediction of a narrow spread in ages will hold for all galaxies.
This is broadly consistent with observations of the Milky Way. From homogenous age dating of 55 Galactic GCs carried out by Salaris02, we know that the majority of the metal-poor GCs ([Fe/H] -1.6) formed approximately Gyr ago with a narrow spread in ages, consistent within the errors with a single formation redshift.
3.4 Spatial Distribution
Figure 4 shows a projection of the final GC locations relative to the most-bound particle in the host galaxy. Filled circles represent those clusters that form before our reionisation cut, and open circles those that form afterwards. The latter population is much more extended than the former: this is to be expected as they form later.
Fig. 5 shows a present day cumulative radial distribution of the Milky Way metal-poor GCs and the AqA2 metal-poor GCs using three different reionisation cuts, correspoding to 105, 210 and 420. The radial distribution of Milky Way GCs and of our model GCs are in agreement, whereas earlier-forming haloes are too centrally-concentrated and later-forming ones too extended. This is very strong evidence that the metal-poor globular clusters do indeed form in the high-redshift haloes that we have identified.
The limited observational evidence on the kinematics of metal-poor globular cluster populations in late-type galaxies is described in Brodie06. There is some indication that different sub-populations may have different rotation properties, as if the galaxies have been built up from mergers of smaller systems. This is consistent with the idea that metal-poor globular clusters formed early on, before their host galaxy.
Overall, there appears to be little net rotation of the observed metal-poor GC population. Our models agree with this result, showing rotation speeds of order 10 km s that are consistent with no net rotation within the sampling errors.
The velocity dispersion of the model globular cluster population is radially biased, having an aniotropy parameter of (here and are the root-mean-square radial and total velocities relative to the galactic centre, respectively). The observational evidence is currently too weak to place strong constraints on .
4 Metal-Rich Globular Clusters
The Milky Way metal-rich GC population is more centrally-concentrated than is the metal-poor one. This is a problem for any formation model that invokes accretion from haloes more massive, and hence later-forming, than those used to define the metal-poor population. As an example, we show in Figure 6 the current-day distribution of metal-rich GCs in a model in which they are the nuclei of stripped dwarf galaxies that have been disrupted by the tidal forces of the Galaxy. In this plot, we estimate the current location of the GC from that of the most bound particle in the dwarf galaxy at the timestep before it lost its identity. In order to match the observed number of metal-rich GCs in the Milky Way, we have included dwarf galaxy haloes whose maximum mass before disruption exceeded but other mass-cuts produce similar profiles. As can be seen from the plot, the radial distribution of such objects is far too extended and cannot possibly represent the metal-rich GC population.
To find a distribution that is more centrally-concentrated than that of the metal-poor GCs we must abandon models that form clusters within sub-haloes and instead form them directly within the galaxy itself. As mentioned in Section 2.3.2, there is observational evidence for the formation of massive star clusters in merging galaxies and so that is the model that we turn to here.
It is not obvious how large a merger is needed in order to generate enough disturbance to trigger GC formation. To match the number of metal-rich GCs in the Milky Way, we assume that a single GC forms whenever a galactic halo of mass or more merges with a larger halo. We further need to assume that the infalling halo must have a mass of at least 1 per cent of the mass of the larger one; otherwise large numbers of GCs would be formed by accretion of small satellites onto the galaxy at late times.
As the cold gas is located at the centre of the haloes, that is where the GC will form. For satellite haloes, we use the location of the most-bound particle as a tracer of the GC position at later times. For the main halo, we assume that the GC forms at a small fraction, 0.1, times the half-mass radius.
There are a lot of ad-hoc assumptions in this model: we don’t know how massive an infalling satellite must be to trigger GC production; we don’t know how many star clusters will form and what their mass will be, and we don’t know the precise location in which they will form. Nevertheless, we present broad-brush results below in order to demonstrate that the model GCs have approximately the right properties and to motivate further study.
4.1 Spatial Distribution
The current-day spatial distribution of the GCs formed in the merging model is shown as the dashed line in Figure 6. The central concentration of GCs within 10 kpc originates from mergers with the main halo. This agrees well with the observed distribution of metal-rich GCs in the Milky Way. In addition, about a quarter of the GCs result from mergers in satellite galaxies that later fell into the main halo. These latter objects have a more extended distribution and do not seem to have analogues in the Milky Way metal-rich GC population.
In Fig. 7, we show the number of metal-rich GCs formed as a function of look-back time. Also shown, for comparison, are the metal-poor GCs.
We know from various observational studies of metal-rich GCs, that the majority formed over a much larger temporal range than metal-poor GCs, probably because they form via processes which require larger dynamical times, (DeAngeli05; Salaris02; Harris91). These observations are broadly consistent with our findings. However, there are seven GCs that have younger ages, between 0.8 and 3.8 Gyr, which we have not shown in Figure 7. These do not have analogues in the Milky Way, although we note that young “intermediate age” globular clusters are found in M31 (SSL09, and references within). It is possible that they are an accident of the formtion history of this particular galactic halo, or that the galactic disk has become so depleted by these late times that it is no longer susceptible to GC production in minor mergers.
In this section we discuss more general aspects of our simulation of the globular clusters starting with the achievements and limitations of the models.
The model proposed in this paper for the formation of metal-poor GCs is based on just two parameters: the temperature threshold for cooling, and the ionising efficiency . The first of these is fixed by cooling physics and is not a variable in the model. The second parameter, the ionising efficiency, is relatively uncertain (especially for the escape fraction) as discussed in Section 2.2.3. In our calculation of we have chosen reasonable values for the uncertain parameters, but they were adjusted to give a good match to the total number of metal-poor GCs observed in the real Milky Way.
Our adopted value () of the ionisation efficiency corresponds to the suppression of GC formation after a redshift of z = 13. (Even if we allow for a possible factor of 2 uncertainty in , this only gives a range of z = 10–15.) This redshift is significantly higher than the estimates of the reionisation redshift of the universe derived from QSO studies (z = 6.4), but if GCs form high-mass stars first, the local intra-cluster medium is going to be far more ionised than the rest of the Universe.
This model, whereby GC formation is truncated when the clusters themselves reionise the remaining gas in the protogalaxy, has one important feature: it naturally predicts that the number of metal-poor GCs is proportional to the baryonic mass of the galaxy. For a constant mass-to-light ratio this is directly equivalent to saying that galaxies of this type will all have the same specific frequency of globular clusters, , as observed. The prediction of proportionality is robust even with the uncertainty in the ionising efficiency calculation; if that calculation becomes more precise we will also be able to predict the absolute values of specific frequency.
Given that our model has involved some adjustment of one parameter to match one observable (the total number of GCs formed), its success can be demonstrated by testing it against other observations. In Section 3 above we show that the model successfully reproduces both the the radial distributions and the ages of the metal-poor GCs. The agreement with these two independent measurements gives us a high degree of confidence in our model.
An important concern for both the metal-poor and the metal-rich models is that the results should not be biased by the resolution of the particular simulation used. We can test this in a very straightforward manner for the current models by repeating the models with a lower-resolution simulation: if the same results are obtained this demonstrates that our models are not affected by mass resolution. Our main analysis uses the AqA2 simulation (mass resolution of ); we have repeated both models with the lower-resolution AqA3 simulation (). In Figure 8 we show the result of the comparison by plotting the radial distributions of both metal-poor and metal-rich GCs produced in both simulations. For both models the agreement is excellent: the total numbers produced agree to within 5 per cent and a Kolmogorov–Smirnov test shows they are drawn from the same distribution. This agreement shows that the results are not biased by the resolution of the simulations used. Note that complete agreement is not expected for a variety of reasons: the extra substructure will cause haloes to form at slightly different times, and the trajectories of the most-bound particles will also be altered. (The lower-resolution AqA3 model of metal-rich GCs extends to larger radii in Figure 8 than the corresponding AqA2 objects only because of two haloes at large radii are just below the mass limit in AqA2, but have slightly higher masses in the low-resolution simulation.)
In our model of the metal-poor GCs we use the most bound particle (at the time of formation) of the dark matter halo hosting each GC to trace the final positions of GCs whose haloes have merged with the central galaxy halo. This may not be a good estimate of the final position if the GCs evolve in a different way to the single dark matter particle; possible processes are dynamical friction, tidal disruption and disk shocking. The time scale for dynamical friction (Chandrasekhar43) is inversely proportional to mass, so the GCs will be more affected by this process than the dark matter particles (which are of order 100 times less massive). For GCs of mass a few times at radii of 10 kpc the dynamical friction time scale is years, so this is unlikely to affect our results. Tidal effects can totally disrupt globular clusters on time scales of to years (Kruijssen08), so this may be an issue, but a much more detailed model of the evolution of the cluster candidates would be required to investigate this effect which we defer to future work. Finally, disk shocking is known to have a disruptive effect on GCs that pass through the disk of a galaxy. This has a time scale of about years for Milky Way GCs (Ostriker72); this will affect the GCs in our model, but as it primarily removes loosely-bound stars from the outer parts of the GCs, it is unlikely to cause any separation of the most-bound particle from the GC.
We have also used the model to estimate the final masses of the GCs. This involves a much higher degree of uncertainty as we have to estimate the efficiencies of both cluster formation and their subsequent evolution. We adjusted the star formation efficiency to give mean masses consistent with Milky Way GCs, but the results then suggest a slightly narrower range of mass than in observed clusters. Interestingly, our results do not provide strong evidence for a power law distribution in the masses of the metal-poor GCs at formation as has been assumed in some studies of their subsequent evolution (Prieto08). In future work, it would be valuable to investigate the evolution of the GCs formed in our model in a similar way.
The approach used above to estimate the GC masses was necessary in part because most of the GCs we analyse do not survive as separate haloes to the present day, but merge with the main halo. This is demonstrated by Fig. 9, which compares the total radial distribution at redshift zero of all the metal-poor GCs to those which have merged with the central halo. The gap between the two curves in the figure indicates the small number of GCs still found in surviving haloes: this is a very small fraction, less than 10 per cent within a radius of 100 kpc. The presence of these surviving haloes does suggest that some GCs associated with the Milky Way will have retained some dark matter, but at these relatively large radii these are may be associated with dwarf satellite galaxies rather than isolated GCs.
If we now consider all the GCs found in surviving dark matter haloes distinct from the central halo in the full simulation (i.e. at larger distances from the main halo), we can measure their (dark matter) masses and positions as shown in Fig. 10. In this figure there are 47 distinct haloes containing 68 GCs and the haloes which contain more than one GC are indicated by different symbols. The overall spread of properties is quite similar to that of dwarf galaxies in the Local Group, both in distance and mass (Mateo98). The objects in the figure can be divided into two groups. The smaller haloes () are nearly all at small radii (100 kpc) and are single GCs: these presumably correspond to small satellite galaxies of the Milky Way that still retain some dark matter. On the other hand, the more massive haloes () tend to lie at large radii (100 kpc) and often contain multiple GCs. These correspond to the larger dwarf galaxies of the local group. We must note, however, that our metal-poor GC formation model does not necessarily apply to the more distant objects as they would not have experienced the same ionisation environment as those forming closer to the main halo.
In a related study of the Aquarius simulations, Gao09 have investigated the formation of the very first stars, systems with temperatures around K that form at redshifts of or higher. In that context they also considered systems formed by line cooling ( K) at later times. These objects are very similar to what we identify as candidate GC formation sites in this paper, although Gao et al. refer to them as ‘first galaxies’ and focus on a comparison of the surviving objects with observed dwarf galaxies around the Milky Way, obtaining similar results to those we describe here.
As we note above, our model for the formation of the metal-rich GCs relies on several assumptions: it is mainly intended to let us determine if any class of model can come close to reproducing the extremely concentrated distribution of these objects observed in the Milky Way. Our merger model is significantly better than any hierarchical model we have considered in reproducing the observed combination of a centrally-concentrated distribution and young ages for these objects. The model does produce a few metal-rich GCs at larger radii and younger ages than observed in the Milky Way, but this is not significant given the many other uncertainties in the model.
6 Summary Future Work
In this study, we have made use of the exquisite resolution of the Aquarius simulations to test plausible metal-poor and metal-rich GC formation models. Here we summarise our main results and indicate our plans for future work.
6.1 Metal Poor GCs
We adopted a relatively simple formation scenario for metal-poor GCs by identifying all haloes which go above the 10K threshold required to cool and form stars. We then measured the ionising contribution from the very first GCs and calculated when the entire Galaxy became ionised. When this occurred, we halted GC formation.
There is some uncertainty in our calculation of the amount of ionisation from the first GCs. We have therefore treated this as a free parameter which we have adjusted so that our model produces a similar total number of metal-poor GCs to those observed in the Milky Way. Even though this has been adjusted, we should emphasise that the values chosen are quite reasonable and that changing the value by as much as a factor of 2 will only vary the redshift when GC formation is suppressed over the range of z = 10-15. Our adopted value of the ionisation contribution (mass of baryons ionised per mass of baryons in the halo forming each GC) is . This results in the galaxy being ionised by a redshift of z = 13 when 173 metal-poor GCs have formed.
We have tested this model for the effects of numerical resolution by repeating the analysis for a lower-resolution simulation. We obtained almost identical results between the two resolution runs, indicating that the model is not biased by the simulation resolution.
Having fixed the number of GCs, we then find that the model successfully predicts two independent properties of the metal-poor GCs: their spatial distribution and formation ages of the GCs. The radial distribution of the model GCs shown in Figure 5 is a very good match to that of the observed Milky Way population. We obtain mean ages of 13.3 Gyrs, consistent with observations by Salaris02. The agreement of our model with these two independent observations provides strong support for our approach.
We can also estimate the final masses of the GCs if we introduce a second variable parameter, the star formation efficiency in the proto-clusters. We have adjusted this (again within reasonable values) to match the observed masses of Galactic metal-poor GCs. Although the mean mass has been set, we note that we predict a narrower range of mass than observed.
6.2 Metal Rich GCs
Our approach to the formation of metal-rich GCs is much more qualitative as the process is much more dependent on gas processes that are not directly described in the current simulations. Our aim was to find classes of model that could reproduce distributions like those of Galactic metal-rich GCs which are much more centrally-concentrated than the metal-poor GCs. This distribution cannot be produced by hierarchical processes because the metal-rich GCs are younger than the metal-poor GCs; in hierarchical processes it is always the oldest objects that are most centrally-concentrated.
We tested two formation mechanisms: cluster formation triggered in the gas disk of the central halo by large merging events and the stripping of infalling haloes (presumed to contain GCs) as they merge with the central halo. We rejected the stripping model because it produces a spatial distribution of GCs in the present day that is far too extended (see Fig. 6).
The only model that produced a sufficiently central distribution was the merger model. In the model we consider that GCs form when galactic haloes of mass larger than M that merge with a larger halo. We then assume that the GC formation takes place in the central gas disk at a radius estimated as 0.1 times the dark matter half-mass radius. Although there are quite a few assumptions in this model, this model produces the correct distribution. The model also predicts a large spread of ages (8 - 13 Gyrs) that is mostly consistent with the observed Galactic age estimates of Salaris02.
As with the metal-poor GCs, we find that this model is not biased by the simulation resolution: when we apply the model to the lower-resolution simulation there is no observable change in the radial distributions within 100 kpc of the central halo.
6.3 Future Work
There are a number of avenues for future work. Similar analysis could be carried out on the remaining five Aquarius haloes. They have slightly different central halo masses and spatial resolutions, so they will allow us to test the same formation mechanisms across subtly different evolutionary environments. We have already carried out similar calculations on the AqF Aquarius halo and found comparable results.
The next major step is to introduce semi-analytic modeling of Bower06. This could give insight into how merger events and halo accretion can alter GC properties. With respect to the metal-rich GCs, the nature of the gas disk during each merger event could be inferred from such models and used to more accurately determine where in the Galaxy these star-forming regions will occur.
Future work will also include GC formation within the Millennium-II (MII) simulation which, although it has a lower resolution (M), can still locate GC formation sites with a minimum of 10 particles per GC. The scale of the MII simulations (box-side length of 100 Mpc/h) will enable us to correlate GC properties with those of their hosts in a wide range of physical environments.
We wish to thank Bill Harris for many valuable comments about the observational data.
The simulations for the Aquarius Project were carried out at the Leibniz Computing Center, Garching, Germany, at the Computing Centre of the Max-Planck-Society in Garching, at the Institute for Computational Cosmology in Durham, and on the ’STELLA’ supercomputer of the LOFAR experiment at the University of Groningen.
This work was supported by travel funding from the Australian Research Council (grant LX0775963). PAT was partially supported by an STFC rolling grant. Brendan Griffen would like to acknowledge the support provided by the University of Queensland via a University of Queensland Postgraduate Scholarship.
Here we estimate , the mean number of ioising photons per baryon, averaging over a Salpeter initial mass function. The calculation is easily extended to other mass functions that can be approximated as a power law above 10 .
The Salpeter IMF is
where is the number of stars in the mass interval , is the total mass of stars, and and are parameters describing the slope and lower mass-cut of the population.
The number of ionising photons per baryon as a function of stellar mass for a typical Population 2 metallicity of 0.001, derived using the Starburst 99 code of Leitherer99, is given in Figure 2 of Tumlinson04. We approximate this a sequence of piecewise, linear fits in :
Tumlinson04 do not give the value of for masses above 120 ; there are few stars of this mass and so their contribution to is relatively small: we simply take to be constant in this regime.
Within each piecewise interval the contribution to is
where and and are appropriate constants taken from Equation 4. This integrates to give
This expression can be simplified by noting that, when summing these expressions over the whole mass range, the first term in the square brackets vanishes whenever is a continuous function.
Putting in values appropriate to a Salpeter IMF gives .