GLOBULAR CLUSTER POPULATIONS: FIRST RESULTS FROM S4G EARLY-TYPE GALAXIES

Dennis Zaritsky; Manuel Aravena; E. Athanassoula; Albert Bosma; Sébastien Comerón; Bruce G. Elmegreen; Santiago Erroz-Ferrer; Dimitri A. Gadotti; Joannah L. Hinz; Luis C. Ho; Benne Holwerda; Johan H. Knapen; Jarkko Laine; Eija Laurikainen; Juan Carlos Muñoz-Mateos; Heikki Salo; Kartik Sheth

doi:10.1088/0004-637X/799/2/159

1. INTRODUCTION

There are a few well-established empirical findings that help guide our developing understanding of galaxy formation and evolution. Among those are some that hint at the importance of environmental processes (Dressler 1980; Postman & Geller 1984), often referred to as the role of nurture in galaxy evolution, and others that highlight the role of halo mass (Kauffmann et al. 2003; Bundy et al. 2006), often referred to as the role of nature. Although numerous studies have attempted to establish the supremacy of one set of influences over another, the two are intertwined due to the relationship in hierarchical growth models between environment and mass (De Lucia et al. 2012).

At the core of the hierarchical growth paradigm is the agglomeration of mass in galaxies by mergers or accretion. Unfortunately, observational markers of the most significant events, such as tidal tails or bridges, that would evince a recent major interaction or merger are generally ambiguous and difficult to establish for large samples (see Tal et al. 2009; Lotz et al. 2011; Adams et al. 2012 for examples), whereas minor events are nearly imperceptible.

What happens to a galaxy as it merges, accretes smaller neighbors, and grows? The scrambling of stars, gas, and dark matter that occurs during such events erases traces of the progenitors and of the processes that occurred prior to and during the event. Much of our intuition is instead guided by simulations, and those are in critical need of empirical confirmation. Although there are specific examples where we catch a galaxy relatively soon after a burst of star formation—E+A or K+A galaxies (Dressler & Gunn 1983; Couch & Sharples 1987)—and where convincing evidence of a merger exists (Zabludoff et al. 1996; Norton et al. 2001; Yang et al. 2004), these are rare objects in the local universe (Zabludoff et al. 1996). For the more common, less dramatic accretion/merger events there is little in the way of tracers of the hierarchical accretion that is such an integral part of our understanding of the growth of structure. We need to identify a population of objects within galaxies that can testify to the combination of different progenitors and to any new star formation that occurs as part of the growth process.

An understanding of globular cluster populations in galaxies has developed over the past two decades that suggests that clusters could be such a tracer and provide new insights and constraints for modelers. There is empirical evidence that all but the least massive galaxies have multiple populations of globular clusters, where metallicity is the distinguishing characteristic among the populations (for a review see Brodie & Strader 2006), although alternative interpretations exist (Richtler 2006; Yoon et al. 2006), and the observational evidence is more complex than previously appreciated (Chies-Santos et al. 2012). The developing interpretive consensus is that these populations reflect both an early epoch of cluster formation and subsequent evolution influenced by major dissipational mergers (Harris & van den Bergh 1981; Schweizer 1987; Ashman & Zepf 1992), dissipationless accretion (Côté et al. 1998), and ongoing cluster formation (Kravtsov & Gnedin 2005). With increasingly sophisticated spectroscopic surveys (such as that of Brodie et al. 2014), the connection between cluster and stellar populations will be tested and used to constrain formation models.

Despite the promise of this field of study, compiling large, homogeneous samples has been challenging. At a minimum, we need to understand whether the number of clusters per galaxy, normalized in some sensible way, varies with galaxy mass, morphology, environment, or stellar population. There is a long history of studies of the cluster number normalized by galaxy luminosity, a quantity defined as the specific frequency of clusters (S_N; see Harris 1991; Brodie & Strader 2006 for reviews), and a more recent focus on cluster number normalized by galaxy stellar mass (T_N; Zepf & Ashman 1993), but uncertainties due to contamination and completeness corrections were often large, and galaxy samples were limited. The state-of-the-art compilation is that presented by Harris et al. (2013), who scoured the literature to obtain such estimates for 422 galaxies, of which 341 are the early types that are the focus here. Despite the care taken in that work, the one potential weakness of literature compilations is the unavoidable heterogeneity in sample selection, image quality and characteristics, image analysis, and cluster population modeling necessary for completeness corrections. This heterogeneity is particularly worrisome for identifying trends across a broad range of galaxy parameters because individual studies have tended to focus on a particular class of galaxies (for example, primarily galaxies in the Virgo cluster or of low mass). As we describe below, there are a number of decisions involved in undertaking a cluster census and, in a relative sense, one is on surer footing if the data quality and analysis are the same across the entire sample. Therefore, a study such as ours at the very least constitutes a complementary approach with which to address these important questions.

Over the last 10 yr or so, terrific progress has been made in obtaining high-fidelity samples of clusters on which to base the cluster census, using both high-angular-resolution images, particularly those provided by the Hubble Space Telescope (HST), and color information from deep photometry to remove contaminants (for examples, see Peng et al. 2006; Strader et al. 2006; Kundu & Whitmore 2001a, 2001b; Rhode & Zepf 2004; Young et al. 2012). However, with a few exceptions (Peng et al. 2006; Strader et al. 2006), these studies only cover a handful of galaxies because of the observational cost of obtaining such data. Those studies that have both high angular resolution and are photometrically sensitive tend to suffer from small fields of view that do not necessarily cover the full radial extent of the globular cluster population, creating yet another source of uncertainty in the final cluster counts (Rhode & Zepf 2003), and generally cover galaxies in a single environment (such as the Virgo cluster; Peng et al. 2006; Strader et al. 2006).

We have chosen to measure T_N for a large sample of galaxies in a manner that is more reminiscent of earlier treatments. The recent emphasis has been on greater and greater precision in the removal of contaminants. That focus has been in large part driven by the need for pure samples on which to base subsequent photometric (Gebhardt & Kissler-Patig 1999; Larsen et al. 2001; Spitler et al. 2008a) and spectroscopic studies (Brodie et al. 2014). By relaxing the criteria on sample purity, because we are not seeking to define samples for follow-up studies, we can accept larger uncertainties in the measurement of the number of clusters, N_CL, in exchange for larger galaxy samples that span a greater range of galaxy properties. If, as expected, N_CL varies by several orders of magnitude among galaxies, then uncertainties as large as a factor of a few may have little effect on broad trends. Ultimately, the final answer to whether such lower precision is scientifically useful depends on the magnitude of the effects present and the size of the sample. A basic weakness of our approach is that we treat the cluster population as a single entity, despite clear evidence from colors (Zepf & Ashman 1993; Ostrov et al. 1993; Gebhardt & Kissler-Patig 1999; Larsen et al. 2001; Kundu & Whitmore 2001a) and kinematics (Strader et al. 2011; Woodley et al. 2010; Pota et al. 2013) that there are multiple populations. Of course, a similar criticism can be levied on studies of stellar populations, and yet basic relationships are still valid for the whole.

