The MASSIVE Survey VI: The Spatial Distribution and Kinematics of Warm Ionized Gas in the Most Massive Local Early-type Galaxies

TL;DR

This study investigates the presence, distribution, and movement of warm ionized gas in 74 massive early-type galaxies, revealing diverse morphologies, kinematic behaviors, and potential origins of the gas.

Contribution

First systematic analysis of warm ionized gas in the most massive local early-type galaxies using integral field spectroscopy, highlighting diversity in gas properties and origins.

Findings

38% of galaxies have detectable warm ionized gas.

Fast rotators are more likely to contain gas than slow rotators.

Gas morphologies and kinematics vary widely among the sample.

Abstract

We present the first systematic investigation of the existence, spatial distribution, and kinematics of warm ionized gas as traced by the [O II] 3727AA emission line in 74 of the most massive galaxies in the local Universe. All of our galaxies have deep integral field spectroscopy from the volume- and magnitude-limited MASSIVE survey of early-type galaxies with stellar mass log(M_*/M_sun) > 11.5 (M_K < -25.3 mag) and distance D < 108 Mpc. Of the 74 galaxies in our sample, we detect warm ionized gas in 28, which yields a global detection fraction of 38\pm6% down to a typical [O II] equivalent width limit of 2AA. MASSIVE fast rotators are more likely to have gas than MASSIVE slow rotators with detection fractions of 80\pm10% and 28\pm6%, respectively. The spatial extents span a wide range of radii (0.6 - 18.2 kpc; 0.1 - 4R_e), and the gas morphologies are diverse, with 17/28 = 61\pm9% being centrally concentrated, 8/28 = 29\pm9% exhibiting clear rotation out to several kpc, and 3/28 = 11\pm6% being extended but patchy. Three out of four fast rotators show kinematic alignment between the stars and gas, whereas the two slow rotators with robust kinematic measurements available exhibit kinematic misalignment. Our inferred warm ionized gas masses are roughly ~10^5M_sun. The emission line ratios and radial equivalent width profiles are generally consistent with excitation of the gas by the old underlying stellar population. We explore different gas origin scenarios for MASSIVE galaxies and find that a variety of physical processes are likely at play, including internal gas recycling, cooling out of the hot gaseous halo, and gas acquired via mergers.

Tables

Table 1. Table 1 Sample of MASSIVE Galaxies with Warm Ionized Gas

Galaxy	$D$	$R_e$	$R_{\mathrm{gas}}$	$M_K$	$M_{\ast }$	$M_h$	$M_{\mathrm{gas}}$	${ℒ}_{H\beta }$	${\sigma }_{\mathrm{stars},\mathrm{central}}$	[O III]/H$\beta$	${\mathrm{PA}}_{\mathrm{gas}}$	${\mathrm{PA}}_{\mathrm{stars}}$	${\mathrm{PA}}_{\mathrm{gas}-\mathrm{stars}}$	FR/SR	Extended	Rotating
$-$	[Mpc]	[kpc]	[kpc]	[mag]	[$\log M_{\odot }$]	[$\log M_{\odot }$]	[$\log M_{\odot }$]	[$\log$ erg s-1]	[km s-1]	$-$	[deg]	[deg]	[deg]	$-$	$-$	$-$
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)
NGC0708	69.0	16.32	4.09	-25.65	11.75	14.5	5.52	39.1	248.03	2.14	356.8$\pm$1.6	126.0$\pm$17.5	50.8$\pm$17.6	SR	Yes	Yes
NGC1167	70.2	9.97	18.21	-25.64	11.74	13.0	5.73	39.3	231.81	5.71	284.3$\pm$15.3	254.0$\pm$3.8	30.3$\pm$15.8	FR	Yes	Yes
NGC1453	56.4	7.63	8.41	-25.67	11.75	13.6	5.82	39.4	341.04	2.59	312.5$\pm$3.1	35.0$\pm$3.2	97.5$\pm$4.4	FR	Yes	Yes
NGC2208	84.1	10.91	8.86	-25.63	11.74	10.0	5.55	39.1	280.34	2.25	241.8$\pm$13.4	95.0$\pm$9.5†	$-$	SR	Yes	Yes
NGC2320	89.4	8.16	5.65	-25.93	11.87	14.1	5.57	39.1	345.86	2.66	148.2$\pm$16.5	147.0$\pm$4.0	1.2$\pm$17.0	FR	Yes	Yes
NGC2783	101.4	18.24	2.85	-25.72	11.78	12.3	5.54	39.1	274.76	$>$1.18	148.5$\pm$5.0	269.5$\pm$11.0	59.0$\pm$12.1	SR	Yes	Yes
NGC7426	80.0	7.26	5.74	-25.74	11.79	13.7	5.71	39.3	306.75	2.66	261.7$\pm$3.1	251.0$\pm$1.2	10.7$\pm$3.3	FR	Yes	Yes
NGC7550	72.7	9.64	2.35	-25.43	11.65	11.8	5.32	38.9	270.73	3.31	238.7$\pm$10.4	155.2$\pm$15.5†	$-$	SR	Yes	Yes
NGC1700	54.4	6.06	12.27	-25.60	11.72	12.6	4.82	38.4	262.74	$>$0.63	$-$	268.0$\pm$1.8	$-$	FR	Yes	No
NGC6482	61.4	4.44	16.46	-25.60	11.72	13.0	5.61	39.2	330.44	1.26	$-$	63.5$\pm$2.0	$-$	SR	Yes	No
NGC7265	82.8	10.93	7.41	-25.93	11.87	14.6	5.08	38.6	256.21	$>$1.20	$-$	259.5$\pm$9.5	$-$	SR	Yes	No
NGC0315	70.3	8.42	1.68	-26.30	12.03	13.4	4.74	38.3	353.32	$>$1.06	$-$	222.0$\pm$7.2	$-$	SR	No	No
NGC0383	71.3	7.09	1.73	-25.81	11.82	14.3	4.60	38.1	304.02	$>$1.56	$-$	140.5$\pm$2.8	$-$	FR	No	No
NGC0410	71.3	11.30	1.79	-25.90	11.86	14.3	4.66	38.2	313.07	$>$0.97	$-$	161.0$\pm$18.5	$-$	SR	No	No
NGC0533	77.9	15.29	1.88	-26.05	11.92	13.8	4.96	38.5	304.81	1.44	$-$	51.2$\pm$5.1†	$-$	SR	No	No
NGC0547	74.0	7.13	1.81	-25.83	11.83	$-$	4.58	38.1	277.65	$>$2.14	$-$	176.5$\pm$21.8	$-$	SR	No	No
NGC1060	67.4	12.98	1.76	-26.00	11.90	14.0	4.19	37.7	325.91	0.88	$-$	74.0$\pm$7.4†	$-$	SR	No	No
NGC3805	99.4	7.34	2.22	-25.69	11.76	14.7	4.58	38.1	280.80	1.25	$-$	61.5$\pm$1.5	$-$	FR	No	No
NGC3842	99.4	10.20	2.11	-25.91	11.86	14.7	5.38	38.9	272.95	1.50	$-$	1.6$\pm$0.2†	$-$	SR	No	No
NGC5129	107.5	10.11	2.32	-25.92	11.86	10.0	4.66	38.2	277.81	$>$1.21	$-$	3.0$\pm$2.5	$-$	FR	No	No
NGC5208	105.0	5.41	1.48	-25.61	11.73	13.0	4.27	37.8	274.84	2.09	$-$	343.0$\pm$1.5	$-$	FR	No	No
NGC5322	34.2	2.42	0.60	-25.51	11.68	13.5	4.08	37.6	258.94	$>$1.13	$-$	267.5$\pm$7.8	$-$	SR	No	No
NGC6223	86.7	7.85	2.05	-25.59	11.72	13.5	5.14	38.7	294.97	$>$0.72	$-$	288.0$\pm$2.2	$-$	FR	No	No
NGC6375	95.8	9.08	2.26	-25.53	11.69	10.0	4.76	38.3	251.26	$>$2.99	$-$	322.0$\pm$7.0	$-$	FR	No	No
NGC7052	69.3	8.98	1.65	-25.67	11.75	10.0	3.63	37.2	314.09	$>$1.26	$-$	61.5$\pm$2.2	$-$	SR	No	No
NGC7436	106.6	12.39	2.48	-26.16	11.97	14.2	4.59	38.1	311.13	$>$1.07	$-$	197.5$\pm$10.2	$-$	SR	No	No
NGC7626	54.0	6.21	1.16	-25.65	11.75	13.9	3.52	37.1	305.02	$>$0.94	$-$	10.5$\pm$1.1†	$-$	SR	No	No
UGC02261	70.6	5.32	1.68	-25.61	11.73	11.8	4.09	37.6	308.33	$>$0.65	$-$	66.0$\pm$2.5	$-$	FR	No	No

