The MASSIVE Survey XIII -- Spatially Resolved Stellar Kinematics in the Central 1 kpc of 20 Massive Elliptical Galaxies with the GMOS-North Integral-Field Spectrograph

TL;DR

This study uses high-quality integral-field spectroscopy to map stellar kinematics in the central regions of 20 massive elliptical galaxies, revealing detailed velocity structures and black hole influences.

Contribution

It provides unprecedented spatially resolved stellar kinematic maps of massive ETGs, combining GMOS-North IFS data with previous wide-field observations for comprehensive analysis.

Findings

Most galaxies show increasing velocity dispersion towards the center.

Positive h4 values suggest central black holes and velocity anisotropy.

Jeans modeling constrains black hole mass and dark matter halo properties.

Abstract

We use observations from the GEMINI-N/GMOS integral-field spectrograph (IFS) to obtain spatially resolved stellar kinematics of the central $\sim 1$ kpc of 20 early-type galaxies (ETGs) with stellar masses greater than $10^{11.7} M_\odot$ in the MASSIVE survey. Together with observations from the wide-field Mitchell IFS at McDonald Observatory in our earlier work, we obtain unprecedentedly detailed kinematic maps of local massive ETGs, covering a scale of $\sim 0.1-30$ kpc. The high ($\sim 120$) signal-to-noise of the GMOS spectra enable us to obtain two-dimensional maps of the line-of-sight velocity, velocity dispersion $\sigma$, as well as the skewness $h_3$ and kurtosis $h_4$ of the stellar velocity distributions. All but one galaxy in the sample have $\sigma(R)$ profiles that increase towards the center, whereas the slope of $\sigma(R)$ at one effective radius ($R_e$) can be of either sign. The $h_4$ is generally positive, with 14 of the 20 galaxies having positive $h_4$ within the GMOS aperture and 18 having positive $h_4$ within $1 R_e$. The positive $h_4$ and rising $\sigma(R)$ towards small radii are indicative of a central black hole and velocity anisotropy. We demonstrate the constraining power of the data on the mass distributions in ETGs by applying Jeans anisotropic modeling (JAM) to NGC~1453, the most regular fast rotator in the sample. Despite the limitations of JAM, we obtain a clear $\chi^2$ minimum in black hole mass, stellar mass-to-light ratio, velocity anisotropy parameters, and the circular velocity of the dark matter halo.

Tables

Table 1. Table 1: Galaxy properties and GMOS observational details for the 20 MASSIVE galaxies.

Galaxy	$D$	$M_K$	PAphot	$R_e$	${\lambda }_e$	${\sigma }_e$	${\lambda }_{1kpc}$	${\sigma }_{1kpc}$	${\gamma }_{1kpc}$	IFU PA	Semester	Exposure
	[Mpc]		[deg]	[′′ (kpc)]		[km s-1]		[km s-1]		[deg]
NGC0057	76.3	-25.75	40.2	17.1 (6.31)	0.028	257 $\pm$ 2	0.025	301 $\pm$ 2	-0.120 $\pm$ 0.006	41	2016B	12 $\times$ 1150s
NGC0315	70.3	-26.30	44.3	27.0 (9.20)	0.063	341 $\pm$ 1	0.027	325 $\pm$ 1	-0.025 $\pm$ 0.007	0	2012B	10 $\times$ 1200s
NGC0410	71.3	-25.90	35.8	21.9 (7.57)	0.048	258 $\pm$ 1	0.052	288 $\pm$ 2	-0.070 $\pm$ 0.006	35	2016B	7 $\times$ 1150s
NGC0545	74.0	-25.83	57.2	27.1 (9.71)	0.081	236 $\pm$ 2	0.034	236 $\pm$ 1	-0.101 $\pm$ 0.005	60	2013B	8 $\times$ 1200s
NGC0547	71.3	-25.90	98.8	30.5 (10.55)	0.081	230 $\pm$ 2	0.024	246 $\pm$ 3	-0.139 $\pm$ 0.008	94	2016B	8 $\times$ 1200s
NGC0741	73.9	-26.06	88.0	27.2 (9.74)	0.050	289 $\pm$ 3	0.037	274 $\pm$ 1	-0.025 $\pm$ 0.010	177	2012B	6 $\times$ 1200s
NGC0777	72.2	-25.94	148.6	16.8 (5.89)	0.060	293 $\pm$ 2	0.027	320 $\pm$ 2	-0.068 $\pm$ 0.004	148	2015B	8 $\times$ 1050s
NGC0890	55.6	-25.50	53.7	24.5 (6.62)	0.136	196 $\pm$ 2	0.027	206 $\pm$ 1	0.013 $\pm$ 0.007	55	2015B	6 $\times$ 850s
NGC1016	95.2	-26.33	42.8	20.5 (9.47)	0.040	279 $\pm$ 1	0.015	300 $\pm$ 2	-0.069 $\pm$ 0.007	40	2013B	6 $\times$ 600s
NGC1060	67.4	-26.00	74.8	19.5 (6.38)	0.048	282 $\pm$ 2	0.034	301 $\pm$ 2	-0.050 $\pm$ 0.008	74	2016B	7 $\times$ 1150s
NGC1129	73.9	-26.14	61.7	45.0 (16.13)	0.124	267 $\pm$ 2	0.350	228 $\pm$ 3	-0.105 $\pm$ 0.031	5	2016B	18 $\times$ 1200s
NGC1453	56.4	-25.67	30.1	21.9 (6.00)	0.204	276 $\pm$ 1	0.199	293 $\pm$ 2	-0.082 $\pm$ 0.011	20	2015B	6 $\times$ 850s
NGC1573	65.0	-25.55	31.7	17.2 (5.43)	0.056	270 $\pm$ 2	0.026	282 $\pm$ 2	-0.057 $\pm$ 0.011	35	2013B	6 $\times$ 600s
NGC1600	63.8	-25.99	8.8	29.6 (9.14)	0.035	299 $\pm$ 1	0.045	337 $\pm$ 1	-0.048 $\pm$ 0.008	15	2014B	9 $\times$ 1230s
NGC1700	54.4	-25.60	90.6	16.9 (4.45)	0.198	227 $\pm$ 1	0.119	227 $\pm$ 1	-0.062 $\pm$ 0.005	90	2015B	13 $\times$ 850s
NGC2258	59.0	-25.66	150.8	20.2 (5.76)	0.071	258 $\pm$ 3	0.034	285 $\pm$ 1	-0.084 $\pm$ 0.007	135	2016B	14 $\times$ 1200s
NGC2274	73.8	-25.69	165.0	18.4 (6.57)	0.073	261 $\pm$ 1	0.042	270 $\pm$ 1	-0.077 $\pm$ 0.007	169	2015B	8 $\times$ 1050s
NGC2340	89.4	-25.90	80*	32.9 (14.27)	0.032	234 $\pm$ 1	0.042	230 $\pm$ 2	-0.092 $\pm$ 0.013	170	2012B	11 $\times$ 1200s
NGC2693	74.4	-25.76	161.3	15.6 (5.63)	0.294	296 $\pm$ 2	0.337	291 $\pm$ 2	-0.073 $\pm$ 0.011	166	2016B	6 $\times$ 1200s
NGC4874	102.0	-26.18	40*	38.8 (19.20)	0.072	260 $\pm$ 1	0.018	270 $\pm$ 2	-0.110 $\pm$ 0.008	145	2015A	18 $\times$ 1220s

Table 2. Table 2: Power law slopes of velocity dispersion profiles for 5 MASSIVE galaxies.

Galaxy	${\gamma }_{1kpc}$	${\gamma }_e$	${\gamma }_{\mathrm{Loubser}}$
NGC0315	$-0.025\pm 0.007$	$-0.041\pm 0.023$	$-0.028\pm 0.011$
NGC0410	$-0.070\pm 0.006$	$-0.216\pm 0.033$	$-0.068\pm 0.009$
NGC0777	$-0.068\pm 0.004$	$-0.110\pm 0.044$	$-0.078\pm 0.010$
NGC1060	$-0.050\pm 0.008$	$-0.081\pm 0.041$	$-0.067\pm 0.010$
NGC1453	$-0.082\pm 0.011$	$-0.039\pm 0.031$	$-0.068\pm 0.011$

Table 3. Table 3: Gaussian components of the MGE fit to the HST WFC3 IR photometry of NGC 1453 (shown in Figure 13 ). All 11 Gaussians share the same center and position angle of 28.5 degrees. From left to right: central surface density I k = L k / 2 ​ π ​ σ k 2 ​ q k ′ subscript 𝐼 𝑘 subscript 𝐿 𝑘 2 𝜋 superscript subscript 𝜎 𝑘 2 superscript subscript 𝑞 𝑘 ′ I_{k}=L_{k}/2\pi\sigma_{k}^{2}q_{k}^{\prime} (calculated using an absolute solar magnitude M ⊙ , F110W = 4.54 subscript 𝑀 direct-product F110W 4.54 M_{\odot{\rm,F110W}}=4.54 ), dispersion σ k subscript 𝜎 𝑘 \sigma_{k} in arcseconds and axis ratio q k ′ superscript subscript 𝑞 𝑘 ′ q_{k}^{\prime} of each Gaussian component.

$I_k\left[L_{\odot }/{\mathrm{pc}}^2\right]$	${\sigma }_k{[}^{\prime \prime }]$	$q_k^{\prime }$
$23280.1$	$0.041$	$0.891$
$9851.1$	$0.270$	$0.882$
$12303.7$	$0.575$	$0.900$
$9798.8$	$0.963$	$0.816$
$6734.01$	$1.595$	$0.800$
$4000.64$	$2.370$	$0.865$
$2235.06$	$3.403$	$0.800$
$1526.44$	$5.743$	$0.835$
$591.08$	$10.292$	$0.800$
$277.19$	$20.778$	$0.821$
$79.34$	$47.289$	$0.897$