We proceed along these lines in an attempt to eventually obtain the largest homogenous sample of cluster population measurements with which to complement the more intensive, focused work that is ongoing on smaller samples and the broader work being carried out with compilations of published results. Our measurement of the cluster populations is based on a statistical excess of point sources in the S⁴G images (Sheth et al. 2010) of nearby galaxies. We will establish that the methodology presented here is sufficiently accurate and precise to be scientifically useful. The S⁴G data have several direct advantages because imaging at IR wavelengths suffers less extinction than that at optical wavelengths, and the IR luminosity is closely tied to older stellar populations. In addition, the S⁴G data is a roughly volume-limited sample of several thousand galaxies. Indirect advantages accrue from ancillary studies that include photometric decomposition (Kim et al 2014; Muñoz-Mateos et al. 2014; H. Salo et al., in preparation), detailed and homogeneous morphological classification (Elmegreen et al. 2011; Holwerda et al. 2014; Buta et al. 2015), radial and vertical disk structures (Laine et al. 2014; Comerón et al. 2011), classification of asymmetric structures (Zaritsky et al. 2013b; Laine et al. 2014), and stellar and gaseous mass estimates (Zaritsky et al. 2014; Querejeta et al. 2014). In Section 2 we describe the sample, how we constructed the cluster candidate catalog, and how we measured the number of clusters and their radial distributions. We discuss the findings in Section 3 and conclude in Section 4.

2. DATA AND MEASUREMENTS

2.1. Constructing the Point Source Catalog and Surface Densities

The parent sample for this study is the S⁴G sample, which currently consists of 2352 galaxies selected as described by Sheth et al. (2010). Although primarily a volume-limited sample, some additional selection criteria, such as the existence of an HI redshift, and the surface brightness limitations of the existing catalogs from which the sample was selected preclude it from being a complete, volume-limited sample. We observed these galaxies using the Spitzer Space Telescope (Werner et al. 2004) and its Infrared Array Camera (Fazio et al. 2004), as described by Sheth et al. (2010). The data from the original S⁴G survey are publicly available through the NASA IRSA website.¹⁵ An extension of S⁴G to cover the missing HI-weak galaxies in this volume that otherwise match our selection criteria was approved for Cycle 10 and is expanded to increase several fold the number of galaxies relevant for this study, which is why we consider the results presented here to be the first of the set of results to come from S⁴G regarding globular cluster populations in early-type galaxies.

From among the galaxies in the completed original sample (DR1), we have so far limited our study. We focus on galaxies with morphologies identified as early-type (−5 ⩽ T-type ⩽1) by Buta et al. (2015) to avoid those galaxies with significant internal structure that could be mistaken for point sources. All but five in our sample have T ⩽ −2. This morphological cut results in the sample for this study being a small fraction of the full S⁴G sample. We further constrain the sample by including only those galaxies within a suitable range of distances. We set the upper end of the distance range to correspond to a distance modulus (DM) of 32.4 (30.1 Mpc), which we found to be the upper limit at which we have enough physical resolution within the central galaxy to result in a reliable globular cluster population radial profile (see below for more discussion of radial profile limits). We set the lower end of the DM range at 30.25 (11.2 Mpc), which we found to be necessary to ensure sufficient background coverage within the images with which the background source density is constrained. We use redshift-independent distances when available in the NASA/IPAC Extragalactic Database (NED); otherwise we use the redshift and adopt H₀ = 72 km s⁻¹ Mpc⁻¹ to derive a Hubble velocity distance estimate.

For this range of distances, the Spitzer images are of sufficiently large angular size to reach well into the background source population—as defined by a lack of a detectable point source surface-density gradient—but they often contain a significant number of bright stars and, occasionally, other large galaxies. To be specific, S⁴G images cover a projected scale of at least 1.5D₂₅, where D₂₅ is the diameter of the galaxy's isophote at the surface brightness level of 25 mag arcsec⁻², except for a subset of 125 archival galaxies. In terms of physical separations, we probe to galactocentric-projected radii >30 kpc in all but one case and, typically, out to between 50 and 100 kpc. Although there are certainly some clusters beyond these radii in many of our galaxies, we will show in Section 2.3 that any uncertainties introduced by the lack of coverage at larger radii are subdominant to other uncertainties. All bright objects are masked in the critical parts of the analysis using the masks developed as part of the S⁴G processing (H. Salo et al., in preparation; Muñoz-Mateos et al. 2014). We remove from our sample those galaxies where another comparably large galaxy is included within the 3.6 μm frame, both because of the increased difficulty in doing the measurement but also because interpreting the distribution of candidate globular clusters is more complicated. In a few other cases, we exclude the galaxy for technical reasons. Our final sample of 97 galaxies is given in Table 1.

Table 1. Globular Cluster Population Properties

Name	DM	T	M_3.6	M_4.5	N₅₀	T_N	Back.	a	b	Q
ESO 357-025	31.13	−2	14.9	15.4	20	6.59 $_{-4.86}^{+6.16}$	−5.38	12.41	−5.84	1
ESO 419-013	31.45	−2	14.5	14.8	98	16.34 $_{-11.32}^{+11.14}$	−5.00	3.91	−2.98	1
ESO 548-023	31.71	−3	14.4	14.8	111	13.44 $_{-12.49}^{+15.01}$	−5.35	14.96	−5.25	0
IC 335	31.19	−3	11.9	12.4	2	0.04 $_{-0.01}^{+0.86}$	−5.49	0.01	−1.45	0
IC 796	32.06	−3	13.0	13.5	153	3.67 $_{-1.24}^{+1.75}$	−4.80	11.71	−5.20	1

Notes. DM refers to the distance modulus. T is the morphological T-Type of the galaxy from the compilation of Buta et al. (2015). The 3.6 and 4.6 μm magnitudes of the galaxies as measured by Muñoz-Mateos et al. (2014) are in Columns 4 and 5. N₅₀ is the number of globular clusters estimated from our best fit model of fixed power law slope within 50 kpc of the galaxy. T_N is the number of these clusters per 10⁹ M_☉ of stellar mass in the galaxy. The quantity Back. represents the log of the measured background level (sources/parsec), and a and b represent the normalization and slope of the power law fit for models where the power law is allowed to float and the background is set to Back.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as: Data Typeset image

Globular cluster candidates at these distances, in images of this angular resolution, will appear as point sources. It is therefore not possible based solely on morphology to separate foreground stars, galaxies of small angular size, and globular clusters. Our measurement will be statistical in nature, based on the expectation that contaminants will not cluster about the galaxy. Although colors, in principle, can help in distinguishing between these classes, in practice (see below) neither the other available Spitzer band at 4.5 μm, nor the available Sloan Digital Sky Survey (SDSS) optical imaging, is of use here. Therefore, we are unable to identify any specific individual sources as likely clusters but only provide an estimate of the excess of 3.6 μm sources correlated with each galaxy and the radial distribution of this excess. We will argue that this excess is due to the globular cluster population and, where possible, validate this claim with comparison to published measurements. However, these characteristics preclude the use of these data, if not complemented with other bands, in studies of cluster subpopulations or for spectroscopic target selection.

Before cataloging all point sources with apparent magnitudes that are consistent with those of clusters at the corresponding distance, we process the images, to aid us in identifying such objects. In addition to applying the masks mentioned previously (masked regions are not considered in the subsequent analysis other than in corrections for completeness), we use the exposure weight maps to exclude areas with substantially less exposure time and, therefore, lower sensitivity. The exact value of the thresholding we use varies for each image but is selected to exclude the image edges. Problems with detections at the image edges are often noticeable as sharp rises or dips in the final radial-density profiles of sources and occur either at image gaps, for cases where multiple images are used to cover the field around a galaxy, or at the largest radii. We guide our selection of the threshold value by adopting the smallest threshold value that eliminates such features.

