THE SAGES LEGACY UNIFYING GLOBULARS AND GALAXIES SURVEY (SLUGGS): SAMPLE DEFINITION, METHODS, AND INITIAL RESULTS

Into this void, we have launched the SAGES Legacy Unifying Globulars and GalaxieS (SLUGGS)⁷ survey. The survey uses a combination of imaging from Subaru/Suprime-Cam and spectroscopy from Keck/DEIMOS to obtain the spatial, kinematical, and chemical properties of both diffuse stellar light and GCs in two dimensions and over wide areas around 25 nearby ETGs spanning a broad range of galaxy properties and environments.

Some aspects of the survey, such as the GC imaging and spectroscopy, build on established observational techniques and apply them with great success using state-of-the-art instrumentation. Other aspects are novel and potentially transformative for the study of nearby galaxies. These involve Stellar Kinematics from Multiple Slits (SKiMS), a technique that we have pioneered to map out the kinematics and metallicities of diffuse starlight in two dimensions (2D) to large radii (e.g., Norris et al. 2008; Proctor et al. 2008, 2009).

All of the spectroscopic aspects of the survey exploit the strong Ca ii triplet absorption line feature in the near-infrared (∼8500 Å), where DEIMOS has both a high efficiency and a high spectral resolution for separating the galaxy and GC lines from the forest of near-infrared sky lines.

The main advantages of the SLUGGS survey are the superior velocity precision (median 12  km s⁻¹) on a very stable, high-throughput instrument, and the large radial coverage (typically ∼2–3 R_e for stars, extending to ∼8 R_e when GCs are included).

In Figure 2 and Table 1 we have summarized other current and planned 2D chemodynamical surveys of ETGs. These surveys use either an Integral Field Unit (IFU) or a fiber-bundle. We briefly discuss these in turn.

Zoom In Zoom Out Reset image size

Building on the success of the original SAURON survey, the ATLAS^3D survey (Cappellari et al. 2011) of 260 nearby ETGs has had a high scientific impact. However, it has a moderate velocity resolution of 105  km s⁻¹ (limiting its ability to map both low-mass galaxies and kinematically cold features in all galaxies) and, perhaps most importantly, a small field-of-view. To compensate for this, galaxies were often targeted with two pointings, giving a median extent for the survey observations of 0.9 R_e. For more luminous and nearby galaxies, comparable to those in SLUGGS, the typical extent is only 0.6 R_e (see Figure 1 in Arnold et al. 2014).

A recently initiated survey comprising all galaxy types using the SAMI instrument (Bryant et al. 2014) will increase the number of ETGs mapped with higher velocity resolution and wider wavelength coverage, but will not reach much beyond 1 R_e. As noted earlier, most of the mass and angular momentum in ETGs is located beyond 1 R_e, where we also expect key signatures of galaxy formation.

A number of surveys are probing beyond 1 R_e. CALIFA (Sánchez et al. 2012) is well underway and, upon completion (expected in 2015), around 80 ETGs will have been observed out to about 2 R_e with a velocity resolution of 69  km s⁻¹ (blue arm).

Perhaps most comparable to SLUGGS is the study of 50 nearby ETGs that uses the Mitchell Spectrograph IFU (formerly VIRUS-P; Greene et al. 2013; Raskutti et al. 2014). This study operates in the optical and so is capable of obtaining both age and metallicity information, as well as kinematics, out to ∼2.5 R_e. The drawbacks are the relatively poor velocity resolution (150 km s⁻¹) and the lower signal-to-noise ratio (S/N) at large radii, which inhibits the 2D coverage. A follow-up survey, MASSIVE, is focusing on galaxies at the extreme high end of the stellar mass range (Ma et al. 2014).