Equations

$$f(v)=\frac{e^{-\frac{y^{2}}{2}}}{\sqrt{2\pi\sigma^{2}}}\bigg{[}1+\sum_{m=3}^{n}h_{m}H_{m}(y)\bigg{]},$$

$$\frac{{\sigma }_{e/8}}{{\sigma }_e}=1.044\pm \mathrm{0.078.}$$

$$\Sigma(x^{\prime},y^{\prime})=\sum_{k=1}^{N}\frac{L_{k}}{2\pi\sigma_{k}^{2}q_{k}^{\prime}}{\rm exp}\Big{[}-\frac{1}{2\sigma_{k}^{2}}\Big{(}x^{\prime 2}+\frac{y^{\prime 2}}{q_{k}^{\prime 2}}\Big{)}\Big{]},$$

$$\Phi_{\rm DM}(r)=\frac{1}{2}V_{c}^{2}{\rm ln}\Big{(}r^{2}+R_{c}^{2}\Big{)},$$

$${\rho }_{\mathrm{DM}}(r)=\frac{V_c^2}{4\pi G}\frac{3R_c^2+r^2}{(R_c^2+r^2{)}^2}.$$

$${\beta }_z=1-\frac{\overline{v_z^2}}{\overline{v_R^2}},$$

Full text

galaxies: elliptical and lenticular, cD — galaxies: evolution — galaxies: formation — galaxies: kinematics and dynamics — galaxies: structure

1 Introduction

As the final product of multiple merger events, massive early-type galaxies (ETGs) in the local universe provide excellent insight into how galaxies evolve. ETGs are complex, multi-component systems, and a full description of their evolutionary processes must take into account the stars, dark matter, supermassive black holes (SMBHs), and any gas present in the galaxies.

Significant recent progress on understanding local ETGs has been made by surveys using integral field spectrographs (IFS), e.g., SAURON (Emsellem et al., 2004), ATLAS3D (Cappellari et al., 2011), SAMI (Croom et al., 2012), CALIFA (Sánchez et al., 2012), and MaNGA (Bundy et al., 2015). These surveys use wide-field IFS to produce two-dimensional maps of the stellar and gas kinematics and investigate fundamental galaxy properties such as the dichotomy between fast and slow rotators, early-type galaxy morphologies, and molecular and ionized gas content. The ETGs targeted by these surveys are predominantly fast-rotating S0 or elliptical galaxies with $M_{*}<10^{11.5}M_{\odot}$. The spatial sampling of these survey is limited by the IFS fiber diameter ($1.6^{\prime\prime}$, $2^{\prime\prime}$, $2.7^{\prime\prime}$ for SAMI, MaNGA, CALIFA, respectively) or lenslet size ($0.94^{\prime\prime}$ for SAURON/ATLAS3D). A few other recent IFS or long-slit studies of a smaller number of ETGs specifically targeted brightest cluster galaxies (BCGs; e.g. Brough et al. 2011; Jimmy et al. 2013; Loubser et al. 2018), and the SLUGGS survey used multi-slits to reach a sky coverage of $\sim 2-4$ effective radii for a subset of ATLAS3D galaxies (Brodie et al., 2014).

We initiated the volume-limited MASSIVE survey (Ma et al. 2014) to investigate the $\sim 100$ most massive galaxies located up to a distance of 108 Mpc in the northern sky. The survey targets a complete sample of ETGs with an absolute $K$-band magnitude brighter than $M_{K}=-25.3$ mag, or a stellar mass greater than $M_{*}\sim 10^{11.5}M_{\odot}$, a parameter space little explored previously. Wide-field kinematics and stellar population studies from this survey have been published in Greene et al. (2015), Pandya et al. (2017), Veale et al. (2017b, a, 2018), Ene et al. (2018), and Greene et al. (2019) based on IFS data from the Mitchell Spectrograph on the 2.7-m Harlan J. Smith Telescope at McDonald Observatory. These two-dimensional kinematics have a spatial resolution of $\sim 4^{\prime\prime}$ (fiber diameter) and cover a field of view (FOV) of $107^{\prime\prime}\times 107^{\prime\prime}$. Paper V (Veale et al. 2017b) found a dramatic increase in the fraction of slow-rotating ETGs with increasing $M_{*}$, reaching $\sim 90$% at $M_{*}\sim 10^{12}M_{\odot}$. Paper VII (Veale et al., 2017a) examined the relationship of galaxy spin, $M_{*}$ and environment and found that galaxy mass, rather than environment, is the primary driver of the apparent kinematic morphology-density relation for local ETGs. The physical processes responsible for building up the present-day stellar masses of massive ETGs must be very efficient at reducing their spin, in any environment. Paper VIII (Veale et al., 2018) investigated the environmental dependence of the stellar velocity dispersion profiles and found the fraction of galaxies with rising outer profiles to increase with $M_{*}$ and environmental density, a trend likely due to total mass variations rather than velocity anisotropies. Paper X (Ene et al. 2018) analyzed substructures in the stellar velocity maps and found kinematic twists and large scale ($R\gtrsim 10^{\prime\prime}$) kinematically distinct components.

In this paper (Part XIII), we present the first results from the high angular-resolution spectroscopic portion of the MASSIVE survey. While the earlier wide-field IFS studies offered insight to galaxies’ global dynamics and assembly histories, kinematics of the innermost regions of galaxies are critical for measuring the masses of the SMBHs and for elucidating the symbiotic relations among black holes, baryons and dark matter near galactic centers. To achieve these goals, we observed the central $5^{\prime\prime}\times 7^{\prime\prime}$ of 20 MASSIVE galaxies using Gemini Multi Object Spectrograph (GMOS, Hook et al. 2004) with $0.2^{\prime\prime}$ lenslets on the 8.1 m Gemini North Telescope. The exposure times were chosen to yield stellar spectra with signal-to-noise ratio (S/N) of $\sim 120$ for each spatial bin. These high-quality spectra enable us to obtain detailed two-dimensional maps of the velocity, velocity dispersion, as well as the skewness and kurtosis of the stellar line-of-sight velocity distribution (LOSVD). Depending on the galaxy, the maps contains from 50 to 300 spatial bins (with an average of 130 bins) and cover a scale from $\sim 0.1$ kpc to a few kpc.

The size of the sample and the IFS spatial coverage in this paper is similar to that of McDermid et al. (2006), who studied the central $8^{\prime\prime}\times 10^{\prime\prime}$ region of 28 ETGs in the SAURON survey using the OASIS spectrograph with a spatial sampling of $0.27^{\prime\prime}$. The fine spatial sampling enabled them to identify two types of kinematically distinct components (KDCs) in lower-mass ETGs: old, kpc-scale KDCs that exist in slow rotating ETGs, and young, small ($\sim 100$ pc scale), almost counter-rotating KDCs in fast rotators. However, their kinematics were obtained from spectra with a lower S/N of 60. Their sample is in the $M_{*}$ range of $10^{10}-10^{11.6}M_{\odot}$, where 68% are fast rotators and many show emission lines. By contrast, the galaxies studied here are mostly slow rotators and few have emission lines. There is no overlapping galaxy with the two samples

The remaining sections of the paper are organized as follows. In Section 2 we present the sample of 20 MASSIVE galaxies with high resolution IFS observations. In Section 3 we describe the observations and the data reduction pipeline. Section 4 provides an overview of how we use the IFS data sets to measure stellar kinematics. Sections 5, 6, and 7 examine the behavior of the stellar kinematics in the galaxies’ nuclei: Section 5 explores the velocity profiles, Section 6 looks at the radial behavior of velocity dispersion, and Section 7 studies the higher order moments $h_{3}$ and $h_{4}$. In Section 8 we showcase how the combined set of small and large scale kinematics can be use to constrain dynamical models of the fast rotator NGC 1453. Section 9 summarizes our main conclusions.

2 Galaxy sample

An in-depth description of the selection of the parent sample of 116 galaxies is given in Ma et al. (2014). In summary, MASSIVE is a volume-limited survey of the most luminous early-type galaxies with $M_{K}\lesssim-25.3$ mag (from 2MASS; Skrutskie et al. 2006) corresponding to $M_{*}\gtrsim 10^{11.5}M_{\odot}$, that are within a distance of 108 Mpc in the northern sky (declination $\delta>-6\arcdeg$).

In this paper we present results for 20 galaxies that were chosen for follow-up observations with high spatial resolution spectroscopy with GMOS. The key properties of these galaxies are summarized in Table 1. The 20 galaxies are located at a distance in the range of 54.4 Mpc to 102.0 Mpc (with a median distance of $\sim 70$ Mpc) and have $-25.50~{}{\rm mag}\geq M_{K}\geq-26.33~{}{\rm mag}$, which corresponds to stellar mass $10^{11.7}M_{\odot}\lesssim M_{*}\lesssim 10^{12}M_{\odot}$. Since they were selected based on the ability to obtain high-$S/N$ GMOS data for dynamical mass modeling, this sample of 20 galaxies tends to be the closer and more massive part of the general MASSIVE sample: they represent $\sim 50\%$ of galaxies within 80 Mpc and $\sim 60\%$ of galaxies more massive than $M_{K}<-25.8~{}{\rm mag}$.

3 GMOS-N observations and data reduction