The basic preprocessing steps include sky subtraction plus modeling and removal of the primary galaxy. We calculate the background sky value by evaluating the median within either the upper or lower quarter of the image, depending on whether the primary galaxy lies in one or the other of these two regions, of unmasked pixels. We subtract this median sky value from the entire image. We then use the IRAF task ELLIPSE to measure the properties of the central galaxy, create an image of that model using BMODEL, and then, by subtraction, obtain an image that is as nearly free of the primary galaxy as possible (see Figures 1 and 2). We choose to use the ELLIPSE task, with free ellipse parameters, rather than the more detailed model-fitting of S⁴G galaxies (H. Salo et al., in preparation) using GALFIT (Peng et al. 2002) or BUDDA (de Souza, et al. 2004) because the ELLIPSE fitter proved to be superior at removing the smooth galaxy light. The GALFIT or BUDDA fits provide a more physically motivated approach as well as estimates of the parameters of physically distinct components in these galaxies (e.g., bulge and disk) and thereby enable one to address a wide range of questions, but for our purpose, which is to remove as much of the galaxy light as possible, the freedom of the model-unconstrained ellipse-fitting results in smaller residuals in the subtracted image.

**Figure 1.** Demonstration of the model subtraction for NGC 1553. We selected this galaxy as an example of a fairly luminous, extended galaxy in our sample. The results are otherwise not unusual. This is the S⁴G produced mosaic in the 3.6 μm passband. The size of the images is 245 across, and north is up. The box in the right panel marks the area displayed in Figure 2.
Download figure:
Standard image High-resolution image

farcm — **Figure 1.** Demonstration of the model subtraction for NGC 1553. We selected this galaxy as an example of a fairly luminous, extended galaxy in our sample. The results are otherwise not unusual. This is the S⁴G produced mosaic in the 3.6 μm passband. The size of the images is 245 across, and north is up. The box in the right panel marks the area displayed in Figure 2.
Download figure:
Standard image High-resolution image

**Figure 2.** Demonstration of the model subtraction near the core of NGC 1553. The quality of the model subtraction is highly variable within a few kiloparsecs of the galaxy center.
Download figure:
Standard image High-resolution image

Examples of the galaxy subtraction both on the scale of the full image and expanded about the target galaxy are shown in Figures 1 and 2, respectively, for NGC 1553. This galaxy is among the larger (in angular extent) galaxies in our sample. The multipolar residual pattern seen in Figure 2 is typical, but stronger in galaxies with disks, and usually spans ∼1 kpc in projected radius from the center. Within the central region, there are limited areas where individual sources can be identified quite close to the center, as well as others where the modeling errors severely increase the local noise. Beyond this inner region, the modeling clearly works quite well. Our procedure for measuring completeness will account for this variation in detection sensitivity.

Once the residual images are available, we run SExtractor (Bertin & Arnouts 1996) to identify point sources, eventually using the stellarity index to reject extended sources. When running SExtractor, one defines the criteria for an acceptable source using a specified minimum required number of pixels detected above a specified threshold, where that threshold is defined in terms of nσ above background, where σ is the background rms. We found that for a number of our images, SExtractor was miscalculating σ because of the odd shape of the image footprints within the overall rectangular "image." We therefore define a flux level that we consider significant and then evaluate the corresponding σ threshold to reach that flux (over a minimum of two pixels). We evaluate this threshold interactively to reach the level of detection shown in Figure 3, guided by the criterion that detections be visually robust sources. We eventually remove most objects near this detection threshold from our catalogs when we set a uniform absolute magnitude limit, but that step is discussed below. Cataloged objects that are clearly spurious, which tend only to be found near the galaxy center, are generally extended and so removed on the basis of that criterion, but we also use the 4.5 μm images to help remove those, as described below.

**Figure 3.** Demonstration of SExtractor object detection near the core of NGC 1553. We detect point sources well within the galaxy. Some spurious sources are also created by the poorer model subtraction at projected separations <3 kpc; a radius of 2 kpc is represented by the large central circle. The spurious sources are subsequently rejected either because they are identified as extended sources or not matched in the 4.5 μm photometry (see text for details). Beyond the innermost central region, the source detection is of similar quality at all radii.
Download figure:
Standard image High-resolution image

We opt to remove only unambiguously extended sources and accept all detections with stellarity >0.1. There are various ways to select between unresolved and extended sources, guided by the concentration or surface brightness of the image. However, because of the limitations of these data (2'' FWHM point-spread function; PSF), we are not in a position to resolve between faint background galaxies and point sources and therefore use only a basic measure of morphology and err on the side of inclusion.

Images of these galaxies at 4.5 μm are also part of S⁴G and therefore exist for all of our targets. However, due to how observations were structured, these images are typically rotated 180° on the sky. This rotation means that although there is overlap between the 3.6 and 4.5 μm imaging at the location of the target galaxy, the outlying regions, which are critical for determining the background source density, do not have overlapping coverage. For this reason, we cannot use the 3.6–4.5 color as a selection criteria (although it is known to be a weak diagnostic for star clusters in any case; Spitler et al. 2008a). However, in the region near the central galaxy, where our model subtraction is more uncertain, we use this additional band to help us discriminate against spurious sources. We run SExtractor in two-image mode, using the 3.6 μm image as the reference to photometer sources in the 4.5 μm image that lie within 100 pixels of the central galaxy. This is the critical region where the model subtraction is most variable. Because the photometry can have large uncertainties in this region, and because we do not want to select against the real source (even if they are not clusters—because we cannot enforce the same selection throughout the full 3.6 μm image), we accept sources within a large color range |3.6–4.5| < 2.5 mag. We also apply the same procedure using a 5 mag color cut but see no appreciable difference in the results (see below for more details).

Optical photometry would aid us either in confirming sources as real or, better still, in differentiating star clusters from other sources. For example, optical–IR colors have been used to discriminate clusters from contaminants and to make further measurement of cluster metallicities (Kissler-Patig et al. 2002; Spitler et al. 2008a), and so combining our images with SDSS images, which are available for a large portion of the sky, is a natural avenue forward. However, the SDSS images turn out to be insufficiently deep. In Figure 4 we show the galaxy-subtracted residual images for NGC 584. It is evident from this comparison that the 3.6 μm images go much deeper than the optical images. Of course, this difference would not be relevant if either the SDSS image was sufficiently deep to detect the clusters or neither image was sufficiently deep. However, we show in Figure 5 that the apparent magnitude distribution of the SDSS data is grossly incomplete at the relevant magnitudes and that the S⁴G sample is well-matched to reach the top few magnitudes of the globular cluster luminosity function (LF). We conclude that we are unable to use the SDSS data to help in our selection and that completeness corrections will not be extreme for the S⁴G cluster counts. A deep set of optical images that cover the footprints of the S⁴G galaxies will provide value in revisiting this question and enable one to examine the subpopulations of clusters in these galaxies.

**Figure 4.** Comparison of SDSS (left panel) and S⁴G (right panel) residual (model-subtracted) images for NGC 584. The point source density is manifestly much higher in the S⁴G images. Here we select a different galaxy than in the previous plots both because an SDSS image of NGC 1553 does not exist and to provide a second, different example of a residual image. Sources can be detected close to the nucleus along the minor axis but are lost in the residuals along the major axis. The completeness simulations will correct for the variable sensitivity using our simulations, under the assumption that the clusters are spherically distributed. A second image, offset to the northwest, exists but is not shown here and provides coverage to much larger radii. North is up, and the angular size of the image portion shown is roughly 5 arcmin.
Download figure:
Standard image High-resolution image

**Figure 5.** Histogram of the apparent magnitude distribution for the combined source distribution of NGC 584 and NGC 1047. We used two galaxies to increase the statistics and ensure the result was not related to the selection of an inappropriate galaxy. Otherwise, these are just the first two NGC galaxies from our sample in a numerically sorted list for which SDSS data are available. The red unshaded histogram shows the S⁴G source distribution, whereas the blue shaded histogram shows the matched SDSS sources. The vertical dotted line and arrow indicate the apparent magnitudes at which we expect to find globular clusters, set at a value that is 2σ brighter than the peak magnitude of the globular cluster luminosity function (LF), and σ is the dispersion adopted for the Gaussian LF (see Section 2.2 for discussion of the cluster LF). For the full range of distance moduli in our sample, the vertical line will range ±1 mag from the plotted position. The SDSS images are grossly incomplete at the relevant magnitudes.
Download figure:
Standard image High-resolution image

