Exploiting a fundamental characteristic of galaxy assembly in the Λ cold dark matter paradigm, the observed spatial biasing and kinematics of metal-poor globular star clusters are used to constrain the local reionization epoch around individual galaxies. Selecting three galaxies located in different environments, the first attempt at constraining the environmental propagation of reionization in the local Universe is carried out. The joint constraint from the three galaxies ( ) agrees remarkably well with the latest Wilkinson Microwave Anisotropy Probe constraint on zreion for a simple instantaneous reionization model. More importantly, the range of zreion values found here is consistent with the global range of zreion estimates from other observations. We furthermore find a 1.7σ indication that reionization completed in low-density environments before the intergalactic medium in high-density environments was reionized. This is consistent with certain theoretical models that predict that reionization was globally prolonged in duration, with neutral hydrogen pockets surviving in high-density environments, even after the surrounding regions were reionized. More generally, this work provides a useful constraint on the formation history of galaxy stellar haloes.
Reionization marks the phase change of the Universe’s intergalactic hydrogen gas from a neutral to an ionized state. Though the general reionization picture follows from the fundamentals of big bang cosmology ( Barkana & Loeb 2001 ), a detailed timeline of how it propagated throughout the Universe is difficult to constrain because most of the direct observational signatures of this event are > rsim 12 billion light-years distant from Earth.
Despite this challenge, observations show that the intergalactic medium was reionized no later than a redshift of z ∼ 5, not before z ∼ 12 and completed on a time-scale of at least Δ z > 0.06 ( Fan et al. 2006 ; Bolton et al. 2010 ; Bowman & Rogers 2010 ; Ouchi et al. 2010 ; Robertson et al. 2010 ; Larson et al. 2011 ). Though the dominant ionizing source contributing to reionization is still an open question, there is general agreement that the ionizing sources first appeared in high-density environments ( Dijkstra, Haiman & Loeb 2004 ; Fan et al. 2006 ; Lidz et al. 2007 ; Wise & Abel 2008 ; Power et al. 2009 ; Baek et al. 2010 ; Bunker et al. 2010 ; Robertson et al. 2010 ; Srbinovsky & Wyithe 2010 ; Willott et al. 2010 ; Yan et al. 2010 ; Bouwens et al. 2011 ; Dopita et al. 2011 ; Lorenzoni et al. 2011 ; Mitra , Choudhury & Ferrara 2011, 2012 ).
A key piece of information that is not yet constrained is the environmental propagation of reionization: did reionization complete in dense environments first or were low-density voids the first locations to ionize? While high-density environments of the Universe probably had a ‘head-start’ and thus contained more ionizing sources, the same locations were also more likely to be shielded by dense clouds of dust and neutral hydrogen (H i ) gas. Furthermore, if the sources were predominately small galaxies, the ultraviolet (UV) background would have been more uniformly distributed than if the ionizing sources were typically active galaxy nuclei of rare, massive protogalaxies. Also, unlike UV radiation, X-rays produced by gas accretion on to black holes can escape local gas and dust absorption and ‘pre-heat’ gas in lower density regions (e.g. Baek et al. 2010 ). Indeed, large-scale radiative transfer models show that the ionizing front can travel vast distances, from massive galaxies in dense regions to low-density locales (e.g. Weinmann et al. 2007 ; Finlator et al. 2009 ; Iliev et al. 2011 ). Understanding how reionization propagated through various environments will help constrain the properties of the ionizing sources and the state of the intergalactic medium during the reionization epoch.
One possible surviving relic of the reionization epoch may be the bimodal distribution of globular star clusters within large galaxies ( Forbes, Brodie & Grillmair 1997 ; van den Bergh 2001 ; Beasley et al. 2002 ; Santos 2003 ; Bekki 2005 ; Rhode, Zepf & Santos 2005 ; Moore et al. 2006 ; Rhode et al. 2007 ; Bekki et al. 2008 ; Griffen et al. 2010 ; Rhode, Windschitl & Young 2010 ): systems of globular clusters (GCs) show two distinct metallicity subpopulations ( Zinn 1985 ; Gebhardt & Kissler-Patig 1999 ; Neilsen & Tsvetanov 1999 ; Kundu & Whitmore 2001 ; Larsen et al. 2001 ; Peng et al. 2006 ; Strader et al. 2006 ; Kundu & Zepf 2007 ; Strader, Beasley & Brodie 2007 ; Spitler, Forbes & Beasley 2008b ). The subpopulations include the ‘metal-poor’ globular clusters (MPGCs) with typical metallicities a few per cent of solar metallicity and the ‘metal-rich’ globular clusters (MRGCs) with metallicities ∼30 per cent solar.
Various theoretical works have invoked reionization as a way to produce the bimodal metallicity distribution in GC systems (e.g. Beasley et al. 2002 ; Santos 2003 ; Bekki 2005 ; Moore et al. 2006 ; Bekki et al. 2008 ; Griffen et al. 2010 ). More specifically, these studies explored models where MPGCs form within small ( Mvir ∼ 10 8 M ⊙ ) dark matter (DM) haloes at z ≳6 ( Bromm & Clarke 2002 ; Boley et al. 2009 ). In this model, the MPGCs formed until cold molecular clouds within these haloes were heated up during the reionization epoch. This suppressed GC formation for a period of time and allowed gases within the host galaxy to enrich in metals. After some time had passed, a second generation of GCs, the MRGCs, started to form during galaxy merger events (e.g. Ashman & Zepf 1992 ; Beasley et al. 2002 ; Griffen et al. 2010 ) and/or the growth of galaxy discs ( Shapiro, Genzel & Förster Schreiber 2010 ) at redshifts z ∼ 2–4 ( Shapiro et al. 2010 ; Spitler 2010 ). Although the global nature of the reionization epoch may provide a convenient way to explain the ubiquitous GC metallicity bimodality observed in nearby galaxies, a number of interesting alternative models have been proposed to explain GC metallicity distributions without reionization (e.g. Côté, Marzke & West 1998 ; Scannapieco, Weisheit & Harlow 2004 ; Pipino, Puzia & Matteucci 2007 ; Hasegawa, Umemura & Kitayama 2009 ; Gray & Scannapieco 2010 ; Muratov & Gnedin 2010 ; Shapiro, Genzel & Förster Schreiber 2010 ).
If reionization truncated GC formation, the typical ages of MPGCs should correlate with the local reionization epoch. This means MPGCs can be used to measure zreion for individual galaxies. Current MPGC age estimates overlap with the expected reionization epoch. Galactic MPGCs are constrained to be older than ∼11 Gyr or z ∼ 2 ( Krauss & Chaboyer 2003 ), with more stringent constraints on individual MPGCs suggestive of a pre-reionization epoch formation (e.g. Hansen et al. 2004 , though exceptions do exist: Hansen et al. 2007 ). Constraints on the ages of individual extragalactic MPGCs are poorer (> rsim 10 Gyr; see references in Brodie & Strader 2006 ), but are consistent with a pre-reionization formation. Also, within the 1σ uncertainties the mean MPGC metallicities of various early-type galaxies are statistically consistent with preliminary measurements for galaxies at redshifts z ∼ 6–8 ( Bouwens et al. 2010 ; Finkelstein et al. 2010 ; Labbé et al. 2010 ). This is shown in Fig. 1 and is consistent with a scenario where MPGCs were in place at redshifts greater than z ∼ 7, well into the expected reionization epoch.
Unfortunately, the limited observational information about extragalactic GCs and large modelling uncertainties mean that direct methods to age-date extragalactic GCs do not provide a useful constraint on the local reionization epoch around individual galaxies.
A novel way to determine when the MPGC formation epoch finished was recently proposed by Diemand, Madau & Moore (2005) and Moore et al. (2006) . The technique is based upon a fundamental property of a Λ cold dark matter (ΛCDM) universe, where the most massive or rarest DM haloes (and anything that formed within them) will tend to be centrally concentrated within their final host haloes and show unique kinematic signatures at the present epoch. Thus by constraining the spatial and kinematic properties of MPGCs today, information about their progenitor haloes can be recovered. Since the rarity of a halo for a given virial mass depends on redshift, this information can be used to age-date MPGCs with enough accuracy to constrain the reionization epoch around the galaxy hosting the GCs ( Diemand et al. 2005 ; Moore et al. 2006 , see also Bekki 2005 ).
This technique was demonstrated in Moore et al. (2006) using the Milky Way’s GC system. They found evidence for a reionization redshift of zreion = 10 ± 2 (updated for a 7-year Wilkinson Microwave Anisotropy Probe , WMAP 7, cosmology; Komatsu et al. 2011 ). In the present work, this technique is applied to two additional galaxies located in denser environments than the Milky Way to help constrain the temporal propagation of reionization through the local Universe.
Section 2 provides a detailed background, describing in detail both the theoretical framework of the ‘Diemand–Moore’ methodology and relevant GC observations. In Section 3 , the observational data are presented. The main mass modelling and zreion measurements are presented in Section 4 . It is followed by a discussion in Section 5 .
2 METHOD AND ASSUMPTIONS
It has been shown that the growth of structure in the Universe can be thought of as a Gaussian mass-density field where mass-density fluctuations form and assemble hierarchically ( Press & Schechter 1974 ; Lacey & Cole 1993 ). The rare mass-density fluctuations in the Universe are biased towards high-density regions of the Universe ( Cole & Kaiser 1989 ; Sheth & Tormen 1999 ). As a consequence, the rarest fluctuations will tend to be the most centrally concentrated material in a galaxy ( Moore et al. 1998 ; White & Springel 2000 ). By convention, the height of the fluctuations is represented by ν, the number of standard deviations above the mean mass-density level at that epoch.
Diemand et al. (2005) used N -body simulations to explicitly show that the location and kinematic properties of mass structures within a galaxy today depend on the rarity of the mass structure or ‘halo’ when it reached its Jeans mass and gravitationally collapsed. They found that the rarest haloes are more likely to accrete on to a galaxy at earlier times, when the gravitational potential of the galaxy was relatively small. Thus any object that formed within a rare halo is more likely to have a smaller infall velocity and will occupy the central regions of the galaxy until the present epoch. This means that an object’s location and orbit or kinematics within a halo today can be used to recover the properties of the halo it formed within ( Diemand et al. 2005 ; Moore et al. 2006 ). Examples of ‘tracer’ objects include a galaxy’s satellites, its stellar halo, remnants of Population III stars and MPGCs.
The ‘Diemand–Moore’ method is a way use the spatial concentration of MPGCs relative to the hosting galaxy’s halo mass profile to constrain the rarity, νσ, of the haloes that the MPGCs formed within. Results from computer simulations are used to translate the spatial bias of the MPGCs into a νσ constraint. The νσ value alone does not constrain the redshift when MPGCs stopped forming, so an additional assumption is made to break the degeneracy between redshift and the typical halo mass for a given ν. By invoking the argument that star formation only occurs efficiently in haloes with virial temperature more than Tvir = 10 4 K, the degeneracy can be broken and the redshift when MPGCs were suppressed by reionization is found. More details are provided in the following sections.
The Diemand–Moore method requires at least three ingredients to constrain the νσ value of MPGC progenitor haloes: MPGC surface density profile from observations, a halo mass model for the galaxy hosting the MPGCs and a theoretical framework (here Diemand et al. 2005 is used) to interpret the spatial bias of the MPGCs, relative to the host galaxy’s mass profile. MPGC kinematic information from line-of-sight velocity measurements is an optional fourth ingredient that can be used for an additional, albeit weaker, constraint on νσ. The kinematics also proves to be a useful tool to constrain the galaxy mass model in a self-consistent manner to the theoretical framework.
The main physical model that the following analysis depends upon is outlined here.
Each MPGC formed at high redshift in a small, rare DM halo (likely together with other MPGCs and a small protogalaxy).
The bath of UV radiation associated with the reionization epoch was such that MPGC formation was ‘instantaneously’ suppressed 1 in the GC formation sites. This allows one make a simple link between the MPGCs and a single νσ class of haloes, thereby constraining zreion exactly.
The halo containing MPGCs was later accreted by the host galaxy halo, i.e. it became a subhalo.
This accretion process effectively disrupts the early accreted, biased, rare DM subhaloes, but the GC tracers remain intact and will maintain the spatial and kinematic properties of the subhaloes even to the present day ( Diemand et al. 2005 ).
Should the second assumption prove to be incorrect, the analysis still provides an interesting characterization of the progenitor subhaloes of MPGCs and perhaps the other mechanism(s) that truncated their formation. Also, the tracers need not be MPGCs and can include e.g. the Milky Way stellar halo ( Diemand et al. 2005 ) and its satellite galaxies ( Moore et al. 2006 , see also Ocvirk & Aubert 2011 ). Indeed, a useful feature of the Diemand–Moore technique is that it does not care about the detailed baryon physics of e.g. star formation, but only requires an observable tracer of subhaloes that accreted relatively early in a galaxy’s formation history.
3 TARGET GALAXIES AND THE OBSERVATIONS
In this section, the observational data are described. To constrain νσ for MPGC progenitors in a galaxy, the MPGC surface density profile must be measured out to large radii (typically tens of arcminutes), where contamination starts to contribute significantly to the profile. To reduce contamination levels, high-quality imaging is needed. To statistically subtract the remaining contamination, spectroscopic information can be used (e.g. Strader et al. 2011 ) or the profile must extend to large enough radii to characterize contamination levels. Furthermore, the galaxy must have relatively ‘well-behaved’ kinematics so that an estimate of the mass model can be derived.
The extragalactic targets described here are the only galaxies that currently meet these demanding requirements. The galaxies also conveniently reside within different galaxy environments, including the ‘high-density’ environment of a cluster galaxy, a ‘medium-density’ group and a ‘low-density’ field galaxy.
3.1 Messier 87
Messier 87 (M87) is the central galaxy associated with the Virgo cluster of galaxies. It is located 16.5 Mpc away and the cluster contains hundreds of member galaxies. Though M87 and the Virgo cluster are frequently treated as the same thing (at least in terms of their ‘common’ DM halo), recent analysis supports early suggestions that the Virgo cluster is not a single massive galaxy cluster halo, but is instead made up of a set of subhaloes in an unrelaxed state ( Binggeli, Tammann & Sandage 1987 ; Strader et al. 2011 ). The M87 galaxy-sized halo is thought to be virialized, but is a distinct subhalo from the rest of the Virgo structure (see discussion in Doherty et al. 2009 ; Strader et al. 2011 ; Romanowsky et al. 2012 ).
From kinematic tracers within ∼100 kpc, the extrapolated virial mass ( Mvir ) of the M87 halo is Mvir ∼ 10 14 M ⊙ ( Strader et al. 2011 ). The surrounding Virgo cluster halo is Mvir ∼ 10 15 M ⊙ (e.g. McLaughlin 1999b ; Strader et al. 2011 ), which means M87 resides in high-density environment and was therefore one of the first locations in the local Universe to host star formation. It has even been proposed that ionizing radiation from sources in Virgo had enough of a head-start to travel to the Local Group and reionize the Milky Way galaxy (e.g. Weinmann et al. 2007 ). For these reasons, M87 and the Virgo cluster are a particularly interesting location to measure the local reionization epoch.
The M87 GC catalogue is presented in Strader et al. (2011) and was constructed from archival Canada–France–Hawaii Telescope (CFHT)/Megacam ( Boulade et al. 2003 ) images. Kinematic information for the GCs come from Strader et al. (2011) , who present new Keck/Deep Imaging Multi-Object Spectrograph (DEIMOS; Faber et al. 2003 ), Keck/Low Resolution Imaging Spectrometer (LRIS; Oke et al. 1995 ) and MMT/Hectospec ( Fabricant et al. 2005 ) line-of-sight velocities. Out of an estimated MPGC population of ∼14000 MPGCs ( Strader et al. 2011 ), 289 have useful velocity information.
The M87 MPGC surface density was derived in Strader et al. (2011) using standard reduction techniques and methodologies. Briefly, GCs were photometrically selected down to i = 22.5 mag over a region extending to radii ∼130 kpc from M87. Hubble Space Telescope imaging and the velocity information were used to accurately constrain the amount of contamination present in the GC catalogue, which was subtracted off the MPGC surface density profile before analysis.
3.2 NGC 1407
The second target of this study is an elliptical galaxy named NGC 1407, whose galaxy group contains ∼30 members ( Brough et al. 2006 ). It is at 21 Mpc distant and has a sizable GC system ( Forbes et al. 2006 ; Harris et al. 2006 ). Romanowsky et al. (2009) found the galaxy group, which is centred on NGC 1407, to have a relatively high mass-to-light ratio for its total virial mass Mvir ∼ 6 × 10 13 M ⊙ . This target provides a probe of reionization in a ‘medium-density’ environment.
Two Subaru/Suprime-Cam ( Miyazaki et al. 2002 ) image data sets of NGC 1407 were analysed. Each covers a ∼34 × 34 arcmin 2 area on the sky. The first is centred on NGC 1407 and is made up nearly 5 h of g , r , i imaging under 0.5–0.6 arcsec conditions. The resulting GC catalogue has been used for an extensive spectroscopic campaign with Keck/DEIMOS ( Romanowsky et al. 2009 ; Foster et al. 2010 ; Pota et al., in preparation) and a photometric study ( Forbes et al. 2011 ). The second Suprime-Cam data set is an archival data set in B , V , I bands with 0.6–0.7 arcsec seeing and totalling ∼2 h in exposure time. In this image set, NGC 1407 was positioned toward one end of the mosaic so the radial coverage extends ∼30 arcmin from NGC 1407 or ∼200 kpc. The GC catalogues are available upon request from the first author.
MPGC surface density profiles were constructed in each data set from a catalogue of GC candidates found in the images [MPGCs were taken to have ( g − i ) 0 ≤ 0.98]. A literature Hubble Space Telescope , Advanced Camera for Surveys catalogue ( Forbes et al. 2006 ) was incorporated into the analysis to improve the surface density profile in the central ∼3 arcmin of the galaxy centre. Standard procedures were followed to derive the GC surface density profile (e.g. Spitler et al. 2008a ). The surface density data extend to ∼210 kpc and a constant surface density is not reached at large radii, suggesting a small number of MPGCs should be found beyond the region covered by the Suprime-Cam imaging. Since contamination was not subtracted from the MPGC surface density profile, the fits to the NGC 1407 MPGCs described in subsequent sections include a constant background level to model the contamination.
An expanded GC line-of-sight velocity catalogue from Romanowsky et al. (2009) and Foster et al. (2010) is used, which incorporates new Keck/DEIMOS observations (Pota et al., in preparation). The total number of NGC 1407 MPGCs with spectroscopic velocities is 167 GCs [excluding suspected ultracompact dwarf (UCD) candidates; see details below], out of an estimated total of ∼4000.
3.3 Milky Way
The Milky Way MPGC spatial and kinematic information come from the new Harris (2010) catalogue. The MPGC sample contains 110 GCs with [Fe/H] < −0.9 and has a mean metallicity of [Fe/H] ∼−1.6. The Milky Way GC analysis is carried out in three-dimensions rather than in projection as for the preceding galaxies. In order to model the halo mass profile of the Milky Way with the MPGC velocity information, heliocentric distances and velocities must be corrected to the Galactic reference frame. The distance correction is straightforward and included in the Harris (2010) catalogue itself, but the velocities in principle require proper motion measurements, and these are in many cases not available with sufficient precision. Instead, after correcting the line-of-sight velocities for the heliocentric motion (e.g. equation 5 of Xue et al. 2008 ), the data are corrected to the true radial velocities on a statistical basis: given the anisotropy profile β( r ) from the best-fitting model (discussed in the next section), the correction from line-of-sight to radial velocity dispersion ( equation 1 of Battaglia et al. 2006 ) is applied to the observed velocities, and the dispersion profile is calculated. This processes is repeated iteratively and ultimately yields a very small correction that affects the final ν values at only the ∼3 per cent level (after excluding from the analysis the GC kinematics inside r = 8 kpc, which are more difficult to correct).
The Diemand–Moore technique uses the spatial bias of MPGCs relative to the hosting galaxy’s DM distribution to constrain the νσ properties of the MPGC progenitor subhaloes. This is accomplished by fitting a modified Navarro, Frenk & White (1996 , NFW) function from Diemand et al. (2005) to the MPGC spatial distribution. The fit also requires as input the scale radius ( rs ) of the host galaxy halo, and thus an accurate mass model of each galaxy is an important component of this work. Although models for these galaxies exist in the literature, they are derived using different techniques and various mass tracers. To remove possible systematic errors resulting from these inhomogeneities, new mass models are derived for each galaxy. Furthermore, to ensure the νσ fits to the surface density profile are derived in a self-consistent manner with the galaxy mass models, MPGC observations are used to derive the mass model.
The mass modelling methodology is outlined in the first subsection below. In subsequent subsections, the individual galaxy mass models are derived. Though some necessary assumptions and compromises are made to yield the final mass models, the halo mass properties and their uncertainties are robust. In the final subsection, Section 4.6 , the νσ measurements for each are converted into zreion constraints.
In the following, the MPGC spatial profile is denoted j ( r ) and Σ( R ), for three-dimensional and the projected cases, respectively. Also, when νσ constraints from the Diemand–Moore technique are discussed, ν actually refers to all subhaloes above the subhalo virial mass corresponding to that ν. For example, the constraint ν= 3.1 for the Milky Way, actually includes all subhaloes equal to and more massive than a ν= 3.1 halo at the redshift of interest.
4.1 Derivation of the galaxy mass models
Here we outline our modelling methods that will be used for specific galaxies in the following subsections. We begin with the three-dimensional number density profile of the dynamical tracer population (MPGCs in this paper). This profile is modelled as a modified NFW profile, as derived in Diemand et al. (2005) :Diemand et al. (2005) , we set α= 1, γ= 1.2 and β ν = 3 + 0.26ν 1.6 . Thus a higher ν value implies a steeper outer slope. The value of ν is related to a characteristic radius rν by
Ideally, we would compare this model to an observed surface density profile, and fit for ν. However, in practice there are not good prior constraints on rs available, i.e. we do not know the DM distribution of the halo which we need to estimate the relative compactness of the MPGC subpopulation. 2 To solve this problem, we will make use of the fact that there are kinematical predictions of the model, so we can use the observed spatial and kinematical distributions to fit simultaneously for ν and rs .
Our basic approach for the kinematics is to take model distributions for mass, MPGC density and anisotropy, solve a spherical Jeans equation (e.g. Mamon & Łokas 2005 ) and derive a velocity dispersion profile for the MPGCs, σ( r ). In the Milky Way, the σ( r ) prediction can be directly compared to observational data, while in external galaxies, an additional step is needed to project to line-of-sight dispersions, σ p ( R ).
Each mass model consists of a stellar component (with a reasonable value for the stellar mass-to-light ratio ϒ * ) and a ΛCDM-based model for the DM distribution. This DM profile has a similar parametric form to equation (1) , but with density and scale parameters (ρ s and rs ) that are highly correlated. To derive these correlations (or modelling priors), we begin with the relation between virial mass and concentration predicted for relaxed haloes at z = 0 by the Bolshoi and BigBolshoi simulations using WMAP 5 cosmological parameters ( Prada et al. 2011 , equation 25). After converting the mass and concentration into scale radius and density, and subtracting off a cosmological baryon fraction of 0.17 to find the DM density alone, the following power law provides an approximation to the average results of the numerical simulations:Duffy et al. 2008 ; Macciò, Dutton & van den Bosch 2008 ) translates to a scatter in ρ s at fixed rs of ±0.24 dex. This scatter ultimately provides the largest systematic uncertainty on the measurements of ν and hence zreion .
These results involve fits of standard NFW density profiles to the simulations (i.e. equation 1 with γ= 1, α= 1, β= 3). However, the Diemand et al. (2005) analyses used a slightly different DM density model, where the central cusp has a power-law slope of γ=−1.2 rather than −1:Equation (3) then becomes
The final ingredient in the models is the velocity dispersion anisotropy profile, which is normally a major source of uncertainty in modelling observations of pressure-supported galaxies, but in the context of the Diemand–Moore models, is uniquely constrained. This anisotropy profile is
A complication for modelling both NGC 1407 and M87 is that there are some correlations between GC luminosities and kinematics that are probably driven by populations of UCDs at the bright end, with distinct dynamics from the ‘normal’ GCs (see Romanowsky et al. 2009 ; Brodie et al. 2011 ; Strader et al. 2011 , for more details). Although the UCDs comprise only a tiny part of the overall ‘GC’ systems, they represent a substantial fraction of the luminosity-based spectroscopic samples. Such objects were therefore eliminated from the kinematic samples by excluding anything that has a large measured size (half-light radius >5 pc), or that has no size measurement but is bright ( i0 < 21.3 and i0 < 20 for NGC 1407 and M87, respectively). These objects have no impact on the measured surface density profiles.
The following subsections present the details of the modelling for each galaxy, and the resulting ν constraints.
4.2 Rarity of M87 MPGC progenitor subhaloes
For M87, the mass models are tuned not only to fit the MPGC subpopulation velocity dispersion but also to have a total vc ≃ 500 km s −1 at r ∼ 5 kpc for consistency with the stellar dynamics results of Murphy et al. (2011) ; this means adjusting ϒ * accordingly. As discussed in Strader et al. (2011) , the NFW models have some difficulty in reproducing the stellar dynamics inferences at r ∼ 15 kpc; the cuspier Diemand et al. (2005) models fare better, but there may still be some tension with the stellar data. These mass models are also completely unable to reproduce the X-ray-based mass-model of Das et al. (2010) while maintaining strongly radial anisotropy and reproducing the observed velocity dispersions (an additional example of the problem is footnoted in Section 4.1 ).
Fig. 2 shows the elements of the modelling procedure, along with comparisons to data as applicable. The model consists of profiles of vc ( r ), β( r ) and j ( r ) (the latter is not shown), which are then input to the Jeans equation to solve for σ r ( r ) and σ θ ( r ), and projected to Σ( R ) (lower left-hand panel in Fig. 2 ) and vrms ( R ) (lower right-hand panel). As a consistency check, the model expectations for β( r ) are compared with an empirical estimate from the line-of-sight velocity distribution of the data in the upper right-hand panel of Fig. 2 . This estimate makes use of a restricted class of spherical constant-β models which predict the line-of-sight velocity kurtosis if the velocity dispersion profile is approximately constant with radius ( Napolitano et al. 2009 ). The kurtosis measurements from the MPGC data can thus be converted to rough estimates of the anisotropy profile.
For M87, a typical-concentration ΛCDM halo with rs = 126 kpc ( Mvir ∼ 5 × 10 13 M ⊙ ) is found, along with the Diemand et al. (2005) anisotropy prediction, to fit the M87 MPGC dispersion data remarkably well. The implied β( r ) is also nicely consistent with the empirical estimates. This model corresponds to ν= 2.19. The predicted velocity dispersions are fairly sensitive to the value of rs , increasing for larger rs (and ν). Using a χ 2 fit to the dispersion data, rs is determined at the ±10 kpc level, and ν at the ±0.08 level.
These results hold for the strong theoretical ΛCDM prior adopted on the halo parameters. Relaxing this prior allows for a much larger range of solutions, as various combinations of halo density and scale radius can all provide similar mass profiles in the region with MPGC observations. Allowing for the predicted 1σ scatter in the halo properties (i.e. varying the normalization of equation 6 by ±0.12 dex), rs is constrained at the ∼50 per cent level, yielding ν= 2.19 ± 0.41 (see Fig. 2 ).
4.3 Rarity of NGC 1407 MPGC progenitor subhaloes
For NGC 1407, several priors are invoked: the mass–concentration relation from ΛCDM; an estimate of the mass at R ∼ 100 arcmin from satellite galaxy dynamics ( vc ∼ 700 ± 100 km s −1 ; see Romanowsky et al. 2009 ) and a measurement of the central stellar velocity dispersion around 1 Re of σ p ≃ 260 km s −1 ( Proctor et al. 2009 ). It is beyond the scope of this paper to carry out the detailed dynamical modelling required to convert this σ p measurement into an estimate of vc , and instead the empirical results from dynamical modelling of various other pressure-supported galaxies that vc /σ p ∼ 1.4–1.7 ( Wolf et al. 2010 ; Dutton et al. 2011 ; Murphy et al. 2011 ; Trujillo-Gomez et al. 2011 ) are used.
Unlike M87, for NGC 1407 there is no model solution that straightforwardly fits all of the above constraints (see Fig. 3 ). In particular, it appears to be very difficult to match the inner MPGC velocity dispersion profile, where the σ p ∼ 200 km s −1 inside ∼4 arcmin (∼25 kpc) is lower than the stellar velocity dispersion. Normally, the MPGC dispersion should be higher because the spatial distribution is more extended. Lowering the dispersion would require tangential orbital anisotropy, which contradicts expectations in the Diemand et al. (2005) model for mild radial anisotropy. It is speculated instead that GC disruption processes near the galactic centre have preferentially depleted the more radial orbits (e.g. Fall & Zhang 2001 ; Vesperini et al. 2003 ). MPGCs inside 3.5 arcmin (∼20 kpc) are therefore excluded from the analysis (it is not clear why this should be necessary for NGC 1407 but not for M87).
Considering now the outer MPGC velocity dispersion data, it also turns out to be difficult to fit the joint constraints from ΛCDM halo concentration expectations and the satellite galaxies dynamics. The mass required by the latter predicts MPGC dispersions higher than observed, which can be alleviated by a low-concentration halo. The range of plausible solutions is bracketed by adopting a halo that has a low concentration at the 1σ level while fitting the satellites constraint, or one with a normal concentration that does not agree with the satellites (whose constraint is then assumed to be unreliable). Alternative modifications to lower the MPGC dispersions would be to use a less cuspy NFW model (which is within the cosmological scatter for halo profiles), or to decrease the adopted distance since the kinematics data will then effectively be probing less far into the halo (although this galaxy is likely to be more distant, not less, based on estimates summarized in NED).
The resulting range in fit parameters is ν∼ 2.6–3.1; again, the systematic uncertainties outweigh the statistical ones (∼±0.15) from fitting the dispersion and density profiles. An additional issue from these solutions is that the low observed kurtosis for the outer MPGCs suggests tangential anisotropy in these regions too, again at odds with the radial model expectations.
4.4 Rarity of Milky Way MPGC progenitor subhaloes
The Milky Way modelling differs slightly from the previous two galaxies by treating 3D rather than projected profiles (of density and velocity dispersion). For the mass model, a ΛCDM motivated halo is included as before, along with stellar disc and bulge components based on models ‘I’ and ‘II’ from Binney & Tremaine (2008) , which bracket a plausible range of parameters for these components. We also tried the disc and bulge model of Widrow, Pym & Dubinski (2008) , which is roughly halfway between the I and II models of Binney & Tremaine (2008) . As shown in the upper left-hand panel of Fig. 4 , since these components influence the central regions only, all choices give the same results for ν. The model I of Binney & Tremaine (2008) was adopted.
The modelling results are illustrated in Fig. 4 . After excluding the region inside r = 2 kpc from the density profile fit because of the presence of a core ( Parmentier & Grebel 2005 ; Bica et al. 2006 ), a good fit to both the density and radial velocity dispersion profiles is found. The best-fitting value for the MPGC subhalo rarity is ν= 3.12 ± 0.34, where the uncertainty is again driven by the scatter in the assumed halo mass–concentration relation. The parameters for the DM haloes in these models are listed in the figure caption. The different DM model derived here is the main reason why a higher ν is found compared to the earlier work by Moore et al. (2006) .
Note that there is a gap in the MPGC velocity data around r ∼ 50 kpc. If the MPGC velocity dispersion in this region were assumed to be as low as for the halo stars (as indicated by the dashed curve in the figure, from Gnedin et al. 2010 ) then the virial mass and ν value would need to be lower. The default halo mass profile is actually similar to that of Gnedin et al. (2010) despite the higher velocity dispersion found here. This is because the high typical anisotropy β∼ 0.65 in the current model elevates the radial dispersion relative to the β∼ 0.4 that the former authors assumed, and which would be appropriate for the DM particles on average, not for the ‘biased’ subset of halo particles that should correspond to the stars and MPGCs.
4.5 Mass modelling summary
In the Milky Way and M87, the MPGC kinematics comfortably fit within the Diemand et al. (2005) theoretical framework and are consistent with other mass tracers. Though the mass models unfortunately suffer from rather large systematic uncertainties from the ΛCDM prior, the quality of the models mean the νσ measurements described below are on firm ground. The situation for NGC 1407, from a mass-modelling perspective, is less clear. There appears to be lingering tension between the models and observations, which merit further study.
4.6 Constraining the reionization epoch
To convert the νσ constraints from the preceding sections into estimates for the truncation redshift of MPGC formation (e.g. zreion ), a degeneracy must be broken between subhalo virial mass and redshift for a given ν. The evolving properties of a given νσ subhalo can be analytically computed using the ( Press & Schechter 1974 ) Gaussian field formalism. The value ν is explicitly defined as ν=δ c /[ D ( z )σ( Mvir )], where the constant δ c is the criterion for spherical collapse, D ( z ) is the linear growth factor and σ( Mvir ) is the typical mass-density fluctuation ‘height’ for a subhalo of a given virial mass Mvir (see e.g. Mo & White 2002 ). For the present work, a top hat filter is used to select the characteristic spatial scale. The camb ( Lewis, Challinor & Lasenby 2000 ) online interface 3 was used to generate the linear power spectrum for a WMAP 7 cosmology ( Komatsu et al. 2011 ). Fig. 5 shows the relationship between the rarity, ν, of a mass-density fluctuation for a given Mvir –redshift pair.
Fig. 5 also shows how a νσ measurement alone will not constrain the redshift when MPGC formation is suppressed. Further information is required to break this degeneracy. Following Moore et al. (2006) , it is assumed that star and MPGC formation can only occur in subhaloes that have virial temperatures greater than Tvir = 10 4 K, where gas is able to rapidly cool and form stars (e.g. Ostriker & Gnedin 1996 ). The region in Fig. 5 below the Tvir = 10 4 K curve can therefore be excluded. The region above the Tvir = 10 4 K curve can be excluded as well, since MPGC formation would have taken place in subhaloes that would ultimately have less or more concentrated spatial distributions than what is observed in the galaxies. For example, for the Milky Way constraint of ν= 3.1, at redshift of z = 10, subhaloes with mass Mvir ∼ 10 8 M ⊙ could still have formed MPGCs efficiently, which would have led to a spatial distribution resembling 2.5σ subhaloes, contradicting the observations. At redshifts higher than the intersection between the Tvir = 10 4 K relation and the νσ curve from MPGC observations, only rarer, e.g. ν > 3.1, subhaloes could form MPGCs, thus the spatial distribution would have been more concentrated than is observed. The intersection between the νσ constraint and the Tvir = 10 4 K curve is therefore taken to be the MPGC truncation redshift or the local zreion of the host galaxy.
The MPGC νσ fits and the corresponding reionization redshifts are given in Table 1 for each galaxy. The weighted mean of the three ν MPGC measurements translates into a reionization redshift of . These constraints are discussed in Section 5 .
|Environmental density||ν||z reion||Min. virial mass ( h−1 M ⊙ )|
|Messier 87||High/cluster||2.19 ± 0.41||1 × 10 8|
|NGC 1407||Medium/group||2.86 ± 0.36||8 × 10 7|
|Milky Way||Low/field||3.12 ± 0.34||7 × 10 7|
|Mean||2.74 ± 0.21||8 × 10 7|
|Environmental density||ν||z reion||Min. virial mass ( h−1 M ⊙ )|
|Messier 87||High/cluster||2.19 ± 0.41||1 × 10 8|
|NGC 1407||Medium/group||2.86 ± 0.36||8 × 10 7|
|Milky Way||Low/field||3.12 ± 0.34||7 × 10 7|
|Mean||2.74 ± 0.21||8 × 10 7|
The reionization redshift ( zreion ) constraints from MPGC observations (see Fig. 6 ) are summarized in Fig. 7 with existing estimates and limits from the literature. The local constraints from MPGC observations are tabulated in Table 1 .
The only existing constraints on zreion from MPGCs are for the Milky Way. As described in Section 1 , our estimate agrees with that derived in Moore et al. (2006) , who used the same technique employed here. Also, the high-resolution DM Milky Way simulation described in Griffen et al. (2010) can only reproduce the observed spatial distribution of MPGCs and their numbers only when MPGC formation was suppressed by reionization at zreion ∼ 13. This value agrees with our own estimate within its 1σ uncertainty.
5.1 Environmental propagation of reionization
The MPGC zreion values and their uncertainties span the redshift range covered by zreion estimates from the literature, with the weighted mean MPGC constraint showing best agreement with the cosmic microwave background (CMB) value and an instantaneous reionization model ( Komatsu et al. 2011 ; Larson et al. 2011 ). It is apparent that the M87 zreion measurement is somewhat of an outlier compared to the other MPGC estimates. Since M87 is located in the high-density environment of a galaxy cluster, it is possible that its lower zreion results from an environmentally dependent reionization epoch. Unfortunately, the relatively large systematic uncertainties on the galaxy halo mass models mean the difference between zreion in the highest and lowest density environments (M87 versus Milky Way zreion ) is only significant at the 1.7σ level.
In addition to a correlation between the MPGC spatial distribution and the local reionization epoch, Moore et al. (2006) predicted that zreion should also correlate with the virial mass-normalized MPGC numbers. This is because galaxies with a lower zreion had more time to produce MPGCs. Both Spitler et al. (2008a) and Spitler & Forbes (2009) found no evidence for a correlation between MPGC virial mass specific frequencies and environment. However, their virial mass estimates were derived from statistical relationships between stellar and virial masses, which are not ideal for use on individual galaxies.
Using the galaxy virial masses derived in Section 4 and MPGC numbers ( NMPGC ) from Section 3 , the following MPGC virial mass specific frequencies ( VMPGC ; Spitler et al. 2008a ) are found: for M87, NGC 1407, and the Milky Way, respectively. The errors are again dominated by systematic uncertainties in the galaxy mass models. M87 has a VMPGC value that is larger than the Milky Way or NGC 1407, at the 2.5σ significance level. The normalized MPGC numbers thus lend further support to the idea that reionization finished in M87 at lower redshifts compared to the other galaxies. In general, massive galaxies at the centres of galaxy clusters (e.g. M87) are known to have enhanced relative number of GCs when normalized to the galaxy’s star light (e.g. Harris 1991 ; McLaughlin 1999a ). Whether this holds when the MPGC numbers are instead normalized to virial masses (which should not be directly influenced by reionization) will need to be explored in a future work.
The above results provide an exciting hint that reionization completed first in low-density environments.
5.2 The theoretical context
A number of large-scale radiative transfer simulations have been run to understand the propagation of reionization through different galaxy environments. From analysis of these simulations, it is found that protogalaxies in high-density environments produce significant numbers of ionizing sources early on and thus are the first locations to reionize ( Iliev et al. 2006 , 2011 ; Finlator et al. 2009 , however, see Weinmann et al. 2007 ).
Such simulations typically do not have a large enough dynamical range to capture both the environmental propagation of reionization over large volumes and the small-scale physics that are important for understanding localized affects. For example, when an inhomogeneous intergalactic medium is incorporated into reionization models, dense pockets of neutral gas can survive even after surrounding regions had completed the reionization process ( Miralda-Escudé, Haehnelt & Rees 2000 ; Ciardi, Stoehr & White 2003 ; Furlanetto & Oh 2005 ; Iliev et al. 2008 ; Choudhury, Haehnelt & Regan 2009 ; Mitra et al. 2011 , 2012 ). These pockets tend to be located in dense environments, where recombination rates can stay high despite their proximity to large reservoirs of ionizing sources. This means the intergalactic medium in dense environments started to reionize first, but only completely reionized after lower density environments had finished.
The trend hinted at by the MPGC observations agrees with this latter theoretical work. Under this scenario, MPGCs in a high-density environment continued forming to lower redshifts in self-shielded pockets, even as low-environmental-density regions were being completely ionized. The progenitor molecular clouds of MPGCs located in more vulnerable environments (like that surrounding the Milky Way) were not surrounded by gas with sufficiently high recombination rates to survive the ionizing front that originated within or passed through the region.
Finally, the range of the MPGC constraints on zreion agrees with the range of literature constraints (see Fig. 7 ). This seems to support models of reionization that are globally extended in duration (e.g. Mitra et al. 2011 , 2012 ).
The spatial and kinematic properties of MPGCs provide an exciting hint that reionization was inhomogeneous. In the local Universe, reionization completed first in low-density environments around zreion ∼ 11–12 (for the Milky Way and large central group elliptical NGC 1407) and finished last in the high-density environment around the Virgo cluster at zreion ∼ 8. While uncertainties on the galaxy halo mass models limit the strength of this conclusion, another property, the relative number of MPGCs, supports this interpretation. The apparent environmental dependence of zreion also aligns with theoretical work that finds that the ionization front propagated from high- to low-density environments, with the high densities finishing reionization last. Furthermore, the apparent agreement between the range of zreion constraints from MPGC observations and those from the literature may provide additional support to models where reionization was globally prolonged in duration.
More observational work, especially to improve the individual galaxy halo mass models, is required to confirm the results found here. It will also be important to model MPGC formation in large-scale, radiative transfer simulations with realistic models of star cluster dynamical evolution. This will help refine the technique utilized here and clarify the relationship between reionization and the truncation of MPGC formation. For example, given that high star formation rates and gas densities were needed to produce such dense star clusters, it is possible that MPGCs trace the end of reionization, where ionizing radiation was able to finally penetrate the densest gas clouds.
Expanding this initial study to other galaxies will yield a reionization map of the local Universe, further aiding our understanding of galaxy formation and evolution over cosmic time.
If reionization is not responsible for the truncation of MPGC formation, the results presented here provide unique constraints on the formation histories of MPGCs. Furthermore, since MPGCs and stellar haloes of the Milky Way share spatial and chemical properties ( Helmi 2008 ; Martell et al. 2011 , see also M31, Huxor et al. 2011 ) they might share a common origin. If this is the case for other galaxies, then MPGC observations can be used to constrain the star formation history of their stellar haloes.
Finally, the power of the Diemand–Moore technique is that it does not care about the detailed star formation and enrichment physics – it only requires an observable tracer object that accreted on to the galaxy at early times ( Diemand et al. 2005 ).
We want to acknowledge the useful comments provided by the anonymous reviewer. We also thank Chris Blake for assistance with the cosmological derivations and Anna Sippel for help with the initial Subaru reductions. LRS was supported by the ARC Discovery Programme grants DP0770233 and DP1094370. This work was supported by the National Science Foundation through grants AST-0808099, AST-0909237 and AST-1109878. DAF thanks the ARC Discovery Programme for support. This paper was based in part on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. Some of the Subaru data were acquired with the time-swap Gemini programme GN-2006B-C-18. Some of the data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The observatory was made possible by the generous financial support of the W.M. Keck Foundation. Observations reported here were obtained at the MMT Observatory, a joint facility of the Smithsonian Institution and the University of Arizona.