Equations

$${\mathcal{M}}_{\mathrm{H}\beta }=28.2\times 10^8{\mathcal{L}}_{\mathrm{H}\beta ,43}{n}_{\mathrm{e},100}^{-1}{M}_{\odot }.$$

Full text

I. Introduction

Massive early-type galaxies typically have very old stellar populations. Given the lack of ongoing star formation, we expect a corresponding lack of an interstellar medium (ISM). However, a significant fraction of the massive early-type population has been observed to contain a complex and multi-phase ISM (e.g., Knapp et al., 1985; Phillips et al., 1986; Knapp et al., 1996; Goudfrooij, 1999; Sarzi et al., 2006, 2010; Davis et al., 2011, 2016). Roughly $10^{7}-10^{9}M_{\odot}$ worth of molecular gas is detected in $\sim 25\%$ of typical early-type galaxies (Knapp & Rupen, 1996; Combes et al., 2007; Young et al., 2011, and references therein). The warm ionized gas traces the dust in these galaxies quite well in terms of both mass fraction and morphology (roughly $10^{3}-10^{5}M_{\odot}$; e.g., Kim, 1989; Buson et al., 1993; Goudfrooij et al., 1994), and is often accompanied by the detection of molecular gas (e.g., Young et al., 2011; Davis et al., 2011; Alatalo et al., 2013). The relationship between the warm ionized and hot X-ray-emitting gas is far less clear (e.g., Sparks et al., 1989; Macchetto et al., 1996; Goudfrooij, 1999), and may well depend on the temperature profile of the X-ray-emitting gas (e.g., McDonald et al., 2010; Werner et al., 2014).

Despite more than thirty years of studying warm ionized gas in early-type galaxies (e.g., Caldwell, 1984; Demoulin-Ulrich et al., 1984; Phillips et al., 1986; Buson et al., 1993; Goudfrooij et al., 1994; Zeilinger et al., 1996; Macchetto et al., 1996; Caon et al., 2000; Martel et al., 2004; Sarzi et al., 2006, 2010), we still do not have a clear understanding of the different origin and excitation mechanisms of the ionized gas. We do not know if the main sources are internal (e.g., stellar mass loss) or external (e.g., mergers). The total ISM masses that we infer are in general lower than predictions from stellar mass loss models, suggesting that the recycled gas is heated or expelled in a timely manner (although our understanding of stellar mass loss is not complete; e.g., see Athey et al., 2002; Temi et al., 2007; Martini et al., 2013; Simonian & Martini, 2016). Given the typically low star formation rates of early-type galaxies, it is interesting to ask where the gas comes from, whether it is long-lived or recently accreted/cooled, and what excites the gas.

One way to proceed is to ask whether the presence of ionized gas correlates with other properties of the galaxies. In particular, early-type galaxies are typically divided into two classes that are thought to be related to their formation history (e.g., Kormendy & Bender, 2009; Cappellari, 2016). Less massive early-type galaxies tend to be flattened, have younger stellar populations, and show regular rotation around the major axis (fast rotators; e.g., Emsellem et al., 2007). In contrast, the most massive early-type galaxies tend to be round and old and have very little net rotation (slow rotators). Slow rotators show a number of indirect signs of past merging activity, including kinematically decoupled components (e.g., Krajnović et al., 2011; Naab et al., 2014).

An excellent data set to address these questions is provided by the integral-field ATLAS3D survey (Cappellari et al., 2011), a volume-limited sample of early-type galaxies with \log(M_{*}/M_{\odot})\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}9.8 ($M_{K}<-21.5$ mag) out to a distance of 42 Mpc. The ATLAS3D survey (and its predecessor SAURON; de Zeeuw et al., 2002) find a wide range of kinematic misalignment angles between the stars and gas in slow rotators, suggesting that in general the gas in these systems has been accreted (Sarzi et al., 2006; Davis et al., 2011; Sarzi et al., 2013). On the other hand, more than half of fast rotators are observed to host warm ionized gas kinematically aligned relative to the stars, which at least suggests that in a significant fraction of fast rotators, the gas may arise from internal processes such as stellar mass loss (Sarzi et al., 2006; Davis et al., 2011), although kinematic alignment can occur in merger remnants as well (van de Voort et al., 2015). Galaxies in richer environments are more likely to show kinematic alignment between the stars and gas compared to more isolated galaxies (Davis et al., 2011). In addition, the most massive fast rotators in ATLAS3D, with \log(M_{*}/M_{\odot})\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10.9, uniformly show kinematic alignment between the gas and stars. However, because of the small volume limit of the ATLAS3D survey, they have poor coverage of the most massive early-type galaxies (only 14% of their galaxies are slow rotators; Emsellem et al., 2011).

The MASSIVE survey (Ma et al., 2014) of early-type galaxies within 108 Mpc (including the Coma cluster) provides an ideal sample to study the warm ISM in galaxies that are a factor of ten more massive than the typical galaxies included in previous systematic surveys such as ATLAS3D, and to probe a wider range of slow rotator properties. Since the galaxies we will consider in this paper all have $\log(M_{*}/M_{\odot})>11.5$ ($M_{K}<-25.3$ mag), our sample is dominated by slow rotators (Veale et al., 2016). Using integral field unit (IFU) spectroscopy obtained as part of the MASSIVE survey (Ma et al., 2014), we will present the first systematic investigation of the existence, spatial distribution, and kinematics of warm ionized gas in the most massive galaxies in the local Universe.

The structure of this paper is as follows. In section II, we present our sample of MASSIVE galaxies. In section III, we describe our methods for detecting warm ionized gas with our IFU observations. In section IV, we present the detection fractions, spatial extent, line ratios and kinematics of the gas. In section V, we consider the implications of our results for the origins of excitation mechanisms of the gas. We summarize our results and discuss future work in section VI.

II. Sample and Data Description

Since the MASSIVE sample is thoroughly described in Ma et al. (2014), here we will only briefly review the selection criteria. The MASSIVE survey targets the most massive local early-type galaxies by requiring $M_{K}<-25.3$ mag (based on the 2MASS Redshift Survey; Huchra et al., 2012) and applying a morphological cut to remove large spiral and interacting galaxies (based on the Hyperleda database; Paturel et al., 2003). The survey is volume-limited out to a distance of 108 Mpc, which results in a factor of ten greater volume compared to previous similar IFU surveys (e.g., ATLAS3D; Cappellari et al., 2011) and allows for the inclusion of many more massive early-type galaxies such as those in the Coma cluster.111We have removed NGC 7681 from our sample even though it hosts warm ionized gas near its nucleus because it is actually a close pair and each individual galaxy is unlikely to make it into the MASSIVE sample (see also Veale et al., 2016). The full MASSIVE sample consists of 116 galaxies, but the sample size for this paper is 74 galaxies (these constitute, with a few exceptions, the brighter two-thirds of the survey with $M_{K}<-25.5$ mag).

This paper is based on integral field spectra obtained with the Mitchell Spectrograph (formerly VIRUS-P) on the 2.7m Harlan J. Smith Telescope at McDonald Observatory. The Mitchell IFU spectrograph has a $107^{\prime\prime}\times 107^{\prime\prime}$ field of view with 246 $4^{\prime\prime}$-diameter fibers per pointing. The spectra cover the wavelength range from $\sim 3500$Å to $\sim 5800$Å. Data reduction details are given in Greene et al. (2012, 2015), but see also Veale et al. (2016) for additional details. We generally have three dithered pointings for a total of 738 spectra per galaxy. We work with the single-fiber IFU spectra as well as a radially and azimuthally binned IFU dataset wherein fibers are grouped together to ensure continuum $S/N\geq 20$ in each spatial bin. Fibers with spectra that have abnormally high median fluxes, exceptionally low continuum S/N, or are near foreground stars are masked. The spatial bins are constructed from the individual fibers (remaining after masking) from the center outward, leaving central fibers that already meet the S/N threshold unbinned. Beginning at the first fiber below the S/N threshold, we divide fibers into circular radial bins and choose azimuthal bins (of equal size within each annulus but not across annuli) that minimize the aspect ratio to avoid thin rings or wedges. The radius of each annulus is then increased, from center outward (reducing the number of azimuthal bins when needed to maintain the proper aspect ratio), until the S/N threshold is reached (for more details, see Veale et al., 2016).