Given the lack of additional data to aid in selecting clusters, we implement a basic 3.6 μm magnitude cut to exclude sources that are clearly too bright to be clusters at the distance of the target galaxy or that are sufficiently faint that we are beginning to reach within the highly incomplete range of our data. Guided by these limits, we constrain ourselves to sources with −11 < M_3.6 < −8, the lower limit arising from defining a cut that is above the cluster detection limit for all of our galaxies.

Constructing the point source catalog is only the first necessary part of our measurement. As a function of both magnitude and location within the image, the completeness of our source counts will vary. To determine completeness, we add artificial point sources over a range of magnitudes that is greater than that defined for the candidate clusters (because measurement uncertainties could move objects within our magnitude limits). We cannot add too many sources without affecting the statistical properties of the image (where artificial sources start to overlap other artificial sources), so we rerun the process 50 times. To obtain better sampling of the regions near the galaxy, we distribute the artificial sources using a Hubble profile with a core radius of 100 pixels. We do assume the clusters are spherically distributed about the galaxy, but this is again an area where guidance from existing work is lacking. Some studies even suggest that the distribution varies for the different cluster populations (Wang et al. 2013), complicating this issue even further. Within apparent magnitude, we distribute the sources uniformly. The images with the artificial sources are then processed and analyzed in the same manner as the original ones. For each radial bin, we bin by 25 pixels and evaluate both the number of real sources and, for each of the 50 realizations, the fraction of artificial sources that is recovered (for sources within the absolute magnitude range of −11 to −8, assuming all sources are at the distance of the primary galaxy). This procedure results in enough sources, even at the smallest radii, that the uncertainties in our incompleteness corrections are subdominant over counting statistics. We then correct the observed number counts and, when converting these numbers to a surface density, we correct for the missed pixels beyond the image boundary by using only the area of the annulus within the image. The correction for masked pixels within the image boundary comes from our artificial source recovery fractions.

2.2. Measuring the Cluster Populations

Using the radially binned, completeness-corrected surface-source-density values, we now estimate the parameters of a power law profile description of the cluster distribution. The data are insufficient to allow for the fitting of the power law and background simultaneously. We adopt two approaches. First, we measure the background projected density level using the surface-density values at radii > 30 kpc and then fit a power law plus background model, with a free power law slope and normalization, for radii between 1 and 15 kpc. Second, we fix the power law slope (the specific value chosen will be discussed below) and vary the power law normalization and background level, still fitting the model for radii between 1 and 15 kpc. The fits are done by minimizing χ². The choice of one model over another hinges conceptually on whether the background at large radii is a sufficiently accurate description of the background at smaller radii. We will select our preferred approach by comparing our results to previously published results for the limited subsample of galaxies where such measurements are available.

The power law description for the radial distribution of globular clusters is historical (see Brodie & Strader 2006) but fails to describe in detail known systems when very large samples of clusters define the radial profile precisely (for examples of truncations, see Rhode & Zepf 2001; Dirsch et al. 2005). Nevertheless, for the quality level of our data and in the interest of homogeneity in analysis, the power law is adequate and preferred. As shown in Figures 6 and 7, the power law model does a satisfactory job of fitting almost all of the galaxies.

Figure 6. — Download figure:
Standard image High-resolution image

Figure 7. — Download figure:
Standard image High-resolution image

The resulting values for the background (logarithm of the counts per bin; Back.), power law slope (b), and normalization (a) for the models where the background is set and the power law index is free are given in Table 1, along with the DM, Hubble T-Type (T), the 3.6 and 4.5 μm apparent magnitudes, and the quality flag (Q). The integrated number of clusters out to a radius of 50 kpc (N₅₀) and the specific frequency relative to the galaxy's stellar mass (T_N), corresponding to N₅₀ in units of number per 10⁹ M_☉ (as introduced by Zepf & Ashman 1993), is given for the models where the power law index is fixed. The data and fits, shown including the background, are given in Figures 6 and 7 for all 97 galaxies, for fixed background and fixed power law index, respectively.

Uncertainties in T_N are based on Poisson statistics in the individual radial bins, propagated through the fitting using Δχ². In cases where the model fit is statistically acceptable we adopt the uncertainties corresponding to 1σ in the model parameters. In cases where the model is statistically unacceptable, for the adopted Poisson uncertainties in the individual bins, we calculate the χ² for which we would have only a 50% chance of rejecting the model fit and rescale all of the uncertainties by the required amount to reach this χ². The justification for this rescaling is that there are systematic uncertainties in the background that are not captured in the statistical uncertainties. We then evaluate the model parameter uncertainties by defining the 1σ ranges using the larger bin uncertainties. We will validate these uncertainty estimates by comparing to literature values of N_CL in Section 2.3.

Evaluating models on the basis of χ² only judges models where data exist, but the data may not exist over an interesting range of parameter space. Another way of judging our profiles is based on how well the data that constrain them span the key radial range of 1–15 kpc. In some cases, the data span only a minority of this range, and so even if the models are statistically valid for these galaxies, they will have large associated uncertainties due to the lack of good constraints on the inner profile. To quantify this distinction in profile quality, we consider the fits to those galaxies for which the data do not reach interior to log r = 3.5 (3.1 kpc) to be of lower quality. We designate these profile fits and those of the few galaxies that show highly irregular profiles (such as NGC 4762) as Q = 0 galaxies. This quality index is included in Table 1.

To quantify the number of clusters in each galaxy, we would in principle integrate our model profile to r = ∞. However, because we only empirically constrain the profile within 15 kpc, we are reticent to extrapolate the best fit profile far beyond this radius, particularly when we let the power law slope float. Unfortunately, 15 kpc is clearly too small an outer radius to adopt if we intend to have a measurement of all of the clusters. In the Milky Way, one would be fairly close to the total number if one counted all of the clusters within 50 kpc, even though there are a few clusters beyond 100 kpc (Harris & Racine 1979). Our aim is to obtain a measurement of the total number of clusters, not the number of clusters within a fixed physical radius, because the latter will introduce a dependence of T_N on the physical scale of the cluster system. We show in the upper left panel of Figure 8 that the number of clusters obtained by integrating to 50 and 100 kpc for models with floating power law slope differs modestly (∼20%), with only a slightly detectable systematic variation with the richness of the cluster system (models with fixed slope will have a fixed ratio between N₅₀ and N₁₀₀). We choose, therefore, to avoid the problem of extrapolating the fitted model to 100 kpc and treat the integral out to 50 kpc as the global number of clusters—noting that this might be a slight systematic underestimation. We will return to this issue when we compare our results to those in the literature.

**Figure 8.** Test of various aspects of our procedure. S $^\prime _{\rm N}$ denotes the "new" version of the specific frequency, and S_N is our standard version. In the upper left panel we consider how S_N changes if we integrate our fitted profiles to 100 kpc instead of the standard 50 (for models with free power law slope). In the upper right panel we calculate the new values by utilizing 100 simulations from which to derive our completion corrections rather than the standard 50. In the lower left panel we consider the effects of fitting to the data only between 2 ⩽ *r_p*(kpc) ⩽ 15 rather than 1 ⩽ *r_p* ⩽ 15. In the lower right panel we consider the effect of increasing our color cut for spurious sources from 2.5 to 5 mag. None of these cause large changes to the bulk of the results, particularly at S_N > 1, where the statistical errors are somewhat smaller.
Download figure:
Standard image High-resolution image