Our galaxies were observed using the GMOS integral field unit (IFU, Allington-Smith et al. 2002) on the 8.1 m Gemini North telescope. The observations were taken in queue mode over six semesters between 2012 – 2016 under the programs GN-2012B-Q-31, GN-2013B-Q-29, GN-2013B-Q-68, GN-2014B-DD-6, GN-2015A-Q-19, GN-2015B-Q-59, and GN-2016B-Q-18. Total exposure times were chosen to ensure a S/N of $\sim 120$ (after spatial binning) and vary between 1 and 6 hours (see Table 1). All galaxies presented here were observed after the GMOS-N upgrade to e2v deep depletion detectors in 2011 (Roth et al., 2012), but before the upgrade to Hamamatsu fully-depleted detectors in 2017 (Scharwächter et al., 2018).

All observations were taken using the two-slit mode of the GMOS-N IFU. This provides a FOV of $5^{\prime\prime}\times 7^{\prime\prime}$ consisting of 1000 hexagonal lenslets with a projected diameter of $0.2^{\prime\prime}$. An additional 500 lenslets observe a $5^{\prime\prime}\times 3.5^{\prime\prime}$ region of the sky displaced by $\sim 1$ arcmin from the science field. The lenslets are coupled to fibers that map the focal plane into two pseudo-slits (each covering half of the FOV) through which light is passed to the rest of the spectrograph. Each pseudo-slit covers the same spectral range and is projected across the full spatial dimension of the detector array (perpendicular to the dispersion direction). In the spectral dimension of the array, the two pseudo-slits are projected with an offset in the central wavelength and thereby avoid overlap if the spectral range is sufficiently narrow. We use the R400-G5305 grating + CaT filter combination to avoid spectral overlap on the detector. This results in a clean wavelength range of 7800 Å – 9330 Å that has good coverage of the CaII triplet and NaI absorption features used for measuring stellar kinematics and stellar populations, respectively.

The detector array consists of three 2k $\times$ 4k e2v deep depletion CCDs placed in a row with $\sim 37$ unbinned pixels gaps in between. To mitigate the loss of spectral information to the chip gaps, we use spectral dithering with a grating central wavelength $\lambda_{c}$ for half of the exposures and $\lambda_{c}+50$ Å for the other half. For most galaxies, a typical value of $\lambda_{c}$ is between 8600 Å and 8700 Å. We carefully choose this value to ensure that the CaII triplet lines do not fall on either of the two chip gaps.

The basic reduction of the raw data frames is performed using the Gemini package within the image reduction and analysis facility (IRAF) software. For an in-depth example of how to reduce GMOS IFU data using IRAF (and potential pitfalls) see Lena (2014). We use the standard GMOS data reduction procedure. The science, flat lamp, and twilight raw frames are bias subtracted. The arc frames are taken in fast-read mode so they are only overscan subtracted. The Gemini calibration unit (GCAL) flat lamp frames (taken before and after the science exposures) are used to identify and extract the trace of each fiber on the detector array and to determine the flat-field response map. The twilight exposure is used to correct for illumination. The GCAL arc lamp frames are used to determine the wavelength solution by fitting a fourth-order Chebyshev polynomial to known CuAr arc lamp lines spanning the wavelength range of the observations. The science spectra are extracted using the fiber traces identified in the flat lamp exposures. The spectra are then corrected for flat-fielding and illumination. Cosmic ray artifacts are removed from the spectra and spectral processing is performed using the arc lamp wavelength solution. We use a custom routine to perform sky subtraction. For each pseudo-slit, we subtract an average spectrum computed from the dedicated sky fibers corresponding to that particular pseudo-slit. The end result of this step is a reduced science frame that contains the one-dimensional spectrum corresponding to each GMOS lenslet.

For all galaxies observed in semester 2013B or later, arc lamp exposures were recorded immediately adjacent to science exposures of each galaxy, with the telescope at the same pointing. For three galaxies observed in semester 2012B (NGC 315, NGC 741, and NGC 2340), arc lamp exposures were recorded in daytime with the telescope parked. The arc-calibrated science frames from 2012B exhibit residual wavelength errors, most readily apparent as a sharp wavelength offset between bright sky lines in the two GMOS pseudo-slits. For each arc-calibrated science frame (prior to sky subtraction), we parameterized the residual wavelength offset $\epsilon_{\lambda}$ in each GMOS lenslet by fitting $\sim 10$ bright sky lines and fitting a polynomial function $\epsilon_{\lambda}(\lambda)$ across the observed wavelength range. To mitigate the low signal-to-noise of sky lines in individual lenslets, we fit $\epsilon_{\lambda}(\lambda)$ as a first-order polynomial in each lenslet independently. The resulting two-dimensional map of $\epsilon_{\lambda}$ in each science frame is interpolated to a wavelength of interest (i.e. centered on the CaII triplet), smoothed using a 20-lenslet boxcar, and applied to a convolution kernel for the stellar template spectrum during kinematic fitting (Section 4). These calibration steps simultaneously measure the wavelength- and lenslet-dependent instrumental resolution $\Delta\lambda$. Since we account for this instrumental term during kinematic fitting and we ultimately measure kinematics from co-added galaxy spectra, the corresponding instrumental kernels for ($\epsilon_{\lambda}$, $\Delta\lambda$) in individual lenslets are co-added as well.

The individual reduced science exposures are not stacked or mosaicked. Instead, we use a suite of custom routines to extract and co-add one-dimensional spectra from multiple exposures. We construct collapsed (along the wavelength direction) images of the galaxy and determine the location of the galactic center by fitting a Moffat profile to the light profiles of the collapsed images. Then, we extract the one-dimensional spectra and tag each of them by the exposure number and that lenslet’s spatial coordinates relative to the galaxy center. We interpolate the extracted spectra to a common wavelength basis and perform a heliocentric velocity correction if the exposures were recorded over multiple nights.

To increase the S/N of the data, we use the Voronoi binning procedure of Cappellari & Copin (2003). The one-dimensional spectra are irregularly spaced in galactocentric coordinates, due to dithers and pointing offsets between individual frames, and the hexagonal shape of the IFU lenslets. Since Voronoi binning requires regularly spaced coordinates, we construct a square grid in (x,y) with a grid spacing of $0.2^{\prime\prime}$, equal to the width of a hexagonal lenslet. We then tag each spectrum to the nearest grid point. We do not consider overlap with multiple grid points – each spectrum is given 100% weight at a single grid position.

The Voronoi binning procedure does not co-add the spectra, but merely defines the bins based on the estimated signal and noise of each point on the regular grid. For each grid point, we estimate the fluxes and the pixel-to-pixel variance of all contributing spectra by using the residuals between the observed spectra and boxcar-smoothed spectra over a 10 pixel window. The signal assigned to each grid point is then the sum of the fluxes of all contributing spectra, while the noise is the quadrature sum of the contributing spectra noise. We use a custom implementation of the binning step that imposes spatial symmetry over four galaxy quadrants. The bin boundaries are then used to create symmetric bins in the remaining three quadrants. The data (i.e. the 1-D spectra) is never folded during this step. We use a Voronoi binning target S/N of 125, which results in high quality spectra that do not sacrifice the spatial resolution of the innermost bins. The resulting S/N per bin values generally scatter around the target, with a typical RMS scatter of $\sim 10\%$, as can be seen in Fig. 1.

In most cases, a Voronoi bin is composed of multiple $0.2^{\prime\prime}\times 0.2^{\prime\prime}$ grid segments, each affiliated with multiple one-dimensional spectra, usually from different exposures. The final step is to co-add all of the one-dimensional spectra in each Voronoi bin and to create a corresponding Gaussian kernel for the instrumental resolution of the co-added spectrum. This is done with a clipped 3-$\sigma$ mean and rescaling for regions overlapping the chip gaps at one of the two wavelength settings.

4 Measuring stellar kinematics

The stellar LOSVD is extracted using the penalized pixel-fitting (pPXF) method of Cappellari & Emsellem (2004), which convolves a set of template stellar spectra with the LOSVD function $f(v)$. The latter is modeled as a Gauss-Hermite series of order up to $n=6$ (Gerhard 1993; van der Marel & Franx 1993):

[TABLE]

where $y=(v-V)/\sigma$, $V$ is the mean velocity, $\sigma$ is the velocity dispersion, and $H_{m}$ is the $m^{\text{th}}$ Hermite polynomial as defined in Appendix A of van der Marel & Franx (1993).

For each binned spectrum we fit all six Gauss-Hermite moments, although for brevity we only show the first four moments in our kinematic plots. We run pPXF with an initial guess of 0 for $V$ and $h_{3}$ through $h_{6}$, and 300 km s*-1* for $\sigma$. As our model continuum fit, we use an additive polynomial of degree zero (i.e. an additive constant) and a multiplicative polynomial of degree three. In cases where the LOSVD is undersampled or the S/N is low, it is important to use the pPXF penalty term to suppress the large uncertainties of the higher order moments (Cappellari & Emsellem, 2004). Since the galaxies in our sample have large velocity dispersions and the data have high S/N, it is not critical that we penalize deviations from a Gaussian solution – hence we set the pPXF keyword BIAS to zero.

Prior to fitting, we center the spectra on the triplet of calcium absorption features by cropping close to the wavelength range 8420 – 8770 Å. We also mask any prominent residual sky lines. That is, we mask small wavelength regions centered on the locations where the sky lines occur. Overall, this corresponds to excluding $\sim 10\%$ or less of the fit region. We present example pPXF fits at the nucleus, in an intermediate and in an outer spatial bin for NGC 1700 in Fig. 2.