This paper represents part of our comprehensive effort to study gas of all phases in the MASSIVE sample, including using CO to trace molecular gas (Davis et al., 2016) as well as hot X-ray gas (Goulding et al., 2016).

III. Methods

Here we describe our methods for detecting warm ionized gas emission in MASSIVE galaxies.

III.1. Stellar Continuum Subtraction

For robust detection of faint lines, it is crucial that we model and subtract the stellar continuum. We must also ensure that each spectrum is at rest with respect to the stars. For each IFU spectrum, we therefore fit the stellar continuum using the penalized pixel-fitting method (pPXF; Cappellari & Emsellem, 2004). We only solve for the first two moments ($V_{\rm stars},\sigma_{\rm stars}$) of the full Gauss-Hermite decomposition (see van der Marel & Franx, 1993; Gerhard, 1993). We adopt the Bruzual & Charlot (2003) library of stellar population templates that span a grid of ages (5 Myr to 12 Gyr) and metallicities ($Z=0.004$, 0.02, and 0.05). To help account for flux calibration and sky subtraction errors, we fit with sixth order additive and sixth order multiplicative polynomials. We found that our H$\beta$, [O II] and [O III] equivalent widths did not change significantly if we instead used first-order additive and twelfth-order multiplicative polynomials. This is an important check because additive polynomials, unlike multiplicative polynomials, can change the equivalent widths of stellar absorption lines (e.g., see Koleva et al., 2009).

Each best-fit composite template spectrum is subtracted from its corresponding observed spectrum. The residual spectrum is used to measure the emission lines. Note that Sarzi et al. (2006) advocate an iterative approach; we do not implement such an algorithm here. We verified that the RMS of the residual spectrum is roughly comparable to the pixel-by-pixel uncertainties of the original spectrum, which means that we can use the RMS of the local continuum around an emission line to define the significance of our fits (see next subsection). All gas and stellar velocities presented in this paper are relative to the galaxy systemic velocity (see also Veale et al., 2016).

III.2. Emission Line Detection Algorithm

Since our spectra span the wavelength range from $\sim 3500$Å to $\sim 5800$Å, we constrain the presence of the [O II] 3727Å , [O III] 4959,5007Å , H$\beta$ and [N I] 5199Å emission lines as follows. We begin by fitting [O II] with a single Gaussian plus zero-slope additive continuum offset. The Gaussian is constrained to have positive amplitude, $1.0\rm\AA<\sigma_{\rm gas}<4.0\rm\AA$ (80 to 320 km s*-1* velocity dispersion), and a maximum central wavelength of $\mu=3727\rm\AA\pm 10\rm\AA$. The constraint on $\sigma_{\rm gas}$ prevents very narrow Gaussians from being fit to noise spikes and very wide Gaussians from being fit to broad continuum residuals. The constraint on $\mu$ ensures that we allow only a reasonable range of gas velocities ($\pm 800$ km s*-1*). Although it is of course possible that emission lines can exist with velocities beyond that which our range of $\mu$ allows, we have confirmed through visual inspection of the spectra that no such cases exist in our sample.

We define the significance of a fit as the ratio of the line amplitude to the RMS of the local continuum residuals (i.e., after best-fit template subtraction), as in Sarzi et al. (2006): $\mathcal{S}\equiv A/N$. [O II] fits with $\mathcal{S}\geq 3$ are considered to be significant. For fits with $\mathcal{S}<3$, we derive the $3\sigma$ upper limit on the integrated flux of [O II] using the known relation between Gaussian amplitude, width, and integrated area while assuming $A=3N$ and $\sigma_{\rm gas}$=200 km s*-1*.

The fit to [O III] is similar to that for [O II], but we fix $V_{\rm gas}$ and $\sigma_{\rm gas}$ to the values derived for [O II]. H$\beta$ is simultaneously fit with the [O III] doublet by fixing its kinematics to those of [O II]. In cases where [O II] is not detected, we allow $V_{\rm gas}$ and $\sigma_{\rm gas}$ to be free parameters of the shared line profile for [O III] and H$\beta$ (with similar constraints as for [O II]). We require $\mathcal{S}\geq 3$ for significant [O III] detections, but $\mathcal{S}\geq 2.5$ for H$\beta$ fits. It is more difficult to constrain the presence of H$\beta$ because of the strong stellar absorption, but after we subtract the expected absorption contribution from stars based on our best-fit pPXF composite templates, excess H$\beta$ emission can be revealed. When necessary, upper limits on [O III] and H$\beta$ are derived similarly to what would be done for [O II].

The [N I] emission line falls in the Mg$b$ wavelength region where the underlying stellar continuum is only poorly modeled by current stellar population templates, as described by Sarzi et al. (2006). Therefore, we require more stringent significance criteria for it, following Sarzi et al. (2006, 2010). Both [O III] and H$\beta$ must be significantly detected in order for the [N I] fit to be trusted. If both [O III] and H$\beta$ are indeed significantly detected, and [N I] itself has $\mathcal{S}\geq 4$, then [N I] is considered to be significantly detected. Otherwise, an upper limit is derived for [N I] as described above. We detect [N I] only in three spectra, and these are the fibers closest to the photometric centers of NGC 0708, NGC 1167, and NGC 1453.

The resolving power of our spectra is $R\approx 830$ (roughly 150 km s*-1*) near [O II], and $R\approx 1060$ (roughly 120 km s*-1*) near H$\beta$, [O III] and [N I]. We subtract these instrumental broadening values in quadrature from our measured line widths. However, the procedure for [O II] is a bit more involved because it is actually an unresolved doublet with a line separation of only 2.8Å.222[N I] is also an unresolved doublet with a line separation of only 2.4Å, and in principle we should also correct its measured line width for blending. However, [N I] is in a poorly modeled part of the stellar continuum and we rarely detect it, as described above. Therefore, we only correct the measured line width of [N I] for instrumental broadening. Before correcting for instrumental broadening of [O II], we apply an empirical deblending calibration to its measured line width, which is derived as follows (see also Salviander et al., 2007). For a grid of input velocity dispersions (50 to 300 km s*-1*), we create and sum two Gaussians with the same amplitude, fixed line separation, and input line width. We fix the amplitude ratio [O II]$\lambda 3728.8$/[O II]$\lambda 3726.0$ to one because this corresponds to electron number density $n_{e}=100$ cm*-3*, which is sensible for the diffuse ISM of early-type galaxies where the actual densities are much lower than the critical density for both doublet lines, $n_{e}\sim 10^{3}$ cm*-3* (see, e.g., Osterbrock, 1989). We then fit a single Gaussian to each of the summed Gaussians and obtain a one-to-one mapping between intrinsic line width and blended line width. We use a quadratic polynomial fit to the empirical mapping to convert measured line widths to intrinsic line widths. If the deblending calibration results in an intrinsic line width that is below our resolution limit, we simply set the line width to the resolution limit; otherwise, we further subtract the instrumental broadening in quadrature.

Figure 1 shows the central fiber spectrum of NGC 1453 in which [O II], [O III], H$\beta$ and [N I] are all significantly detected. Figure 2 shows example two-dimensional maps of various quantities for NGC 1453.

III.3. Sensitivity Limits

In Appendix A, we discuss simulations that assess our integrated flux completeness limits and uncertainties for [O II], H$\beta$ and [O III]. Here we will only briefly comment on the results of those simulations. Uncertainties on the measurements of detected emission lines (notably the integrated flux) were assumed to be $10\%$ in all cases based on the simulations. We can reliably detect the [O II] emission line down to a minimum integrated flux of $\sim 3\times 10^{-16}$ erg s*-1* cm*-2*.

We reach an EW limit of $\sim 2$Å and $\sim 0.5$Å for [O II] and [O III], respectively, in the central regions of galaxies where we work with individual fibers. In the binned data, we reach nearly an order of magnitude deeper in EW sensitivity, although stacking only works well when the gas fills much of a bin. Given that the stellar continuum level is three to four times fainter near [O II] than near [O III], we have comparable detection limits for the two lines in terms of integrated flux. The higher EW detection limit of [O II] is due to the lower continuum S/N at blue wavelengths. While we use [O II] as our primary detection line, none of our conclusions would change had we used [O III] instead.