**Figure 8.** Test of various aspects of our procedure. S $^\prime _{\rm N}$ denotes the "new" version of the specific frequency, and S_N is our standard version. In the upper left panel we consider how S_N changes if we integrate our fitted profiles to 100 kpc instead of the standard 50 (for models with free power law slope). In the upper right panel we calculate the new values by utilizing 100 simulations from which to derive our completion corrections rather than the standard 50. In the lower left panel we consider the effects of fitting to the data only between 2 ⩽ *r_p*(kpc) ⩽ 15 rather than 1 ⩽ *r_p* ⩽ 15. In the lower right panel we consider the effect of increasing our color cut for spurious sources from 2.5 to 5 mag. None of these cause large changes to the bulk of the results, particularly at S_N > 1, where the statistical errors are somewhat smaller.
Download figure:
Standard image High-resolution image

To obtain the number of clusters in each galaxy, we use the integrated profile to r = 50 kpc and then correct that number for clusters outside of the magnitude range of our detected candidate clusters, assuming a Gaussian LF. Standard parameters for the peak and width of the LF are M_V ∼ −7.4 and σ_V = 1.4 for early types and σ_V = 1.2 for later types (Brodie & Strader 2006). Adopting an estimate for the V−3.6 color of ∼2.4 (Barmby & Jalilian 2012) results in an estimate of the location of the LF peak at M_3.6 = −9.8. The dispersion of the cluster LF is not well-measured at 3.6 μm, but given the decreased sensitivity to age and metallicity differences at this wavelength, we adopt the lower range of σ estimates in the V band, σ_3.6 = 1.2. We adopt the same LF for all galaxies. There is little variation in the globular cluster LF with galaxy luminosity (Strader et al. 2006). Variations of these parameters, within reason, tend to change N_CL by tens of percent, rather than by factors of a few, which is what we shall conclude is the actual uncertainty in our measurement (see Section 2.3).

The specific frequency can then be defined consistently with previous definitions as the number of clusters per parent galaxy luminosity, although here we would be using the 3.6 μm luminosity rather than the historical B or V luminosity. However, because we expect the ratio of clusters to total stellar mass to be the more physically interesting measurement, we go further along this path by using the Spitzer magnitudes and their calibration to stellar mass (Eskew et al. 2012) to calculate a mass-dependent specific frequency. We produce a measurement of the ratio of the number of clusters to stellar mass rather than to the luminosity in one photometric band, normalized so that the values are in a similar numerical range (T_N is in units of clusters per 10⁹M_☉, following the suggestion of Zepf & Ashman 1993).

The Eskew et al. (2012) method for estimating stellar masses comes from a region-by-region comparison of reconstructed star formation histories in the Large Magellanic Cloud (LMC) based on analysis of optical color–magnitude diagrams (Harris & Zaritsky 2009) and Spitzer 3.6 and 4.5 μm photometry from Meixner et al. (2006). Differences in stellar populations among the regions provide an estimate for the robustness of the mass estimates against such variations and result in an uncertainty estimate of 30% (presumably lower in whole galaxies, which should average over the more extreme star formation histories seen in localized regions of the LMC). The fitting formulae presented by Eskew et al. (2012) have been confirmed, apart from uncertainties related to differences in the adopted stellar initial mass function (IMF), by an independent analysis of SDSS spectroscopy (Cybulski et al. 2014) and a detailed analysis of S⁴G photometry (Querejeta et al. 2014).

Before proceeding to discuss the results, we expand a little on our set of choices. First, we have shown in Figure 8 that the choice of integrating our model profile to 50 or 100 kpc produces little change (upper row, left panel). We opt to remain conservative in our extrapolation. Second, we show in Figure 8 a comparison of S_N obtained for the standard parameter choices and that obtained using 100 rather than 50 simulations (upper row, right panel). There is no systematic difference between the two values and only some slight variance at extremely low S_N (<1), where statistical uncertainties dominate. We will continue with the results drawn from 50 random simulations for the completeness corrections. Third, we show in the bottom row, left panel, a comparison of our standard S_N measurement with that obtained using an inner cutoff of 2 kpc in our model-fitting. This panel shows the largest scatter among the four panels, but the scatter is again mostly for S_N < 1, where the statistical error bars dominate. Given that for most systems the smaller cutoff radius results in no significant change in our S_N measurement and that the extra data at small radii can help reduce fitting uncertainties, we proceed with an inner cutoff of 1 kpc in our fitting. Finally, we explore changing the color cutoff used to remove spurious sources from 2.5 to 5 mag in the lower row, right panel. Here we find almost no detectable difference between the resulting S_N measurements, indicating that our cutoff of 2.5 mag is not resulting in the rejection of real sources.

2.3. Comparison to Literature and the Validation of Our Adopted Methodology

To close this section, we compare, where we can, our estimates of N_CL to those in the literature. This comparison is not straightforward. Serious differences can arise in the results from the range of completeness corrections, both in terms of detection of sources in an image and then in terms of correcting for the entire LF of clusters. Differences in the adopted distance can further complicate matters. It is often the case that the completeness corrections are not particularly explicit and certainly differ at least in detail (such as in the peak and width of the Gaussian LF) with ours. Furthermore, some of the best data with which to identify clusters, from the HST, suffer from the small field of view, and therefore spatial completeness corrections must be made. In certain studies (Rhode & Zepf 2004; Spitler et al. 2008b), ground-based data are combined with HST to provide both the improved resolution in the core of the galaxies and the larger coverage beyond. For all of these reasons, it is difficult to construct a homogeneous comparison sample from the literature. These concerns echo those expressed previously regarding literature compilations in general, but we will show below that these concerns appear to be at a level of precision below that which we or the literature studies achieve.

We use the results in the compilation of literature values by Harris et al. (2013) for the comparison shown in Figure 9. We include in the comparison galaxies to which we had assigned a quality flag of 0 in order to maximize the sample size. Figure 9 consists of two panels where we compare the results of the approach in which we fix the background to the measurement obtained at large radii (left panel) and where we fix the power law slope (right panel). The line plotted in both panels is the 1:1 line, not a fit. It does a satisfactory job of describing the mean trend, although there is clearly significant scatter about the line in the left panel, beyond that described by the internal uncertainties.

We find that the fixed slope modeling has two significant advantages over the fixed background modeling. First, the asymmetric tail of points at low N_CL that is seen in the left panel of the figure is absent in the right panel. The tail in the fixed background measurements is presumably due to systems that are poor in globular clusters but whose background level is slightly overestimated by the measurement at large radius. In such cases, we conclude that there are no clusters, whereas the more precise literature samples are able to recover the actual, but small, number of clusters. Second, the scatter about the line is visibly reduced using the fixed slope approach. Again, we suspect that systematic errors in the background estimation are what is driving the majority of the increased scatter in the left panel, although systems with poorly constrained power law slopes contribute to scatter because of the necessary extrapolation from 15 to 50 kpc. Therefore, both because of background systematic errors and our inability to constrain the power law slope beyond ∼15 kpc, we conclude that the fixed slope approach is the more robust. With the exception of a few outliers, which mostly have a quality flag value of 0, the internal uncertainty estimates do a credible job of representing the scatter between our measurements and those in the literature, when we adopt the fixed slope approach.

For the fixed slope approach, the scatter about the 1:1 line is 0.29 dex for the Q = 1 sample. If we attribute all of the scatter to our measurements, this scatter suggests that the uncertainty in our measurements is roughly a factor of two. Although this is obviously a large number in absolute terms, relative to the range in N_CL of several orders of magnitude the uncertainty is relatively modest. The mean offset from the 1:1 line is only −0.046 dex (we underestimate the literature values by 10%), confirming that the use of N₅₀ as a proxy for N_CL is valid and that our 50 kpc integration limit produces a systematic error that is significantly smaller than the random errors. This quantitative agreement with the literature values in the mean also supports our decision not to include the correction for the background bias described by Harris (1986) due to "contamination" of the background estimates by the cluster population itself. We empirically determined the correction to be small in our data relative to our uncertainties. Hereafter, we refer to N₅₀ as N_CL.

