THE SLUGGS SURVEY: WIDE-FIELD STELLAR KINEMATICS OF EARLY-TYPE GALAXIES

The relict signatures of the mass assembly and growth of galaxies are evident in their present-day internal dynamical structure. Within the ΛCDM cosmological framework, galaxies grow through gas infall and in-situ star formation, along with hierarchical merging and accretion of smaller galaxies (e.g., White & Rees 1978). The relative dominance of these modes is conventionally thought to connect to morphological type, with late-type (spiral) galaxies dominated by the in-situ mode, and early-type galaxies (ETGs; ellipticals and lenticulars) by the merging mode.

Merging naturally produces stellar spheroids through violent relaxation, but four decades of detailed observations have shown that ETGs are more than simple featureless spheroids. While remarkably uniform in some respects, e.g., following tight scaling relations (Djorgovski & Davis 1987; Dressler et al. 1987), ETGs also exhibit a diversity of structural features and dynamical behaviors. These include wide variations in central rotation (Bertola & Capaccioli 1975; Illingworth 1977; Davies et al. 1983), isophote shape (Bender et al. 1989), surface brightness profiles (e.g., Lauer et al. 1995; Kormendy et al. 2009), and radio and X-ray luminosity (e.g., Bender et al. 1987, 1989). Because many of these properties are correlated (Bender et al. 1989; Nieto et al. 1991b), the ellipticals can be sub-divided into two fairly distinct varieties: boxy-nonrotating and disky-rotating, where the latter are also closely connected to the lenticulars (Kormendy & Bender 1996).

This bimodal ETG classification scheme has recently been refined through the landmark ATLAS^3D survey (Cappellari et al. 2011a), which used data from the SAURON integral field spectrograph for detailed, homogeneous analyses of kinematics, dynamics, and stellar populations. It was determined that ∼80%–90% of ETGs belonged to the family of centrally fast/regular rotators, characterized by rapid rotation and oblate axisymmetric shapes that reflect underlying disk-like components (Emsellem et al. 2011; Krajnović et al. 2011, 2013; see also Rix & White 1990; Kuntschner et al. 2010). The remaining ∼10%–20% of ETGs were classified as centrally slow/non-regular rotators, which show a variety of complex features such as kinematically distinct cores.

The ATLAS^3D findings motivated a proposed paradigm shift away from traditional morphological classifications and toward a kinematically motivated grouping of S0s with fast-rotator ellipticals (Emsellem et al. 2011; Cappellari et al. 2011b). These galaxies are in turn thought to have close structural and evolutionary connections with spiral galaxies, with remnant disks formed through dissipative processes, while the slow rotators have lost their disks through more active merger histories (Khochfar et al. 2011).

Although the ATLAS^3D results are remarkable, the conclusions about galaxy types are not unequivocal. It is well established from simulations that fast rotators are a very generic category that may include a variety of formational pathways (e.g., Burkert et al. 2008; Naab et al. 2014; Moody et al. 2014), while the finer details of the kinematic data are suggestive of distinct sub-types (Krajnović et al. 2008; Hoffman et al. 2009). Furthermore, the SAURON data were spatially limited. The survey extent in galactocentric radius was to ∼0.9 R_e typically, and to only ∼0.6 R_e for the more luminous galaxies (M_K < −23.5), where R_e is the effective radius enclosing half of the projected light (see Figure 1).

Zoom In Zoom Out Reset image size

This limitation leads to several concerns and questions. First, how certain is it that the defining characteristics of ETGs are found on these scales, and can be used to transform our understanding of galaxies? Second, given the critical importance of kinematics for recognizing the pervasive presence of disk-like structures in ETGs, how well can the properties of the disk and spheroid subcomponents be determined with the available data? Third, how much of the evolutionary story of ETGs is missed without observations of their outer regions, which may undergo different assembly processes than their inner parts?

The last point is particularly relevant in light of the emerging picture of two-phase assembly for ETGs, which is motivated both by observational evidence of dramatic size-growth (e.g., Daddi et al. 2005; Cimatti et al. 2008; van Dokkum et al. 2008) and by cosmologically based simulations. The first phase involves the formation of a compact (∼1 kpc) in-situ bulge at redshifts z > 2 through violent dissipative events such as gas rich major-mergers or a fragmenting turbulent disk fed by cold flows (e.g., Noguchi 1999; Elmegreen et al. 2008; Dekel et al. 2009; Ceverino et al. 2010). The second phase is slower, taking place from z ∼ 3 down to z = 0 (Oser et al. 2012), and consists of mostly gas-poor minor mergers which drive an outward growth of the galaxy (Naab et al. 2009; Bezanson et al. 2009; Oser et al. 2010; van Dokkum et al. 2010; Newman et al. 2012).

In this view, nearby ETGs are layered structures where the central regions are the remnants of the primordial in-situ bulges, and the outer extended halos consist of predominantly accreted material from disrupted satellites. At intermediate radii there should be a transition zone comprised of both accreted material and stars scattered outward from the inner bulge (Zolotov et al. 2009). Tantalizing observational evidence for this two-component structure has been found through photometry of nearby ETGs (Forbes et al. 2011; Huang et al. 2013a, 2013b; Petty et al. 2013), but kinematical and dynamical information is needed for a clearer picture.

Wide-field kinematics observations have been obtained in ETGs using integrated-light stellar kinematics (e.g., Graham et al. 1998; Rix et al. 1999; Kronawitter et al. 2000; Halliday et al. 2001), as well as planetary nebulae (PNe) and globular clusters (GCs) as proxies for the field stars (e.g., Coccato et al. 2009; Arnold et al. 2011; Pota et al. 2013). The latter methods have turned up kinematical transitions at large radii that imply deviations from the central region trends, along with hints that fast and slow rotators follow similar total angular momentum scaling relations (Romanowsky & Fall 2012). However, the results from discrete velocities tend to be somewhat noisy, and there are uncertainties about how these tracers connect to the underlying stellar populations. The gold standard for galaxy kinematics is direct stellar-light spectroscopy, including two-dimensional (2D) information both to map out any off-axis rotation and to provide strong constraints on the internal dynamics (Cappellari 2008). Occasional forays have been made in this direction but have been too restricted either in radius or in 2D coverage (e.g., Statler & Smecker-Hane 1999; Weijmans et al. 2009; Murphy et al. 2011).

A remarkable new approach in this context is to use general-purpose wide-field multi-slit spectrographs to obtain spatially well-sampled 2D spectra out to several R_e. This "Stellar Kinematics from Multiple Slits" (SKiMS) technique was pioneered by Norris et al. (2008) and Proctor et al. (2008), and reached fruition in Proctor et al. (2009) by using the Keck/DEIMOS spectrograph (Faber et al. 2003). DEIMOS SKiMS was applied to a total of eight ETGs through a subsequent series of papers (Arnold et al. 2011; Foster et al. 2011, 2013), and here is expanded to a sample of 22 (including reanalysis of the previously observed systems using new observations and an updated methodology). The large spatial scales of these data, relative to ATLAS^3D, are illustrated by Figure 1, with galaxy parameters reported in Table 1. This is the first systematic study of wide-field stellar kinematics for a large sample of ETGs, which represents a critical step forward in understanding the details of stellar halo assembly.

Table 1. Galaxy Information

Galaxy	R.A.	Decl.	Re	Rmax	$\lambda _{R_e}$	Rot$_{R_e}$	MK	Morph.	P.A.		Vsys	Dist.	Env.
NGC	(h:m:s)	(d:m:s)	(arcsec)	(arcsec)			(mag)		(deg)		(km s−1)	(Mpc)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)
0720	01:53:00.5	−13:44:19	33.9	68	0.12	S	−24.6	E5	142.3	0.43	1745	23.4	F
0821	02:08:21.1	+10:59:42	39.8	78	0.27	F	−24.0	E6	31.2	0.35	1718	23.4	F
1023	02:40:24.0	+39:03:48	47.9	97	0.39	F	−24.0	S0	83.3	0.63	602	11.1	G
1400	03:39:30.8	−18:41:17	29.3	76	0.27	F	−24.3	E1/S0	36.1	0.11	558	26.8	G
1407	03:40:11.8	−18:34:48	63.4	142	0.08	S	−25.4	E0	58.3	0.05	1779	26.8	G
2768	09:11:37.5	+60:02:14	63.1	115	0.25	F	−24.7	E6/S0	91.6	0.57	1353	21.8	G
2974	09:42:33.3	−03:41:57	38.0	57	0.66	F	−23.6	E4	44.2	0.37	1887	20.9	F
3115	10:05:14.0	−07:43:07	34.4	145	0.58	F	−24.0	S0	43.5	0.49	663	9.4	F
3377	10:47:42.3	+13:59:09	35.5	107	0.52	F	−22.8	E5–6	46.3	0.33	690	10.9	G
3608	11:16:58.9	+18:08:55	29.5	54	0.04	S	−23.7	E1–2	82.0	0.20	1226	22.3	G
4111	12:07:03.1	+43:03:57	12.0	45	0.62	F	−23.3	S0	150.3	0.79	792	14.6	G
4278	12:20:06.8	+29:16:51	31.6	89	0.18	F	−23.8	E1–2	39.5	0.09	620	15.6	G
4365	12:24:28.3	+07:19:04	52.5	112	0.09	S	−25.2	E3	40.9	0.24	1243	23.3	G
4374	12:25:03.7	+12:53:13	52.5	136	0.02	S	−25.1	E1	128.8	0.05	1017	18.5	C
4473	12:29:48.9	+13:25:46	26.9	84	0.23	F	−23.8	E5	92.2	0.43	2260	15.3	C
4486	12:30:49.4	+12:23:28	81.3	151	0.02	S	−25.4	E0/cD	151.3	0.16	1284	17.7	C
4494	12:31:24.1	+25:46:31	49.0	119	0.21	F	−24.1	E1–2	176.3	0.14	1342	16.6	G
4526	12:34:03.1	+07:41:58	44.7	99	0.45	F	−24.6	S0	113.7	0.76	617	16.4	C
4649	12:43:40.0	+11:33:10	66.1	178	0.13	F	−25.5	E2/S0	91.3	0.16	1110	17.3	C
4697	12:48:35.9	−05:48:03	61.7	146	0.32	F	−23.9	E6	67.2	0.32	1252	11.4	G
5846	15:06:29.3	+01:36:20	58.9	106	0.03	S	−25.0	E0–1/S0	53.3	0.08	1712	24.2	G
7457	23:00:59.9	+30:08:42	36.3	59	0.52	F	−22.4	S0	124.8	0.47	844	12.9	F