Since we will compare our results to those of the SAURON (de Zeeuw et al., 2002) and ATLAS3D (Cappellari et al., 2011) surveys, it is useful to mention their sensitivities. They both use the SAURON instrument (Bacon et al., 2001) and utilize [O III] for line detection ([O II] is not in their observing band). The SAURON survey quotes limits in [O III] and H$\beta$ of $\sim 0.3$Å and $\sim 0.2$Å, respectively (Sarzi et al., 2006).

IV. Warm Ionized Gas Properties

IV.1. Detection Fractions

With the emission line measurements in hand, we now address whether a galaxy hosts warm ionized gas using [O II] as our primary tracer. This is complicated by false positive [O II] detections that can crop up in the outer parts of an IFU dataset (typically at the \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10\% level). We thus take two approaches to identify which galaxies host warm ionized gas. First, for the purpose of assigning a detection flag to each galaxy, we apply an additional [O II]-local continuum $S/N\geq 3$ cut. This additional $S/N$ cut removes most fibers within an IFU dataset except those within $\sim 15$ arcsec of the photometric center, which generally have [O II]-local continuum S/N\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10 and [O III]-local continuum S/N\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}100. The much higher local continuum $S/N$ ratios of these central fibers drastically reduces the chance of identifying noisy residual spikes as false positive [O II] detections.

If, after applying this continuum $S/N$ cut, there are $\geq 3$ fibers with [O II] detections remaining in the IFU data, then the galaxy is flagged as hosting warm ionized gas. If there are no fibers with [O II] detections, then the galaxy is flagged as hosting no detectable warm ionized gas. In cases with only one or two fibers with detections, we find that only fibers close to the center and with $\mathcal{S}>3.5$ are reliable (in such cases we flag the galaxy as hosting a nuclear warm ionized gas reservoir). Although one may worry that these criteria discard galaxies that only have off-center [O II] emission, we verified through visual inspection that there are no such cases in our sample.

The second complementary approach involves using the binned IFU spectra (see section II and Veale et al., 2016), which afford us significantly greater continuum $S/N$ at large galactocentric radii. In our outermost bins, we typically have [O II]-local continuum $S/N>4$ and [O III]-local continuum $S/N>30$. If the gas is uniformly distributed rather than patchy, binning can increase the significance of a detection. In every case where [O II] is detected in our binned spectra with S/N $\sim 20$ at large galactocentric radii, we also detect [O II] in the central fibers. When individual fibers at large galactocentric radii do not reveal significant emission, sometimes the corresponding binned IFU spectra do reveal real emission.

We find that $28/74\approx 38\pm 6\%$ of the MASSIVE galaxies in our sample host a warm ionized medium.333Unless otherwise noted, the uncertainties that accompany all detection fractions are computed assuming a binomial distribution: $\delta f=\sqrt{f(1-f)/N}$, where $f$ and $N$ are the detection fraction and sample size, respectively. The properties of these 28 galaxies are given in Table 1, ordered according to whether the gas is extended and rotating, extended and non-rotating, or centrally concentrated. Two-dimensional maps of various quantities for these galaxies are shown in Appendix B using the same order as in Table 1 (NGC 1453 is shown separately in Figure 2). Most of these galaxies host warm ionized gas that is detected only in the central fibers: $17/28=61\pm 9\%$. In $3/28\approx 11\pm 6\%$ galaxies, the warm ionized medium is spatially extended but patchy and non-rotating, and in $8/28\approx 29\pm 9\%$ of the MASSIVE galaxies, the warm ionized gas is spatially extended and shows regular rotation. If we split the full MASSIVE sample by stellar kinematic classification, then we find that $12/15=80\pm 10\%$ of MASSIVE fast rotators host warm ionized gas compared to only $16/58=28\pm 6\%$ MASSIVE slow rotators. A two-proportion z-test confirms that the detection fractions for MASSIVE slow and fast rotators are indeed statistically significantly different, with a z-score of 3.7 corresponding to a p-value $\ll 0.01$. The distance and stellar mass distributions of MASSIVE fast and slow rotators are not significantly different from each other, as confirmed by non-parametric statistical tests. This means that the significantly different detection fractions for the two kinematic classes cannot be attributed to an underlying bias in their distance and stellar mass distributions.

We are clearly more likely to detect warm ionized gas in MASSIVE fast rotators than slow rotators. We now ask whether there are other defining characteristics of the galaxies with warm ionized gas since that may provide indirect clues about the origin of the gas. For example, if cooling flows represent an important gas origin scenario for MASSIVE galaxies, then we may expect different parent populations of halo masses (e.g., Kereš et al., 2005). If instead the dominant factor is whether a galaxy can recycle and retain the mass lost by its old stellar population, then we might expect different underlying distributions of stellar velocity dispersion. We thus investigate the following physical properties: stellar mass, halo mass, and central stellar velocity dispersion, the distributions of which are shown in Figure 3. We used three different non-parametric tests (two-sample Kolmogorov-Smirnov, Mann-Whitney-Wilcoxon, and Anderson-Darling) to check for statistically significant differences between the distributions of the three physical properties for the MASSIVE sample with and without gas. All three tests returned p-values $>0.01$ in every case, which suggests no statistically significant differences between the underlying parent populations.

[FIGURE:]

IV.2. Spatial Extent

Because of the roughly $2\times 2$ arcmin2 field of view of the Mitchell Spectrograph, we are able to study the spatial extent of the gas out to generally $2-4\times R_{e}$. We use the binned IFU datasets described in section II to determine via visual inspection the radially furthest bin in which an [O II] emission line exists.

Of the 28 MASSIVE galaxies with warm ionized gas, 11 galaxies have spatially extended gas and the remaining 17 galaxies have [O II] emission only in a few fibers near their photometric centers. For the 11 spatially extended galaxies, we adopted the midpoint of the furthest bin with reliable [O II] emission as the maximum radial extent of the warm ionized gas ($R_{\rm gas}$), and for the 17 compact cases, we assumed an upper limit on $R_{\rm gas}$ of 5 arcsec (recall that the fiber radii are 2 arcsec, but there are usually three dithers). The adopted $R_{\rm gas}$ values are given in Table 1. The maximum physical extents span a wide range from only 0.6 kpc up to 18.2 kpc, with a median of 2.3 kpc; likewise, the $R_{\rm gas}/R_{e}$ fractions run from 0.12 up to 3.71, with a median of 0.25.

The furthest radial bin with detectable [O II] emission contains on the order of $\sim 100$ stacked individual fibers. Since we are binning so many fibers, a few of which can have detected emission (occasionally false positives), we ask whether the individual fiber detections might dominate the signal detected in the binned spectrum. In fact, the individual fibers constitute $<10\%$ of the integrated flux of the line detected in the binned spectrum. Thus, the binned detection suggests that there is extended low-surface-brightness emission filling a large fraction of the bin.

IV.3. Comparison to SAURON and ATLAS3D

Our $\sim 39\pm 6\%$ global detection fraction is lower than the $36/48=75\pm 6\%$ detection fraction for SAURON’s less massive early-type galaxies (Sarzi et al., 2006). While the Sarzi et al. (2006) work has not been extended to the full ATLAS3D sample, Davis et al. (2011) report that $67\pm 3\%$ of fast rotators and $86\pm 6\%$ of slow rotators in ATLAS3D have detectable warm ionized gas emission. Similarly, Nyland et al. (2016) report a warm ionized gas detection fraction of $54/101\approx 53\pm 5\%$ for the subset of ATLAS3D galaxies with radio observations. Both of these quoted ATLAS3D detection fractions are higher than ours. In this section, we will explore the detectability and resolvability of ATLAS3D and SAURON gas reservoirs if they instead existed in the more luminous and distant MASSIVE galaxies.444The median distances of the SAURON and ATLAS3D samples are 16 Mpc and 24 Mpc, respectively (Sarzi et al., 2006; Cappellari et al., 2011), substantially lower than the median distance of 80 Mpc for the MASSIVE sample (Ma et al., 2014). We compare with an ATLAS3D catalog of ionized gas measurements (M. Sarzi, private communication). We only consider ATLAS3D galaxies that have detected warm ionized gas by requiring that the spatially integrated [O III] EW $>0.02$Å as in Davis et al. (2011); there are 183 such objects.