The results described above are valid for any choice of power law exponent between the plausible range of −2 to −2.5 in projected surface density valid for normal ellipticals (Brodie & Strader 2006). However, we have gone beyond that in producing the results in Figure 9 in that we have chosen the value of the fixed power law index to minimize the scatter in the right panel. We obtain the result by recalculating N_CL for different choices of index, stepping in units of 0.1, and evaluating the scatter produced in the analog of the right panel of Figure 9. We find that we prefer a power law slope of −2.4 and adopt this slope universally in this study. The three-dimensional density will have a power law dependence that is one unit steeper, so the integral easily converges.

We adopt a constant power law index, despite previous findings that the index depends on M_* (Brodie & Strader 2006). If the radial slope does depend on mass, then fixing the slope could introduce an artificial trend in T_N with M_*. The previous claims are that the slope steepens as one progresses to lower-mass systems. Such an effect could result in our underestimating N_CL in high-mass galaxies and overestimating in low-mass ones. However, when we examine the correspondence between our data and the literature values, the best fit slope obtained using the ordinary least squares (OLS) bisector method presented by Isobe et al. (1990) to account for uncertainties in both axes is 1.14 ± 0.24 and therefore consistent with the 1:1 line. We also find no significant correlation between our measured power law slopes and M_* for log(M_*/M_☉) > 10 (below that mass, our fits with free power law index are unreliable). These results do not demonstrate the absence of such a relationship, only that it is too weak to detect with a sample of this size and scatter in N_CL that is a factor of two.

Similarly, we use the consistency of the literature comparison with the 1:1 line to argue that other simplifications we have adopted, such as the constancy of the peak magnitude of the cluster LF and the fixed integration to 50 kpc regardless of galaxy size or mass, are not affecting our measurements at a level that is noticeable, given the uncertainties. Such a statement could simply indicate that we have inferior measurements—after all, the inability to measure an effect is not a desirable attribute—but, as shown in Figure 9, our estimated uncertainties are comparable or only slightly larger in general than those claimed by other investigators. Furthermore, in defense of our approach we cite the homogeneity of the sample and our algorithmic approach, which are critical in comparisons across galaxy luminosity, type, and mass.

The reader may still wonder why our measurements appear to be robust to variations in the radial power law slope and the peak magnitude of the cluster LF. The former likely arises because of the robustness of the integral under the model fits. Even with the wrong power law slope, the fit is likely to represent the mean surface-density cluster values reasonably well over the fitted radial range, and when we integrate over the radius to arrive at a total cluster population we are relatively insensitive to the slope of the fit, as long as we evaluate the integral over a similar radial. As for the peak magnitude, our images are complete to below the peak magnitude, even with the anticipated possibility of a varying peak magnitude, and therefore our completeness corrections are never going to be larger than a factor of two and clearly are often much better. Because the observational scatter is a factor of two, variations due to improper corrections related to a variable peak LF magnitude are not easily detected.

3. DISCUSSION

It is evident from the surface-density profiles (Figures 6 and 7) that for ∼50% of the galaxies there exists a clustered source population, which consists presumably of globular clusters; that a power law is an adequate description for the radial distribution of these sources, given the current state of the observations for most galaxies; and that there is a range in the properties of that population (numbers and radial extent) among the galaxies in our sample. Our argument that these sources are globular clusters is circumstantial because we lack the resolution to confirm their nature. We base our conclusion on the following: (1) they evidently cluster about the parent in those cases where an excess is seen; (2) an excess is not found in every case, even when the data extend to small radii, demonstrating that the sources do not spuriously arise from residuals in the parent galaxy subtraction process; (3) the sources have absolute magnitudes consistent with those expected for globular clusters; (4) the radial distribution of sources, where that is well-measured, lie in the range of r⁻² to r⁻⁴, which is consistent with previous measurements (Ashman & Zepf 1998; Brodie & Strader 2006); (5) the number of excess sources, ranging from a few up to a thousand, are in the range expected based on previous cluster population studies (Harris 1991; Brodie & Strader 2006); and (6) in the few cases where we can compare our measurements to higher-fidelity measurements in the literature, we reproduce prior results. From now on, we will refer to these excess sources as clusters, although it must be understood that these populations could be partially contaminated by other sources that also cluster about the parent galaxy. In general, one possible such source would be star-forming knots (Thilker et al. 2005; Gil de Paz et al. 2005; Zaritsky & Christlein 2007; Herbert-Fort et al. 2012) and intermediate-age versions of such structures, although for our early-type galaxies those should not be a major source of contamination. Individual stars are generally insufficiently luminous to match our magnitude criteria (− 11 < M_3.6 < −8), although some especially luminous stars (Blum et al. 2006)—again, not as likely in early-type galaxies—could fall in this range.

3.1. Numbers of Clusters

We show the number of clusters, N_CL, as a function of parent galaxy stellar mass, M_* in Figure 10. Some basic qualitative results can be drawn quickly from the figure: (1) there is a general increase in the N_CL that roughly follows the rise in stellar mass; (2) despite this mean trend, for M_* < 10¹⁰ M_☉, the variation in N_CL approaches a factor of 100 at a given M_* and may reflect different galaxy characteristics; (3) for M_* < 10¹⁰M_☉, the behavior of N_CL can be characterized as flat; and (4) neither the trend at large M_* nor the scatter at low M_* correlate strongly with morphology (within this admittedly early-type sample). We expand on each result below.

A linear fit to the data using orthogonal regression (Isobe et al. 1990) results in a best fit slope of 0.94 ± 0.19, supporting the qualitative suggestion that at least roughly N_CL∝ M_*. If we only use the data for log (M_*/M_☉) > 10, which avoids the more ambiguous lower-mass galaxies, then the slope is 1.39 ± 0.36 (1.24 ± 0.19 if using the bisector method). Slopes >1 suggest that T_N is increasing with M_*. This result is consistent with that of Harris et al. (2013) for this same mass range. The implication of such a result is that more massive galaxies are more capable of producing clusters. Such a result has repercussions for merger models in which clusters may form (Schweizer 1987; Ashman & Zepf 1992) and for models in which merging is dissipationless. However, as we will discuss in the next section, such interpretations are premature until we understand the nature of the stellar IMF and its effects on our estimates of M_*.

Among galaxies with log (M_*/M_☉) < 10, T_N exhibits large scatter, with values that can differ by as much as two orders of magnitude (Figure 11). These galaxies have limited populations of clusters, and therefore statistical errors in combination with subtle systematic errors may be responsible for the large scatter. The trend for at least some of these systems to have proportionally larger N_CL than expected for their mass is clearly visible here and in the study of Harris et al. (2013). Furthermore, our comparison to the literature values of N_CL suggests that we are obtaining reliable estimates of N_CL, with plausible internal uncertainty estimates, even for galaxies with low N_CL. Unfortunately, the object-by-object comparison to the literature is limited to a few galaxies with similarly low N_CL, and therefore we may simply have been fortunate in those galaxies that overlapped existing studies. However, other studies find significant scatter in specific frequency as well among low-mass galaxies. For example, Miller et al. (1998) find S_N values that range from 0 to 23 for a set of dwarf elliptical galaxies. The quantitative values are not directly comparable to our results because their normalization is done with respect to an optical luminosity, but the range in values is greater than an order of magnitude. Also, Strader et al. (2006) suggest that there are two families of galaxies among the dwarf ellipticals, one with S_N ∼ 2 and another with S_N ∼ 5–20. Again, the numbers are not directly comparable to ours because of the use of bluer bands and our conversion to stellar mass, but the suggestion of two families is consistent with the visual impression of our Figure 10. The existence of two populations is not evident in Harris et al. (2013).