Barth et al. (2002) tested the robustness of using the Ca triplet spectral region for velocity dispersion measurements. Similar to Dressler (1984), they found little sensitivity in the measurements to the choice of template stars. We performed our own template mismatch tests using two different sets of template stars chosen from the CaT Library of 706 stars in Cenarro et al. (2001). The first template set contains the identical 15 K and G stars as in Table 2 of Barth et al. (2002). The second template set contains all $\sim 360$ G and K stars in the CaT Library. For the latter, we ran pPXF for all the bins in one of our 20 galaxies and examined the 40 stars that were assigned the highest weights. Only one of the 40 stars is in common with the 15 stars in Barth et al. (2002). Despite the difference in the two template choices, we find that the rms difference in the kinematic moments measured from the two sets of templates is $\sim 5$ km s*-1* for $V$ and $\sigma$, and $\sim 0.01$ for $h_{3}$ and $h_{4}$, well within the measurement errors of $\sim 10$ km s*-1* and $\sim 0.02$, respectively. We therefore confirm that our kinematic measurements are robust to template choices. The results reported in this paper use the 15 template stars in Barth et al. (2002).

We note that we performed similar tests in Veale et al. (2017b) for the Mitchell IFU data in the wavelength range $3650-5850$ Å. The results there also indicated that our kinematic measurements are relatively insensitive to the choice of input template library.

The stellar spectra in the CaT Library cover the wavelength range 8348 – 9020 Å with a spectral resolution of 1.5 Å FWHM. To match the spectral resolution of our binned spectra, we convolve the templates with a Gaussian distribution with appropriate dispersion. We determine the instrumental resolution of our data by using the known sky lines prior to sky subtraction, as described in Section 3, or arc lines from the CuAr calibration lamp. We first fit a Gaussian to each individual sky or arc line. Then we fit a low-order polynomial to the full width at half-maximum (FWHM) of the lines versus wavelength and determine the corresponding FWHM for our fit region centered at $\sim 8600$ Å. We determine this best-fitting FWHM for each individual lenslet of each exposure. Typical values for the FWHM are $\sim 2.5$ Å, with variations of $\sim 0.3$ Å with lenslet position. This corresponds to an instrumental resolution of $\sim 37$ km s*-1* at 8600 Å with a sampling of 0.67 Å per pixel. Finally, we generate a Gaussian kernel for each binned spectrum by averaging the kernels of each individual lenslet assigned to that bin. While it is possible for the line spread function to deviate mildly from a Gaussian shape (e.g., slightly flat-topped for the KCWI spectrograph; Morrissey et al. 2018; van Dokkum et al. 2019), in practice the uncertainty in the LSF introduces non-negligible bias in the recovery of the kinematic moments only when the measured velocity velocity dispersion is smaller than the instrumental dispersion (Cappellari, 2017). Typical velocity dispersions for MASSIVE galaxies are $\gtrsim 250$ km s*-1*, much higher than the spectral resolution of GMOS ($\sim 40$ km s*-1*). Uncertainties in the LSF are therefore subdominant to other sources of systematic error.

The error bars on the kinematic moments are obtained through a bootstrap approach. This choice is motivated by the fact that each Voronoi bin contains tens (inner bins) to hundreds (outer bins) of individual lenslet spectra from different science exposures that have different noise properties. A bootstrap trial consists of taking the individual lenslet spectra that belong to a Voronoi bin and drawing a sample with replacement, then using this sample to generate a new co-added spectrum, which we run through pPXF to determine the best-fitting Gauss-Hermite moments. We repeat this process 100 times. Finally, the error for each kinematic moment is computed as the standard deviation of the pPXF results for the 100 bootstrap trials. Using this bootstrap approach, we find errors that are $50\%-100\%$ larger than the errors from a simpler Monte Carlo approach.

The Monte Carlo approach assumes the noise of the co-added galaxy spectrum to be representative of the noise for all the individual lenslets spectra that make up the co-added spectrum. Since the individual lenslets that make up a bin often come from exposures that were taken several days (or months) apart, it is very likely that there exist significant differences in the noise properties from one lenslet spectrum to the next (caused by variations in nighttime conditions, instrument performance, or spatially within the galaxy). This variation is captured only to a lesser extent by the noise estimates that are used in the Monte Carlo approach, since these are generated from the final spectrum after co-adding all the input lenslet spectra. We believe that the errors from the bootstrap approach are instead more realistic since this approach doesn’t make any assumptions about the noise properties of the data and incorporates any systematic differences between spectra extracted from different science frames and lenslets prior to co-adding.

5 Velocity features

Fig. 3 (left panel) shows the radial velocity profiles measured from our GMOS IFS data along the photometric major axis of the central region of each of the 20 galaxies. Sixteen of the 20 galaxies have low rotation velocities with $|V|\lesssim 30$ km s*-1*. Three galaxies exhibit fast central rotation, with $|V|$ reaching values of $\sim 100-150$ km s*-1*: NGC 1129 (light blue), NGC 1453 (orange), and NGC 2693 (magenta). The rotation speeds at $R\sim 1$ to 1.5 kpc ($3^{\prime\prime}$ to $4^{\prime\prime}$) of these three galaxies in Fig. 3 are maintained at a similar level out to one effective radius and beyond in our wide-field Mitchell IFS data, which measured a maximum $|V|$ of 70 km s*-1*, 95 km s*-1*, and 150 km s*-1* for NGC 1129, NGC 1453, and NGC 2693, respectively (Table 1, Ene et al. 2018).

To quantify the importance of rotation relative to dispersion in the central region probed by the GMOS IFS, we compute the spin parameter $\lambda$ measured within a circular aperture of radius 1 kpc, $\lambda_{\rm 1kpc}$, where $\lambda(<R)\equiv\langle R|V|\rangle/\langle R\sqrt{V^{2}+\sigma^{2}}\rangle$, and the brackets refer to luminosity-weighted averages. Individual values of $\lambda_{\rm 1kpc}$ are listed in Table 1 (column 8) and plotted in the right panel of Fig. 3. When compared to the larger-scale spin measured within one effective radius, $\lambda_{e}$ (column 6 of Table 1), Fig. 3 (right panel) shows that the galaxies with higher $|V|$ featured in the left panel also have larger $\lambda_{\rm 1kpc}$ as well as larger $\lambda_{e}$. The majority of our sample galaxies, however, has both low central and low global spins and clusters in a distinct part of the $\lambda_{\rm 1kpc}$-$\lambda_{e}$ parameter space in Fig. 3.

The fourth galaxy highlighted in the left panel of Fig. 3, NGC 1700 (dark blue), shows a prominent KDC in our GMOS data. This component rotates in exactly the opposite direction of the main body rotation (see Fig. 4). The existence of this counter rotating core (CRC) was first pointed out by Franx et al. (1989) and later examined through detailed stellar kinematics up to four effective radii of Statler et al. (1996). Both studies used long-slit spectroscopic observations and neither achieved the high angular resolutions and high $S/N$ presented here.

Previous works have also uncovered kinematically distinct features at the center of ETGs. Ricci et al. (2016) studied the stellar kinematics of the circumnuclear regions of 10 massive, nearby ETGs and found a KDC with an extent of $\sim 200$ pc, while McDermid et al. (2006) found that some fast rotating SAURON ETGs have central small scale counter-rotating KDCs and that these KDCs have stellar populations that are younger than those of the main galaxy body. This is also the case for the NGC 1700 KDC, which Kleineberg et al. (2011) found to have a distinct stellar population, much younger than the main body and suggest that it formed as the result of a merger between the main galaxy and a small stellar companion on a retrograde orbit. On larger scales ($R\sim 10^{\prime\prime}$), two KDCs in MASSIVE galaxies NGC 507 and NGC 5322 were detected in our wide-field Mitchell IFS data (Ene et al., 2018). The rotation axis of the inner component of NGC 507 was misaligned by $\sim 105\arcdeg$ from the rotation axis of the outer component, while for NGC 5322 the two components are misaligned by almost $180\arcdeg$. We noted in Ene et al. (2018) that more MASSIVE galaxies may contain KDCs on smaller scales not resolvable by the $4.1^{\prime\prime}$ diameter fibers of the Mitchell IFS. NGC 1700 is one such example.

Kinematic features such as central KDCs and any misalignment between the kinematic and the photometric axes provide useful clues about the merger history of massive ellipticals. Visual inspections of the 20 GMOS velocity maps in Appendix A suggest that the kinematic axis of the central $\sim 8^{\prime\prime}$ region – when rotation is detected – is not necessarily aligned with the photometric axis (which is typically the long axis of the GMOS FOV). The relative angle ranges from well aligned (e.g., the 3 galaxies with high core rotations) to almost maximally misaligned (e.g., NGC 1700). In addition, the GMOS data now enable us to quantify how the kinematic features in the central regions of massive galaxies are connected to and aligned with the large-scale kinematic features measured from our Mitchell data in Ene et al. (2018). A detailed analysis of these properties will be presented in a separate work.

6 Velocity dispersion

6.1 Central velocity dispersion

The high spatial sampling of the GMOS IFS enables us to measure the stellar velocity dispersion $\sigma$ within apertures of different sizes. We measure the central stellar velocity dispersion by taking a luminosity-weighted average of $\sigma$ within a circular aperture of radius 1 kpc. The resulting values are listed in column 9 of Table 1. We find that 6 of the 20 galaxies have high central velocity dispersion $\sigma_{\rm 1kpc}\gtrsim 300$ km s*-1*, 4 have $\sigma_{\rm 1kpc}\lesssim 230$ km s*-1*, and the rest is in between. NGC 890 has the lowest value with $\sigma_{\rm 1kpc}=206$ km s*-1*.