First, if the ATLAS3D gas reservoirs were relocated to 80 Mpc (the median distance of the MASSIVE sample), would they be spatially resolved with our observations? Consider first that the 183 ATLAS3D galaxies with detected warm ionized gas span a range of physical radii from $0.2-4.4$ kpc, with a median of $1.6$ kpc. With a few exceptions on the smaller end, these values are consistent with the physical radii of MASSIVE compact (unresolved) gas reservoirs, whose upper limits on $R_{\rm gas}$ of five arcsec correspond to physical radii of $0.6-2.5$ kpc, with a median of $1.8$ kpc. Indeed, we find that $121/183\approx 66\%$ of ATLAS3D gas reservoirs would have a radial gas extent of less than five arcsec if they were relocated to 80 Mpc, which means that the majority of the ATLAS3D sample would be classified as spatially unresolved in our observations.

Second, would our observations be sensitive enough to detect emission from ATLAS3D gas reservoirs, if they were relocated to the median MASSIVE distance of 80 Mpc? Only $101/183\approx 55\%$ of ATLAS3D gas reservoirs shifted to 80 Mpc would have [O III] EW above 0.5Å within a single Mitchell fiber. If we crudely convert the ATLAS3D [O III] EW to [O II] EW by multiplying the former by a factor of six, then $123/183\approx 67\%$ of ATLAS3D gas reservoirs relocated to 80 Mpc have [O II] EW above 2Å within a single Mitchell fiber.555The factor of six assumes [O II]/[O III]$=2$, which is typical for our MASSIVE galaxies as well as stacked SDSS spectra of extended LINER-like objects (without reddening corrections; Johansson et al., 2016). We also assume that the stellar continuum is a factor of three fainter near [O II] as compared to near [O III] (see subsection III.3). This exercise suggests that the approximately factor of two lower warm gas detection fraction for MASSIVE galaxies, as compared to both the SAURON and ATLAS3D samples, might simply be attributable to sensitivity differences. Therefore, while the warm ionized gas detection fraction might indeed drop as a function of mass, we do not have good enough statistics to robustly constrain that expected trend.

We also note that the warm ionized gas detection fractions of fast rotators and slow rotators in ATLAS3D are comparable, which is not the case for MASSIVE. In addition, the slow rotator detection fraction is much higher for ATLAS3D ($\sim 86\%$) than for MASSIVE ($\sim 28\%$). These differences between the two surveys appear to be significant even though there are relatively few fast rotators in MASSIVE and relatively few slow rotators in ATLAS3D. While the distributions of [O III] EW for slow and fast rotators in ATLAS3D are statistically indistinguishable (according to three non-parametric tests), there is a longer tail toward high EWs for the fast rotators. Since we are not as sensitive, it is possible that a large fraction of ATLAS3D slow rotators with relatively low [O III] EW would not be detected by us. Alternatively, it may be that the slow rotators have patchier gas distributions than the fast rotators, and thus our larger spaxels hurt us disproportionately.

IV.4. Kinematics: Gas versus Stars

A comparison of the kinematics of the stars and the warm ionized gas can reveal clues about the origin of the gas. For this reason, we run the fit_kinematic_pa routine described in Krajnović et al. (2006) to determine the kinematic position angle (PA) of the warm ionized gas for the eight galaxies in which the gas is rotating and extended.666The fit_kinematic_pa Python routine is publicly available at http://www.ascl.net/1601.016. We define the kinematic gas PA as the angle between north and the line that bisects the gas velocity field, measured counterclockwise from north to the redshifted emission side. In order to alleviate the impact of false positive fibers at large galactocentric radius and explore how well the PA can be constrained, we run fit_kinematic_pa three times, applying different fiber selection criteria each time and then taking the mean and standard deviation of the resulting values. The resulting kinematic gas PAs are given in Table 1. Note that it is difficult to constrain the kinematic position angle of the gas in NGC 1167 because rotation is observed in the gas only out to $\sim 1R_{\rm e}$ and because NGC 1167 harbors a nearly face-on disk (see also subsection V.3).

We have also measured kinematic stellar PAs for MASSIVE galaxies using kinemetry (Krajnović et al., 2006). These values are listed in Table 1; additional details and measurements for the full MASSIVE sample will be presented in Ene et al. (in preparation). To ensure consistency with the kinematic gas PAs, our kinematic stellar PAs are also measured counterclockwise from north to the redshifted emission side. MASSIVE is dominated by slow rotators, and in some cases, there simply is not a preferred axis of rotation. Of the four slow rotators with extended and rotating warm ionized gas, NGC 2208 and NGC 7550 do not have measurable kinematic stellar PAs and we therefore do not consider them further in this section. The kinematic and photometric stellar PAs of the two other slow rotators, NGC 0708 and NGC 2783, are misaligned. In contrast, all four of the MASSIVE fast rotators that we consider here (NGC 1167, NGC 1453, NGC 2320, and NGC 7426) have measurable stellar kinematic PAs, and their kinematic and photometric stellar PA values are well-aligned.

In Figure 4, we show the kinematic misalignment angle between the stars and warm ionized gas for MASSIVE slow and fast rotators as compared to ATLAS3D slow and fast rotators. Kinematic classifications of MASSIVE galaxies into slow and fast rotators are presented in Veale et al. (2016) whereas the kinematic misalignment angles and kinematic classifications of ATLAS3D galaxies are taken from Krajnović et al. (2011), Davis et al. (2011) and Emsellem et al. (2011). Misalignment angles for MASSIVE galaxies are defined as $\rm PA_{\rm gas-stars}\equiv\rm PA_{\rm gas}-\rm PA_{\rm stars}$, where $\rm PA_{\rm gas}$ is the kinematic position angle of the warm ionized gas and $\rm PA_{\rm stars}$ is the position angle of the stars (kinematic when available, else photometric). The misalignment angles for all ATLAS3D galaxies are based on kinematic stellar PAs. All misalignment angles are defined to lie between [math] and $180$ degrees.

In general, misalignment between stars and gas points towards an external origin of the gas, while alignment between gas disks and stars suggests that the gas originated from the stars. However, see Serra et al. (2014) and Lagos et al. (2015) for the expected misalignment when gas cools out of the halo, and van de Voort et al. (2015) for the timescales on which externally accreted gas is expected to become aligned with the stars due to gravitational torques. The MASSIVE fast rotators are predominantly aligned ($3/4=75\pm 22\%$), which is consistent with the results from ATLAS3D that more massive fast rotators tend to host the warm ionized gas kinematically aligned with the stars (Davis et al., 2011). In contrast, the two MASSIVE slow rotators are consistent with misalignment between the stars and gas, which is in agreement with the interpretation that slow rotators predominantly get their gas via external means (Sarzi et al., 2006; Davis et al., 2011). Our total number of rotating and spatially extended gas reservoirs is too small to say anything definitive about the distribution of misalignment angles. Nevertheless, the MASSIVE fast rotator NGC 1453 is especially interesting because its gas and stellar PAs are nearly perpendicular to each other; we will return to this galaxy in section V.

IV.5. Line Ratios and Radial Profiles

We now turn to the sources of ionization for the warm ionized gas in MASSIVE galaxies. There are five ionizing sources typically considered (see Sarzi et al., 2010, and references therein): photoionization by an active galactic nucleus (AGN; e.g., Ho et al., 1997; Ho, 2008), photoionization by evolved stars (e.g., Binette et al., 1994; Stasińska et al., 2008), photoionization by a young stellar component (e.g., Graves et al., 2007; Kaviraj et al., 2007), excitation by shocks (e.g., Koski & Osterbrock, 1976; Dopita & Sutherland, 1995; Allen et al., 2008), and thermal conduction via the hot X-ray-emitting gas (e.g., Sparks et al., 1989; Mathews & Brighenti, 2003; Sarzi et al., 2013).

Since we do not have access to standard excitation diagrams that make use of the [N II]$\lambda\lambda 6548,6584$/H$\alpha$ ratio (e.g., Baldwin et al., 1981; Veilleux & Osterbrock, 1987), it is difficult to distinguish between the first four sources of ionization for our sample. Nevertheless, we can speculate using two alternative methods: (1) use non-standard excitation diagrams to find extraordinary cases where a single source of ionization can be unambiguously identified (e.g., see Sarzi et al., 2006, 2010; Annibali et al., 2010; Yan & Blanton, 2012), and (2) check whether the radial profile of the ionized gas flux follows the radial profile of the stellar continuum flux as probed by the emission line EW. As we obtain more X-ray observations of MASSIVE galaxies using Chandra (PI: Goulding), we will be in a position to also examine trends between X-ray properties and the presence of warm ionized gas.