Notes. (1) Galaxy NGC number. (2) Right ascension and (3) declination in J2000 coordinates taken from the NASA/IPAC Extragalactic database (http://ned.ipac.caltech.edu/) (NED). (4) Effective (half-light) radius in units of arcseconds. Effective radii for ATLAS^3D galaxies were taken from Cappellari et al. (2011a), while the remaining values (for NGC 720, NGC 1400, NGC 1407, NGC 3115) were derived using the same methodology. (5) Maximum galactocentric extent of SKiMS data, in arcseconds (see Section 4.2 for definition). (6) Dimensionless proxy for stellar specific angular momentum within 1 R_e (see Emsellem et al. 2007) with the majority of values from Emsellem et al. (2011) and Cappellari et al. (2007, NGC 720). For galaxies without available IFU data (i.e., NGC 1400, NGC 1407 and NGC 3115), we estimated $\lambda _{R_e}$ using SKiMS and literature kinematic measurements. (7) Central rotator designation of fast (F) or slow (S) with values from Emsellem et al. (2011) or derived from available kinematic data (for NGC 720, NGC 1400, NGC 1407, NGC 3115). (8) Extinction corrected absolute K-band magnitudes derived from the 2MASS extended source catalog (Jarrett et al. 2000). Values taken from Emsellem et al. (2011) or measured using the same methodology for galaxies not in the ATLAS^3D sample. (9) Morphology from NED, combining the RSA and RC3 classifications. See Brodie et al. (2014) for the Hubble stage parameters. (10) Photometric position angles and (11) ellipticities. Values taken from Krajnović et al. (2011), Cappellari et al. (2007, NGC 720), Spolaor et al. (2008, NGC 1400, NGC 1407), and Capaccioli et al. (1987, NGC 3115). (12) Systemic velocity in km s⁻¹ from ATLAS^3D (Cappellari et al. 2011a) where available, and otherwise from NED. (13) Distances measured in Mpc from Cappellari et al. (2011a). For NGC 720 and NGC 3115 we have used Tonry et al. (2001) surface brightness fluctuation (SBF) distances with a −0.06 mag correction. For NGC 1400 and NGC 1407 we assume that they are at the same distance and assign an average SBF distance of 26.8 Mpc. (14) Local environment type: field (F), group (G), or cluster (C). For the environmental densities, see Brodie et al. (2014), where we note also that there are some slight differences in the values adopted for R_e, M_K, P.A., , and distance—particularly for the non-ATLAS^3D galaxies.

Download table as:  ASCIITypeset image

These 22 galaxies were observed as part of the ongoing SAGES Legacy Unifying Globulars and GalaxieS Survey (SLUGGS⁸) of 25 ETGs in the local universe (D < 30 Mpc; Brodie et al. 2014). The survey combines stellar and GC system photometry, kinematics, and metallicities to explore the halo regions of ETGs. A key goal of the survey is to target a representative sample of galaxies, spanning a range of luminosities, morphologies, and environments. The diversity in shape and intrinsic brightness of our sample is immediately obvious from Figure 2, which displays color images of each galaxy on a common spatial scale. The environment of our 22 galaxies includes 5 field galaxies (23%), 12 in groups (54%), and 5 in clusters (23%). Further details of the SLUGGS survey are described in Brodie et al. (2014). Analyses of GC metallicities for a subset of the SLUGGS sample have been presented in Usher et al. (2012), of GC kinematics in Pota et al. (2013), and of stellar metallicities in Pastorello et al. (2014).

Zoom In Zoom Out Reset image size

This rich stellar kinematics data set will be analyzed in increasing depth through a series of papers, and here we take a first look at some critical questions motivated above. Do the kinematical trends established with ATLAS^3D extend to larger radii, and support their revamped galaxy classification schemes? Or do novel patterns emerge that require a reassessment of the situation? Do the wide-field data provide better constraints on bulge–disk decomposition? Are kinematical transitions detected that support the two-phase assembly paradigm?

We describe the observational methodology and data reduction process in Section 2. In Section 3, we describe our kinematic measurements of V, σ, h₃, and h₄, and perform repeatability tests and literature comparisons. In Section 4 we introduce a new 2D smoothing algorithm in order to produce wide-field two-dimensional kinematic maps for our sample. We explore some general kinematic trends of the sample in Section 5 and summarize in Section 6.

The SLUGGS observational strategy is designed around obtaining precise radial velocities of individual GCs around nearby ETGs (see Pota et al. 2013 and Brodie et al. 2014 for details). This is accomplished using custom-made multi-slit masks and the DEIMOS spectrograph at the Nasmyth focus of the 10 meter Keck II telescope. DEIMOS is wide-field (167 × 5'), highly multiplexing (∼100 slitlets per slitmask), and equipped with a flexure compensation system, making it well suited for precision kinematic studies. We make optimal use of this instrument by extracting integrated galaxy-stellar-light spectra from the background regions of slits placed on GCs, thereby obtaining both types of spectra simultaneously. This stellar-light aspect of the SLUGGS survey is called Spectroscopic Mapping of EArly-type Galaxies to their Outer Limits (SMEAGOL). The ensuing SKiMS measurements are discussed in Section 4, while the GC measurements are presented in Pota et al. (2013), in references therein, and in additional papers in preparation.

2.1. Observations

DEIMOS has a roughly rectangular field-of-view with a few missing regions (see top of Figure 3). For comparison, the typical R_e value for galaxies in our sample is less than an arcminute. A single slitmask contains ∼100 slits, each with a minimum length of 4 arcsec. Additional filler slits are sometimes placed near the galaxy when there are no suitable GC candidates, and also at large radius to get sky spectra that are virtually uncontaminated by galaxy light. Normally, at least four slitmasks are observed per galaxy, and are arranged in a spoke pattern centered on the galaxy to ensure good azimuthal coverage (see Figure 3 for an example).

Zoom In Zoom Out Reset image size

We collected over a thousand spectra for 22 galaxies over the course of 7 years and 19 observing runs. Our baseline observing strategy is to integrate on a mask for two hours in a series of four exposures, though conditions frequently cause deviations from this plan. The resulting spectra are very deep and permit kinematic information to be extracted from individual slits at surface brightnesses fainter than μ_i = 23 mag arcsec⁻² (see Figure 4).

Zoom In Zoom Out Reset image size

All of the spectra discussed here were obtained using a 1200 line mm⁻¹ grating centered at 7800 Å, and a 5500 Å order-blocking filter. This setup is optimized to target the near-infrared (NIR) Calcium ii triplet, with prominent lines at 8498, 8542, and 8662 Å. In this wavelength range the typical spectral dispersion of 0.32 Å pixel⁻¹ translates to a velocity scale of 11 km s⁻¹ pixel⁻¹. With a slit width of 1 arcsec, the average FWHM resolution is 1.6 Å (see Figure 5), or equivalently ∼55 km s⁻¹. This translates into an instrument velocity resolution of σ_inst ∼ 24 km s⁻¹. Although our full wavelength coverage of ∼6600–9100 Å sometimes includes the Hα line, it is not used in our kinematic analysis.

Zoom In Zoom Out Reset image size

2.2. Data Reduction and Sky Subtraction

The DEIMOS spectra were all reduced using a modified version of the spec2D pipeline (Cooper et al. 2012). The pipeline uses internal quartz-lamp flatfields to rectify each slit, measure a slit function, and make a fringing correction (that is minimal and straightforward to correct using the spectral flat-fields, owing to the accuracy of the flexure compensation system; Wirth et al. 2004). Wavelength solutions are computed from NeArXeNe arclamp spectra taken through each mask, with residuals of 0.04 Å or better (equivalent to 1.4 km s⁻¹) across the whole field. The background spectrum in each slit is modeled as a b-spline and the separate exposures are co-added together before a final one-dimensional spectrum is produced. We modified the pipeline to additionally return the associated pixel-by-pixel inverse variance array of the sky for use during kinematic fitting.

Each background spectrum extracted from a slit is composed of two primary components: the integrated galaxy light that we are ultimately interested in, and foreground sky that must be subtracted out. An accurate measure of the latter component is crucial given the strong and variable emission of the night sky in the NIR. Fortunately, the large field-of-view of DEIMOS relative to the luminous size of the galaxies under study means that there are numerous sky-slits toward the edges of each mask that are essentially free of galaxy light. We select a number of these slits in order to build up a suitable catalog of sky spectra for each mask. Specific care is taken to avoid slits with vignetting or other reduction issues that become more common toward the edges. Also, for reasons that are explained below, we select sky spectra from a variety of positions across the short axis of the mask.

The next step would usually involve constructing a single high-S/N master sky spectrum by co-adding several of the selected sky spectra (see, e.g., Proctor et al. 2009; Foster et al. 2011). However, the application of the conventional approach was complicated by the discovery of a spatially dependent resolution in DEIMOS. Slits located at the bottom of a slitmask have slightly higher resolution (smaller FWHM) than slits placed near the top (see Figure 5). Consequently, the spectrum from an arbitrarily placed slit may have a slightly different resolution than the constructed master sky spectrum.

Ultimately, the effect is subtle, with the FWHM changing by less than a pixel between the top and bottom of the mask. However, we wanted to make every attempt at optimizing the sky subtraction given the potential for large residuals from strong near-infrared sky lines. We opted to use pPXF (Penalized Pixel-Fitting method; Cappellari & Emsellem 2004) to construct a unique sky spectrum for each slit, modeled as a weighted combination of the input sky spectra. This approach is suitable for sky-dominated spectra, which is the case for the vast majority of our data, and should result in a very accurate sky subtraction (see, e.g., Weijmans et al. 2009). An example of a sky-subtracted science spectrum can be see in Figure 6.

Zoom In Zoom Out Reset image size

The line-of-sight velocity distribution (LOSVD) of each spectrum is parameterized as a truncated Gauss-Hermite series (van der Marel & Franx 1993; Gerhard 1993); including V (mean velocity), σ (velocity dispersion), h₃ (a measure of skewness), and h₄ (a measure of kurtosis). These parameters are measured in pixel space using the pPXF software. The best-fitting solution is determined by finding the weighted combination of template stars, which, when convolved by the parameterized LOSVD, most closely matches the input spectrum in a maximum likelihood sense. We include a fourth-order additive legendre polynomial to account for deviations in continuum shape between our spectra and the template stars.

A potential pitfall in making kinematic measurements on low S/N spectra is fitting noise in the spectrum with anomalous h₃ and h₄ values. To help guard against this, pPXF uses regularization in the form of a penalty term that biases the fit toward a Gaussian when the S/N of a spectrum is insufficient to constrain the higher-order moments of the LOSVD. From the Bayesian perspective, this is equivalent to imposing a prior on the h₃ and h₄ parameters that encapsulates the notion that higher-order moments are less likely to be constrained as S/N decreases. A non-zero penalty requires a minimum decrease in the variance of the fit residuals in order to accept a larger amplitude h₃ or h₄ value (see Cappellari & Emsellem 2004). While this penalization results in a bias toward Gaussian LOSVDs, it also reduces the likelihood of measuring spurious higher-order moments.

The appropriate penalty value is at the very least resolution and S/N-dependent, and can be determined via Monte Carlo simulation. Following the method outlined in Cappellari et al. (2011a), we measured the requisite penalty at a range of S/N levels, and fit these results with a third-order polynomial. The coefficients are c₀ = 0.119, c₁ = 5.76 × 10⁻³, c₂ = −1.26 × 10⁻⁵, and c₃ = 9.37 × 10⁻⁹. The penalty value for a given science spectrum is determined using its measured S/N and the above polynomial.

We restrict the pPXF fit to the spectral region immediately surrounding the NIR Calcium ii triplet between 8480 Å and 8750 Å. The latter bound was increased on a galaxy by galaxy basis when a high systemic velocity redshifted the wings of the reddest CaT feature anywhere near the wavelength cutoff (e.g., NGC 4473). No skylines are explicitly masked out since those pixels should, in principle, be appropriately weighted by the inverse variance arrays. However, we do allow for outlying pixel values to be clipped, which can occur when skylines are improperly subtracted.

Template mismatch can be an issue with this method, so an appropriate suite of input template stars is crucial. Previous SKiMS analyses (e.g., Proctor et al. 2009) used 13 template stars observed with DEIMOS that included both giant and main sequence stars at several different metallicities. Here, we use an expanded set of 42 stars carefully selected from the Cenarro et al. (2001) library of 706 stars. The primary differences in the new template catalog are the inclusion of main sequence early-type stars to allow for the presence of Paschen lines, and M-type dwarfs to account for any molecular absorption features. The resolution of the selected template stars is 1.5 Å, which is well matched to the spectral resolution of our spectra (1.5–1.7 Å).

pPXF requires reasonable initial guesses of the V and σ fitting parameters in order to efficiently find the best fitting LOSVDs. This is especially important for lower S/N spectra where noise artifacts (e.g., large skyline residuals) can create likelihood barriers in parameter space. Thus, we devised the following automated method for finding reasonable initial parameter values. First, we create a grid of 120 (V, σ) pairs that span a realistic range of values for the galaxies under study. Second, we run pPXF in two moment mode (i.e., not fitting h₃ or h₄) using each of these pairs as initial guesses. Third, the resulting ensemble of solutions and reduced chi-squared $\tilde{\chi}^{2}$ values are used to determine the best initial guesses of V and σ. Finally, these values and the appropriate penalty terms are used to run pPXF in four moment mode to derive the V, σ, h₃, and h₄ values that maximize the likelihood function.

After determining the best fitting LOSVD for each spectrum, a Monte Carlo method is used to estimate the uncertainty of each kinematic parameter (V, σ, h₃, h₄). First, we use the best-fit model LOSVD to create a realization of each spectrum, and then add noise using the inverse variance array. Each model spectrum is subsequently refitted using pPXF with the penalty set to 0 to derive conservative error estimates as suggested by Cappellari & Emsellem (2004). This process is repeated 100 times per spectrum, with the error on each kinematic parameter taken as the standard deviation of the returned values. The best fitting solutions generally have $\tilde{\chi }^2\le 1$, indicating that the inverse variance arrays are realistic estimators of the pixel by pixel uncertainties of the spectra. In those cases where $\tilde{\chi }^2>1$, we attempt to derive conservative Monte Carlo errors by multiplicatively scaling the inverse variance array such that $\tilde{\chi }^2=1$. Finally, the measured velocities are corrected for heliocentric motion, and each individual fitted spectrum is inspected by eye for quality control. The right panel of Figure 6 shows an example of a measured LOSVD and its 99% uncertainty envelope.

3.1. Measurement Repeatability

Repeated spectroscopic measurements can be useful in assessing the accuracy of Monte Carlo error estimates and in identifying unaccounted-for systematic error terms. We define two spectra to be co-spatial (i.e., repeated) if they are within 3 arcsec of each other. This search radius is small enough to ensure that only spectra from separate masks are paired. Figure 7 shows the difference between repeat measurements, normalized by their Monte Carlo errors added in quadrature, versus the average S/N per resolution element of each spectrum. If the Monte Carlo errors capture all sources of uncertainty in the measurements, then the standard deviation of normalized differences should be nearly unity (i.e., $\tilde{\chi }^2=1$). Indeed, this is case for the vast majority of the measurements, reflecting the superb measurement precision of the DEIMOS spectrograph.

Zoom In Zoom Out Reset image size

However, it is also clear that the distribution widens with increasing S/N, indicating that the uncertainties are underestimated for the brightest spectra. This behavior is expected when there are small systematic error terms that only become non-negligible when the Monte Carlo error approaches zero. After some experimentation, we measured systematic errors in V and σ of 5 km s⁻¹ and 6 km s⁻¹, respectively. These values are added in quadrature to the previously derived Monte Carlo error estimates. The uncertainties on h₃ and h₄ were already conservative, leading to no apparent need for additional systematic error terms.

3.2. Literature Comparisons

While the majority of the DEIMOS spectra are at radii beyond 1 R_e, there are still a fair number of spectra that are within the ATLAS^3D/SAURON observational footprint (Cappellari et al. 2011a). The galaxies in common between the two surveys span a range of systemic velocities, velocity dispersions, h₃, and h₄ values, providing a valuable opportunity to cross check measured kinematic quantities. Figure 8 shows a comparison between measurements from each survey that are spatially coincident to within 3 arcsec. The dashed black lines in each panel are linear fits to the measurements that take x- and y-axis error bars into account while the solid black lines show the 1:1 relation.

Zoom In Zoom Out Reset image size

A few things are immediately obvious from this figure. For instance, the general agreement between the two surveys is very good, with the measurements showing a high degree of correlation, and the median offsets for each parameter being very small. The exception is the velocity dispersion, where the SAURON values are systematically higher than the DEIMOS values by 13 km s⁻¹. We note that Murphy et al. (2011) found higher velocity dispersions than SAURON for NGC 4486 (M87), and that the mismatches were markedly reduced for a few galaxies after the ATLAS^3D improved reanalysis of the original SAURON data (Emsellem et al. 2004).

Before discussing some specific galaxies, we should point out some important caveats regarding this comparison. First of all, most of the matched spectra are near the outer edge of the SAURON field-of-view, and therefore represent some of the noisier measurements from that survey. Co-adding such data may lead to higher on average velocity dispersions being measured. Conversely, these data comprise the high S/N tail of our DEIMOS observations. Second, the spectral resolution and dispersions of each instrument are quite different. Operating in low resolution mode, SAURON has an instrumental velocity resolution of σ_inst = 108 km s⁻¹ and a dispersion of 66 km s⁻¹ pixel⁻¹. The corresponding values for our DEIMOS measurements are 24 km s⁻¹ and 11 km s⁻¹ pixel⁻¹. Third, and perhaps most importantly, these two surveys utilize very different spectral regions, with SAURON probing shorter wavelengths (4800–5380 Å) than our survey using the Calcium Triplet lines (∼8500 Å). These spectral windows are quite different, each having benefits and drawbacks over the other. LOSVD fitting at bluer wavelengths benefits from more absorption features and fewer skylines than in the near-infrared. However, this range is also more complicated, can have emission lines to contend with, and is sensitive to template mismatch issues that do not affect the near-infrared CaT (Barth et al. 2002). Given that the SAURON and DEIMOS spectral regions probe slightly different stellar populations, this could in principle be the cause of the systematic offset seen in velocity dispersion. This effect is explored in Appendix C, but we conclude that it is very unlikely to be responsible for the offsets we see.

In some cases, the discrepancy could be due to sampling issues in the SAURON spectra. For example, in NGC 7457 the matched SAURON velocity dispersion values range between 30 and 70 km s⁻¹ (i.e., at or below their spectral resolution) while the corresponding DEIMOS values cluster tightly around 40 km s⁻¹. However, this cannot be the case for all galaxies.

A possible contribution to the systematic offset may be the different values of the h₄ parameter in each survey. The h₄ parameter measures the kurtosis or how peaked the velocity profile is. The velocity profiles in the radial range of overlap of our data often have h₄ parameters close to zero, whereas those from the SAURON instrument are slightly positive. If the velocity profile is more peaked than a Gaussian, then a fit with h₄ = 0 will tend to underestimate the true width of the profile and hence derive a smaller velocity dispersion.

The primary focus of this paper is to present the large-scale kinematic structure of ETGs, and to highlight general trends with radius. To do this, we will be combining SAURON data and other literature data sets with our own in order to construct 2D kinematic maps (see the next section). Since we are concerned with general trends, we remove any mean systematic offset between the data sets so that the final maps do not have unphysical step functions in them.

In this section, we outline the methodology used to create 2D maps of V, σ, h₃, and h₄, for each galaxy in our sample. These 2D kinematic maps are meant to provide an immediate sense of large-scale kinematic behavior and to highlight structures and general trends in the data. While these types of maps are most commonly used to present the data products of IFU observations, they are also well suited for presenting other spatially resolved data sets, such as ours. However, the spatial sampling density resulting from slitmask spectra is quite sparse compared to that of an IFU. This results in large gaps in spatial coverage that require smoothing and/or interpolation to fill.

Previous approaches toward solving this problem include using a kinemetric model (Proctor et al. 2009), and kriging fitting (Arnold et al. 2011; Foster et al. 2013). Both of these methods have their pros and cons, e.g., the former requires a specific parameterization, while the latter is prone to over-smoothing or ringing. Here we introduce a new methodology that requires no parameterization, preserves much of the spatial information, and minimally smooths the data to suppress noisy measurements.

4.1. Smoothing Algorithm

We account for the empty regions between spectra by generating a Voronoi tessellation that assigns a spatial Voronoi bin to each stellar spectrum. For example, the locations of the stellar spectra for NGC 821 are plotted in Figure 9(b), with the corresponding Voronoi bins visible as single colored polygons in Figure 9(c). We point out that the raw velocity map displayed in Figure 9(c) is noisy enough to obscure the general kinematic trends in the data, and is not very useful as is. In principle, we could have partially ameliorated this problem by co-adding several noisy spectra together into a single higher S/N spectrum before fitting the LOSVD. However, while this approach is suitable for IFU measurements with ample spatial sampling, it would necessarily destroy valuable spatial information in our data. Additionally, many of our large-radius spectra are separated by rather large distances on the sky, and may be unsuitable for direct co-addition since they sample kinematically distinct regions of the galaxy.

Zoom In Zoom Out Reset image size

A better approach is to smooth the noisy data just enough to bring out the underlying trends, but not so much as to destroy spatial information. A common technique is to use a Gaussian smoothing kernel that de-weights more distant, and thus less correlated, points. However, the varying spatial density of measurements makes a single-width Gaussian kernel unsuitable, since it would lead to maps with zones of under or over smoothing. The adaptive Gaussian kernel method used in Coccato et al. (2009) attempts to solve this problem, but still relies on an arbitrarily shaped Gaussian weighting function and an assumed kinematic model for the galaxy.

4.1.1. Weighting Function

Here, we use an empirical approach that uses the spatial correlation information contained in each data set to optimally smooth the associated kinematic maps. For each separate data set (e.g., h₃ values of NGC 821 or σ values of NGC 1023), we calculate the semi-variance γ, as defined in Equation (1) below. This is measured by finding all pairs of values that are separated by a particular distance Δ to within the specified tolerance δ, and then measuring half the variance of their differences.

$$\begin{eqnarray} \gamma (\Delta) &=& \frac{1}{2\,N_S} \sum _{(i,j)\in S}^{N_S} [Z(\boldsymbol {x}_i) - Z(\boldsymbol {x}_j)]^2 \nonumber\\ S &=& \lbrace (i,j) : -\delta < \Vert \boldsymbol {x}_i - \boldsymbol {x}_j \Vert - \Delta \le \delta ; i > j \rbrace . \end{eqnarray} \tag{ 1 }$$

Here, $\boldsymbol {x}_k$ is the position of the k-th spectrum and $Z(\boldsymbol {x}_k)$ is the value of interest at that location.

For each data set we use the inverted form of the semi-variance as an empirical weighting function in the smoothing process outlined below.

4.1.2. Grouping Measurements

With the weighting function in hand, the next step is to combine the measurements together to suppress noise while not washing out legitimate large-scale structures. Our adopted technique is an extension of the one-dimensional moving-average algorithm, commonly used to filter time-series data, to two dimensions. The basic idea is to estimate the local mean of a data set for a series of spatial locations by taking the appropriately weighted average of nearby measurements. A key point here is that individual measurements can belong to many unique groupings of data, each centered at a different spatial location.

The groupings (or clusters) of measurements are selected using a custom 2D search algorithm and the requirement that clusters meet certain connectedness criteria. Each group is assigned a unique identifier, based on the associated measurements. In each cluster, the semi-variance is calculated for each grouped measurement using (1) its distance to the centroid of the cluster and (2) the polynomial fit for the associated semi-variogram. These values are then combined with the Monte Carlo error estimates to derive a final weighting factor for each measurement. Next, these weights are used to compute a weighted sum of measurements for each cluster. At this stage, there are many unique clusters that overlap heavily with other groupings. The final step involves combining the overlapping clusters together, associating the results with the Voronoi bins associated with each individual measurement, and creating the final map.

However, we have not yet specified how many points should go into a cluster—a crucial issue since this will control the degree of smoothness in the final map. Importantly, the number of points is not constant, and is determined separately for each grouping of measurements. We require that each cluster contains just enough points such that the resultant error of the weighted sum is within some tolerance (e.g., ∼10–15 km s⁻¹ for velocity maps). This feature makes the smoothing method naturally adaptive because groupings in high spatial-measurement-density regions will be physically small and localized, and just the opposite in low-density regions. In basic terms, the smoothing is optimized so that the centers of maps are not over-smoothed, and the edges of maps are not under-smoothed.

Figures 9(c) and (e) show maps for the raw and smoothed versions of the NGC 821 velocity data. The latter map reveals far more structure than the former, including a falling rotational amplitude along the major-axis. For comparison, Figure 9(d) shows the map that results from a flat weighting kernel, which is similar to taking a straight average in each cluster.

4.2. Final Maps

With the kinematics data and maps in hand, we carry out several different analyses to characterize the 22 galaxies in our sample, and to connect to the science questions raised in Section 1. We examine variations of angular momentum with radius in Section 5.1, bulge–disk structure in Section 5.2, and implications for galaxy formation in Section 5.3.

5.1. Angular Momentum

The principal, novel parameter used in the ATLAS^3D survey is the stellar specific angular momentum proxy, λ_R,

$$\begin{equation} \lambda _R \equiv \frac{\langle R |V|\rangle }{\langle R\sqrt{V^2+\sigma ^2}\rangle }, \end{equation} \tag{ 2 }$$

calculated cumulatively inside either 0.5 or 1.0 R_e. This parameter is closely related to traditional V/σ measurements, and has limiting values of λ_R = 0 to λ_R = 1 for minimum and maximum rotation support, respectively. In Emsellem et al. (2007), a fixed value of λ_R = 0.1 was proposed as the dividing line between slow and fast rotators. This divide was refined in Emsellem et al. (2011) by including a link with ellipticity:

$$\begin{equation} \lambda _R = \lambda _0 \sqrt{\epsilon } \,, \end{equation} \tag{ 3 }$$

where λ₀ ≃ 0.3 and the exact value depends on whether 0.5 or 1.0 R_e apertures are used.

The 0.5–1.0 R_e apertures were motivated not only by the instrumentation constraints but also by the finding that these observations reached beyond most of the kinematically distinct components (KDCs) and into the region where λ_R apparently plateaued. However, the asymptotic values of λ_R including larger radius data were not explored, and it should be noted that while 1 R_e is a natural scale for the stellar luminosities and stellar masses of galaxies, the natural scale for total masses and angular momenta is much farther out (see the discussion in Romanowsky & Fall 2012).

Our goal is therefore to investigate further the scale dependence of stellar specific angular momentum trends, and the robustness of the ATLAS^3D conclusions to the inclusion of larger-radius data. Here we will take a first look at these issues, and defer a full examination to another paper (C. Foster et al., in preparation). Also, in the interest of focusing on variations with radius, here we will work with differential rather than cumulative values (which require careful luminosity weighting). This means that we will calculate a local λ as in Equation (2) but in annuli rather than over filled apertures.

Our λ calculations proceed as follows. Each smoothed V and σ map is subdivided into a fine 2D grid of pixels. These pixels are binned into annuli with widths of 5 arcsec, and with fixed ellipticities from Table 1. The value of λ within each annulus is calculated from Equation (2), using the true galactocentric radius R rather than the elliptical radius R_m, and assuming constant flux weighting per pixel. Profiles of λ with R_m are generated using an oversampled moving window approach, where each annulus is advanced outward 1 arcsec relative to the previous one. For stability, these profiles are computed only out to a radius R_max (see Section 4.2 for definition).

In Figure 10, we show the resulting radial dependence of the λ parameter for our sample. Galaxies are color-coded by their central fast/slow rotator designation. A diversity of λ profiles is seen when the kinematics are probed beyond 1 R_e. Similar results for a sample of galaxies were found by Coccato et al. (2009) from PN kinematics, and here we directly confirm that this complexity exists in the stellar kinematics (as also seen by Proctor et al. 2009 in a much smaller galaxy sample; variations in central kinematic profiles were also seen by Emsellem et al. 2011 but are difficult to compare since these were cumulative rather than differential λ profiles).

Zoom In Zoom Out Reset image size

The seven centrally slow rotators continue to have low angular momentum to large radii, indicating relatively slow rotation in the bulk of the stellar mass for these galaxies (note that although λ_R increases up to 0.2 for some of them, there is an accompanying increase of ellipticity, which suggests a continued slow-rotator classification as in Equation (3)). Slow rotating galaxies tend to be the most massive galaxies and located at the high density cores of galaxy groups and clusters. They also reveal dynamically interesting features such as kinematically distinct cores (KDCs; e.g., NGC 3608, NGC 4365), minor-axis rotation (e.g., NGC 4365, NGC 4374, NGC 4486), and kinematic twisting (e.g., NGC 1407).

The centrally fast rotators, on the other hand, are more varied. Half of the galaxies have λ profiles that continue to increase with radius (NGC 720, NGC 1023, NGC 2768, NGC 2974, NGC 4111, NGC 4526, NGC 4649, NGC 7457), while half have λ nearly constant or decreasing (NGC 821, NGC 1400, NGC 3115, NGC 3377, NGC 4278, NGC 4473, NGC 4494, NGC 4697). Some decreasing cases (NGC 821, NGC 4278, NGC 4473) even reach the domain of formally slow rotation.

We next consider a more simplified view of the λ variations by introducing a λ-gradient, Δλ, which is the change in local λ between 0.5 and 2.0 R_e (or the farthest extent of the data, in a few cases where this is less than 2.0 R_e). In Figure 11, we plot Δλ versus the local λ at 1.0 R_e for the 22 galaxies in our sample. Galaxies with Δλ > 0 have rotation becoming more dominant with radius, and less dominant for Δλ < 0. The point styles reflect galaxy luminosities, ellipticities, physical sizes, and morphological types, all of which could, in principle, relate to the gradients in specific angular momentum.

Zoom In Zoom Out Reset image size

We do not discern any clear connecting patterns between Δλ and the galaxy parameters including environment – except for possibly the galaxy types. The seven fast rotators with the most strongly increasing angular momentum are all S0 or borderline E/S0 (if one considers that NGC 2974 was classified as an S0 or later by all six observers in Naim et al. 1995, and even exhibits spiral arms). The remaining nine fast rotators, with more constant or declining λ, include seven ellipticals, one E/S0, and one S0.

This intriguing segregation suggests that the standard morphological classifications may still be relevant: although ATLAS^3D has established that the Es and S0s have similar kinematical properties inside 1 R_e, at larger radii their kinematics may diverge systematically. Whether these galaxies are fundamentally the same or different will require more detailed structural–dynamical analyses (including the links between radial variations in rotation and ellipticity), and ideally a larger sample of galaxies. However, we will look at possible interpretations in the next subsections.

One notable case here is NGC 4649, which has the second highest luminosity in the entire ATLAS^3D sample of 260 ETGs. It is somewhat peculiar in its central regions (at 0.35 R_e), where, unlike all of the other most massive ETGs, it is (barely) a fast rotator (λ_R = 0.13). Its behavior is even more deviant in its outer regions, where it becomes remarkably fast-rotating (λ ∼ 0.5 at 2 R_e). We will briefly consider the nature and origin of this galaxy in the following subsections.

In summary, we can now provide some insight into the radial dependence of specific angular momentum beyond one effective radius. We find that the central fast and slow rotator classifications remain valid at larger radii in many, but not all, cases. This suggests that centrally fast rotators are not necessarily a homogeneous class of galaxies, and that the elliptical versus lenticular galaxy morphological distinction may still be relevant.

Following the initial submission of this paper, there appeared a related survey of extended kinematics in ETGs using the Mitchell Spectrograph (formerly VIRUS-P) integral field unit (Raskutti et al. 2014). Although detailed comparisons with our results are difficult owing to the low spatial resolution of the VIRUS-P survey at both small and large radii, the qualitative conclusions appear the same: Raskutti et al. found that the slow and fast rotator classifications generally remain valid out to the limit of their data, but with a few galaxies showing drops in rotation and extended kinematically distinct components.

We note that even our stellar kinematics data are still restricted in radius to ∼2.5 R_e, and we cannot test how the results may change much farther out (where we recall that ∼5 R_e is the characteristic scale length for total stellar angular momentum). Indeed, there are already two cases known from GC and PN kinematics where fast rotation at ∼2.5 R_e dives toward slow rotation at larger radii (NGC 1023 and NGC 3115; Noordermeer et al. 2008; Arnold et al. 2011).

5.2. Embedded Disks

One key finding from the previous subsection is the dramatic decrease of local angular momentum with radius in a significant fraction of the galaxies. Coccato et al. (2009) suggested a natural interpretation, that the central regions host a significant stellar disk which fades toward larger radii. Hence the central kinematics may reflect a combination of disk and bulge components, while the outer kinematics are "pure bulge."

This interpretation leads us to one of the recurring themes in extragalactic astronomy, of decomposing galaxies into their constituent parts, primarily bulges and disks. This means determining the relative fractions of these subcomponents in luminosity, as well as their subcomponent density profiles and flattening. There have been many techniques and software packages developed to carry out bulge–disk decomposition based on photometry, but the results are notoriously model-dependent and uncertain.

In this context, kinematics offers a ray of hope, as it may be used in combination with photometry to alleviate the degeneracies in deprojection and decomposition (Cinzano & van der Marel 1994; Scorza & Bender 1995; Romanowsky & Kochanek 1997; Magorrian 1999; Cappellari et al. 2013a). The potential here is evident by the clear disk-like features found in the ATLAS^3D/SAURON stellar kinematic maps, when the surface brightness maps appeared much more ambiguous about the underlying galaxy structure. In practice, however, photometric–kinematic decomposition has involved fairly arduous modeling (e.g., based on detailed LOSVD shapes), and the ATLAS^3D survey has so far applied only a conventional photometric approach (Krajnović et al. 2013).

We suggest instead a new approach where the availability of wide-field stellar velocity maps may permit relatively painless and well-constrained basic bulge–disk decompositions. The starting assumption is that the disk and bulge components each have fairly uniform rotational properties, and that the overall variation in rotation with radius then reflects the transition between disk and bulge dominance. The combined photometric and rotation profiles may then be fitted with a disk–bulge model, as has been shown in recent pioneering work using PN and GC data to decompose S0 galaxies (Cortesi et al. 2011, 2013; Forbes et al. 2012). The same idea could also be used for fast-rotator ellipticals, where kinematic transitions are found at smaller radii than in the S0s, and our SKiMS data could be sufficient to constrain the decompositions.

As an illustration, we consider the case of NGC 3377, a low-luminosity elliptical with a strong decline in λ beyond ∼1 R_e. The decline in rotation amplitude alone is even stronger, as shown in the left panel of Figure 12, which also combines PN kinematics data in order to verify the behavior at larger radii (Coccato et al. 2009). At a glance, we may suppose that the disk-to-bulge transition occurs at ∼1.5 R_e.

Zoom In Zoom Out Reset image size

The photometric data for NGC 3377 are shown in the right panel, as a 1D surface brightness profile (Goudfrooij et al. 1994). We have carried out a three-component fit, where, in addition to standard disk and bulge components (Sérsic 1968 indices n = 1 and n = 4), we have included a cusp of "extra light" based on the decomposition of Hopkins et al. (2009). This third central component is generic to many ellipticals (where it may trace ancient starbursts) and here impacts the overall decomposition results for the rest of the galaxy. Our fitted parameters are R_e = 35 for the cusp, a disk scale radius of R_d = 139, and R_e = 517 for the bulge. The fractional luminosities in the different components are 19%, 20%, and 61%, respectively. These contrast with previous disk-fraction inferences of 5% and 56%–82% from Scorza & Bender (1995) and Krajnović et al. (2013), respectively.⁹

After fitting the photometry, we assign simple constant v_rot values to each component, constructing a total v_rot profile that we compare to the data. We find that values of 160 km s⁻¹, 170 km s⁻¹, and 25 km s⁻¹ (cusp, disk, bulge) fit the rotation profile remarkably well (note that the maximum circular velocity was found from dynamical modeling to be 212 km s⁻¹; Cappellari et al. 2013b). We also found it difficult to reproduce this profile without a decent photometric fit that included all three components, and we further note that the ATLAS^3D decomposition, with R_e = 7'' for the bulge and R_d = 22'' for the disk, is not consistent with these kinematic data.

This simplistic 1D model fit is meant to be schematic only, and not to fully explore the parameter uncertainties, nor to make use of all of the available data. In particular, a robust fit should be carried out in two dimensions for both the photometry and the kinematics, and would also incorporate stellar population differences (see Johnston et al. 2012; Coccato et al. 2013; Fabricius et al. 2014). Figure 13 illustrates some of the available 2D photometric information, where the galaxy appears flattened and disky in the inner regions, and round and boxy in the outer regions. This morphology supports the notion of a distinct disk component embedded in a more extended bulge, and the isophote shape transition radius seems consistent with the kinematic transition (for more discussion of how ellipticity correlates with rotational support at fixed anisotropy, see Emsellem et al. 2011, Appendix B). A similar morphology has also been found for another case from our sample with a declining rotation profile (i.e., NGC 4697; Scorza & Bender 1995; Huang et al. 2013b).

Zoom In Zoom Out Reset image size

As another consistency check on the embedded-disk scenario, we consider the velocity dispersion maps (Figure 14). We expect that along the major axis, a cold disk will produce a lower overall velocity dispersion, relative to the minor axis, which should be dominated by the hotter bulge component. Unambiguous instances of this behavior are indeed observed in the maps for three of the declining rotators (NGC 821, NGC 3115, NGC 3377).

Our final set of clues comes from the Gauss-Hermite h₃ moment. An anti-correlation between h₃ and V/σ implies an extended tail of stars in the direction of the galactic systemic velocity, which in turn suggests a fast rotating disk superimposed on a more slowly rotating bulge (Rix & White 1990; Scorza & Bender 1995). Such an anti-correlation has been frequently observed in the central regions of ETGs (e.g., Bender et al. 1994; Krajnović et al. 2008), and now we have the opportunity to see what happens in their outer regions.

In our sample, all of the fast rotator galaxies with embedded disk-like structures identified by their rapid rotation in the previous section have anti-correlated V/σ − h₃ values in the same region. The same is true for the massive galaxy NGC 4649 as well as for the slow rotator NGC 1407, which has a rotating structure within 1 R_e that is misaligned with the photometric major-axis and clearly exhibits anti-correlated V/σ − h₃. These observations support the embedded disk picture.

Intriguingly, the anti-correlation extends beyond the regions of obvious fast rotating or disk-like embedded components. This includes not only along the major axis at larger radii, but also off-axis and well away from the disk. The implication, along with the observation of weak kinematic twists, may be that the central disk is embedded in a much puffier, hot-disk component (see Kormendy & Bender 2012)—which seems consistent with the ATLAS^3D finding that these galaxies remain oblate–axisymmetric out to ∼2.5–3 R_e, based on kinematic alignments (Krajnović et al. 2011). These components could in turn be embedded in a more extended slow-rotating spheroid that would produce a second rotation drop, given observations at even larger radii (see, e.g., Arnold et al. 2011). We thus speculate that ETGs may be kinematically multi-component systems, with multiple layers of embedded structure.

In summary, our results reinforce the long-standing picture that many ellipticals as well as lenticulars host embedded stellar disks. These two morphological classes may simply differ in the relative masses and length scales of the disks, and represent a continuous family that may even connect smoothly to spirals (e.g., Scorza & Bender 1995; Cappellari et al. 2011b). Whether or not the galaxies appear kinematically different would then depend on the spectroscopic field of view relative to the disk scale length, which is nearly always larger than the SAURON aperture, but in many cases smaller than our DEIMOS aperture.

If this simple picture is correct, it could validate the claims from ATLAS^3D that Es and S0s comprise a unified family of fast rotators. However, just because two galaxies have disks and bulges does not necessarily mean that they are closely related: it could be a case of parallel evolution, as in the wings of birds and bats whose homology is driven by common physical principles. Resolution of this question will require detailed examination of the properties of the galaxy subcomponents and of their scaling relations (see, e.g., Figure 34(b) in Romanowsky & Fall 2012).

5.3. Two-phase Assembly

We have presented wide-field (to ∼2.5 R_e) kinematic maps of V, σ, h₃, and h₄ for 22 nearby early-type galaxies. Our sample consists of 16 centrally fast rotators and 6 centrally slow rotators—close to the observed fractions measured in the volume limited sample from ATLAS^3D (Emsellem et al. 2011). These galaxies span a range of sizes (R_e = 0.8–8.6 kpc), ellipticities (0.04–0.58), luminosities (−22.4 > M_K > −25.5), morphologies (S0–E0), and environments (field, group, and cluster). Given this diversity of properties, it is perhaps unsurprising that these galaxies also exhibit a range of large-scale kinematic behaviors (see Figure 14).

The wealth of structure evident in these large-scale kinematic maps will require detailed dynamical modeling to fully understand and interpret within the context of ETG assembly. Here, our goal was to present the data and our analysis method. We also highlight some general trends. For example, we find that the centrally slow rotator ETGs, which reside almost exclusively in the high density cores of groups and clusters, remain slowly rotating in their halos. On the other hand, the centrally fast rotating galaxies exhibit different behaviors at large radius with about half of our sample continuing to rotate rapidly at large radii while the other half have declining angular momentum profiles so that they rotate slowly in their halos. In other words, observations of central kinematics miss a significant part of the overall structural and dynamical picture of ETGs.

The variations of rotation with radius appear to correlate with the traditional galaxy morphological types (E and S0), and we suggest that it would be premature to retire such classifications. The rotational gradients are probably connected to the relative dominance of embedded stellar disks, and we provide a conceptual demonstration of how wide-field kinematics data might be used to dramatically improve bulge–disk decompositions.

The dramatic cases of rapidly declining rotation profiles suggest galaxies with distinct or decoupled inner and outer regions. Such transitions may be fossil signatures of two-phase galaxy formation, although they have yet to be produced in detail by cosmological simulations.

Future papers in this series will further characterize and explore the interesting structures shown here, while the data itself will become public at the conclusion of the survey.

We thank A. Burkert for useful discussions, E. Ramirez-Ruiz for computational support, R. van den Bosch, M. Cappellari, E. Emsellem, D. Krajnović, and T. Lauer for helpful comments, and the referee for a positive report. The data presented herein were obtained at the W. M. Keck Observatory, operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration, and made possible by the generous financial support of the W. M. Keck Foundation. Based in part on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. We thank the NSF for a Graduate Student Research Fellowship and Sigma-Xi for their financial support. We thank the National Science Foundation for funding via AST-0909237 and AST-1211995. We thank the ARC for funding via DP130100388.

In Figure 14, the data for each galaxy are laid out as follows. At left, the galaxy's name is listed above the isophotal contour map of the respective I-band DSS image—oriented so that the photometric major-axis is horizontal.¹⁰ The associated compass, which points north by east, illustrates the rotation of the data (see Table 1 for adopted position angles). The panel below the contour map shows plots of the locations of each measurement used to construct the kinematic maps and includes the 1 R_e contour and the convex hull (bounding polygon) of data locations. This hull is slightly larger than R_max, the maximum radial extent of the data, which we define as the largest intermediate-axis radius at which at least 75% of the area within an elliptical annulus lies inside the outer measurement hull.

In the smoothed kinematic maps shown at right, each measurement is assigned a Voronoi bin (see Section 4.1.2) and plotted as a function of X and Y position in the galaxy, normalized by R_e. Panels are labeled according to measurement type (V, σ, h₃, or h₄), and Voronoi bins are colored according to value using the color bar shown at the bottom left. The minimum/maximum of the range is indicated in the bottom right corner of each panel.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this Appendix, we list a (non-exhaustive) selection of salient properties for our sample galaxies along with a qualitative description of their respective combined stellar and GC (when available) kinematics. In all cases, morphological classifications are from the RC3 as reported in NED.

The following literature kinematic data (IFU and long slit) are included along with our Keck/DEIMOS measurements in order to provide coverage of the galaxy centers and improve spatial sampling: NGC 720 (Cappellari et al. 2007); NGC 821 (Weijmans et al. 2009; Pinkney et al. 2003; Proctor et al. 2005; Forestell & Gebhardt 2010; Cappellari et al. 2011a); NGC 1023 (Cappellari et al. 2011a; Fabricius et al. 2012); NGC 1400 (Spolaor et al. 2008); NGC 1407 (Spolaor et al. 2008); NGC 2768 (Cappellari et al. 2011a); NGC 2974 (Cappellari et al. 2011a); NGC 3115 (Norris et al. 2006); NGC 3377 (Coccato et al. 2009; Cappellari et al. 2011a); NGC 3608 (Cappellari et al. 2011a); NGC 4111 (Cappellari et al. 2011a); NGC 4278 (Cappellari et al. 2011a); NGC 4365 (Cappellari et al. 2011a); NGC 4374 (Cappellari et al. 2011a); NGC 4473 (Cappellari et al. 2011a); NGC 4486 (Cappellari et al. 2011a; Murphy et al. 2011); NGC 4494 (Coccato et al. 2009; Cappellari et al. 2011a); NGC 4526 (Cappellari et al. 2011a); NGC 4649 (Pinkney et al. 2003; Cappellari et al. 2011a); NGC 4697 (Cappellari et al. 2011a); NGC 5846 (Cappellari et al. 2011a); NGC 7457 (Cappellari et al. 2011a).

NGC 720. This is a highly elongated E5 galaxy possibly residing in a triaxial halo (Buote et al. 2002). Cappellari et al. (2007) found clear but moderate central rotation within 1 R_e, and similarities to both fast and slow rotators. Our kinematic data show that this slow but significant rotation continues essentially constant out to the edge of the SLUGGS data around 2 R_e.

NGC 821. This E6 galaxy is a classic example of a disky elliptical (Lauer 1985), although Kormendy & Ho (2013) argued for a reclassification as an S0. Its characteristic disky isophotes disappear beyond 1 R_e (about 30 arcsec; Bender et al. 1988; Nieto et al. 1991a; Goudfrooij et al. 1994). Rotation is prominent inside of this radius (Dressler & Sandage 1983; Bender et al. 1994; Pinkney et al. 2003; Emsellem et al. 2004), but diminishes farther out (see Figure 14; Proctor et al. 2009; Coccato et al. 2009; Weijmans 2010). Our data show that both the velocity dispersion and the h₄ maps exhibit large-scale structure with distinctly different amplitudes along major and minor axes. Only the metal-poor GCs show significant rotation and that is along the photometric minor axis in the inner parts (Pota et al. 2013), consistently with the PN data (Coccato et al. 2009).

NGC 1023. This is a barred lenticular galaxy (SB0). Our measured velocity dispersion profile declines strongly with radius, with σ values eventually falling below the resolution limit. Over the same radial range, the rotational amplitude gradually increases with radius before plateauing near 1 R_e (to ∼200 km s⁻¹) and eventually dropping beyond 3 R_e (see also, e.g., Noordermeer et al. 2008; Coccato et al. 2009; Cortesi et al. 2011).

NGC 1400. This galaxy (E1/S0) is the second brightest galaxy of the NGC 1407 group. The combined velocity map indicates a moderately sized embedded disk on the scales of the effective radius within an aligned and moderately rotating bulge. Additional evidence for this interpretation comes from the clear anti-correlation between V/σ and h₃ that is present only in this region. Its GC system kinematics are consistent with the stellar kinematics (Pota et al. 2013).

NGC 1407. This is a giant elliptical galaxy (E0) at the center of the Eridanus A group. Weak rotation at all radii with anti-correlated h₃ within a half-light radius are evident. Overall, the velocity dispersion falls with radius, though a notable upward "annular" spike is apparent around 1 R_e (first detected by Proctor et al. 2009). There is a hint of a spatially corresponding structure in h₄. Pota et al. (2013) found that the metal-rich GCs rotate along the major axis while metal-poor GCs only show significant rotation in the inner parts along the minor axis. This rotation declines and progressively aligns with the major axis in the outskirts. The metal-rich GC velocity dispersion is consistent with that of the galaxy stars. Both GC subpopulations have tangentially biased orbits (Romanowsky et al. 2009).

NGC 2768. This is a well studied flattened (E6/S0) galaxy. Its cylindrical rotational field within ∼40 arcsec (Emsellem et al. 2004) and sustained overall rapid rotation are obvious in our combined maps and have already been well documented (Fried & Illingworth 1994; Proctor et al. 2009). V and h₃ are globally anti-correlated as expected from an S0 with bulge and disk components of similar scale size and relative brightness (B/T = 0.7; Forbes et al. 2012). According to Pota et al. (2013), the metal-rich GC subpopulation rotates along with the stars, while the metal-poor GCs do not show significant rotation.

NGC 2974. Though classified as an E4, this galaxy exhibits some spiral structure in addition to more subtle substructure in the form of shells (Tal et al. 2009). It rotates rapidly within 1 R_e with corresponding anti-correlated h₃ values (Cinzano & van der Marel 1994; Emsellem et al. 2004). At larger radius the rotation amplitude plateaus out to 3 R_e along the major-axis, and also remains anti-correlated with h₃.

NGC 3115. The kinematics for this nearby S0 is described in Arnold et al. (2011), who also studied the GC kinematics. A rapidly rotating thin disk, apparent out to about 2 R_e, is embedded within a bulge that also exhibits fairly rapid rotation. Anti-correlated velocity and h₃ is apparent out to about 2 R_e in the bulge. As evidenced by the dispersion map, these two structures are both embedded in a more spheroidal and slowly rotating structure that only becomes evident at larger radius. Arnold et al. (2011) also found rapid metal-rich GC rotation marginally consistent with the stellar rotation in the inner part that decreases steeply farther out. The metal-poor GCs rotate with slightly lower amplitude and also exhibit an outer decline beyond 4 R_e.

NGC 3377. The combined stellar kinematic map of this elongated elliptical (E5–6) galaxy presents several key characteristics of an embedded disk within a more spheroidal component. The inner regions are dusty, have disky isophotes (e.g., Bender et al. 1988), rotate rapidly, and have depressed velocity dispersions compared to larger radii. The outer regions are boxier and rotate more slowly but remain aligned with the inner regions within about 4 R_e (see also Coccato et al. 2009). Pota et al. (2013) reported significant rotation along the major axis with outer twist of the metal-rich GCs in agreement with stars and long-slit data. The metal-poor GCs do not show significant rotation. The velocity dispersion of the metal-rich GCs is consistent with that of the stars while that of the metal-poor GCs is higher.

NGC 3608. This slightly flattened low luminosity elliptical galaxy (E1–2) shows only mild outer rotation along the major axis. It is one of the first known galaxies with a central counter-rotating kinematically decoupled core (Jedrzejewski & Schechter 1988). As such, it is one of the few non-regular rotators in the ATLAS^3D sample (Krajnović et al. 2011).

NGC 4111. This is the most disk-like galaxy (S0) in our sample and has the smallest effective radius (i.e., R_e <1 kpc). As a result, our kinematic data extend out to 4 R_e. In the inner parts, long slit data from Loyer et al. (1998) show a bump in the major axis rotation at a galactocentric radius of roughly 5 arcsec with otherwise constant amplitude farther out around 150 km s⁻¹. Our map shows that this rotation is sustained up to 4 R_e.

NGC 4278. This well-studied small round elliptical (E1–2) contains a distorted H i disk (Knapp et al. 1978). Our combined kinematic maps indicate that the stars are rotating only within 0.5 R_e, with very weak to no rotation farther out. The velocity dispersion falls rapidly with radius. Pota et al. (2013) did not find evidence that the metal-rich GCs rotate but concluded that the metal-poor GCs rotate mildly in the outer regions along a direction between the major and minor axes. The velocity dispersion of both GC subpopulations agrees broadly with the stellar kinematics.

NGC 4365. This giant elliptical (E3) galaxy is well known for its kinematically decoupled core rotating along the major axis and otherwise minor axis "rolling" rotation (Davies et al. 2001), which can be explained as a signature of long-axis tube orbits in a triaxial potential (Statler et al. 2004; van den Bosch et al. 2008). Such a configuration could arise through a major merger involving a spiral galaxy, with the central "normal" rotating component originating from newly formed stars, and the outer rolling component from the progenitor's cold disk (Hoffman et al. 2010). Our large-scale combined kinematic maps reveal that the minor-axis rolling continues to large radii and hence involves the bulk of the stellar mass. The galaxy also appears to be in the process of further halo growth via an interaction with NGC 4342 in the same group (Blom et al. 2014). Pota et al. (2013) found that the two "normal" GC subpopulations have identical velocity dispersion radial profiles. Blom et al. (2012) reported that the unusual intermediate metallicity subpopulation (green) of GCs rotates like the stars along the minor axis while both metal-rich and metal-poor GCs rotate along the major axis.

NGC 4374. This Virgo round elliptical (E1) galaxy harbors a rapidly rotating nuclear gas disk (Bower et al. 1998). Napolitano et al. (2011) used bright PNe as tracers to measure the stellar kinematics out to ∼350 arcsec. They found a steeply declining dispersion profile in the inner ∼30 arcsec beyond which the profile flattens out at a value of ∼230 km s⁻¹. As might be expected for such a slow rotator (Emsellem et al. 2011), we detect only mild to no rotation in the outer parts of the combined kinematic map, in agreement with Coccato et al. (2009).

NGC 4473. This elongated elliptical (E5) galaxy presents disky isophotes (Bender et al. 1988). It is an interesting case of a so-called "2σ" galaxy (Krajnović et al. 2011) which has off-center velocity dispersion peaks. The favored explanation for this feature is the presence of two counter-rotating disks (Cappellari et al. 2007). The inner rotation quickly declines with radius, eventually giving way to two-component retrograde major-axis and minor-axis rotation in the outskirts. This outer kinematic configuration is a tell-tale sign of triaxiality, and the galaxy overall may be considered a case of a "kinematically distinct halo" (Foster et al. 2013).

NGC 4486. Also known as M87, this well known bright central galaxy of Virgo is morphologically classified as an E0–1 and cD galaxy. Our kinematic map shows that the velocity dispersion falls with radius, though the values remain high at all radii probed (≳ 200 km s⁻¹). Some weak amplitude ∼20 km s⁻¹ off-axis rotation is apparent beyond 0.5 R_e (see also Murphy et al. 2011). Its GC kinematics are well-studied (e.g., Côté et al. 2001; Strader et al. 2011) and only weak rotation along the major axis is detected in both subpopulations. The overall GC velocity dispersion is only mildly declining at all probed radii. Helped by the sheer number of measured GC velocities, Romanowsky et al. (2012) found significant substructures in phase space indicating that this otherwise apparently quiet galaxy is actively accreting other systems.

NGC 4494. This "ordinary" elliptical (E1–2) galaxy contains an inner dust ring (Lauer et al. 2005). The central regions contain a kinematically distinct core (Bender et al. 1994) and double rotation amplitude maxima (Krajnović et al. 2011). Our map shows that rotation is sustained in the outskirts, along with a 10 deg twist in the kinematic position angle around 1.5 R_e (Proctor et al. 2009). In agreement with our velocity dispersion map, Napolitano et al. (2009) found a radially mild decrease in velocity dispersion using combined PN and long-slit data. The GC kinematics were studied in Foster et al. (2011), who found significant rotation for the metal-poor GCs consistent with the stellar motions, but little to no rotation in the metal-rich GCs. On the other hand, both GC subpopulations have a velocity dispersion profile consistent with that of the stars.

NGC 4526. This is an edge-on barred lenticular (SB0) galaxy with indications of possible spiral structure (Michard & Marchal 1994). Our combined velocity map exhibits a steep decrease to very low values below our resolution. The central regions of both the σ and h₄ maps suggest embedded components along the major and minor axes. There is a cold central disk (Bureau et al. 2002) in the inner parts that transitions into a hotter co-rotating bulge.

NGC 4649. Also known as M60, this massive Virgo E2/S0 galaxy is of comparable size to NGC 4486. The velocity maps show a clear contribution from the nearby spiral galaxy NGC 4647 to the northwest although there is no obvious sign of interaction. In contrast to M87, it exhibits rapid rotation with radially declining σ (also see Teodorescu et al. 2011 for a similar result using PNe). Kinematic maps for higher-order moments show clearly anti-correlated V − h₃ and radially dropping h₄. Das et al. (2011) used long-slit and PN tracers to model the mass of the system and concluded that the stellar orbits are centrally isotropic to mildly radial. Bridges et al. (2006) found no significant rotation of the GC system around NGC 4649, and a relatively flat velocity dispersion radial profile out to 3.5 R_e.

NGC 4697. This elongated elliptical (E6) galaxy exhibits fast rotation in the inner parts (Emsellem et al. 2011). Our extended kinematic map shows that this fast rotation is due to an embedded disk with a clear rotation amplitude decline with radius. The inner parts show a clear V − σ anti-correlation that vanishes in the outskirts. The velocity dispersion map is radially declining. The declining outer rotation and velocity dispersion radial profile was detected using PN tracers in Méndez et al. (2009). The PN data also imply mild radial anisotropy (de Lorenzi et al. 2008) and evidence has been found that they may in fact belong to two separate kinematic populations (Sambhus et al. 2006).

NGC 5846. This massive early-type (E0–1/S0) galaxy has quite circular isophotes (Michard & Marchal 1994). Our combined kinematic map shows that there is no obvious rotation in the central regions nor in the outskirts. There are clear contributions from the nearby galaxy NGC 5846A to the south in all four maps. Pota et al. (2013) demonstrated that neither GC subpopulation shows significant rotation apart from a marginal detection in the metal-poor GCs around 150–300 arcsec. They concluded that this may be due to "contaminating" GCs associated with NGC 5846A. The velocity dispersion radial profile of the metal-rich GCs is relatively flat, in contrast to that of the declining PNe (Das et al. 2008).

NGC 7457. This small S0 hosts a small bar (Michard & Marchal 1994) and a pseudo-bulge (Pinkney et al. 2003). As with other lenticular galaxies, it is classified as a fast rotator based on its inner kinematics (Emsellem et al. 2011). Our combined kinematic map shows global rotation (also see Cortesi et al. 2013) and anti-correlated h₃. The velocity dispersion values are universally low. The kinematics of this poor GC system have been presented by Pota et al. (2013), who found that the GCs are globally rapidly rotating, in agreement with the stars.

The conventional wisdom is that ETGs are primarily composed of old stars. However, accretion can bring in small amounts of gas that reinvigorate star formation for short periods of time, thus contaminating the predominantly old stellar population with some number of young stars. Given a dissipational origin, this young stellar component may then exhibit its own distinct kinematic structure; e.g., it could be both colder and more rapidly rotating than the preexisting stars. Since young stellar populations emit more of their light in the blue region of the spectrum, we were interested to know whether adding a young stellar component led to measured LOSVD shapes that were systematically different between blue and red spectral regions. Specifically, will a small population of young stars yield different kinematic measurements at blue wavelengths (4800–5380 Å) than at red wavelengths (8480–8750 Å)?

To test for this differential sensitivity to young stars, we created a set of composite spectra using the MILES single stellar population model spectra (Vazdekis et al. 2010). Each composite spectrum contains a particular fraction of younger (3 Gyr) stars ranging between 0% and 100%, with the balance made up of old (10 Gyr) stars. To keep the test simple, we kept both sub-populations at solar metallicity. To simulate the effect of a moderately rotating component, we convolved the 3 Gyr old spectrum with a Gaussian LOSVD with V = 100 km s⁻¹ and σ = 75 km s⁻¹. Conversely, the 10 Gyr population is meant to represent a dynamically hotter and slowly rotating spheroid, and was thus convolved with a Gaussian LOSVD with V = 0 km s⁻¹ and σ = 150 km s⁻¹.

If the young stellar population preferentially affects the blue spectral measurements, then we should detect a systematic offset in the derived kinematic parameters between the red and blue spectra. Specifically, we expect that the LOSVD parameters derived from the blue spectra will be biased toward the kinematics of the young stellar population. We used pPXF to measure the kinematic parameters for each composite spectrum separately in the blue and red wavelength ranges. The results are shown in Figure 15. Filled histograms (at left) show the LOSVDs for composite spectra with a young star fraction of 0, 25, 50, 75, and 100%. We note that young stellar populations are brighter than older populations, which results in a natural overrepresentation of younger stars in the LOSVD. The overplotted open histograms show the resulting LOSVDs if this luminosity difference were unaccounted for. Importantly, the relative brightnesses of young and old stellar populations is similar in the chosen blue and red spectral ranges. This means that red and blue spectra are roughly equally effected by young (3 Gyr) stellar populations, and thus no significant differential sensitivity to young stars is expected.

Zoom In Zoom Out Reset image size

The right panels illustrate this point a bit more clearly. In each panel we plot the difference in the measured kinematic quantity between the red and blue spectral regions versus its value as measured in just the blue spectral region. The curve follows our ensemble of composite spectra from a young-star fraction of 0% to 100% as shown by the color bar at the top left. When the young star fraction is 0, the spectrum contains only old stars and has a purely Gaussian LOSVD with zero mean velocity and a velocity dispersion of 150 km s⁻¹. As the young star fraction increases, the composite LOSVD becomes non-Gaussian (non-zero h₃ and h₄) until the population is purely made of young stars and the LOSVD is again Gaussian. Note that this effect is very non-linear in that a small amount of young stars is able to appreciably shift all four kinematic parameters—a consequence of the extremely bright intrinsic luminosity of young stellar-populations.

If the red and blue spectra were equally affected by the young stellar population, then each curve would be exactly horizontal. While this is not strictly the case, the deviations are small enough to be considered negligible given the typical uncertainties of our stellar kinematic measurements. In other words, there is relatively little differential sensitivity to stellar population effects between these wavelength regions. This simple test suggests that the discrepancies in measured kinematic values between the SAURON and SLUGGS data sets (see Section 3.2) are not caused by a differential sensitivity to young stellar populations in the different wavelength ranges used in each survey.