As another measure of central $\sigma$, we take a luminosity-weighted average of the GMOS $\sigma$ within an aperture of $4.1^{\prime\prime}$ diameter, which is the size of one Mitchell IFS fiber. This GMOS $\sigma$ can then be compared directly to the Mitchell $\sigma$ measured from the central Mitchell fiber. Fig. 5 shows an excellent overall agreement, with a median fractional difference of $({\rm GMOS}~{}\sigma-\sigma_{c})/\sigma_{c}=-0.01\pm 0.02$.

The observed rms scatter in the difference of the two velocity dispersion measurements is $\sim 12$ km s*-1*, while the estimated scatter computed from the individual measurement errors is $\sim 6$ km s*-1*. Since the GMOS $\sigma$ is an average of the velocity dispersion over many spatial bins, each of which is measured from a high S/N spectrum, the statistical errors on the GMOS $\sigma$ in Fig. 5 are very small ($\sim 1-2$ km s*-1*) and are likely subdominant to various systematic errors. Using the measurement errors on individual GMOS bins ($\sim 10$ km s*-1*) instead of the luminosity-averaged statistical erros, we find that the scatter in the difference of the velocity dispersions is in good agreement with the scatter estimated from the errors. Keeping this in mind and considering that the GMOS and Mitchell $\sigma$ are measured using data from different instruments and telescopes over different spectral regions using different template libraries, the agreement shown in Fig. 5 is nonetheless reassuring. A similar comparison for the higher-order moments $h_{3}$ and $h_{4}$ finds that the observed rsm scatter is $\sim 0.025$ while the scatter estimated from the measurement errors is $\sim 0.015$. Once again, accounting for the fact that the luminosity-averaged statistical errors are underestimating the observed errors ($\sim 0.02$), we recover a much better agreement between the observed and estimated scatter.

We note that measuring the GMOS $\sigma$ as a luminosity-weighted average over the dispersions from individual bins within some spatial region is not identical to measuring $\sigma$ from a single co-added spectrum within the same spatial region (e.g., a Mitchell fiber sized aperture). To quantify this difference, we use the GMOS data to generate a co-added spectrum (where we sum all the individual lenslet spectra that fall within a $4.1^{\prime\prime}$ diameter aperture) and derive the pPXF best-fitting kinematic moments for the co-added spectrum. Although the $\sigma$ values computed through these two methods (luminosity-weighted average $\sigma$ and $\sigma$ from a co-added spectrum) are not identical, we found the differences to be small, on the order of $\lesssim 5\%$. A few outliers correspond to galaxies with large velocity gradients in the region probed by GMOS, which the $4.1^{\prime\prime}$ diameter Mitchell fibers do not resolve. The unresolved high velocity leads to an increase in the recovered value of Mitchell velocity dispersion. If we instead compute $v_{rms}=\sqrt{V^{2}+\sigma^{2}}$, we find better agreement with the Mitchell value. Veale et al. (2017b) also found agreement between $\sigma$ computed from these two methods.

6.2 Comparison with literature

To compare the GMOS measurements with literature values of $\sigma$ compiled in the HyperLeda catalog (Paturel et al., 2003), we compute the GMOS $\sigma$ as a luminosity-weighted average within the aperture used by HyperLeda (0.595 kpc radius). The left panel of Fig. 6 shows good agreement between GMOS and HyperLeda values, with a median fractional difference of $({\rm GMOS}~{}\sigma-{\rm HL}~{}\sigma)/{\rm HL}~{}\sigma=0.01\pm 0.02$. The HyperLeda $\sigma$ values are compiled from heterogeneous measurements in the literature and usually have large error bars.

Deviations of more than $10\%$ occur for three galaxies: NGC 315 ($325\pm 1$ km s*-1* versus $296\pm 10$ km s*-1* in HyperLeda), NGC 1129 ($237\pm 3$ km s*-1* versus $326\pm 15$ km s*-1* in HyperLeda), and NGC 2693 ($299\pm 3$ km s*-1* versus $339\pm 10$ km s*-1* in HyperLeda). We note that the HyperLeda catalogue contains a large range of values for NGC 315 (260 km s*-1* to 360 km s*-1*) and NGC 2693 (270 km s*-1* to 400 km s*-1*).

For 18 galaxies in our sample, we can also compare the GMOS $\sigma$ with those reported in the HET catalogue (van den Bosch et al., 2015). The HET catalogue provides stellar kinematic parameters for 1022 galaxies measured from long-slit spectra in the wavelength range 4200 – 7400 Å taken with the Marcario Low Resolution Spectrograph (Hill et al., 1998) on the Hobby-Eberly Telescope (HET) at McDonald Observatory. The HET values are measured using a $3.5^{\prime\prime}\times 2^{\prime\prime}$ aperture and the typical S/N inside the aperture is greater than 100. The right panel of Fig. 6 compares the HET $\sigma$ with our GMOS $\sigma$ measured within an equivalent circular aperture with radius $1.5^{\prime\prime}$. The median fractional difference is $({\rm GMOS}~{}\sigma-{\rm HET}~{}\sigma)/{\rm HET}~{}\sigma=-0.02\pm 0.02$. We also note that for the 3 outliers with HyperLeda mentioned in the previous paragraph, our values agree better with HET: $334\pm 4$ km s*-1* for NGC 315, $231\pm 4$ km s*-1* for NGC 1129, and $331\pm 6$ km s*-1* for NGC 2693.

6.3 Aperture correction

Aperture correction relations are used frequently to transform velocity dispersions measured with fiber-fed spectrographs such as SDSS to values measured within a uniform aperture of a standard physical size. These correction relations are important for systematic studies of galaxies at different distances and particularly so in the case of high redshift studies where the typical IFS fiber size covers a significant fraction of the small apparent size of the galaxies. A typical application consists of correcting the velocity dispersion of galaxies to apertures with sizes related to the effective radius, e.g., $R_{e}$, $R_{e}/2$, or $R_{e}/8$. Our GMOS and Mitchell IFS data together span a wide radial coverage (from $\sim 0.5^{\prime\prime}$ to beyond $50^{\prime\prime}$) of each galaxy and can be used to derive an aperture correction relation between velocity dispersions measured within two commonly used radii, $R_{e}/8$ and $R_{e}$.

For most of our galaxies, the $5^{\prime\prime}\times 7^{\prime\prime}$ coverage of the GMOS FOV corresponds to a radial extent of $\sim 1-1.5$ kpc from the center of each galaxy, allowing us to measure the value of $\sigma$ within $R_{e}/8$. Here we use $R_{e}$ measured from deep $K$-band photometric data taken with WIRCam on the Canada-France-Hawaii Telescope as part of the MASSIVE Survey (Quenneville et al. in prep). To determine $R_{e}$, the photometry package ARCHANGEL (Schombert, 2007) is used to fit elliptical isophotes to the stacked image of each galaxy. A curve of growth is then constructed from the aperture luminosity for each isophote as a function of radius. The total luminosity and half-light radius $R_{e}$ are then measured from the curve of growth. The values of $R_{e}$ for our sample of 20 galaxies are given in column 5 of Table 1. They range from $16^{\prime\prime}$ to $45^{\prime\prime}$, with the average $R_{e}$ being $25^{\prime\prime}$ ($\sim 9$ kpc).

In Fig. 7 we present the aperture correction relation between the velocity dispersion at $R_{e}/8$ (measured from GMOS data) and the velocity dispersion at $R_{e}$ (measured from Mitchell data). We find the mean correction to be only 4.4% but the scatter is large:

[TABLE]

The fact that the average of the ratio of $\sigma$ at the two radii is close to unity by no means implies that the radial profiles of $\sigma$ are flat between $R_{e}/8$ and $R_{e}$. The diverse radial profiles of $\sigma$ are evident in Figs. 16-25 of Appendix A and are reflected in the the large scatter in Eqn. 2.

Four galaxies in our sample were identified in Veale et al. (2018) as having rising outer velocity dispersion profiles (green points in Fig. 7). These galaxies have higher $\sigma$ at $R_{e}$ than at $R_{e}/8$. When these four galaxies are excluded, the ratio increases to $1.062\pm 0.076$.

The fact that the ratio of the two $\sigma$ measurements is so close to unity could also be caused by luminosity-weighting, which biases $\sigma_{e}$ to higher values by assigning more weight to the inner bins. On average, the luminosity-weighted $\sigma_{e}$ is $\sim 5-10$ km s*-1*  higher than the corresponding arithmetic average. Using the arithmetic average value for $\sigma_{e}$ for all 20 galaxies leads to a ratio of $1.064\pm 0.084$.

The aperture correction reported here is in agreement with our earlier measurement of $\sigma_{e/8}/{\sigma_{e}}=1.062\pm 0.079$ using the Mitchell data for the 41 brightest MASSIVE galaxies (Veale et al., 2017b). For 40 ETGs with lower mass than our sample, the SAURON IFU data gave a higher mean aperture correction but a comparably large scatter: $1.147\pm 0.083$ (Cappellari et al., 2006). Other aperture correction studies for nearby ETGs found a ratio of 1.087 using a compilation of kinematic and photometric data from literature for 51 E and S0 galaxies (Jorgensen et al., 1995), and a ratio of 1.133 from long-slit spectroscopic data for 35 ETGs in the Coma cluster (Mehlert et al., 2003), but neither work reported error bars.

6.4 Velocity dispersion radial profiles

Fig. 8 shows the velocity dispersion $\sigma$ as a function of radius from $\sim 0.3^{\prime\prime}$ to $\sim 50^{\prime\prime}$ from our GMOS (purple circles) and Mitchell (green squares) IFS data. For more details, see the individual kinematic maps and radial profiles (unfolded along the IFU PA) of the four velocity moments $V,\sigma,h_{3}$, and $h_{4}$ in Appendix A. In the inner few arcsecs, the $\sigma$ profiles for most galaxies in our sample decrease from peak values close to the galaxy centers to smaller values further out.