The first clue about the sources of ionization comes from emission line ratios.777Due to the limited wavelength coverage of our spectra, we can only measure the [O III]/H$\beta$ and [N I]/H$\beta$ ratios. Since we can only compute robust line ratios for the central parts of our galaxies, we will use only the value in the central fiber for each galaxy. However, we should note that our 4 arcsec diameter central fiber typically covers the inner $\sim 1.5$ kpc of our galaxies. In such a large physical aperture, even in cases where an AGN plays some role in the photoionization of the gas, there are likely to be contributions from other sources (e.g., Eracleous et al., 2010). Therefore, aside from relatively luminous Seyfert galaxies such as NGC 1167 (see below), we cannot say much about optical AGN activity from these data. In standard excitation diagrams (e.g., Baldwin et al., 1981; Veilleux & Osterbrock, 1987; Kewley et al., 2006), the [O III]/H$\beta$ ratio is combined with another strong line, typically [N II]/H$\alpha$, to distinguish between H II regions, shocks, and an AGN. Star-forming galaxies follow a path toward higher [O III]/H$\beta$ and lower [N II]/H$\alpha$ values as metallicity decreases and the ionization parameter increases (e.g., see Kennicutt, R., 1998; Kewley & Dopita, 2002; Kauffmann et al., 2003; Tremonti et al., 2004; Kewley et al., 2006; Moustakas & Kennicutt, 2006). If we assume that the high stellar metallicity that characterizes our galaxies (Greene et al., 2015) also applies to the gas, then we can more easily map the [O III]/H$\beta$ ratio to sources of ionization (see discussion in Section 6.2 of Nyland et al., 2016). High values of [O III]/H$\beta$ (\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}3) are likely to point towards AGN activity, while very low values (\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}0.5) would point towards star formation.

In our sample, NGC 1167 has [O III]/H$\beta\approx 6$, which is indicative of AGN activity (Groves et al., 2004a, b). The remaining galaxies have intermediate [O III]/H$\beta$ measurements (or lower limits when H$\beta$ is not detected). These intermediate [O III]/H$\beta$ values are a defining feature of low-ionization nuclear emission regions (LINERs; e.g., see Heckman, 1980; Kewley & Dopita, 2002), which may arise from excitation by evolved stars, shocks, or a low-luminosity AGN (e.g., Leitherer et al., 1999; Dopita et al., 2000; Ho, 2008; Allen et al., 2008; Eracleous et al., 2010; Groves & Allen, 2010). In several of our galaxies, non-central fibers tend to also exhibit these so-called “extended LINER-like” excitation ratios (i.e., low ionization emission regions or LIERs; Belfiore et al., 2016; Johansson et al., 2016). In principle, the [N I]/H$\beta$ ratio can help us distinguish between excitation grids for fast shocks and H II regions (see detailed discussion by Sarzi et al., 2010), but the [N I] line is detected in the centers of only three galaxies. One of those galaxies is the aforementioned NGC 1167, which is already consistent with central excitation by an AGN based on its high [O III]/H$\beta$ ratio alone. For the remaining two galaxies (NGC 0708 and NGC 1453), the combination of [O III]/H$\beta$ and [N I]/H$\beta$ ratios is more consistent with excitation grids for fast shocks than H II regions, which again broadly suggests LINER-like emission.

Given the limitations of the line ratio approach for our dataset, we now turn to the spatial regime. If the source of ionization was centrally concentrated, then we would expect an upturn in emission line EW toward the center of the galaxy along with central excitation ratios suggestive of, say, an AGN. In the simplest case where a uniformly distributed old stellar population is responsible for ionizing the gas, we would expect a roughly constant radial EW profile. In Figure 5, we show the radial profiles of logarithmic [O II] EW as well as the best-fit straight lines for the 11 MASSIVE galaxies with extended warm ionized gas. Note that the radial scale of the EW profile is different for every galaxy because both $R_{e}$ and $R_{\rm gas}$ are different. A straight line model for $\log_{10}$(EWOII/[Å]) as a function of galactocentric radius does not always capture the complexity of the data, but at least half of the cases shown are consistent with a slope of zero.

In several cases, the data unambiguously reveal a non-constant trend in the equivalent width radial profile. An example is NGC 1167 in which the [O II] EW rises non-linearly and rapidly toward the center of the galaxy, suggesting that an AGN dominates as the source of ionizing photons near the center of this galaxy (this upturn toward the center is also seen in the radial EW profile of H$\alpha$ by Gomes et al., 2016b, in the CALIFA survey). A similar phenomenon can be seen for NGC 6482, but it is less clear due to the different radial scale of the EW profile. In contrast, NGC 7265 presents an increasing [O II] equivalent width radial profile. Since the observed EW depends on both the surface brightness of the starlight and the absorption of the photoionizing continuum by the gas (which itself depends on the geometry of the gas, dust, etc.), the expected radial dependence of the EW in different heating scenarios is not obvious (e.g., Yan & Blanton, 2012; Papaderos et al., 2013; Gomes et al., 2016b, a). In a simple geometric model in which the stars are spherically distributed and the gas is in a disk, the EW would actually rise with radius (Sarzi et al., 2010), but generally in the SAURON galaxies the EWs are flat or fall slightly with radius.

Taking post-asymptotic giant branch (pAGB) stars as the main source of photoionization, stellar population synthesis models predict an H$\alpha$ EW range of $\sim 0.5-3$Å (e.g., Gomes et al., 2016b; Belfiore et al., 2016). Indeed, the stacked SDSS spectrum of extended LINER-like objects presented by Johansson et al. (2016), which are presumably powered by pAGB stars, has an H$\alpha$ EW of 1.3Å. Their empirical spectrum also has an H$\alpha$/[O II] EW ratio of $\approx 0.3$, which corresponds to a predicted [O II] EW of $\approx 4.3$Å (see their Table 2 and Figure 3). This predicted [O II] EW, which is generally attributed solely to pAGB stars, is well above our own [O II] EW detection limit of 2Å. Thus, there is apparently a real difference between the emission line EWs in MASSIVE galaxies and the average galaxy considered by Johansson et al. (2016). There are many possible sources for this discrepancy: the typical line EW due to pAGB stars might be different in the most massive galaxies (e.g., because of a higher metal content or different mean stellar ages), the total amount of warm ISM available to be ionized might genuinely be lower, or the line ratios could simply be different as a function of galaxy mass (again due to differences in metallicity or the ionizing spectrum from pAGB stars).

IV.6. Warm Ionized Gas Masses

For all 28 MASSIVE galaxies that host warm ionized gas, we estimate gas masses as follows. Within $R_{\rm gas}$, we take all fibers which have $\mathcal{S}\geq 3$ for [O II] and assume that [O II]/H$\beta=6$ (based on the mean of the central fibers of NGC 1453, in which we have the greatest number of H$\beta$ detections) to compute an H$\beta$ luminosity. To compute the total warm ionized gas mass as traced by the H$\beta$ recombination line, we use equation (1) of Nesvadba et al. (2011), which we rewrite here for completeness:

[TABLE]

Here, $\mathcal{L}_{\rm H\beta,43}$ is the H$\beta$ luminosity in units of $10^{43}$ erg s*-1* and $n_{\rm e,100}$ is the electron number density in units of 100 cm*-3*. Due to the limited wavelength coverage of our spectra, we do not have independent constraints on $n_{e}$ for each galaxy, so we adopt the common value of $n_{e}=100$ cm*-3* based on observations of other early-type galaxies (see Mathews & Brighenti, 2003). The derived warm ionized gas masses and corresponding $\mathcal{L}_{\rm H\beta}$ are given in Table 1. There is roughly an order of magnitude systematic uncertainty in the gas masses due to the unknown gas density and line ratios, but the derived masses are quite consistent at this level with previous work (see also Goudfrooij, 1999).

Figure 6 shows the maximum radial extent of the warm ionized gas as a function of the warm ionized gas mass and the halo mass. The warm ionized gas masses span $10^{4}-10^{5}M_{\odot}$. As expected, cases of compact centralized gas detections have lower warm ionized gas masses than spatially extended gas reservoirs. The warm ionized gas tends to have the greatest radial extent in galaxies with a patchy, non-rotating spatial gas distribution.