Due to start in 2014 and collect data for six years, the MaNGA survey (PI Kevin Bundy, http://www.sdss3.org/future/manga.php) will map an impressive 10,000 galaxies, roughly half of which will be ETGs. It will cover UV to near-IR wavelengths (3600–10,000 Å) with a velocity resolution of 64 km s⁻¹. Although it will reach only about 1.5 R_e for the primary sample, a secondary sample of ∼1500 will be observed to 2.5 R_e.

Partially overlapping with SLUGGS both in theme and in galaxy sample is the Next Generation Virgo Cluster Survey (NGVS), a multi-band photometric study of the Virgo cluster with an emphasis on GCs (Ferrarese et al. 2012). However, NGVS does not incorporate a full, systematic spectroscopic component, which is the main focus of our survey, along with an exploration of galaxies in a wide range of environments.

The SLUGGS survey has been designed to address the following science themes: What are the basic, global chemodynamical properties of ETGs? (including their masses, angular momenta, orbital structures, and metallicities over a wide range in radius). What is the distribution of dark matter in ETGs? How are the outer regions of ETGs assembled? How does ETG assembly depend on mass, environment, and other variables? Do the observations agree with theoretical models of galaxy formation?

SLUGGS began as a pilot project in 2006 November, using DEIMOS to measure velocities of GCs around the giant elliptical NGC 1407 (Romanowsky et al. 2009). Since then it has evolved into a full survey, leading to over 30 papers to date (many of them discussed below), with the first large installments of data coming in Usher et al. (2012) on GC metallicities, Pota et al. (2013) on GC kinematics, Arnold et al. (2014) on stellar kinematics, and Pastorello et al. (2014) on stellar metallicities.

As a preview of the nature and quality of the SLUGGS results, Figure 3 shows the observed 2D stellar velocity field of a galaxy, zooming out to show different types of data on three spatial scales. The inner regions of the data come from the well-known SAURON integral field spectrograph, the intermediate regions are from SKiMS, and the outer regions come from GCs. The SKiMS data show initially the same pattern as SAURON out to 1 R_e, with a disklike rotation pattern that is well-aligned with the isophotes. However, by 2 R_e, there are indications of a kinematic twist that becomes even more prominent in the GC data at ∼5 R_e where there are additional asymmetries that suggest unrelaxed post-merger material. Thus at large radii, we see the emergence of new information about the galaxy structure and formation that is not apparent from the central regions alone.

Zoom In Zoom Out Reset image size

This paper provides a general overview of SLUGGS, including the galaxy sample definition (Section 2) and summaries of our observing and data reduction techniques (Section 3), and modeling methods (Section 4). Some initial results are presented in Section 5 and conclusions are provided in Section 6. An Appendix provides further details about the galaxy sample. More information and news can be found on the survey Web site.⁸

3.1.1. Background

3.1.2. Imaging Acquisition and Reduction

Imaging is carried out in three filters for each galaxy: g, r and i (except for a few galaxies with adequate archival data available). Our nominal target for each filter is S/N ∼20 at one magnitude fainter than the turnover of the GC luminosity function (occurring at M_i ∼ −8). The total exposure times depend on the band, galaxy distance, and observing conditions, ranging from ∼300 s for the r-band in the nearest galaxies in the sample, up to ∼4000 s for the g-band in the most distant. We use a five-point dither pattern to fill in the chip gaps. We also take very brief, ∼10 s exposures in order to obtain unsaturated images of the galaxy centers and of any bright UCDs.

Our reduction of Suprime-Cam imaging is largely based on the SDFRED package (Yagi et al. 2002; Ouchi et al. 2004), but incorporates custom modifications to improve the sky subtraction and alignment of images. High-quality astrometric solutions across the Suprime-Cam field of view are essential for spectroscopic follow up. In order of preference, astrometric calibration is linked to the SDSS DR7 (Abazajian et al. 2009), USNO-B2 (Monet et al. 2003), or 2MASS (Skrutskie et al. 2006), with typical rms values of ∼015. An example of a reduced Suprime-Cam image is provided in Figure 7.

Zoom In Zoom Out Reset image size

3.1.3. Globular Cluster Photometry and Selection

Object detection and photometry are performed using the DAOPHOT suite (Stetson 1992). All photometry is aperture photometry, using optimal extraction radii and aperture corrections calculated for each set of images. Photometric calibration is achieved using standards taken during the run, or by reference to point sources in SDSS DR7. Corrections for foreground dust extinction are applied using standard maps (Schlegel et al. 1998).

Significantly extended objects are removed as probable background galaxies (slightly resolved objects could be UCDs and are not excluded). The GC sequence is defined in g−i versus r−i color–color space (see example in Figure 8 and cf. Sinnott et al. 2010). This sequence is fairly well-separated from the foreground-star sequence, except at the bluest end (where near-infrared imaging would help in principle; Muñoz et al. 2014). We thereby obtain a sample of GC candidates that is relatively contaminant-free, as is important both for spectroscopic target selection, and for general photometric studies of the GC systems. The bright magnitude GC selection limit is set by the expected range of UCDs (up to M_i ∼ −13; Brodie et al. 2011); the faint limit varies according to the actual depth of the imaging.

Zoom In Zoom Out Reset image size

Detection, secure identification, and photometric analysis of GCs from ground-based imaging become challenging near the centers of bright galaxies. Fortunately, almost all of the SLUGGS galaxies have archival imaging available from HST, allowing for improved GC detection and photometry closer to the galactic centers. The HST analysis is similar to that for ground-based imaging, although only one color is typically available (see Spitler et al. 2006; Strader et al. 2006). The principal difference is that the vast majority of the GCs in our survey are resolvable by HST, which allows for another step in contaminant reduction (eliminating the foreground stars in the central regions of the galaxies). We also use information from HST to fine-tune the parameters used in the Suprime-Cam GC selection.

Once all of the data have been processed, full photometric catalogs and imaging data for the SLUGGS galaxies will be made publicly available on the project Web site.⁹

3.1.4. Galaxy Surface Photometry

Our wide-field imaging is also suitable for studying the host galaxy, although this is currently the least advanced aspect of our survey. Suprime-Cam is not optimally suited for mapping out the extended low surface brightness regions of galaxies, owing to limitations with flat-fielding and scattered light that are typical for most telescopes, and require special instrument configurations or heroic processing efforts to circumvent (e.g., Tal et al. 2009; Janowiecki et al. 2010; Paudel et al. 2013). However, our survey imaging still reaches fainter surface brightness levels than SDSS, which is important for studying the interesting outer regions of the galaxies.

Using the standard IRAF task ELLIPSE, we measure the isophotes in the annuli, thus building up a 2D description of the galaxy light (see Blom et al. 2012a). The resulting data include the radial profiles of surface brightness, ellipticity, position angle, and higher order Fourier components (such as boxiness and diskiness). These profiles are used throughout the survey in various contexts, such as in dynamical models, derivations of galaxy color profiles, and in comparisons of GC system properties (spatial and kinematical) to the diffuse galaxy light.

The extended stellar light can also be inspected for substructures produced by galaxy interactions and mergers. Although our sample was designed to avoid galaxies where such features are dominant, we have serendipitously discovered several halo substructures in our Suprime-Cam data.

Figure 9 shows one example, presented here for the first time. An S0 galaxy, NGC 4111, is seen to be linked by a faint tidal stream to a smaller spiral galaxy, UGC 7094 (=PGC 38375). The UGC 7094 stream has a maximum surface brightness of μ_g ∼ 26.4 mag arcsec⁻², spans an apparent length of ∼75 kpc, has a slightly redder color than its parent galaxy, (g − i)₀ ∼ 0.93 versus ∼0.81, and a similar luminosity (L_g ∼ 2 × 10⁸ L_{g, ☉} versus ∼3 × 10⁸ L_{g, ☉}).

Zoom In Zoom Out Reset image size

The stream thus appears to represent an ongoing minor merger, with a stellar mass ratio of ∼1:20 (based on i-band luminosities). Another S0 galaxy toward the north east, NGC 4117, shows distorted outer features (not shown in the figure), suggesting that we may be witnessing multiple interaction events happening within a small group of galaxies as it falls into the Virgo cluster.

The Suprime-Cam imaging also allows us to do far more than characterize substructures photometrically; we can identify embedded GCs for spectroscopic follow up with DEIMOS, allowing the chemodynamics of the features to be explored. So far, we have used this approach for stellar streams around two galaxies in the survey, M87 and NGC 4365 (Romanowsky et al. 2012; Blom et al. 2014). These data allowed us to identify, for the first time, NGC 4342 as the origin of a stream in the NGC 4365 system.

The number of discrete spectroscopic tracers in substructure can be increased beyond GCs by acquiring supplementary narrow-band imaging with Suprime-Cam to identify PNe, which we then follow up with DEIMOS (Foster et al. 2014).

3.2.1. DEIMOS Mask Design and Observations

3.2.2. Globular Cluster Spectroscopy

3.2.3. Stellar Spectroscopy

The spec2d reduction pipeline generates the subtracted background spectra (loosely referred to as "sky"). However, because our point-like targets are projected onto the host galaxy, these background spectra are a combination of pure sky and galaxy background light. We use the SKiMS technique (Norris et al. 2008; Proctor et al. 2009; Foster et al. 2009; Schuberth et al. 2012; Richtler et al. 2014) to extract the galaxy spectrum from these background spectra. Briefly, we use a specially defined "sky index" to identify pure-sky slits at large galactocentric distances (typically 6–7 R_e), where the contribution of galaxy light is insignificant.

In earlier SLUGGS work, the pure-sky spectra were co-added and normalized into a single master sky spectrum that was then scaled for non-local sky subtraction for each slit (Proctor et al. 2009; Foster et al. 2009, 2011). In order to minimize the impact of small variations of the spectral resolution across the spatial axis of the DEIMOS mask, in recent SLUGGS work we have used the penalized Pixel Cross-correlation Fitting (pPXF) routine by Cappellari & Emsellem (2004) to derive a linear combination of pure-sky spectra that best matches the sky component in each slit (see Arnold et al. 2014).

This sky spectrum, individually optimized for each slit, is then subtracted from the corresponding slit, and the remaining signal is the galaxy stellar-light spectrum. Because these SKiMS spectra are spread in two dimensions around the galaxy, in accordance with the GC targets, we can thereby effectively use DEIMOS as a sparsely sampled integral-field unit.

This method yields a final continuum level in the sky index definition region accurate to <1% of the noise in the skyline residuals (see Foster et al. 2009). As an example, Figure 10 illustrates the sky index definition and sky-subtracted spectrum for the galaxy light of NGC 4494. Using this procedure, we obtain galaxy spectra that typically extend to 2–3 R_e, with some reaching 4 R_e.

Zoom In Zoom Out Reset image size

Stellar kinematics are measured using pPXF. Briefly, the code uses the fully propagated variance array to proportionally "penalize" pixels in the spectrum. Template spectra are convolved with Gauss–Hermite polynomials and compared to the observed spectra. The routine identifies the set of kinematic parameters and templates that minimizes residuals between the observed and combined template spectrum. We simultaneously fit the first four velocity moments (i.e., recession velocity, velocity dispersion, higher-order Gauss–Hermite moments h₃ and h₄). Uncertainties are measured using Monte Carlo methods. Our kinematic measurements are described in detail in Arnold et al. (2014). Figure 10 shows an example SKiMS spectrum together with its best fit pPXF output template spectrum.

For the highest S/N spectra, we measure the strength of the CaT index for the diffuse stellar light to large radii using the technique described in Foster et al. (2009). The index is the weighted sum of the equivalent width of all three CaT lines. A correction is applied to account for the velocity dispersion broadening. CaT indices are then converted into metallicity as described in Section 4.2.3. In this manner, we are able to generate metallicity maps out to typically ∼1.5–2 R_e, with some reaching 3 R_e (not quite as far as the kinematics measurements owing to the greater S/N requirements; for more detail, see Foster et al. 2009, 2011; Pastorello et al. 2014).

4.1.1. GC Numbers and Spatial Distributions

Accurate GC number counts, including reliable uncertainties, are important for interpretations of their formation histories (e.g., Bekki et al. 2008; Coenda et al. 2009), estimates of stellar halo masses (Brodie & Strader 2006), and for analyses of their correlations with interesting galaxy properties such as supermassive black hole masses and dark matter halos (e.g., Spitler & Forbes 2009; Burkert & Tremaine 2010; Harris & Harris 2011; Snyder et al. 2011; Rhode 2012; Harris et al. 2013, 2014; Wu & Kroupa 2013; Forte et al. 2014; Durrell et al. 2014). One of the goals of the SLUGGS survey is to provide an updated set of reliable GC numbers (including the red and blue subpopulations separately), and then to explore what impact these have on the scaling relations. As described below, some of the literature values of N_GC are in error by factors of three or more, and we anticipate that our more accurate estimates will clarify these fundamental relations.

Although every care is taken to ensure the GC catalogs from photometric data contain very little contamination (from foreground stars and background galaxies; see Section 3.1.3), a certain amount is unavoidable and will impact measurements of the true GC number counts. Fortunately, the contaminants contribute a fairly constant surface number density over the regions covered by the GC systems. This surface density is typically estimated as a free parameter in fits to the GC surface density profiles, and can be cross-checked via spectroscopic contamination rates (Strader et al. 2011; Usher et al. 2013).

To derive total GC numbers around a galaxy, the surface density profile is integrated with the radius. Correction must also be made for the fraction of GCs that were too faint to be included in the catalogs. Since the functional form of the GC system luminosity function is predictable (Jordán et al. 2007), the missed GCs can be accounted for by measuring the brighter parts of the luminosity function and extrapolating to faint magnitudes (which can be done securely because the SLUGGS images are designed to reach past the turnover luminosity). Corrections are thus made for completeness and contamination as a function of magnitude.

As an illustration of the importance of reliable surface density distributions for deriving accurate estimates of N_GC, Figure 11 shows the surface density and cumulative number profiles of blue and red GCs with radius around NGC 1407, based on Suprime-Cam imaging. The cumulative profiles flatten out at a radius of ∼14 arcmin (∼110 kpc), corresponding to the edge of the Suprime-Cam image. Thus, with very high quality, wide-field imaging, one can just barely get a reliable count of total GC numbers around a high-mass elliptical galaxy at 27 Mpc.

Zoom In Zoom Out Reset image size

For poorer-quality data, narrower fields, or more nearby massive galaxies, the total GC counts can become precarious. For NGC 1407, a single HST/Advanced Camera for Surveys (ACS) pointing picks up only ∼10% of the total GCs (Forbes et al. 2006), and ground-based imagers with ∼10–20 arcmin fields would subtend only half of them. Naturally such limitations are obvious, and attention has been given in the literature to making reasonable extrapolations while accounting for the associated uncertainties, but the results until now have not been encouraging. As a case in point, the giant elliptical NGC 4365 was analyzed with a central ACS pointing supplemented by parallel WFPC2 imaging, and estimated to have 3246 ± 598 GCs by Peng et al. (2008). Our subsequent analysis of Suprime-Cam imaging found a lower-limit of 6450 ± 110 GCs (Blom et al. 2012b), which is higher than the previous result by at least 5σ.

Besides total numbers, our Subaru data allow accurate surface density distributions to be derived for the SLUGGS galaxies out to large radii. We can improve on past procedures for determining important basic parameters for these galaxies, such as R_e, by fitting Sérsic models both to the galaxies and to their GC systems (see, e.g., Blom et al. 2012b; Usher et al. 2013; Kartha et al. 2014, for GCs).

On a related note, the analysis of the radial "extent" of GC systems (e.g., Rhode et al. 2007) is also sensitive to limitations with the imaging data. For example, analysis of the SLUGGS data for NGC 720 by Kartha et al. (2014) revealed that its GC system extent had previously been underestimated by a factor of three.

4.1.2. GC Kinematics

The kinematics of a GC system can be characterized by the rotation amplitude v_rot along a position angle θ₀, and by the rotation-subtracted velocity dispersion σ_p. These parameters quantify the amounts of coherent rotation and random motion of the system, and they can be compared to the same quantities for stellar kinematics. We study these quantities using a maximum-likelihood approach (Kissler-Patig & Gebhardt 1998) that we have extended to allow for flattened galaxies (see details in Pota et al. 2013 and discussion in Romanowsky et al. 2009 about problems with other methods). The approach is equivalent to a simplified version of kinemetry used for galaxy stellar-light kinematics (see Section 4.2.2).

One particular kinematic measurement of interest is the projected rms velocity v_rms which, as a metric of specific kinetic energy, can be connected through dynamical modeling to estimates of mass. Figure 12 shows v_rms profiles for two example SLUGGS galaxies, one high-mass (NGC 1407) and one low-mass (NGC 3377), using SLUGGS data from Pota et al. (2013). For each galaxy, the data are separated into profiles for the blue and red GC subpopulations. The case of NGC 3377 illustrates the generic result that the blue GCs have higher rms velocities than the red GCs, as expected from the larger spatial extent of the blue GC system density profile (which must be supported by hotter kinematics in the same gravitational potential). In some cases like NGC 1407, the radial behaviors of the v_rms profiles are also different, which probably reflects very different orbital types between the blue and red GCs.

Zoom In Zoom Out Reset image size

Interpreting these kinematic data through dynamics, in terms of orbital properties and mass profiles, will be discussed further in Section 4.3. It should be noted however that the observed v_rms profiles are not generally smooth, which is a complication that should ultimately be addressed by allowing for strong orbital variations in the models. In addition, some of the observed fluctuations may be produced by substructures, thereby providing valuable clues about recent merging events (cf. Romanowsky et al. 2012; Coccato et al. 2013; Schauer et al. 2014; Murphy et al. 2014).

4.1.3. GC Metallicities

The metallicities of the GCs are derived using the CaT index (Section 3.2.2), whose strength has previously been shown to be sensitive to changes in mean metallicity (Armandroff & Zinn 1988; Diaz et al. 1989; Cenarro et al. 2001; Foster et al. 2009, 2010). CaT is independent of age for stellar populations older than 3 Gyr (Vazdekis et al. 2003). Although the index is sensitive to a bottom heavy initial mass function (IMF; Vazdekis et al. 2003), our CaT measurements are insensitive to horizontal branch morphology and to α element enhancement (Brodie et al. 2012).

We have compared CaT metallicities with literature values derived using traditional Lick Index analyses, finding excellent agreement between the two techniques (with a reduced χ² value of 0.59; Usher et al. 2012). We use the single stellar population models of Vazdekis et al. (2003) to derive the following relation for converting the CaT index strengths into metallicities:

$$\begin{equation} [Z/{\rm H}] = (0.461 \pm 0.013) \, {\rm CaT} + (-3.750 \pm 0.080) \,. \end{equation} \tag{ 1 }$$

Further details of the CaT metallicity measurement techniques are given in Usher et al. (2012).

4.2.1. Stellar Density Distributions

4.2.2. Stellar Kinematics

Once the SKiMS-based stellar kinematics measurements are obtained (Section 3.2.3), the overall kinematic structure of the galaxy can be visualized through 2D maps. Several methods have been used for preparing these maps, including Voronoi binning with semi-variogram-based smoothing (Arnold et al. 2014), and kriging interpolation (Krige 1951; Matheron 1963; Cressie 1988; Foster et al. 2013). An example of the latter is shown by Figure 13.

Zoom In Zoom Out Reset image size

As with the GCs (Section 4.1.2), the full 2D kinematical information for the stellar light can be reduced to one-dimensional profiles through modeling. We adopt a "kinemetry" approach, originally designed by Krajnović et al. (2006) for IFU data-cubes, and extended to sparsely sampled data by Proctor et al. (2009) and Foster et al. (2011). Kinemetry is analogous to ellipse-based galaxy surface photometry and assumes that the rotation and the dispersion are stratified on ellipses. These profiles generally include the rotation amplitude, position angle, and "flattening" or ellipticity (which is often, but not always, similar to the photometric ellipticity; cf. Krajnović et al. 2008). The kinemetry results for the survey will be presented in Foster et al. (2015).

4.2.3. Stellar Metallicities

As discussed in Section 3.2.3, CaT indices can be measured for the integrated starlight out to radii of typically ∼1.5–2 R_e. These values are converted to total metallicities [Z/H] using Equation (1), but with an empirical correction that takes into account the observed dependence of CaT-derived metallicities on galaxy mass. Pastorello et al. (2014) suggested that this dependence may be a reflection of the index's sensitivity to the slope of the IMF, and derived the correction by comparing our values to metallicities from overlapping SAURON data (Kuntschner et al. 2010).

The uncertainties are evaluated using Monte Carlo resampling of the spectra, in a method similar to the one adopted for GC metallicities (Section 4.1.3).

In order to obtain smooth 2D stellar metallicity maps for the galaxies in our sample, we interpolate between our data points with the kriging technique (Pastorello et al. 2014). An example of a 2D stellar metallicity map derived in this way is shown in Figure 14 for NGC 4111.

Zoom In Zoom Out Reset image size

From these 2D maps we extract azimuthally averaged radial metallicity profiles, from which we measure the inner (R < 1 R_e) and outer (R > 1 R_e) metallicity gradients. These gradients are powerful tools for constraining galaxy formation models and feedback processes (e.g., Mihos et al. 2013; Hirschmann et al. 2015). For example, galaxies with more active merger histories should have fairly flat metallicity profiles at large radii (White 1980), while galaxies dominated by in situ star formation should have steep gradients, with their chemical enrichment being dependent on the distance from the center of the gravitational potential.

One of the main science drivers for SLUGGS is to use stars and GCs as dynamical tracers in the outer regions of ETGs and thereby to map out their dark matter halos, since H i rotation curves are not generally available in these galaxies. Although this approach, along with the use of PNe, has been pursued for decades, robust conclusions have remained elusive. This is in part because of modeling uncertainties and degeneracies (e.g., de Lorenzi et al. 2009), but also because the lack of a large sample of galaxies with homogeneous data has limited the statistical weight of the conclusions.

Furthermore, the galaxies previously well-studied with GC dynamics were limited to the high-luminosity ellipticals at the centers of galaxy groups and clusters, while it is only with the advent of more efficient spectroscopy that the realm of ∼L* early-types can be investigated properly (see Figure 6). This is an important area of parameter space not only because such galaxies are more "ordinary" and allow for fairer comparisons with ∼L* spirals such as the Milky Way (e.g., Trujillo-Gomez et al. 2011), but also because of reports based on PNe that the mass profiles of these galaxies are systematically different from those of the brightest ellipticals, with much less dominant dark matter halos (Romanowsky et al. 2003; Napolitano et al. 2011; Deason et al. 2012).

In order to focus the discussion about the limitations and opportunities for dynamical modeling, we construct a set of toy galaxy models and present the circular velocity and projected velocity dispersion profiles in Figure 15. These models consist of spherical ETGs, with either 4 L* or L* luminosities, typical GC systems, and either standard Λ cold dark matter (ΛCDM) halos (Navarro et al. 1997)¹⁰ or dark matter halos with constant-density cores,¹¹ as might arise either from baryonic effects or from alternative types of nonbaryonic dark matter (see discussion in Newman et al. 2013a).

Zoom In Zoom Out Reset image size

The top panel of Figure 15 illustrates some aspects of the mass distribution that may be learned through dynamical constraints that extend from the centers of galaxies to their outer parts, as in SLUGGS. Within the central few kpc of the galaxy, the circular velocity profile is hardly affected by large variations in the dark matter halo properties. This is because the stars dominate the central regions, which is an advantage for constraining the stellar mass (and IMF) using detailed data and models in these regions, but a disadvantage for learning about the dark matter distribution (e.g., Treu et al. 2010; Tortora et al. 2010; Cappellari et al. 2012).

It is clearly important to employ data extending to larger radii to constrain the dark matter halos. Here our SKiMS data have great potential, owing to their 2D coverage out to ∼3 R_e (∼10 kpc), where the dark matter fraction is expected to be around 50%. These data, combined with central kinematics from SAURON/ATLAS^3D, may provide a unique opportunity to measure the slope of the central dark matter profile—which has so far been determined for dwarf, spiral, and brightest cluster galaxies (e.g., Oh et al. 2011; Chemin et al. 2011; Walker & Peñarrubia 2011; Newman et al. 2013a).

The GC kinematics data extend even farther out, to ∼8 R_e (tens of kpc), which is close to the expected dark matter halo scale radius. This is an important region to probe since it lies beyond the galaxy centers where myriad baryonic effects could have altered the dark matter profiles. At these large radii, it becomes possible to test the primordial mass–concentration relation that is a robust prediction of ΛCDM theory, at least for a given set of cosmological parameters (e.g., Dutton & Macciò 2014).

While the SLUGGS stellar and GC kinematics data sets are each expected to provide powerful constraints on their own, the full potential of the survey will be realized when all of the data are integrated into "global" dynamical models of the galaxies with large baselines in log-radius. In particular, each galaxy has between two and four tracer populations (stars and one, two or three GC subpopulations, depending on how many distinct subpopulations are present), each of them providing independent constraints on the mass profile for a distinct radial regime. This is an established approach for dwarf spheroidal galaxies with multiple stellar populations (e.g., Walker & Peñarrubia 2011; Amorisco et al. 2013), and can now be applied to GC systems with large data sets (Schuberth et al. 2010; Napolitano et al. 2014; Agnello et al. 2014b; V. Pota et al., in preparation).

Like GCs, PNe are also used as large-radius dynamical tracers. Here the approach is somewhat simpler because the spatial distribution of PNe can be assumed to follow the field starlight (see discussion in Courteau et al. 2014). However, GCs have the potent advantage of providing additional leverage through their subpopulations, with the blue GCs in particular extending to larger radii than the PNe (and red GCs; see bottom panels of Figure 15). Given the extensive overlap between SLUGGS and PN-studied galaxies (see Section 2) we plan to try incorporating joint GC and PN modeling constraints (cf. Forbes et al. 2012).

Most of our sample galaxies furthermore have detectable X-ray emitting hot gas halos (e.g., O'Sullivan et al. 2001; Boroson et al. 2011), some of which are suitable for mass analysis. We have started to compare optically based galaxy masses with those from X-rays and hence to test the assumptions inherent in both methods. So far we have found in two galaxies that the mass estimates in the outer regions (∼20–30 kpc) as derived from the hot gas were much greater than those from the GCs (and PNe; Romanowsky et al. 2009; Napolitano et al. 2014)—which presumably reflects departures from hydrostatic equilibrium. More generally, the mass estimates from the SLUGGS galaxies should prove useful in understanding the properties of the interstellar medium (e.g., Humphrey et al. 2013; Kim & Fabbiano 2013).

Given the above science goals and available data, there are a number of challenges for dynamical modeling. These include proper treatment of galaxy flattening, permitting and constraining orbital anisotropy, and uncertainties in the stellar mass (from IMF variations in particular). To tackle these issues, we are pursuing a layered modeling approach, ranging from simple mass estimators that may be applied quickly to many galaxies (Strader et al. 2011; Agnello et al. 2014a) to detailed, slower Schwarzschild-type orbit models (Romanowsky & Kochanek 2001; Gebhardt & Thomas 2009). We are also using techniques of intermediate complexity such as spherical Jeans models (Romanowsky et al. 2009; Napolitano et al. 2014), flattened Jeans models (M. Cappellari 2008, and in preparation), the projected virial theorem (Agnello et al. 2014b), and distribution function methods (Deason et al. 2012; Mamon et al. 2013). It should be noted that many of these mass-modeling methods are designed to simultaneously constrain the orbital anisotropy, e.g., by making use of higher-order information in the line-of-sight velocity distributions.

So far, the SLUGGS results support the dominance of a two-phase mode of galaxy assembly, as described in Section 1. The principal pieces of evidence include: broken radial metallicity gradients in the GC subpopulations (Arnold et al. 2011; Forbes et al. 2011), distinct core and halo kinematic behaviors of both the GCs and the underlying starlight (Arnold et al. 2011, 2014; Pota et al. 2013; Foster et al. 2013), and the presence of stellar streams, kinematical and positional fluctuations, and cold substructures in phase space (Romanowsky et al. 2012; Blom et al. 2014; see also Figure 9 and Section 4.1.2, as well as references in Section 1 and other GC-based results in the literature from Côté et al. 2003; Schuberth et al. 2010; Woodley & Harris 2011; D'Abrusco et al. 2014). Rotation is found to be relatively weak in the outer regions, at odds with general expectations for major mergers, but consonant with multiple, minor mergers in a cosmological context (Vitvitska et al. 2002; Bekki et al. 2008; Romanowsky & Fall 2012; Pota et al. 2013). Note that one rough prediction from cosmological accretion rates is that around half of nearby massive galaxies should have residual substructure in their GC systems from minor mergers (Figure 1).¹³

Despite this general agreement between observation and theory, there is a lingering puzzle about orbital anisotropy in the galaxy halos. Almost any model of spheroidal galaxy formation (monolithic collapse, major mergers, multiple minor mergers) predicts radially biased orbits as a relic of the infall of gas and stars (e.g., Dekel et al. 2005; McMillan et al. 2007; Oñorbe et al. 2007; Sales et al. 2007; Prieto & Gnedin 2008). In the context of two-phase assembly, the radial bias should be prevalent especially in the most massive ellipticals, where the ratio of accreted to in situ stars and GCs is the highest (Wu et al. 2014). However, there are long-standing indications that the orbits of GCs around both the Milky Way and giant ellipticals are near-isotropic or even tangentially biased (Frenk & White 1980; Romanowsky & Kochanek 2001; Côté et al. 2001; Hwang et al. 2008)—a result which has only strengthened with the advent of data from SLUGGS (Romanowsky et al. 2009; Deason et al. 2012; Spitler et al. 2012; Pota et al. 2013; Agnello et al. 2014b; K. A. Woodley et al., in preparation).

This observation has at times been attributed to a natural selection effect wherein GCs on near-radial orbits have been disrupted from passing close to the galactic center. However, this effect has not been shown to operate in the outer halos where low anisotropy is found (cf. Baumgardt 1998; Fall & Zhang 2001; Vesperini et al. 2003; Brockamp et al. 2014). The only other explanation so far proposed is preferential orbit dragging during galaxy formation by slow accretion and adiabatic contraction (Agnello et al. 2014b).

For over two decades it has been realized that most massive galaxies possess GC systems with bimodal color distributions (Zepf & Ashman 1993; Ostrov et al. 1993). Since the vast majority of these clusters have been shown to be old (e.g., Strader et al. 2005), the color bimodality implies a corresponding bimodality in metallicity, as is well-established for the Milky Way, where the metal-rich "bulge" and metal-poor "halo" clusters also have distinctly different kinematics and spatial distributions (Harris 2001). This observation is significant because GC formation is widely assumed to accompany the major star forming episodes in a galaxy's history, and therefore the presence of metallicity subpopulations suggests multiple epochs or phases of galaxy assembly that can be explored via their GC tracer populations (e.g., Section 5.4).

From a kinematic point of view, literature studies (e.g., Côté et al. 2003; Richtler et al. 2004; Schuberth et al. 2010; Lee et al. 2010; Woodley et al. 2010a) and more recently SLUGGS data have been used to establish the kinematic distinctiveness of the subpopulations (Strader et al. 2011; Pota et al. 2013), and to identify cases where more than the standard two subpopulations co-exist, due to a more complex or more recent accretion history (Blom et al. 2012a, 2012b; Agnello et al. 2014b). Bimodality is also evident from the different ellipticities of the red and blue GC systems (Park & Lee 2013).

The GC bimodality paradigm has been disputed in favor of a highly nonlinear relation between metallicity and optical color that transforms a unimodal metallicity distribution into a bimodal color distribution (e.g., Yoon et al. 2006, 2011a, 2011b; Blakeslee et al. 2012). This scenario has been tested with near-infrared colors, which should be more linear, and in the few cases with adequate S/N, bimodality was confirmed (Spitler et al. 2008; Chies-Santos et al. 2012a; Cantiello et al. 2014).

The more definitive approach is to measure metallicities via spectroscopy, which has so far not been done with sufficiently large and homogeneous GC data sets for drawing general conclusions (though see Park et al. 2012). We have now addressed this issue in a series of papers using SLUGGS data (Foster et al. 2009, 2011; Alves-Brito et al. 2011; Brodie et al. 2012; Usher et al. 2012). We found that, in old GC populations, our CaT index is insensitive to alpha-element enhancement or horizontal branch morphology and can be reliably used to infer GC metallicities. Using CaT measurements for over 900 GCs in 11 galaxies, we established directly that the GC metallicity distributions were typically bimodal—with some cases unclear owing to low-number statistics or to the likely presence of more than two metallicity subpopulations.

Here we note that, as with all aspects of GCs, it is important to exclude the highest-luminosity objects from the analysis since they are heavily contaminated by UCDs—which are likely to be chemically and dynamically distinct (Section 5.3). The natural focus on bright objects in earlier work may account for uncertain conclusions about bimodality (e.g., Cohen et al. 1998). We have also found that the color–metallicity relation is a broken linear function that may vary with host galaxy mass (Usher et al. 2014).

As a summary of these results, Figure 16 (top panel) shows "stacked" GC spectra for eight SLUGGS galaxies. The blue stack contains 390 spectra from GCs with 0.72 < (g − i) < 0.88  has a S/N = 373 Å⁻¹, and has a CaT-based metallicity of [Z/H] = −1.2. The red stack contains 353 spectra from GCs with 1.00 < (g − i) < 1.16, has S/N = 419 Å⁻¹, and has a CaT-based metallicity of [Z/H] = −0.4. The two stacks differ not only in their CaT line strengths but also in the strengths of the interspersed, weaker iron lines (marked in the figure), confirming that the CaT is a direct tracer of metallicity. We are now further exploring the dependencies of other weak metal lines (such as Na, Mg, Ti) on metallicity, the IMF, and GC mass (Usher et al. 2014).

Zoom In Zoom Out Reset image size

The bottom panel of Figure 16 shows the spectroscopic metallicity distribution for all 502 GCs with high quality data (Δ[Z/H] <0.3 dex) from a total of 11 galaxies from Usher et al. (2012, 2013). The individual galaxies show metallicity bimodality, which persists in the stacked distribution (no shifts were applied before stacking since the metallicity peaks for each galaxy were all at similar locations). Using the Gaussian mixture modeling code of Muratov & Gnedin (2010), the χ² fit, peak separation, and kurtosis of the distribution all strongly prefer bimodality to unimodality, with p < 0.001, p = 0.158, and p = 0.001, respectively. The code fitted peaks of [Z/H] = −1.47 ± 0.07 and −0.45 ± 0.05 dex. Very similar results are obtained using the full sample of 908 GCs.

In summary, GC bimodality is supported by spectroscopic metallicities and the spatial distributions and kinematics of GC systems. Although nonlinear relations between color and metallicity might in principle play a minor role in shaping GC color distributions, and the GC color/metallicity distributions of galaxies may occasionally not distinctly reveal the underlying GC subpopulations (e.g., due to a more complex accretion history), the evidence strongly supports the result that the widely observed color bimodality is a true reflection of distinct subpopulations with inherently different mean metallicities.

UCDs, compact spheroidal objects too large to fit within the general classification of GCs, were predicted two decades ago to exist as the remnant nuclei of tidally stripped galaxies (Bassino et al. 1994), and were subsequently confirmed by observations in the Virgo and Fornax clusters (Hilker et al. 1999; Drinkwater et al. 2000). Since that time there have been sustained observational efforts to understand whether these objects actually are stripped nuclei (generally from dwarf ellipticals, but possibly also from compact ellipticals) or are extended star clusters—or both (e.g., Da Rocha et al. 2011; Norris & Kannappan 2011). Accumulating a large catalog of UCDs to begin studying them as a population is contingent on ascertaining their distances—usually through spectroscopy. SLUGGS thus provides a key opportunity to advance the UCD inventory by selecting suitable candidates as targets along with normal GCs in our DEIMOS masks.

Some results to date from SLUGGS and associated work include the first UCD found around a non-cluster galaxy (Hau et al. 2009); the discovery of a new class of UCD, fainter than any UCDs previously studied (Strader et al. 2011; Brodie et al. 2011); the extension of objects into the so-called zone of avoidance in the size–luminosity plane for galaxies, UCDs and star clusters (Forbes et al. 2013); as well as the discoveries of the densest known GC (Foster et al. 2011) and the densest known galaxy (Strader et al. 2013), which heralded a link between the UCD and cE populations (Norris et al. 2014; see Figures 17 and 18).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These results have eroded the classical distinction between star clusters and galaxies, since objects are now found to exist with a continuum of structural properties intermediate to these two populations. Deciphering the evolutionary pathways becomes the next challenge after mapping out parameter space. The UCD population appears to consist of both massive star clusters and stripped galaxies, with a ratio that changes across the size–luminosity diagram (e.g., even objects with classical GC-like sizes of ∼4 pc could originate from stripped nuclei; Pfeffer & Baumgardt 2013).

Given the central place of GCs in our efforts to understand the assembly of galaxies in the SLUGGS survey, it is worthwhile to revisit current thinking on the formation of GCs themselves. We take as our starting place the review of Brodie & Strader (2006), while noting that the earlier reviews of Ashman & Zepf (1998) and Harris (2001) still offer useful insights.

Brodie & Strader (2006) discussed the universal appearance of at least two subpopulations of GCs in all but the least massive galaxies, arguing that the metal-poor GCs likely formed at very high redshift (z ∼ 10–15; Elmegreen et al. 2012), while the metal-rich GCs originated either in the dissipational merging that formed the main field star population of the host galaxy (at z ≳ 2 for L* ETGs; e.g., Kruijssen et al. 2012), or/and in clumpy star formation regions in turbulent galaxies (Kravtsov & Gnedin 2005; Shapiro et al. 2010). This scenario was consistent with the ages, metallicities, numbers, and spatial distributions of the GCs in each subpopulation (where the spatial distributions along with the kinematics are a reflection of the GC orbits). These are four fundamental observational constraints that must all be reproduced for any model to be wholly successful in explaining the formation of GC systems (and ultimately of the host galaxies themselves).¹⁴

Subsequent work has clarified and confirmed aspects of this picture while other parts have held up less well. As discussed in Section 5.2, in-depth spectroscopic work has confirmed GC bimodality, which has become a recognized benchmark for modeling the formation of GC systems. As the general field of galaxy formation has shifted over the years from a focus on major mergers to a two-phase assembly paradigm—with early in situ star formation and later satellite accretion—the context and implications of GC modeling have shifted as well (e.g., Beasley et al. 2002; Prieto & Gnedin 2008; Muratov & Gnedin 2010; Tonini 2013; Katz & Ricotti 2014; Li & Gnedin 2014). Where it had been difficult to reproduce distinct GC metallicity peaks from major mergers, the "new" minor-merger dominated paradigm lends itself more naturally to bimodality, with the metal-rich and metal-poor populations forming predominantly in the host and accreted galaxies, respectively.

However, while the models have now succeeded relatively well in reproducing GC metallicity distributions, they have to some extent neglected the other critical constraints of ages and orbital distributions, which have also been firmed up over the years. The compilations of spectroscopic age estimates by Puzia et al. (2005) and Strader et al. (2005) have been reinforced by subsequent findings of generally old GCs in ETGs (e.g., Pierce et al. 2006a, 2006b; Cenarro et al. 2007; Norris et al. 2008; Woodley et al. 2010b; Chies-Santos et al. 2012b; Colucci et al. 2013).¹⁵ However, there are also important novelties from studies of the Milky Way, whose GC system is assumed to be broadly analogous to those of ETGs. Here, while isochrone fitting of GCs suggests ages ranging from ∼11 to 13.5 Gyr (formation epochs of z ∼ 2.5–20; Dotter et al. 2011; Leaman et al. 2013), more precise ages are now claimed based on the white dwarf cooling sequence (Hansen et al. 2007, 2013; Bedin et al. 2009). From studies of one metal-rich and two metal-poor Milky Way GCs, it is inferred that they formed at z ∼ 2 and 3, respectively. These conclusions fit in nicely with the existing picture for the metal-rich GCs—but not for "cosmological" (very high z) formation for the metal-poor GCs.

The small number statistics here are obviously problematic for drawing general conclusions about the overall Milky Way GC system, much less the GCs of ETGs. Moreover, since galaxies began forming at z > 7.5 (Finkelstein et al. 2013), some 13 Gyr ago, it is likely that at least some metal-poor GCs also formed this early. However, if most metal-poor GCs did share a z ∼ 3 formation epoch, then this would present a severe challenge to any model that forms them mainly in association with low-mass dark matter subhalos or dwarf galaxies that are subsequently assembled into larger galaxies. There is a well-established method for connecting the collapse epochs of subhalos with their spatial and orbital distributions at z = 0, with the implication for metal-poor GCs that they must have formed very early, at z ∼ 9–11, based on their tightly bound orbits (Diemand et al. 2005; Moore et al. 2006; Spitler et al. 2012; Corbett Moran et al. 2014). This issue has been noticed in cosmologically motivated models that form GCs at later epochs, where the spatial distribution was found to be much too extended compared to observations (Prieto & Gnedin 2008; a similar problem is implied for the models of Tonini 2013).

Besides the ages issue, there is also the anisotropy puzzle discussed in Section 5.1, wherein the metal-poor GCs are not observed to follow radially biased orbits as expected in accretion scenarios. The overall implication seems to be that a large fraction of these GCs formed in situ within the host galaxy, in a burst of metal-poor star formation that is missing from current models.

This inference seems to mesh nicely with observational evidence in the Milky Way for inner and outer components of the metal-poor stellar halo (Carollo et al. 2007), which have been explained through cosmological hydrodynamical simulations as in situ and accreted components (Zolotov et al. 2009; Font et al. 2011; Tissera et al. 2014). A dual origin for the metal poor GCs is further supported by observed correlations between ages, metallicities, and orbits (Marín-Franch et al. 2009; Forbes & Bridges 2010; Keller et al. 2012), and between metallicity and galactocentric radius, with steep central gradients that transition to constant mean metallicities in the outer parts—both in the Milky Way and in ETGs (Harris 2001; Forbes et al. 2011).

This metallicity gradient behavior has also been observed for metal-rich GCs (Forbes et al. 2011), and taken together, the evidence from GCs suggests three phases of galaxy formation: early in situ metal-poor starbursts, later metal-rich central starbursts, and finally accretion of mixed-metallicity satellites.

The central metal-rich starburst phase may be described more specifically in the context of the current understanding of massive galaxy formation, where the bulk of star formation does not occur in major mergers but in situ, where cold streams from the cosmological web feed violently unstable disks at z ∼ 2 with high star formation rates (e.g., Parry et al. 2009; Dekel et al. 2009; Hopkins & Hernquist 2010). These disks are thought to host super-giant molecular clouds that are natural sites for GC formation (Shapiro et al. 2010), providing a tidy explanation for the origin of many metal-rich GCs.

The appeal of this "wild-disk" metal-rich GC formation scenario is enhanced by contrasting it with the more traditional merger-formation scenario, as in Muratov & Gnedin (2010). Here most of the metal-rich GCs date from z < 1, with very few from z > 2, in marked contrast with the data on GC ages (Dotter et al. 2011; but see also Griffen et al. 2010, who were able to reproduce older GCs through mergers). The merger model does appear successful in at least reproducing the "young" branch of the Milky Way GC age–metallicity relation, but such a relation is probably generic to scenarios of protracted GC formation epochs, e.g., within dwarf galaxies that are later accreted.¹⁶

These developments from the last half-decade represent emerging insights into the mysterious origins of GCs, and of galaxies in general, while the full picture still remains elusive. We expect that substantial further progress will be made through considering the full chemodynamical constraints from the entire SLUGGS sample.