In the transition region at $3^{\prime\prime}$ – $4^{\prime\prime}$, the GMOS values are well matched to the Mitchell values, with the exception of a few cases such as NGC 1129 and NGC 2693. These galaxies represent the outliers mentioned at the end of Section 6.1 with high GMOS velocity gradients that cannot be resolved by the $4.1^{\prime\prime}$ diameter Mitchell fiber(s) and for which the velocity dispersion of the inner few Mitchell fibers should be compared to $v_{rms}=\sqrt{V^{2}+\sigma^{2}}$ measured from the GMOS data. It should also be noted that the exact location of the innermost Mitchell data point in each radial plot in Fig. 8 is uncertain to within a few arcsecs due to the large diameter of the Mitchell fiber.

For most galaxies $\sigma$ continues to decrease in the region probed by Mitchell data. We find that $\sigma$ at large radii is typically smaller than near the center by $\sim 100$ – 150 km s*-1*. Four galaxies – NGC 545, NGC 741, NGC 1129 and NGC 4874 – have rising $\sigma$ profiles towards larger radii, with NGC 1129 and NGC 4874 showing the most prominent rise. The one outlier is NGC 890, whose $\sigma$ has a peak value at $R\sim 3^{\prime\prime}$ and decreases gently both towards smaller and larger radii.

We quantify $\sigma(R)$ in the central region probed by GMOS using a single power law: $\sigma(R)\propto R^{\gamma_{\rm 1kpc}}$. The best-fit logarithmic slope $\gamma_{\rm 1kpc}$ for each galaxy is given in column 10 of Table 1, and each best-fitting profile is shown as a blue straight line over the radial range used in the fit in Fig. 8. We use only data points beyond $R\approx 0.3^{\prime\prime}$ for the fit to mitigate the effects of seeing that can flatten the shape of $\sigma(R)$ in the innermost part. The $\sigma$ profiles from $0.3^{\prime\prime}$ to $\sim 4^{\prime\prime}$ are reasonably approximated by a single power law, with the exception of NGC 315 and NGC 1700, which show a break at $R\sim 2^{\prime\prime}$. For these two galaxies we restrict the fit to $R<2^{\prime\prime}$. Our aim here is to find a simple form to approximate the shape of the GMOS $\sigma(R)$ outside of $R\sim 0.3^{\prime\prime}$, and for this purpose, a single power-law provides a reasonable fit to the data. All galaxies have $\chi^{2}$ per degree of freedom (DOF) $\sim 3$ or less, with the exception of NGC 1129 which has $\chi^{2}$ per DOF $\sim 7$. We do not think it is worthwhile to attempt to find a closer fit to the diverse shapes of $\sigma(R)$ using a more complicated functional form with PSF convolution. Instead, we will perform full dynamical modeling using the complete 2D measurements of $\sigma$ as well as other kinematic moments and will report the results in later papers.

We find that 19 out of the 20 galaxies have negative $\gamma_{\rm 1kpc}$, i.e., $\sigma$ increases from $\sim 3^{\prime\prime}$ inward. To compare the central behavior of $\sigma(R)$ with that at larger radii, we compare $\gamma_{\rm 1kpc}$ with the logarithmic slope of $\sigma(R)$ at the effective radius. For the latter, we measure $\gamma_{e}$ using the wide-field Mitchell data to fit a broken power law to $\sigma(R)$ (Veale et al., 2018), and then use the asymptotic logarithmic slopes $\gamma_{1}$ and $\gamma_{2}$ to compute the local logarithmic slope at $R_{e}$. The two slopes $\gamma_{\rm 1kpc}$ and $\gamma_{e}$ are plotted in Fig. 9. In contrast to the negative $\gamma_{\rm 1kpc}$, $\gamma_{e}$ is spread out from $-0.25$ to $+0.25$. Furthermore, more than half of the galaxies have $\gamma_{\rm 1kpc}<\gamma_{e}$, i.e., $\sigma(R)$ changes logarithmic slopes in the inner $\sim 3^{\prime\prime}$ and rises more steeply towards the galaxy’s center.

Loubser et al. (2018) measured the slopes of the velocity dispersion profiles using long-slit spectroscopic observations of a combined sample of brightest cluster galaxies (BCGs) and brightest group galaxies (BGGs). The radial coverage of their long-slit data extends up to $10$ kpc for BGGs and $15$ kpc for BCGs. For the five GMOS galaxies in common with their BGG sub-sample, Loubser et al. (2018) found negative slopes. The values are shown in Table 2, along with our own measurements of $\gamma_{\rm 1kpc}$ and $\gamma_{e}$. The Loubser et al. (2018) values of the velocity dispersion slope are in good agreement with our own $\gamma_{\rm 1kpc}$ for all five galaxies and in moderate agreement with $\gamma_{e}$, except for NGC 410.

An increase in the line-of-sight $\sigma$ towards small radii can be accounted for by the presence of either a central mass concentration (e.g., a black hole) or radial anisotropy at small radii (Binney & Mamon, 1982; Gerhard et al., 1998). This mass-velocity anisotropy degeneracy can be broken somewhat by using information about the LOSVD kurtosis $h_{4}$, which we will discuss in Section  7. Our dynamical modeling of the full set of kinematic moments ($V,\sigma,h_{3},h_{4}$) will determine which combination of velocity anisotropy and mass profiles would best fit the data presented in this paper. Prior such modeling of a handful of massive elliptical galaxies have found tangential velocity anisotropy and massive black holes in the central regions and radial velocity anisotropy at larger radii (Gebhardt et al. 2000, 2003; Shen & Gebhardt 2010; Gebhardt et al. 2011; McConnell et al. 2012; Thomas et al. 2014, 2016).

7 Higher-Order Velocity Moments

7.1 Skewness $h_{3}$

In Fig. 10 of Veale et al. (2017b), we examined the behavior of the two odd velocity moments $V$ and $h_{3}$ on scales of $\sim 1$ kpc up to $\sim 20$ kpc measured from the Mitchell IFS data for a sample of MASSIVE galaxies. For the fast rotators ($\lambda_{e}\gtrsim 0.2$), we found a clear anti-correlation between the spatially-resolved $h_{3}$ and $V/\sigma$ within each galaxy. The slope of the anti-correlation, $\Delta h_{3}/\Delta(V/\sigma)$, ranged between $-0.1$ and $-0.2$. The slow rotators ($\lambda_{e}\lesssim 0.2$), on the other hand, showed positive, negative or no correlations between the two odd moments. This result is consistent with the interpretation that anti-correlations between $h_{3}$ and $V/\sigma$ within a galaxy are associated with its internal disc kinematics (Bender et al., 1994), and only fast-rotating ETGs exhibit such features.

The high-resolution GMOS IFS data now enable us to extend this analysis to the core regions of 20 ETGs in the MASSIVE survey. As shown in Fig. 10 (top panel), we find a clear anti-correlation between $h_{3}$ and $V/\sigma$ for the four galaxies (NGC 1129, NGC 1453, NGC 1700, NGC 2693) with noticeable core rotations in Fig. 3. Since the rotations continue from $\sim 1$ kpc to $1R_{e}$ and beyond (see Section 5), these four galaxies have relatively high spin parameters (measured at both 1 kpc and $1R_{e}$): $\lambda_{\rm 1kpc}=0.350,0.199,0.119,0.337$ and $\lambda_{e}=0.124,0.204,0.198$ and $0.294$, respectively.

The slope of the spatially-resolved $h_{3}$ and $V/\sigma$ measured within each galaxy is plotted against the galaxy spin parameter for all 20 MASSIVE galaxies in the bottom panel of Fig. 10. For the faster rotators with $\lambda_{\rm 1kpc}$ above 0.1, the anti-correlation slopes are between $-0.1$ and $-0.25$, similar to the range measured on larger scales in Veale et al. (2017b). The large spread in the slopes with both positive and negative values for the remaining galaxies is also similar to that seen on large scales for slow rotators in Veale et al. (2017b).

As described in Section 5, NGC 1700 has a counter-rotating core that rotates in exactly the opposite direction from the main body rotation. Despite its relatively low central rotation ($\sim 40$ km s*-1*; Fig. 3), the $h_{3}$ and $V/\sigma$ within the core region of NGC 1700 show a steep anti-correlation with a slope of $-0.25$. Hints of $h_{3}-V/\sigma$ anti-correlation for galaxies with counter-rotating cores or kinematically decoupled components have also been found in lower-mass ETGs in the ATLAS3D survey (Krajnović et al., 2011).

7.2 Kurtosis $h_{4}$

In Fig. 11 we compare the luminosity-weighted average $h_{4}$ measured within $1$ kpc using the GMOS data to the luminosity-weighted average $h_{4}$ measured within one effective radius using the Mitchell data. On both small and large scales, we find that the average $h_{4}$ is predominantly positive: 14/20 galaxies have $h_{4,1kpc}>0$ and 18/20 galaxies have $h_{4,e}>0$. We have previously found that MASSIVE galaxies show predominantly positive large scale average $h_{4}$ that is in $\sim 50\%$ of cases accompanied by a rising outer velocity dispersion profile (Veale et al. 2017b, 2018).