V. Discussion

We will attempt to use the morphology, kinematics, and ionization of the gas to infer some constraints on its origins. We consider origins both internal (stellar mass loss or cooling from a hot halo) and external (merging or gas accretion).

V.1. Gas Morphology

Among our sample of 28 MASSIVE galaxies with detected warm ionized gas, there exist a wide variety of gas morphologies. It is instructive to consider whether the morphology of the gas can help distinguish between various possible gas origin scenarios. Motivated by the work of Macchetto et al. (1996), we can roughly classify all of our different gas reservoirs into three broad categories: (1) smooth, extended, and rotating disks, (2) extended and filamentary/patchy gas, and (3) unresolved/small reservoirs. The majority of our detections are of the unresolved/compact variety ($17/28\approx 61\pm 9\%$). In our spatially resolved cases, we are detecting the most extreme gas reservoirs in terms of spatial extent, and filamentary emission is less frequently detected ($3/28\approx 11\pm 6\%$) than smooth rotating emission ($8/28\approx 29\pm 9\%$). Our detection fractions for each of the three morphological categories are roughly in agreement with the similarly constructed categories of Macchetto et al. (1996).

There is little that we can say about the unresolved gas, except that it likely resides in the small disks that are seen in 75% of early-type galaxies (e.g., see van Dokkum & Franx, 1995; Ferrarese & Ford, 1999; Martini et al., 2003). These dust disks can have semi-major axes up to several hundred parsecs (see Table 2 of van Dokkum & Franx, 1995). Given that our 2 arcsec radius fibers correspond to $\sim 0.75$ kpc at 80 Mpc and that we cannot resolve $R_{\rm gas}$ below $1-3\times 0.75$ kpc (based on how many central fibers the gas is detected in), it is entirely possible that in our 18 galaxies with compact gas reservoirs, the ionized gas emission comes from such a central dust disk. It is also possible that there is extended but non-detected gas in these galaxies (see subsection IV.3).

V.2. Origin of the Gas

We will now step through the different possible gas origin scenarios by presenting emblematic examples of each.

V.2.1 External accretion versus internal origins

Misalignment angles are often used as clues to the origin of the gas. While alignment may point to internal or external processes, misalignment is a clear sign of external accretion (e.g., Sarzi et al., 2006; Davis et al., 2011). In contrast to external accretion, stellar mass loss is a prominent internal source of gas. It is difficult to say for sure which of our galaxies have gas reservoirs consistent with an internal recycling origin. Although kinematic alignment between the stars and gas is a clue in favor of stellar mass loss, such alignment can also occur in merger remnants within a few dynamical times due to gravitational torques exerted on the gas by the stars (e.g., van de Voort et al., 2015).

An extreme example of external accretion is the MASSIVE fast rotator NGC 1453 in which the warm ionized gas is kinematically misaligned with respect to the stars by $\sim 98$ deg. Buson et al. (1993) and Zeilinger et al. (1996) also found that the warm ionized gas is misaligned relative to the stars in NGC 1453, though with a larger misalignment angle of $\approx 128$ deg. It is unlikely that stellar mass loss from the old stellar population is the primary source of the warm ionized medium in NGC 1453 because such recycled gas would tend to remain in the potential well of the stars, leading to photometric and kinematic alignment between the stars and gas.

V.2.2 Cooling flows

In NGC 6482, the gas is patchy and its velocity structure may hint at cooling flows as the dominant origin scenario (e.g., see Fabian, 1994). NGC 6482 is known to be a fossil group system; it is a single elliptical galaxy with a group-like X-ray halo and in which all other group members are at least two optical magnitudes fainter (see Mulchaey, 2000; Khosroshahi et al., 2004). This galaxy is also known to harbor a patchy dust reservoir (Alamo-Martínez et al., 2012). We do not have the spatial resolution to study the detailed morphology of the gas, but if the origin of the warm ionized gas is indeed via cooling flows, then we would expect to see a connection with the hot X-ray-emitting halo gas. For example, McDonald et al. (2011) find a clear relationship between H$\alpha$ filaments and the X-ray properties of “cool core” clusters (negative radial temperature gradients; e.g., see Hudson et al., 2010). Future comparisons of deep narrow-band and X-ray imaging could help us test the cooling flow origin hypothesis in MASSIVE galaxies like NGC 6482. Also powerful would be comparisons to known “cool-core” cluster galaxies in MASSIVE which have warm ionized gas detections (NGC 0315, NGC 0533, NGC 0547 and NGC 0708) using deep high-resolution images that clearly reveal the morphology of the gas.

V.3. Special Cases

The MASSIVE fast rotator NGC 1167 represents an interesting case of both gas rotation near the center of the galaxy and a patchy distribution of gas at larger galactocentric radius. The galaxy itself is likely a face-on disk, and exhibits extended UV emission in images taken with the Galaxy Evolution Explorer (GALEX; Martin et al., 2005). Interestingly, Struve et al. (2010) detected a $\sim 160$ kpc disk of neutral hydrogen in NGC 1167 and suggested satellite accretion as its origin. In total, the warm ionized gas in NGC 1167 reaches roughly 18 kpc, but rotation is evident in the gas only within about 10 kpc. Since gas is detected in the outermost bins for NGC 1167 and there might be more detectable emission outside the field of view of the Mitchell Spectrograph, our $R_{\rm gas}$ can be considered a lower limit for this galaxy. We see evidence of an AGN based on the [O III]/H$\beta$ line ratio, and our detection of [O II] in binned fibers at large galactocentric radius hints at excitation by H II regions. As part of the CALIFA survey, Gomes et al. (2016b) found clear evidence of localized H$\alpha$ clumps at large galactocentric radii (see their Figure 1), which they attribute to low-level ongoing star formation. They also comment on the patchiness of the gas by classifying NGC 1167 as having a “perturbed” rotational pattern (see their Table 1). This is one clear case where residual star formation may accompany and indeed ionize the gas (an assertion further supported by the discovery of $3.3\times 10^{8}M_{\odot}$ worth of molecular gas, and a mid-infrared SFR of $\sim 0.3M_{\odot}$ yr -1; O’Sullivan et al., 2015; Davis et al., 2016).

Another very interesting case is NGC 7265. It is unclear what the source of the gas is in NGC 7265 since it is redshifted with respect to the systemic velocity and has a very asymmetric spatial extent with no blueshifted component. The gas velocities are consistent with the systemic velocity near the center of the galaxy but extend up to 250 km s*-1* at a radius of 7.4 kpc (see Figure 8). It is tempting to think we are observing an outflow driven by an AGN (e.g., Ciotti et al., 2016), in analogy with the ATLAS3D galaxy NGC 1266 (Alatalo et al., 2011; Davis et al., 2012; Nyland et al., 2013; Alatalo et al., 2015). However, with only redshifted emission, we cannot rule out an inflow, rotation, or simply an ongoing merger. We do detect a compact radio source (Nyland et al. in preparation) and a soft nuclear X-ray source (Goulding et al., 2016) in the center. Since we do not detect H$\beta$ we cannot say much about the nuclear line ratios. If this is an outflow, however, it is difficult to understand why we observe only redshifted emission.

VI. Future Work and Summary

We have investigated the existence and spatial distribution of warm ionized gas in a volume- and magnitude-limited sample of massive early-type galaxies with $\log(M_{*}/M_{\odot})>11.5$ ($M_{K}<-25.3$ mag) and distance $D<108$ Mpc using deep wide-field integral field spectroscopy. The main conclusions of this paper are:

We find that $28/74\approx 38\pm 6\%$ of MASSIVE early-type galaxies host warm ionized gas down to an [O II] EW limit of $\sim 2$Å. In eight of these cases, the warm ionized gas exists in the form of a large rotating disk, whereas in three cases, the gas is extended but patchy/non-rotating. The MASSIVE galaxies NGC 1167, NGC 1700, and NGC 6482 host the most extended warm ionized gas reservoirs, out to roughly 18, 12 and 16 kpc, respectively. The remaining 17 galaxies host detectable warm ionized gas only near their centers. Our detection fraction is lower than the $\sim 70\%$ measured by the SAURON and ATLAS3D surveys of nearer but less massive early-type galaxies. This difference is consistent with arising from our shallower surface brightness sensitivity. 2. 2.

