GLOBULAR CLUSTER POPULATIONS: RESULTS INCLUDING S⁴G LATE-TYPE GALAXIES

As in Paper I, the parent sample is the S⁴G sample, which currently consists of 2352 galaxies (Sheth et al. 2010). It is primarily a volume-limited sample, but additional selection criteria, such as the existence of an H i redshift, preclude it from being a complete, volume-limited sample. We provide and analyze images of these galaxies obtained with the Spitzer Space Telescope (Werner et al. 2004) and its Infrared Array Camera (IRAC; Fazio et al. 2004). The data are publicly available through the archive on the NASA IRSA website.²¹

To complement the early-type sample of Paper I ($-5\;\leqslant$ T-type $\leqslant \;1$), we focus now on galaxies morphologically classified as late type ($4\;\leqslant$ T-type $\leqslant \;10$) by Buta et al. (2015) and select only those that are nearly edge-on (inclination ≥80°) to minimize the area on the sky in which the internal structure in the parent galaxy could be mistaken for clusters and to maximize the area in which the clusters will be detected and not hidden by the parent galaxy. The inclinations are based on ellipticities measured in the S⁴G pipeline (Muñoz-Mateos et al. 2015). The inclination cut, in particular, results in a sample that is a small fraction of the full sample. We further constrain the sample by including only those galaxies within a suitable range of distances, again following the criteria defined in Paper I. We set the upper end of the distance range to correspond to a distance modulus of 32.4 (30.1 Mpc), where we have enough physical resolution to adequately sample the globular cluster population radial profile. We set the lower end of the distance modulus range at 30.25 (11.2 Mpc), where we ensure sufficient background coverage within the images with which to constrain the background source density. Finally, we remove from the sample any galaxies that have bright nearby neighbors that would compromise the analysis and the one galaxy (UGC 1839) in the S⁴G sample that satisfies all of these criteria but for which there is not a well-determined magnitude in the S⁴G catalog. Because of these restrictive constraints, we re-examined all edge-on galaxies in the full sample for additional candidates. We found two (IC 5269 and NGC 1145) that are classified as early types (T-types −1 and 0, respectively) but appear to be disk dominated and satisfy our distance and isolation criteria. We add these two galaxies to our sample and list the final 73 galaxies in Table 1. 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.

Table 1.  Globular Cluster Population Properties

Name	DM	T	m3.6	m4.5	N50	TN	Background	Q
ESO 107-016	32.26	8	15.3	15.7	129	192${}_{-164}^{+320}$	−4.64	1
ESO 146-014	31.64	8	15.2	15.8	27	35${}_{-32}^{+99}$	−5.26	1
ESO 249-035	31.77	8	17.0	17.5	75	627${}_{-607}^{+962}$	−4.92	0
ESO 292-014	32.24	7	13.1	13.5	403	68${}_{-36}^{+49}$	−5.14	1
ESO 346-001	32.10	7	13.4	13.7	38	11${}_{-9}^{+49}$	−4.78	1
ESO 356-018	31.59	9	15.0	15.5	81	137${}_{-67}^{+96}$	−5.07	1

Note. 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. (2015) 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. Background is the logarithm of the surface number density of background objects. The quality flag Q is defined by whether the data extend sufficiently over the radial range of interest to provide a robust constraint on the fitted models (Q = 1 is good, Q = 0 is poor).

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:  DataTypeset image

To produce our measurements of the cluster population size, we follow the procedure described in detail in Paper I and briefly summarize here. We use the bright object masks and exposure weight maps developed as part of the S⁴G processing (Muñoz-Mateos et al. 2015; Salo et al. 2015). 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 adopt the smallest threshold value that eliminates such features. We perform basic preprocessing steps including sky subtraction plus modeling and removal of the primary galaxy. We calculate the background sky value by evaluating the median of the unmasked pixels 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. 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. 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 100 (simply the first NGC galaxy on our list). The multipolar residual pattern seen in Figure 2 is typical. Except along the disk major axis, we can identify individual sources quite close to the center. Once the residual images are available, we run SExtractor (Bertin & Arnouts 1996) on unmasked versions to identify point sources, eventually using the stellarity index to reject clearly extended sources. For the inner regions, where the model subtraction is most important, we require consistent detection and photometry of candidate clusters in the 3.6 and 4.5 μm images to help reject spurious sources. We apply this criterion only in the central regions because the two images generally do not overlap at large radii. Candidate clusters are defined as objects within an absolute magnitude range given by $-11\lt {M}_{3.6}\lt -8$, defined at the upper end by the known globular cluster luminosity function and at the lower end by requiring high completeness over our entire sample. We apply source count corrections as a function of both magnitude and location within the image, determined by adding 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). Details of parameter choices in this and previous steps are given in Paper I. We made no modifications to the procedure between this paper and our previous work and present validation of these parameter choices in Paper I.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 of insufficient quality to allow for the fitting of the power law and background simultaneously. We follow the preferred approach from Paper I, where we fix the power-law slope at −2.4 and vary the power-law normalization and background level, still fitting the model for radii between 1 and 15 kpc by minimizing χ². The power-law exponent was set in Paper I by maximizing the concordance between our results and those of Harris et al. (2013) for galaxies in common. The specific radial range was set at the lower end by the minimum radius at which the model subtraction provided low-enough residuals for cluster detection and at the upper end by the maximum radius at which there was often signal above the background in the cluster surface density profiles. As shown in Figure 3, the power-law plus constant-background model does a satisfactory job of fitting almost all 73 galaxies.