As mentioned in Section 6, $h_{4}$ can help break the degeneracy between the effects of mass and velocity anisotropy on the behavior of the line-of-sight velocity dispersion. For example, Veale et al. (2018) found a statistically significant positive correlation between the radial gradient in $h_{4}$ and the velocity dispersion outer logarithmic slope $\gamma_{\rm outer}$ for a sample of 90 MASSIVE galaxies. They argue that this trend is more likely due to mass profile variations than to velocity anisotropy, since radial anisotropy at large radii would imply that a positive $h_{4}$ gradient comes along with a more negative $\sigma$ gradient. Their arguments applied to the behavior of $\sigma$ and $h_{4}$ (and the implications for mass and velocity anisotropy) at large radii, where the dominant contribution comes from the dark matter halo mass.

By contrast, the GMOS data studied here probe the central regions of galaxies in which stars and possibly black holes dominate the mass. In Fig. 12 we plot the $h_{4}$ radial gradient versus the velocity dispersion logarithmic slope $\gamma_{\rm 1kpc}$ for the core region of the galaxies in our sample. We do not find any significant correlation between the $h_{4}$ radial gradient and the velocity dispersion slope for the sample as a whole. Even after excluding the outlier NGC 890, the $p$-value of the significance of correlation for a linear fit to the data (with slope $0.34\pm 0.19$) is $0.17$. However, we note that 14 out of 19 galaxies with negative $\gamma_{\rm 1kpc}$ have positive $h_{4}$ gradients (and positive $h_{4}$ values in general, as shown in Fig. 11).

For galaxies without central black holes, the information about $h_{4}$ can be used to infer the velocity dispersion anisotropy: at small radii close to the galaxy center, various types of models predict that tangential anisotropy produces a peaked LOSVD ($h_{4}>0$), while radial anisotropy gives a flat-top LOSVD ($h_{4}<0$) (Bender et al. 1994 and references therein). However, using information about $h_{4}$ alone to estimate the effect of a black hole on the LOSVD and/or the anisotropy profile is not as straightforward. For example, Baes et al. (2005) studied a two-parameter family of isotropic models with a central black hole. They found that the $h_{4}$ profile is significantly affected by the presence of a central black hole for the isotropic models with shallow central density profiles (i.e. the central regions show positive $h_{4}$ peaks), but didn’t see any effects for models with steep density cusps. Additionally, there are several cases of SMBHs with measured negative $h_{4}$ where Schwarzschild orbit modeling predicts tangential anisotropy in the galaxy center (Pinkney et al. 2003; Gebhardt et al. 2003; McConnell et al. 2012). This behavior was also reproduced by merger simulations of dynamical systems with central SMBHs that have produced cases with central tangential anisotropy and negative $h_{4}$, particularly in the systems with higher SMBH masses (Rantala et al., 2018).

Comparing to the isotropic models of Baes et al. (2005), the fact that our galaxies do not exhibit central positive $h_{4}$ peaks suggests that, if the data resolve the black hole SOI, then the galaxies are unlikely to be isotropic. The simulations of Rantala et al. (2018) further suggest that central tangential anisotropy could explain the low observed $h_{4}$, but detailed dynamical modeling is needed to make any robust statements about the behavior of the velocity anisotropy in the center of galaxies.

8 Jeans modeling

The stellar kinematics over the large radial coverage of the combined GMOS and Mitchell datasets provide powerful constraints on the mass components in ETGs. Here we apply the computationally simple method of Jeans Anisotropic Modeling (JAM; Cappellari 2008) to illustrate how the two datasets can be used to infer the central black hole mass, the stellar mass-to-light ratio and the dark matter halo content of the regular fast rotator NGC 1453 in our sample.

We choose NGC 1453 due to its regular photometric properties and almost purely elliptical isophotes (Fig. 13; Goullaud et al. 2018). In addition to being a regular fast-rotator, Ene et al. (2018) showed that the photometric and kinematic position angles of NGC 1453 are almost perfectly aligned, suggesting that the galaxy may be modeled as an axisymmetric system. As an interesting note, NGC 1453 also contains a large ($R\sim 20^{\prime\prime}$) warm ionized gas disc that is rotating almost perpendicularly to the stellar kinematic axis, suggesting that the gas accreted through external processes (Pandya et al., 2017).

The rest of our sample is either slowly rotating and/or shows kinematic and photometric twists. The orbit superposition model (Schwarzschild, 1979) is likely a more suitable method for determining the dynamical masses of these galaxies. Orbit modeling results are beyond the scope of the present paper and will be reported in future work.

8.1 Mass Model

We assume that the mass in each galaxy consists of three parts: a central black hole, a stellar mass component and a dark matter halo. To model the stellar component we use the Hubble Space Telescope Wide Field Camera 3 (WFC3) IR photometry in the filter F110W from Goullaud et al. (2018). The surface brightness is fitted using the Multi-Gaussian Expansion (MGE) method (Emsellem et al. 1994; Cappellari 2002) with a sum of 2-dimensional Gaussian components that share a common center and position angle:

[TABLE]

where $x^{\prime}$ and $y^{\prime}$ are projected coordinates measured from the galaxy center, with $x^{\prime}$ being along the photometric major axis. The subscript $k$ labels the individual Gaussian components, $L_{k}$, $\sigma_{k}$, $q_{k}^{\prime}$ are the luminosity, dispersion and projected axis ratio of each Gaussian respectively. To obtain the MGE fit, the model predictions need to be convolved with the point spread function (PSF), which is taken from Goullaud et al. (2018), and is also expressed as an MGE using 5 nearly-circular Gaussian components (with axis ratios $>0.98$). Our best-fit MGE to the surface brightness of NGC 1453 consists of 11 Gaussian components, which are summarized in Table 3 and plotted in Fig. 13. The small fitting residuals (bottom panel) demonstrate that the MGE model agrees very well with the data.

Given an inclination and stellar mass-to-light ratio $\Upsilon_{\rm F110W}$, the 2-dimensional MGE fit can be deprojected and converted into a realistic 3-dimensional stellar mass density profile. Throughout this work we assume an edge-on orientation (inclination of $90^{\circ}$) and a spatially constant $\Upsilon_{\rm F110W}$.

We parametrize the dark matter halo using a two-parameter logarithmic potential:

[TABLE]

where $R_{c}$ is the characteristic radius of the halo potential and $V_{c}$ is the circular speed at $r\rightarrow\infty$ (such a DM description was, e.g., favored by orbit-based models of M87 in Murphy et al. 2011). The logarithmic potential corresponds to a density profile given by (Binney & Tremaine 2008):

[TABLE]

Before being passed to JAM, the dark matter density profile is fitted using a 1-dimensional MGE model.

8.2 Modeling velocity anisotropy

We augment the four parameters describing the mass profile with a description for the anisotropy,

[TABLE]

which quantifies the flattening of the velocity ellipsoid along the minor axis.

Previous works with JAM have typically assumed a globally constant $\beta_{z}$, although there have been some attempts at introducing moderate spatial variation, by having $\beta_{z}$ vary between the different Gaussian components of the MGE (e.g., Cappellari et al. 2015). The large radial span of our combined two data sets motivates including at least some spatial dependence. Perhaps most importantly, our model should distinguish between the core and the outskirts of the galaxy. The core-like part of the light profile of massive galaxies is commonly interpreted as being the result of black hole scouring (Faber et al. 1997), which in turn predicts a bias towards tangential orbits in the centre and more radial orbits in the outskirts (Quinlan & Hernquist 1997; Milosavljević & Merritt 2001). In order to at least partially replicate orbit-type variation, we assign separate anisotropies $\beta_{z}^{G}$ to the Gaussians with $\sigma_{k}<1^{\prime\prime}$ and the Gaussians with $\sigma_{k}>1^{\prime\prime}$. The choice of $1^{\prime\prime}$ is motivated by the fact that the light profile of NGC 1453 starts to fall off more steeply once $R\gtrsim 1^{\prime\prime}$ (see bottom panel of Fig. 13).

8.3 Fitting the kinematic data

Our mass model is specified by six parameters: $M_{\rm BH}$, $\Upsilon_{\rm F110W}$, $R_{c}$, $V_{c}$, $\beta_{z}^{G}(\sigma_{k}<1^{\prime\prime})$ and $\beta_{z}^{G}(\sigma_{k}>1^{\prime\prime})$. JAM places constraints on these parameters by comparing the model predictions for $v_{\rm rms}=\sqrt{V^{2}+\sigma^{2}}$ (including instrumental PSF convolution) to the observed $v_{\rm rms}$ for each Voronoi bin. We use the GMOS PSF (modeled as a Gaussian with dispersion $0.297^{\prime\prime}$) since the PSF is most important near the center, at the locations of the GMOS data points. We compute the PSF by fitting a Gaussian function to point sources from the GMOS acquisition images, which are usually taken before a group of 3 – 4 science exposures. The final PSF estimate is a weighted average of the Gaussian FWHMs, weighed by the total exposure time of the group of science exposures that each acquisition image precedes. In the fit, we exclude the few Mitchell kinematic points with $R<5^{\prime\prime}$ since this region is better measured by GMOS.

We find the JAM parameters by first running a broad regular-grid search, followed by a best-fit estimation using Bayesian inference. We calculate the posterior probability distribution using an implementation (Cappellari 2013) of the Adaptive Metropolis algorithm of Haario et al. (2001). The best-fitting model (with 3$\sigma$ errors) is $\Upsilon_{\rm F110W}=2.06\pm 0.13\ \Upsilon_{\rm\odot,F110W}$, $M_{\rm BH}=3.29\pm 0.75\times 10^{9}\ M_{\odot}$, $V_{c}=364\pm 134\ {\rm km/s}$, $R_{c}=7.2\pm 8.3\ {\rm kpc}$, $\beta_{z}^{G}(\sigma_{k}<1^{\prime\prime})=-0.58\pm 0.62$ and $\beta_{z}^{G}(\sigma_{k}>1^{\prime\prime})=0.15\pm 0.04$.