The warm ionized gas detection fraction does not appear to depend on stellar mass, halo mass, or central stellar velocity dispersion. However, there is a statistically significant ($3.7\sigma$) dependence in terms of net angular momentum: $\sim 80\%$ of MASSIVE fast rotators host warm ionized gas compared to only $\sim 30\%$ of MASSIVE slow rotators. The much lower detection fraction for MASSIVE slow rotators compared to ATLAS3D slow rotators ($\sim 85\%$) can at least partially be attributed to sensitivity differences between the two surveys. 3. 3.

Of the four fast rotators with rotating and spatially extended gas, three show kinematic alignment between the stars and gas, which is consistent with previous findings by ATLAS3D that more massive fast rotators tend to have the gas and stars kinematically aligned. On the other hand, both slow rotators with measurable kinematic stellar and gas position angles show misalignment, which agrees with the idea that slow rotators generally get their gas through external means. 4. 4.

The most prominent gas excitation mechanisms are likely photoionization from evolved stars and fast shocks on extended kpc scales with possible additional contributions from low-luminosity AGN on nuclear scales. Using radial profiles of the [O II] EW and non-standard excitation diagrams, we cannot easily pinpoint the sources of ionization for the gas. The one exception is NGC 1167, whose central line ratios and radial [O II] EW profile suggest photoionization by an AGN for gas near the center. 5. 5.

It is likely that a variety of physical processes are responsible for the origin and ionization of the gas in MASSIVE early-type galaxies. We see examples of gas with large misalignment angles that were likely acquired via external accretion, a $\sim$16 kpc-scale patchy gas reservoir likely cooling from a hot halo, and even a possible AGN-driven outflow.

One exciting possibility opened up by our observations is to use the extended rotating disks we have uncovered to measure enclosed dynamical masses. These dynamical masses will be complementary to those that we are deriving from modeling the stars (e.g., Thomas et al., 2016), and will give us the opportunity to test the stellar dynamical models as well as examine the universality of the initial mass function (e.g., Treu et al., 2010; Cappellari et al., 2012; Conroy & van Dokkum, 2012). Such dynamical modeling will also allow us to test whether what we are calling warm ionized gas disks are actually so-called “red geysers,” which are thought to be large-scale AGN-driven outflows (Cheung et al., 2016).

Also of considerable interest is connecting our work on the warm ionized medium of MASSIVE early-type galaxies to other phases of their ISM. The relationship between the warm ionized and hot X-ray-emitting gas is still unclear (e.g., Forman et al., 1985; Heckman et al., 1989; Macchetto et al., 1996; Goudfrooij, 1999; Sarzi et al., 2013; Goulding et al., 2016), and we have approved Chandra time (PI: Goulding) to increase the MASSIVE sample with hot gas measurements. Similarly, the occurrence rate of warm ionized gas (and dust) as a function of radio continuum detection may reveal relationships to accreting supermassive black holes and the presence of central dust disks; this will be explored in a future study. Finally, finding direct relationships between the warm ionized and cold gas phases (both molecular and atomic) is important for addressing the fate of the gas in these massive quiescent systems; our studies of CO emission from these galaxies probe fundamental questions along these lines (Davis et al., 2016, and Davis et al. in preparation).

Acknowledgements

We thank Marc Sarzi for providing us with a catalog of ATLAS3D warm ionized gas measurements and for helpful suggestions. We also thank Duncan Forbes, Brent Groves, Yifei Jin, Davor Krajnovic and Freeke van de Voort for helpful comments. We thank the referee for a thorough and useful report that significantly improved the paper. VP thanks Bill Mathews, Joel Primack and X. Prochaska for helpful discussions. The MASSIVE survey is supported in part by NSF AST-1411945, NSF AST-1411642, HST-GO-14210, and HST-AR-14573.

Appendix A Integrated Flux Completeness Limits from Simulations

In this appendix, we present simulations that help to quantify the integrated flux completeness limits of our emission line detection algorithm. Our IFU spectra span a wide range of full continuum S/N ratios that vary with galactocentric radius and differ from galaxy to galaxy. The individual fibers closest to the photometric center have full continuum S/N \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}20 (reaching above 100 for some galaxies), and there is a gradual drop with radius to values as low as $\sim 1$ for the outermost individual fibers (but the binned spectra at these large galactocentric radii always have full continuum S/N \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}20). Therefore, we created a large set of synthetic spectra assuming three different continuum S/N ratio regimes: low ($S/N=10$), intermediate ($S/N=40$), and high ($S/N=70$). The synthetic spectra were constructed as follows. We began with the central fiber spectrum of a random galaxy in our sample (NGC 2274) that does not have emission lines, and degraded its continuum $S/N$ to one of the aforementioned fixed values by scaling its error spectrum and adding Gaussian random noise per pixel. By using the central fiber spectrum of a real galaxy rather than a best-fit composite pPXF stellar template (which naturally does not have emission lines), we are to some degree able to include systematic observational errors in our simulations (as well as features in wavelength regions that current stellar population models are not able to reproduce; e.g., the [N I] 5200Å region).

We then superimposed Gaussian emission lines assuming a range of input significance values: $\mathcal{S}=0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0,10.0$. With an assumed line significance level and a measure of the local continuum RMS level, we can solve for the input amplitude of a line ($A=N\times\mathcal{S}$). For a given combination of degraded continuum $S/N$ and input line significance $\mathcal{S}$, we created 100 realizations. Since we will use our simulations to only probe the detectability of integrated flux of emission lines, the velocity and velocity dispersion of simulated emission lines are all fixed to 0 and 200 km s*-1*, respectively. Our simulations will test the detectability of the [O II], [O III] 4959,5007, and H$\beta$ emission lines – we do not simulate the detectability of the [N I] line because we adopt quite stringent significance criteria for it anyway (see subsection III.2). We note that our simulations are not suitable for studying systematic uncertainties due to stellar continuum modeling and subtraction on a galaxy-by-galaxy basis; rather, our simulations quantify the detectability of emission lines based on a sufficiently large and representative grid of local continuum S/N and line significance values.

With the synthetic spectra in hand, we re-did our entire analysis on the simulated dataset in the same way as for the observations (running pPXF, subtracting the best-fit composite stellar template, fitting the residual emission lines, etc). Figure 7 shows our ability to recover the [O II], H$\beta$ and [O III] 5007 emission lines as a function of input integrated flux and degraded continuum S/N ratio. We find the sensible result that as the continuum S/N decreases, an emission line must become stronger for it to be significantly detected. In the highest continuum S/N regime, we can get down to about $3\times 10^{-16}$ erg s*-1* cm*-2* integrated flux for [O II], though fainter lines in the observed spectra are detected and can be visually confirmed. The integrated flux detectability thresholds for [O III] are a bit higher than those of [O II] probably for two reasons: (1) [O III] 5007 is in general fainter than [O II], and (2) in our red galaxies, the stellar continuum is much stronger in the [O III] rather than [O II] wavelength region. H$\beta$ has the highest completeness threshold because it is significantly fainter, is dependent on proper continuum subtraction, and is considered detected only if [O II] is also already detected.

The fact that our simulated [O II] integrated fluxes get down to about $3\times 10^{-16}$ erg s*-1* cm*-2* agrees well with the typical measurements from the observations. Considering NGC 1453 as a prototypical example, the integrated flux of [O II] in the central fiber is $\sim 7\times 10^{-15}$ erg s*-1* cm*-2*. In contrast, the integrated flux of [O II] in the fiber near $R_{\rm gas}$ is only $\sim 1.1\times 10^{-17}$ erg s*-1* cm*-2*. Upper limits on the integrated flux for [O II] non-detections are all less than $3\times 10^{-16}$ erg s*-1* cm*-2*, with a median value of $\sim 1.7\times 10^{-16}$ erg s*-1* cm*-2*. The integrated fluxes are recovered to within $5\%$ error (in the median) relative to the input value, and the line widths are recovered to within $3\%$ (in the median). We therefore assign all integrated flux measurements in the real data a conservative relative uncertainty of $10\%$.

Appendix B Two-dimensional Maps of MASSIVE Galaxies with Warm Ionized Gas

In this appendix, we show two-dimensional maps of various quantities for MASSIVE galaxies in which we detect warm ionized gas. The ordering of the galaxies in Figure 8 is the same as in Table 1: extended and rotating gas, extended but patchy gas, and unresolved compact gas (NGC 1453 is shown separately in Figure 2).