We also find that there is marked change in the behavior of N_CL with M_* below a galaxy stellar mass of ∼10¹⁰ M_☉, consistent with previous studies such as Georgiev et al. (2010) and Harris et al. (2013). As we discussed when reviewing the increased scatter at low M_*, these systems are more susceptible to systematic errors. Nevertheless, there is little if any dependence of N_CL with stellar mass for log(M_*/M_☉) < 10, suggesting that some of the low-stellar-mass galaxies have the largest T_N. It appears that large variations in T_N, perhaps driven by environment or star formation history (for examples of related phenomenon; see Georgiev et al. 2010), are present in galaxies with low stellar masses and that those differences average out as one examines galaxies with increasing M_*. Alternatively, these systems may also have lost much of their baryonic material (presumably not through tidal stripping, which would affect stars and clusters, but through stellar feedback which could have removed material destined to form stars but not clusters). As such, the range of T_N may help constrain models of stellar feedback (Dekel & Silk 1986). These speculations motivate more study of low-mass galaxies to understand the drivers of cluster formation. Of course, increased verification of N_CL measurements in low M_* galaxies would need to be part of such a program.

At the high-mass end, we do not find the sharp rise in T_N found by Peng et al. (2008). There are various possible reasons—some trivial, some interesting—for this disagreement. First, we have only four galaxies in this mass range, so we may simply have been unfortunate (likewise, Peng et al. 2008 have a small number in this mass range). Potentially more interesting is that the Peng et al. (2008) sample consists exclusively of Virgo galaxies, and the increased number of clusters that they find around the massive galaxies may be the result of stripping of clusters from lower-mass galaxies in the cluster environment.

Finally, morphology, at least over this limited range of early-type galaxies, appears to play at most a minor role in determining T_N. In Figure 12 we compare the distribution of T_N for three different subtypes of galaxies. No strong difference stands out, with all three types having predominantly more galaxies with low T_N and all having at least a few galaxies with moderate T_N. Larger samples will be needed to tease out the hints of differences that appear among the panels, if they exist, but such differences would complicate the analysis. For example, Spitler et al. (2008b) found stronger correlation between T_N and M_* than what we find, but they include spirals in their sample. Judging from their Figure 15, the inclusion of spirals, which appear as a class to have lower T_N than the earlier types, helps drive the correlations because they also tend to have lower M_*. It will therefore be imperative in future studies that include a broad range of morphological types to allow for a dependence of T_N on both M_* and morphology.

3.2. Effect of IMF Variations

When considering the behavior of N_CL with M_*, we have exclusively concerned ourselves so far with uncertainties in the ordinate of Figure 10. However, the stellar masses are also vulnerable. Although there are systematic uncertainties in the overall mass normalization due to uncertainty in the true IMF (see Meidt et al. 2014 for discussion), changing the mass normalization from that advocated by Eskew et al. (2012) simply shifts all the points the same distance along the abscissa and leaves the N_CL − M_* slope unchanged. However, if the IMF variation is correlated to a galaxy's total stellar mass, such behavior will alter the slope of the N_CL − M_* relationship. All of the existing mass estimators (such as that from Bell & de Jong 2001; Eskew et al. 2012; Meidt et al. 2014) are predicated on a universal IMF, so no existing simple mass estimator resolves this issue.

We approach the problem in an orthogonal manner, by postulating that the cluster specific frequency should not depend on M_*. Although there is no known reason why this must be the case, it is a natural possibility. Because our current data indicate that it depends weakly on M_* (see Figure 10 and previous discussion), we now ask what IMF behavior would completely remove this trend. If the true stellar mass is $M^\prime _*$ , and M_* is the mass inferred using a Salpeter IMF, M_Salpeter, which is the case for the Eskew et al. (2012) estimator that we use, then $\log ({M}^\prime _*/{M}_{{\rm Salpeter}}) \propto 0.39 \log ({M}_{{\rm Salpeter}}/{M}_\odot)$ (where the measured M_* ≡ M_Salpeter) is required to set the slope of the relationship between log N_CL and log (M_*/M_☉) to 1 (using our previous fit to the relation between these quantities for galaxies with log(M_*/M_☉) > 10).

In Figure 13 we show the relationship between M_*, $M_*^{\prime }$ , and M_Salpeter as measured by Conroy & van Dokkum (2012) for their sample of galaxies on the basis of detailed analysis of integrated spectra. We evaluate M_* using the I-band luminosities of these galaxies and their tabulated values of M/L_I. To obtain the ratio of the true stellar mass to that derived with a Salpeter IMF for each galaxy, we use their tabulated values of M/L_K for each galaxy and the M/L_K value for the Milky Way (they use the ratio of these quantities to normalize out age and metallicity differences among their sample galaxies). We then correct this ratio by 1.6, to account for the fact that the MW galaxy in their modeling has a Kroupa IMF rather than a Salpeter IMF; however, this simply results in a renormalization of the coordinate axis in Figure 13. Last, in the same figure we plot the relationship discussed previously, $\log ({M}^\prime _*/{M}_{{\rm Salpeter}}) \propto 0.39 \log ({M}_{{\rm Salpeter}}/{M}_\odot)$ , normalized to best match the Conroy & van Dokkum (2012) data.

The quantitatively excellent match of the observed trend in the Conroy & van Dokkum (2012) sample and our deduced IMF behavior from the proposition that T_N is independent of M_* help bolster the arguments on both sides—that the IMF variations are indeed real and that T_N is independent of M_*. To be precise, the result of the orthogonal regression for the Conroy & van Dokkum (2012) sample shown in the figure (again, using the formulas from Isobe et al. 1990) has a slope of 0.37 ± 0.08, well within 1σ of our slope estimate of 0.39 derived from the proposition that T_N equals a constant.

We pause briefly here to stress that the results we present are only a weak indirect argument for IMF variations. A dependence of T_N on stellar mass could arise from a variety of physical factors or even from systematic errors in our estimates of N_CL. Our principal point here is that the magnitude of such effects is consistent with the magnitude of the effect arising from IMF variations that are suggested by completely different lines of evidence. As such, IMF variations cannot be ignored when considering trends in T_N.

The most pertinent physical measurement of the specific frequency is perhaps the ratio of mass in globular clusters to the total stellar mass of the galaxy (rather than the number of clusters). This approach has been advocated by McLaughlin (1999) and Harris et al. (2013), and that ratio is defined as epsilon ^b (they also advocate using the baryonic mass rather than the stellar mass and perhaps even relating this all to the halo mass). In the standard approach of a universal IMF and cluster LF, this distinction would not result in any difference in what we have discussed prior to this section. Even if the IMF varies, but it varies in concert for both the stellar populations within galaxies and their globular clusters, then a specific frequency defined as the ratio of the masses of clusters to total stellar mass would be independent of the adopted IMF and any variations of that IMF with M_*.

In general, this type of discussion assumes that the IMF within the clusters themselves is universal. There is, however, evidence that the IMF among clusters can also vary (Zaritsky et al. 2012, 2013a), although the driver of that variation has not yet been identified, so we cannot attempt to model its potential impact on the results described here. Such variations would affect analyses that are based on the stellar mass within clusters (see Harris et al. 2013 for a discussion relating stellar mass in clusters to galaxy total mass).