The $\chi^{2}$ of a large suite of JAM models shows a clear minimum for each of the six mass model parameters (Fig. 14). The constraints on the two dark matter halo parameters are weaker primarily because at $r\ll R_{c}$, $\rho_{\rm DM}\propto V_{c}^{2}/R_{c}^{2}$. As a result, there is significant degeneracy between the two halo parameters. The negative $\beta_{z}^{G}(\sigma_{k}<1^{\prime\prime})$ and positive $\beta_{z}^{G}(\sigma_{k}>1^{\prime\prime})$ of our best-fit model indicate a radially increasing $\beta_{z}$ profile (the $\beta_{z}$ profile that can be obtained from the best-fit $\beta_{z}^{G}$ is not a step function, but smoothly increasing from $\beta_{z}\approx-0.25$ at small radii to $\beta_{z}\approx 0.15$ outside of the core) 111We also ran JAM with the two-component $\beta_{z}$ replaced by a single, constant $\beta_{z}$. The resulting best-fit model ($\Upsilon_{\rm F110W}=2.25\ \Upsilon_{\rm\odot,F110W}$, $M_{\rm BH}=2.27\times 10^{9}\ M_{\odot}$, $V_{c}=404\ {\rm km/s}$, $R_{c}=11.6\ {\rm kpc}$, $\beta_{z}=0.11$) produces a worse fit, with $\chi^{2}=292.7$.. JAM’s $\beta_{z}$ cannot be directly related to tangential/radial anisotropy (defined in spherical coordinates), but our data clearly favors a velocity anisotropy that is different in the inner and outer parts of the galaxy. Figure 15 shows our best-fitting $v_{rms}$ predictions, demonstrating that the trends in both the GMOS ($R<5^{\prime\prime}$) and Mitchell ($R>5^{\prime\prime}$) data sets can be well reproduced by JAM and a model without a central black hole fails to reproduce the observed kinematics at small radii.

While JAM is plausibly a suitable method for NGC 1453 and other fast rotators, the utility of JAM is limited by its assumptions of axisymmetry, cylindrically aligned velocity ellipsoids and by the issue that its solutions could be unphysical. These assumptions are not well motivated in slowly rotating, triaxial ETGs, which constitute the majority of our galaxies. Furthermore, the MGE-based anisotropy profiles, which are globally constant or with limited spatial variation, may be too simplistic for representing a galaxy’s true anisotropy structure. Future papers in the MASSIVE survey series will present results from full-scale dynamical mass modeling using stellar orbit libraries that take advantage of the full set of kinematic moments from GMOS and Mitchell observations (Liepold et al. in preparation; Quenneville et al. in preparation).

9 Conclusions

We have presented the first results from the high spatial-resolution spectroscopic component of the MASSIVE survey. The spatially-resolved spectroscopic observations were obtained with the Gemini GMOS-N IFS for the central $5^{\prime\prime}\times 7^{\prime\prime}$ (about 2 kpc across) regions of 20 ETGs in the MASSIVE survey. These galaxies have $M_{K}\leq-25.5$ ($M_{*}\gtrsim 10^{11.7}M_{\odot}$) and are located at a median distance of $\sim 70$ Mpc. We measured the stellar kinematics from high-$S/N$ ($\sim 120$) spectra using a Gauss-Hermite parameterization of the LOSVD. We obtained two-dimensional maps of the first four kinematic moments ($V$, $\sigma$, $h_{3}$, and $h_{4}$) and compared them to the large-scale stellar kinematics obtained from the Mitchell IFS (Veale et al. 2017b, a, 2018). The two IFS datasets together cover a length scale of $\sim 0.3^{\prime\prime}-100^{\prime\prime}$, or $\sim 0.1-30$ kpc. Our main findings are as follows:

• •

The velocity maps of most galaxies in the sample show some level of rotation in the central regions (upper left panel in Figs. 16-25 of Appendix A). Three galaxies (NGC 1129, NGC 1453, and NGC 2693) exhibit fast rotations with $|V|$ rising up to $\sim 100-150$ km s*-1* within 1 kpc, and a fourth galaxy (NGC 1700) has a counter-rotationg core that rotates in the opposite direction of the main body rotation (Fig. 3 and Fig. 4). These 4 galaxies also have noticeable rotation out to one effective radius, $\lambda_{e}\gtrsim 0.1$. The rest of the galaxies have central rotations (if detected) with $|V|\lesssim 30$ km s*-1*, as well as low velocity rotation up to $R_{e}$ ($\lambda_{\rm 1kpc},\lambda_{e}\lesssim 0.1$). The kinematic axis of the central rotation is not necessarily aligned with the photometric axis or the large-scale kinematic axis, indicating a diverse merger history within this sample of 20 high-mass ETGs.

• •

The velocity dispersion $\sigma$ within $1$ kpc reaches $\sim 300$ km s*-1* or beyond in 6 of the 20 galaxies and is greater than 250 km s*-1* for 14 galaxies. We measure the luminosity-weighted average $\sigma$ within an aperture of radius $R_{e}/8$ versus $R_{e}$ and find an aperture correction relation of $\sigma_{e/8}/\sigma_{e}=1.044\pm 0.078$ (Fig. 7).

• •

For all but one galaxy, the velocity dispersion profiles $\sigma(R)$ in the radial range of $\sim 0.3^{\prime\prime}-3^{\prime\prime}$ are well fit by a single power law form with negative logarithmic slope (Fig. 8). A rising $\sigma(R)$ towards smaller radii is indicative of the presence of a central SMBH but can also be caused by central radial anisotropy in the stellar velocities. On large scales, by contrast, the logarithmic slope of $\sigma(R)$ at $1R_{e}$ ranges from $-0.25$ to $+0.25$, reflecting varying degrees of contributions from dark matter mass as well as velocity anisotropies. We will use the observed higher-order velocity moment $h_{4}$ to break this mass-anisotropy degeneracy in future dynamical mass modeling work.

• •

For the galaxies with clear rotation in the inner $\sim 1$ kpc, the spatially-resolved skewness $h_{3}$ in this region is anti-correlated with the velocity, i.e., the spatial bins with higher $V$ values have more negative $h_{3}$. The slope of this correlation, $\Delta h_{3}/\Delta(V/\sigma)$, ranges between $-0.1$ and $-0.25$ for the galaxies with relatively high spins in our sample ($\lambda_{\rm 1kpc}\gtrsim 0.1$). The $h_{3}$-$V$ anti-correlation within a galaxy indicates disc-like kinematics in both the core and main body of fast-rotating massive ETGs.

• •

The kurtosis $h_{4}$ is generally positive, with 14 of the 20 galaxies having a positive average $h_{4}$ within $1$ kpc, and 18 of them having positive average $h_{4}$ within $1R_{e}$ (Fig. 11). Most galaxies also show a rising or flat radial profile $h_{4}(R)$. In the absence of a central mass, peaked LOSVDs ($h_{4}>0$) at small radii are often associated with tangential anisotropy, which tends to lower the central line-of-sight $\sigma$. Yet 19 out of the 20 galaxies in this study are observed to have $\sigma(R)$ that increase towards the center. The presence of central black holes is a possible explanation. We will perform full-scale stellar orbit modeling to determine the velocity structures, stellar and dark matter mass distributions, and black hole masses in these galaxies.

• •

Using the GMOS and Mitchell stellar kinematics, we applied the Jeans modeling method to measure the mass distributions in NGC 1453, the most regular fast rotator in our sample. To partially account for spatial variations in the stellar velocity anisotropy, we allow the anisotropy parameter $\beta_{z}$ in the inner and outer parts of the galaxy to be different. The JAM results show that our kinematics point towards a nonzero central black hole mass and a spatially varying velocity anisotropy. The best-fit mass model (with 3$\sigma$ errors) has a black hole mass of $3.29\pm 0.75\times 10^{9}M_{\odot}$, a stellar mass-to-light ratio of $\Upsilon=2.06\pm 0.13\ \Upsilon_{\odot}$ (in WFC3 F110W band), and a circular velocity of $V_{c}=364\pm 134$ km s*-1* for the dark matter halo.

High spatial-resolution kinematics that resolve the sphere of influence of the SMBH are a key requirement for any attempt at black hole mass determination through dynamical modeling. Combining the small scale results presented here with the large-scale kinematics of previous MASSIVE papers will allow us to model the mass distributions of nearby massive ETGs and study their assembly histories. New SMBH mass measurements will help refine the various correlation relations between black holes and their host galaxy properties. In particular, the high-mass range sampled by our galaxy sample may be relevant for the exploration of the $M_{\bullet}-\sigma_{e}$ saturation at high $\sigma_{e}$ (e.g., Lauer et al. 2007; McConnell et al. 2011, 2012; McConnell & Ma 2013; Kormendy & Ho 2013).

This work is based on observations obtained at the Gemini Observatory processed using the Gemini IRAF package, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), Ministério da Ciência, Tecnologia e Inovação (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). We thank Richard McDermid for his help in the early stages of acquiring and preparing the GMOS data. The MASSIVE survey is supported in part by NSF AST-1411945, AST-1411642, AST-1815417 and AST-1817100, HST GO-14210, GO-15265 and AR-14573, and the Heising-Simons Foundation. J.L.W. is supported in part by NSF grant AST-1814799. Gemini North is located on Maunakea, a sacred place for indigenous Hawaiians who have honored it before and since the construction of modern astronomical facilities. We hold great privilege and responsibility in using the Maunakea summit.

Gemini:Gillett (GMOS)

Appendix A Kinematic moments plots