Zoom In Zoom Out Reset image size

We present the adopted and resulting values for the distance modulus (DM), Hubble T-type (T), the 3.6 and 4.5 μm apparent magnitudes, the integrated number of clusters from a radius of 5 kpc out to a radius of 50 kpc (N₅₀), 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), the logarithm of the surface number density of background sources, and a quality flag (Q, described below) in Table 1. The choice of the inner boundary in the integration affects the total numbers of clusters inferred, but not the relative numbers among galaxies because of our choice of a fixed power-law slope. The specific choice of a 5 kpc integration lower limit is inherited from Paper I, where the sample consists of large elliptical galaxies. Perhaps a more appropriate choice now is a radius of 3.1 kpc, which is the smallest radius for which we have measurements in all of the galaxies we will consider (see below). Changing the lower bound of the integration from 5 to 3.1 kpc results in an increase of all the globular cluster specific frequencies by 35%, but we retain the 5 kpc inner bound for consistency with Paper I. We will return to this issue when we discuss globular cluster specific frequencies in an absolute context.

We base our uncertainty estimates for T_N on Poisson statistics in the individual radial bins, propagated through the fitting using Δχ². In cases where the model fit is statistically acceptable (defined as having a probability >0.33 of the data being drawn from the model), we adopt the uncertainties corresponding to the 1σ range in the model parameters. In cases where the model is statistically unacceptable, for the adopted Poisson uncertainties in the individual bins, we calculate the value of χ² at which the probability of the data being drawn from the model is 0.5. We increase the uncertainties in each bin, all by the same scaling factor, so that the resulting χ² now has this value. Rescaling is required for the majority of the galaxies because the internal (Poisson) statistics are often underestimates of the full uncertainty. This approach corresponds to using the scatter about the fit to estimate the true uncertainty. In Paper I we externally validated the determined uncertainties by comparing our measurements of the number of clusters, N_CL, to those in the literature for the subset of galaxies for which this comparison is possible. Unfortunately, we cannot expand on that comparison here because none of the galaxies in this study are also in the Harris et al. (2013) study.

We noted in Paper I that χ² only judges models where data exist, but that data may not exist over the critical range of parameter space. We then quantified how well the data constrain the model over the key radial range of 1–15 kpc by applying a binary flag, Q = 0 for galaxies in which the data do not reach interior to $\mathrm{log}r=3.5$ (3.1 kpc) and 1 otherwise. This quality index for our new sample of galaxies is included in Table 1. In Paper I and here, we opt to use only Q = 1 profiles for our subsequent discussion and, by doing so, remove 33 galaxies from our sample for subsequent discussion. This choice was supported in Paper I by the superior matches between our measurements and those of Harris et al. (2013) for the Q = 1 sample.

As in Paper I, we quote the integrated number of clusters within 50 kpc based on the best-fit profile. The question of how far out in radius to integrate the profiles had greater bearing in Paper I because there we explored fitting different power-law slopes to different galaxies. Ultimately, we decided that the data were insufficiently constraining to allow this freedom and settled on the fixed power-law slope. Because we adopt that constant power-law slope here, the decision to limit our cluster counts to radii $\leqslant$ 50 kpc does not affect the relative number of clusters we measure among galaxies. We integrate the power-law 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 luminosity function that is the same for all galaxies. We adopt the same standard parameters for the peak and width of the luminosity function (LF) as in Paper I: M_V ∼ −7.4 and σ_V = 1.4 for early types and σ_V = 1.2 for later types (Brodie & Strader 2006). For a V − 3.6 color of ∼2.4 (Barmby & Jalilian 2012) the location of the LF peak lies at M_3.6 =  −9.8. The dispersion of the cluster LF is not well measured at 3.6 μm, so we adopt the lower range of σ estimates in the V band, σ_3.6 = 1.2. There is little variation in the globular cluster luminosity function with galaxy luminosity (Strader et al. 2006). Variations of these parameters, within reason, tend to change the numbers of clusters by tens of percent, rather than by factors of a few, which is what we concluded in Paper I to be the actual uncertainty of our measurements. Nevertheless, one of the greatest sources of systematic uncertainty comes from our adoption of universality in the globular cluster population radial profile shape and luminosity function. A comparison to external studies (see Paper I) that make different assumptions provides guidance on the magnitude of this uncertainty and supports our adopted uncertainties.

