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. The ETGs are complex, multicomponent 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 (IFSs), e.g., SAURON (Emsellem et al. 2004), ATLAS 3D (Cappellari et al. 2011), SAMI (Croom et al. 2012), CALIFA (Sánchez et al. 2012), andMaNGA (Bundy et al. 2015). These surveys use wide-field IFSs 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, ETG morphologies, and molecular and ionized gas content. The ETGs targeted by these surveys are predominantly fast-rotating S0 or elliptical galaxies with MM 10 11.5 * < . The spatial sampling of these surveys is limited by the IFS fiber diameter (1 6, 2″, and 2 7 for SAMI, MaNGA, and CALIFA, respectively) or lenslet size (0 94 for SAURON/ATLAS 3D ).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 multislits to reach a sky coverage of ∼2-4 effective radii for a subset of ATLAS 3D galaxies (Brodie et al. 2014).

We initiated the volume-limited MASSIVE survey (Ma et al. 2014) to investigate the ∼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 * ∼10 11.5 M e , a parameter space little explored previously. Wide-field kinematics and stellar population studies from this survey have been published in Greene et al. ) 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10″) kinematically distinct components (KDCs).

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 into galaxies' global dynamics and assembly histories, the kinematics of the innermost regions of galaxies are critical for measuring the masses of the SMBHs and elucidating the symbiotic relations among black holes, baryons, and dark matter near galactic centers. To achieve these goals, we observed the central 5″×7″ of 20 MASSIVE galaxies using the Gemini Multi Object Spectrograph (GMOS; Hook et al. 2004) with 0 2 lenslets on the 8.1 m Gemini North Telescope. The exposure times were chosen to yield stellar spectra with a signal-to-noise ratio (S/N) of ∼120 for each spatial bin. These high-quality spectra allow us to obtain detailed two-dimensional maps of the velocity and velocity dispersion, as well as the skewness and kurtosis of the stellar line-of-sight velocity distribution (LOSVD). Depending on the galaxy, the maps contain from 50 to 300 spatial bins (with an average of 130 bins) and cover a scale from ∼0.1 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″×10″ region of 28 ETGs in the SAURON survey using the OASIS spectrograph with a spatial sampling of 0 27. The fine spatial sampling allowed them to identify two types of KDCs in lower-mass ETGs: old, kpc-scale KDCs that exist in slow-rotating ETGs, and young, small (∼100 pc scale), almost counterrotating 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 e , 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-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 used to constrain the dynamical models of the fast rotator NGC1453. Section 9 summarizes our main conclusions.

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 ETGs with M K -25.3 mag (from 2MASS; Skrutskie et al. 2006) corresponding to M * 10 11.5 M e that are within a distance of 108 Mpc in the northern sky (decl. δ > -6°).

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.T h e2 0g a l a x i e sa r e located at a distance in the range of 54.4-102.0 Mpc (with a median distance of ∼70 Mpc) and have -25.50magM K  -26.33mag, which corresponds to stellar mass 10 11.7 M e  M * 10 12 M e . Since they were selected based on the ability to obtain high-S/NG M O Sd a t af o rd y n a m i c a lm a s s modeling, this sample of 20 galaxies tends to be the closer and more massive part of the general MASSIVE sample: they represent ∼50% of galaxies within 80 Mpc and ∼60% of galaxies more massive than M K < -25.8mag.

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 and 2016 under the programs

, and GN-2016B-Q-18. Total exposure times were chosen to ensure an S/N of ∼120 (after spatial binning) and vary between 1 and 6 hr (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 an FOV of 5″×7″ consisting of 1000 hexagonal lenslets with a projected diameter of 0 2. An additional 500 lenslets observe a 5″×3 5 region of the sky displaced by ∼1′ 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 Ca II triplet and Na I absorption features used for measuring stellar kinematics and stellar populations, respectively.

The detector array consists of three 2k×4k e2v deep depletion CCDs placed in a row with ∼37 unbinned pixel gaps in between. To mitigate the loss of spectral information to the chip gaps, we use spectral dithering with a grating central wavelength λ c for half of the exposures and λ c +50 Åfor the other half. For most galaxies, a typical value of λ c is between 8600 and 8700 Å.W e carefully choose this value to ensure that the Ca II 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 IRAF software. For an indepth 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 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 flatfielding 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 onedimensional spectrum corresponding to each GMOS lenslet. For all galaxies observed in semester 2013B or later, arclamp 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 arccalibrated 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 ò λ in each GMOS lenslet by fitting ∼10 bright sky lines and a polynomial function ò λ (λ) across the observed wavelength range. To mitigate the low S/N of the sky lines in individual lenslets, we fit ò λ (λ) as a first-order polynomial in each lenslet independently. The resulting two-dimensional map of ò λ in each science frame is interpolated to a wavelength of interest (i.e., centered on the Ca II 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 Δλ. 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 (ò λ , Δλ) 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, 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 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 one-dimensional spectra) are 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 ∼10%, as can be seen in Figure 1.

In most cases, a Voronoi bin is composed of multiple 0 2×0 2 grid segments, each affiliated with multiple onedimensional spectra, usually from different exposures. The final step is to co-add all of the one-dimensional spectra in each Voronoi bin and create a corresponding Gaussian kernel for the instrumental resolution of the co-added spectrum. This is done with a clipped 3σ mean and rescaling for regions overlapping the chip gaps at one of the two wavelength settings.

Measuring Stellar Kinematics

The stellar LOSVD is extracted using the penalized pixelfitting (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),

where y=(v -V )/σ, V is the mean velocity, σ is the velocity dispersion, and H m is the mth 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 zero for V and h 3 through h 6 and 300 km s -1 for σ. 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 ∼10% or less of the fit region. We present example pPXF in a central, intermediate, and outer spatial bin for NGC1700 in Figure 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 same 15 K and G stars as in Table 2 of Barth et al. (2002). The second template set contains all ∼360 G and K stars in the CaT Library. For the latter, we ran pPXF for all of 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 ∼5 km s -1 for V and σ and ∼0.01 for h 3 and h 4 , well within the measurement errors of ∼10 km s -1 and ∼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 FWHM of the lines versus wavelength and determine the corresponding FWHM for our fit region centered at ∼8600 Å. We determine this best-fitting FWHM for each individual lenslet of each exposure. Typical values for the FWHM are ∼2.5 Å, with variations of ∼0.3 Åwith lenslet position. This corresponds to an instrumental resolution of ∼37 km s -1 at 8600 Åwith a sampling of 0.67 Åpixel -1 . 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 (LSF) 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 nonnegligible bias in the recovery of the kinematic moments only when the measured velocity dispersion is smaller than the instrumental dispersion (Cappellari 2017). Typical velocity dispersions for MASSIVE galaxies are 250 km s -1 , much higher than the spectral resolution of GMOS (∼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 circular velocity of V c =364±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 • -σ e saturation at high σ e (e.g., Lauer et al. 2007;McConnell et al. 2011McConnell et al. , 2012;;McConnell & Ma 2013;Kormendy & Ho 2013).