These subtleties become important when trying to place the results in the more global context of galaxy formation. For example, the massive (log(M_*/M_☉) > 10) portion of the sample shown in Figure 11 has a mean value of T_N of 11.8 that, when combined with the mean cluster mass of 1.2 × 10⁵ M_☉ obtained from our adopted LF and the Eskew et al. (2012) mass calibration, suggests a cluster formation efficiency by stellar mass of 0.14% for an assumed universal IMF. Uncertainties of a factor of two in these assumptions are well within our expectations. McLaughlin (1999) measures epsilon ^b = 0.26% ± 0.05%, about a factor of two larger than our estimate. Whether this difference is due to differences in the galaxy sample, differences in stellar mass estimates, cluster counts, or simple statistical noise (see Figure 11) is not evident. As such, uncertainties of this magnitude impact arguments relating to whether the cluster population scales more directly to stellar or total mass (e.g., Strader et al. 2006; Georgiev et al. 2010; Harris et al. 2013).

Interestingly, the high efficiency suggested by the richest low-mass systems, about a factor of 10 larger than that for the higher-mass galaxies, is comparable to the efficiency derived for cluster formation in dwarf galaxies at high redshift (Elmegreen et al. 2012). This agreement suggests that the cluster populations in these galaxies may be undisturbed, with little cluster destruction beyond cluster "infant mortality." In more massive galaxies, cluster destruction is expected to be significant (Gnedin & Ostriker 1997), with destruction fractions that can be as large as 90% (allowing for a reconciliation of our efficiencies in low- and high-mass galaxies). Unfortunately, any quantitative prediction of the destroyed cluster fraction depends on various assumptions, including the characteristics of the original population.

4. SUMMARY

The S⁴G images (Sheth et al. 2010) provide another opportunity to explore the properties of globular cluster populations in galaxies. In particular, the images are sufficiently deep to reach well into the globular cluster LF and have large fields of view, enabling reasonably complete surveys for clusters in a sample of thousands of galaxies that surveys the local volume. Although the data provide neither the angular resolution nor multiwavelength coverage of the latest best examples of similar efforts, the strength of this program is in its large, homogeneous sample of galaxies. The results obtained here on the basic properties of globular cluster populations, number, and specific frequency as a function of various galaxy properties complement the higher-fidelity surveys undertaken with the HST and deep ground-based surveys. Our empirical findings are consistent with those based on compilations of existing data (e.g., Harris et al. 2013), but given the nature of the literature samples, where different categories of galaxies are the focus of independent studies, we advocate that quantitative measures and comparison be based on homogeneous data and analysis.

We present our initial foray into this topic with the S⁴G data, concentrating here on early-type galaxies in the core survey. This subsample consists of 97 galaxies drawn from the over 2300 galaxies in the S⁴G data set. We avoided late-type galaxies to minimize confusion between HII regions, luminous stars, and globular clusters, but an extension of this work to edge-on late-type galaxies should be straightforward. A large sample of well-studied edge-on galaxies exists within S⁴G (Comerón et al. 2011). In addition, we have begun a program to enlarge the complement of early-type galaxies in S⁴G, galaxies that our initial selection criteria discriminated against; therefore an extension of this work to include the larger sample of early-type galaxies will be forthcoming.

Using this initial set of galaxies, we established the following.

1.
Accurate measurements of the number of clusters, N_CL, in each galaxy can best be obtained by fitting a radial surface-density profile to all detected point sources with absolute magnitudes between −11 and −8, if sources are assumed to be at the distance of the galaxy, i.e., a power law with an exponent of −2.4. The normalization and background level (i.e., surface density of contaminants) are fit. Comparison to independent measurements of N_CL demonstrates that our measurements have an uncertainty of a factor of ∼2 and are well-described by the internal uncertainties we calculate.
2.
Over the full range of galaxy stellar masses, M_*, N_CL rises nearly proportionally to M_* (log(N_CL)∝(0.94 ± 0.19)log (M_*/M_☉)). However, the behavior of N_CL appears to change significantly below log(M_*/M_☉) ∼ 10. This result confirms previous findings (Georgiev et al. 2010; Harris et al. 2013). If we only fit to galaxies with log(M_*/M_☉) > 10, then the relationship steepens to log(N_CL)∝(1.39 ± 0.36)log (M_*/M_☉). Because of the large uncertainty, we have only marginal evidence that this slope is larger than 1 (although a bisector fit of the data results in a similar result, slope =1.24 ± 0.19). Definitive evidence of a slope larger than 1 would demonstrate that more massive galaxies, or their progenitors, are somehow more efficient producers of globular clusters and cannot therefore simply be the result of dissipationless mergers of lower-mass early-type galaxies.
3.
Because determining if T_N varies with galaxy stellar mass is critical for cluster formation models, we took care to quantify the effect of a varying stellar IMF on the determination of the galaxy stellar mass. In particular, we demonstrate that if the IMF of these galaxies depends on galaxy mass (as has been suggested by a number of recent studies; Conroy & van Dokkum 2012; Cappellari et al. 2012), the apparent, but weak, increase in T_N with M_* is entirely removed. In fact, when we postulate that T_N is constant with M_*, we recover the exact dependence of the IMF with stellar mass observed in the data of Conroy & van Dokkum (2012). We conclude that T_N is independent of mass for early-type galaxies with log (M_*/M_☉) > 10. It is possible that the most massive ellipticals, or central dominant ellipticals, deviate from this trend (McLaughlin 1999; Peng et al. 2008); we have insufficient data to constrain T_N for that subclass. However, the lack of a dependence of T_N on M_* in general is consistent with simple predictions from models of dissipationless mergers and of dissipational mergers where any increased star formation during the merger stage produces a corresponding number of globular clusters. With respect to the latter scenario, we note that globular cluster candidates consistent in age with the star burst phase in E+A galaxies have been detected (Yang et al. 2004).
4.
At log(M_*/M_☉) < 10, the scatter in N_CL is large (two orders of magnitude) and so may suggest that we are detecting strong differences in cluster formation and destruction efficiencies among this population of galaxies. This result confirms previous studies of the specific frequency among low-mass galaxies (Miller et al. 1998; Strader et al. 2006; Georgiev et al. 2010). Our data also hint at the possibility of two T_N families of galaxies at low M_*, also supporting previous speculations (Miller et al. 1998; Strader et al. 2006). Progress in understanding the drivers of cluster formation may come principally from a focused effort in studying the cluster populations of low-mass galaxies. However, because these populations are inherently poor in numbers, large samples will be needed to counter the larger statistical uncertainties associated with the cluster populations of individual galaxies.

Using the S4G images and the larger sample of early-type galaxies currently being collected, we expect this work to expand to a significantly larger sample of galaxies that both produces more statistically significant results for the cluster populations of early types and a comparison sample of the cluster populations of late-type galaxies. Ultimately, the power of having measurements of T_N for the S⁴G sample will also lie in the detailed ancillary analysis products such as photometric decomposition, detailed morphologies, radial and vertical disk studies, and measurements of stellar and gas masses. In all, we expect this series of papers to provide a broad but comprehensive outline of the basic properties of globular cluster systems in galaxies.

D.Z. acknowledges financial support from NASA ADAP NNX12AE27G and NYU CCPP for its hospitality during long-term visits. L.C.H. acknowledges support from the Kavli Foundation, Peking University, and the Chinese Academy of Science through grant XDB09030102 (Emergence of Cosmological Structures) from the Strategic Priority Research Program. E.A. and A.B. acknowledge financial support from the CNES (Centre National d'Etudes Spatiales, France) and from the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme FP7/2007-2013/ under REA grant Agreement PITN-GA-2011-289313 to the DAGAL network. The authors thank Ryan Leaman and Mike Beasley for valuable comments. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Finally, we thank our conscientious referee for valued comments and suggestions.

GLOBULAR CLUSTER POPULATIONS: FIRST RESULTS FROM S⁴G EARLY-TYPE GALAXIES

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

ABSTRACT

1. INTRODUCTION

2. DATA AND MEASUREMENTS