To convert the number of clusters to a specific frequency, we use the stellar mass of the parent galaxy as calculated from the Spitzer magnitudes and their calibration to stellar mass (Eskew et al. 2012). More detailed stellar mass and stellar population modeling of the S⁴G galaxies exists in a set of studies (Meidt et al. 2012, 2014; Querejeta et al. 2015; Röck et al. 2015), but using the Eskew et al. (2012) calibration, which is coarser but consistent with those other studies, provides a direct, self-consistent mass estimate for the clusters and enables straightforward extension to other galaxies beyond those in S⁴G.

We use the mean trend between T_N and M_* to address two simple questions: What fraction of stars end up in long-lived, massive clusters? Where are most clusters found in the current universe? Using the results we present in Figure 4, we adopt the following expression for T_N:

$$\begin{eqnarray}{T}_{{\rm{N}}}=\left\{\begin{array}{ll}{10}^{6.7}{M}_{*}^{-0.56} & \mathrm{if}\;{10}^{8.5}\leqslant \;({M}_{*}/{M}_{\odot })\lt {10}^{10.5}\\ 8.3 & \mathrm{if}\;{10}^{10.5}\leqslant \;({M}_{*}/{M}_{\odot })\leqslant {10}^{11}\\ {10}^{-6.11}{M}_{*}^{0.63} & \mathrm{if}\;({M}_{*}/{M}_{\odot })\gt {10}^{11}\end{array}\right.\end{eqnarray} \tag{ 1 }$$

which simply follows the fit shown in Figure 4 for $8.5\lt \mathrm{log}({M}_{*}/{M}_{\odot })\lt 10.5$, is flat for intermediate M_*, and then rises for the most massive galaxies, as found by Peng et al. (2008). Although we choose a T_N slope to match the Peng et al. (2008) data in this mass range, we normalize their relation to achieve continuity with ours at lower masses. For an estimated mean cluster mass of 1.2 × 10⁵M_⊙, obtained from our adopted cluster luminosity function and the Eskew et al. (2012) mass calibration, we express the cluster stellar mass fraction, the fraction of stars in a galaxy that are in clusters versus the field, as 1.2 × 10⁻⁴T_N, which when combined with the expression above for T_N suggests a moderately increasing cluster stellar mass fraction with decreasing parent galaxy mass over the stellar mass range that we best constrain, $8.5\lt \mathrm{log}({M}_{*}/{M}_{\odot })\lt 10.5$.

The resulting values of the mass fraction range from 10⁻³ at M_* = 10^10.5M_⊙ to 0.013 at M_* = 10^8.5M_⊙, or roughly a factor of 10 increase over the range of parent galaxy stellar mass explored best here. We discussed in Paper I how interpreting the cluster stellar mass fraction as representative of the cluster formation efficiency for the lowest-mass galaxies in this range matches theoretical expectations in some models for the cluster formation efficiency in dwarf galaxies at high redshift (Elmegreen et al. 2012), which further suggested that the cluster populations in these galaxies may be dynamically undisturbed to the present day. We caution that global changes to the inferred cluster mass fractions will also arise with changes to the adopted inner integration boundary of the cluster radial distribution profiles used to calculate the total number of clusters in each galaxy. Any such change will affect all galaxies in equal proportions, thereby leaving comparisons among galaxies unaffected.

The decreasing cluster stellar mass fraction, as we move toward higher parent-galaxy stellar masses, could reflect either a true decrease in cluster formation efficiency or a greater rate of cluster disruption. Theoretical models show that cluster destruction can be extremely effective in massive galaxies (Gnedin & Ostriker 1997), resulting in the destruction of as many as 90% of the original clusters. If the destroyed fraction is this large at the upper end of our stellar mass range and if the rate of dynamical evolution decreases in lower-mass galaxies, then dynamical evolution of the cluster populations would be a straightforward explanation of the observed trend. The problem, however, involves a complex interplay of various factors and must be treated carefully. A recent study by Gnedin et al. (2014) that explores the interplay between the evolution of the globular cluster population and the growth of nuclear star clusters, and eventually central black holes, is one example of how such modeling can proceed. For higher parent-galaxy stellar masses, ∼10¹¹M_⊙ and above, they predict an increasing fraction of stellar mass in globular clusters, consistent with the Peng et al. (2008) results that we have spliced onto Figure 4; unfortunately they do not model lower-mass parent galaxies. Mieske et al. (2014) specifically show how the qualitative behavior we observe in T_N versus M_* can arise from tidal disruption.

Although the cluster stellar mass fraction, as defined by present-day clusters, is highest in lower-mass galaxies, addressing whether these galaxies also contain the bulk of clusters is complicated by having to account for the relative numbers of galaxies of different stellar masses. The cluster distribution function among galaxies of different M_* is given by the combination of T_N and ϕ(M_*), where ϕ(M_*) represents the volume density of galaxies with M_*. Taking parameterized expressions for each of these into account, adopting the stellar mass function measured from the GAMA survey by Baldry et al. (2012), we calculate the globular cluster distribution function shown in Figure 6. This figure clearly shows that the bulk of today's clusters are found in galaxies with M_* ∼ 10^10.8M_⊙, which are relatively massive galaxies. The greater numbers of less-massive galaxies combined with the somewhat larger T_N values for those galaxies were insufficient to counterbalance the fact that the more-massive galaxies also have a larger absolute number of clusters. To enable further calculations with this globular cluster distribution function, G(M_*), we provide the following fitting function, which is also shown in Figure 6:

$$\begin{eqnarray}G({M}_{*}) & = & 0.9\mathrm{exp}({-({M}_{*}/{10}^{11.3})}^{2.3})\\ & & \times ({(\mathrm{log}{M}_{*}-8.4)}^{2.8}+6.64).\end{eqnarray} \tag{ 2 }$$

The parameters (normalizations, exponents, and constants) in Equation (2) where determined by minimizing χ² for the selected functional form. We did not explore a wide range of functional forms, so the equation is simply meant to be a convenient fitting function over the range of stellar masses plotted in Figure 6, and we do not ascribe physical meaning to the function or the fitted parameters.

Zoom In Zoom Out Reset image size

There are two caveats to this result. First, we do not measure the cluster populations of the most-massive galaxies ourselves. We have adopted a functional form for T_N at the highest masses that is consistent with the Peng et al. (2008) data. While this approach may not be appropriate for galaxies outside of clusters, the most-massive galaxies tend to be found mostly in clusters. Furthermore, as discussed in Paper I, the possibility that the stellar initial mass function varies systematically among early-type galaxies could account for the entire observed rise in T_N at these masses. Fortunately for our calculation, these massive galaxies are exceedingly rare, so significant variations in the behavior of T_N in this mass regime have modest effects on the globular cluster distribution function. To demonstrate that the effects are minor, we recalculate G(M_*), adopting instead a constant T_N for M_* > 10^10.5M_⊙, and show the difference with our previous calculation in Figure 6. Second, we have not extrapolated our T_N relation to M_* < 10^8.5M_⊙ because we have little data at those masses. The Baldry et al. (2012) galaxy stellar mass function also does not extend below M_* = 10⁸M_⊙. However, if we extend the T_N relation down to M_* = 10⁸M_⊙, then we do find the cluster numbers rising at low M_*, although slowly. If T_N is grossly larger for low-mass galaxies, as hinted at by the three galaxies in our sample at lower M_*, then a significant number of clusters could be hosted by such galaxies.

The scatter in our measurements of T_N is large, although not manifestly larger than that plausibly attributable to measurement uncertainties (Figure 4). However, we might still be able to uncover physical sources of scatter in T_N if those contribute a comparably large level of apparent scatter. To investigate, we correlate various characteristics of our galaxies with deviations, ΔT_N, from the mean relation between T_N and M_*, within the stellar mass range for which that relationship is best defined, $8.5\lt \mathrm{log}({M}_{*}/{M}_{\odot })\lt 10.5$.

The first characteristic we explore is the galaxy's large-scale environment as measured using the cosmic shear field (Courtois & Tully 2015; Courtois et al. 2015). We have previously used this measurement to investigate the nature of extremely massive and gas-rich galaxies within the S⁴G sample (Courtois et al. 2015). This environment measurement is a coarse one, indicating only whether the galaxy lies in a void, filament, sheet, or knot (designations 0, 1, 2, and 3, respectively). We find no significant correlation between ΔT_N and environment, although the scatter among T_N values appears to be smaller in the densest environments, the knots (Figure 7). This visual impression is supported by the results of a statistical F-test to determine the likelihood that any two sets of data were drawn from a Gaussian distribution with the same variance. Comparing the filament and sheet galaxies with those in the knots results in probabilities that they were drawn from the same parent sample of 0.01 and 7 × 10⁻⁴, respectively, suggesting that the variances of the underlying populations are actually different. The one strong caveat for this test is that it is highly sensitive to non-Gaussianity in the underlying distributions. This result merits attention with larger samples, but for now we conclude that we find no clear sign of large-scale environmental effects on the specific frequency of globular clusters for galaxies in this intermediate-mass range.

Zoom In Zoom Out Reset image size

Next we explore the relations between ΔT_N and measurements of the galaxies' baryonic and dark-matter content. The rationale for such explorations rests with questions regarding the relative efficiency of cluster formation (and destruction) with other measurements related to the star-formation efficiency, such as the condensed baryon fraction (the fraction of the universally apportioned baryons within a dark-matter halo that cool sufficiently below the halo virial temperature that they can be detected either as stars and cool gas, f_C) and the star-formation efficiency as quantified by the ratio of stellar mass to cool gas mass, M_G, or as the stellar mass to dark-matter mass, M_D. We have the necessary measurements of M_G and rotation speeds, v_C, for a subset of our galaxies, mostly among the late types because we are utilizing the H i compilation of Courtois et al. (2011). Both M_G and M_D are calculated as described in detail by Zaritsky et al. (2014). We present these data in Figure 8 and the statistical significance of the resulting correlations between ΔT_N and these quantities (calculated by determining the probability that a random set of galaxies would have the same or larger Spearman rank correlation coefficient) in Table 2.

Zoom In Zoom Out Reset image size

There is no correlation detected between f_C and ΔT_N, and there marginal correlations (at or slightly below 2σ significance) between ΔT_N and the two measurements of star-formation efficiency. With either measurement of star-formation efficiency, the putative trend suggests that as the overall star-formation efficiency increases (either relative to cold gas content or dark matter) the relative efficiency of cluster formation also increases. If a galaxy (of a given current-day stellar mass) has been more efficient at turning its baryons into stars, it has also been more efficient at forming (or retaining) its clusters.

It is always difficult to interpret correlations, particularly among parameters that contain a measurement in common. In this case the difficulty is that all of the quantities under discussion depend on M_*. To examine the sense of the effect expected in terms of parameter correlations should M_* be incorrectly measured, consider that if M_* is incorrectly overestimated (underestimated), then T_N will be underestimated (overestimated) and the star-formation efficiencies will be overestimated (underestimated). The result of this behavior is that one would find cluster efficiency going in the opposite sense as star-formation efficiency, opposite in sense to that observed. This argument suggests that our results are not due to correlated errors among the parameters arising from their common use of M_*. In fact, errors in M_* may be weakening a stronger underlying correlation of the sense we observe.

Even so, the statistical significance of the measured correlations is modest. To test these correlations further, we examine the behavior of our early-type galaxies. For most of these galaxies, we have no H i measurement and hence also no measurement of v_C, which is why these galaxies are absent from Figure 8. If we presume that the gas content is a small fraction of the stellar content, we can place these galaxies at the right-hand side of a revised version of the rightmost panel of Figure 8, which we do in Figure 9. Consistent with the expectation from the suggested correlation, most of these galaxies have ΔT_N > 0. We conclude that the existence of the two correlations relative to different estimators of star-formation efficiency plus the properties of the early-type galaxies all support the suggestion that T_N rises relative to the mean for galaxies that have converted a larger fraction of their baryons into stars.

Zoom In Zoom Out Reset image size

Finally, we address the question of whether a more fundamental specific frequency is obtained when one normalizes relative to stellar mass or total mass. The latter normalization could be appropriate not necessarily because cluster formation somehow involves dark matter but because it depends on a quantity that scales more closely to dark matter than to stellar mass. One likely such candidate quantity is the baryonic mass of a galaxy. Harris et al. (2013) explored this question and found that using a dynamical estimator of halo mass did indeed lead to a tighter relationship between specific frequency and galaxy properties. With our sample we can only do this test with the limited subset of galaxies for which we have a measured v_C. We find that our correlation is stronger when using M_* (a Spearman rank correlation probability of arising randomly of 0.029) rather than halo mass (comparable probability of 0.11). We conclude that we do not find support in our data for expressing specific frequency as a function of total mass, but we note that we have a limited sample of only 29 galaxies currently with which we can do this test. This question merits further investigation.