Quasar-mode Feedback in Nearby Type 1 Quasars: Ubiquitous Kiloparsec-scale Outflows and Correlations with Black Hole Properties

The present low-redshift sample ($z\lt 0.3$) was selected from the larger Quasar and ULIRG Evolution Study sample (QUEST; Schweitzer et al. 2006; Veilleux et al. 2006, 2009a, 2009b; Netzer et al. 2007). The QUEST sample consists of local ultraluminous infrared galaxies (ULIRGs) from the 1 Jy sample (Kim & Sanders 1998) and optically selected Palomar-Green quasars (Schmidt & Green 1983). The properties of the Gemini subsample of QUEST are listed in Table 1. The heterogeneous selection of this subsample was due to the combination of three archival data sets with seven observations designed to meet the current science goals. For the sources observed for this program, selection was based on a combination of observability, proximity (to maximize spatial resolution), and diversity in galaxy properties.

Table 1.  Sample Properties

Galaxy	Other Name	z	log(${L}_{\mathrm{bol}}/{L}_{\odot }$)	AGN Frac.	log(${M}_{\mathrm{BH}}/{M}_{\odot }$)	σ ($\mathrm{km}\,{{\rm{s}}}^{-1}$)	log[${P}_{\mathrm{radio}}/(W\,{\mathrm{Hz}}^{-1})$])
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
I Zw 1	PG0050+124, F00509+1225	0.0608	12.07	${0.93}_{-0.10}^{+0.07}$	${8.53}_{-0.5}^{+0.5}$	188 ± 36	22.7
F05189−2524	⋯	0.04288	12.22	${0.71}_{-0.21}^{+0.29}$	${8.32}_{-0.5}^{+0.5}$	201 ± 64	23.1
F07599+6508	⋯	0.1483	12.58	${0.75}_{-0.39}^{+0.25}$	${8.59}_{-0.5}^{+0.5}$	⋯	24.4
Mrk 231	F12540+5708	0.0422	12.60	${0.71}_{-0.07}^{+0.07}$	${8.58}_{-0.5}^{+0.5}$	233 ± 113	24.1
F13218+0552	⋯	0.2047	12.68	${0.83}_{-0.21}^{+0.17}$	${8.55}_{-0.5}^{+0.5}$	⋯	23.8
F13342+3932	⋯	0.1797	12.49	${0.69}_{-0.24}^{+0.23}$	${9.12}_{-0.5}^{+0.5}$	⋯	23.8
PG1411+442	⋯	0.0898	11.78	${1.00}_{-0.00}^{+0.00}$	${8.54}_{-0.48}^{+0.46}$	216 ± 31	22.3
PG1613+658	Mrk 876	0.12925	12.29	${0.82}_{-0.09}^{+0.11}$	${8.34}_{-0.52}^{+0.47}$	⋯	23.1
PG1700+518	F17002+5153	0.2902	13.12	${0.80}_{-0.19}^{+0.16}$	${8.79}_{-0.45}^{+0.45}$	⋯	24.8
F21219−1757	⋯	0.1127	12.17	${0.78}_{-0.13}^{+0.15}$	${8.61}_{-0.5}^{+0.5}$	121 ± 11	23.7

Note. Column 3: redshift measured from the current data. Column 4: Galaxy bolometric luminosity (Veilleux et al. 2009b). Column 5: AGN fraction, or fraction of the bolometric luminosity of the galaxy due to an AGN (Veilleux et al. 2009b). The error bars reflect the full range of values derived from six mid-infrared diagnostics. Column 6: black hole mass, taken from reverberation mapping (Bentz & Katz 2015) where possible or HST photometric measurements otherwise (Veilleux et al. 2009a). Errors from reverberation mapping have the uncertainty in the virial coefficient f added in quadrature ($\delta f=0.44;$ Woo et al. 2010). Photometric error bars are from the scatter in the M_BH–L_H relationship from which the masses were determined (Marconi & Hunt 2003). Column 7: stellar velocity dispersions (Dasyra et al. 2006a, 2007; Grier et al. 2013; Rothberg et al. 2013). The large error bars for F05189−2524 and Mrk 231 are due to conflicting optical and infrared measurements (Dasyra et al. 2006a; Rothberg et al. 2013); the mean of the two measurements is shown, and the error bar is half the difference between the measurements. Column 8: 1.4 or 1.5 GHz radio power from Barvainis & Lonsdale (1997) [PG1411+442], Carilli et al. (1998) [Mrk 231], Condon et al. (1990) [F05189−2524], Condon et al. (1998) [F07599+6508, F13218+0552, F13342+3932, F21219−1757], or Kukula et al. (1998) [I Zw 1, PG1613+658, PG1700+518].

Download table as:  ASCIITypeset image

All of the quasars in the sample are traditional Type 1s except for F05189−2524. However, this quasar has a broad-line region (BLR) observed in the near-infrared (Veilleux et al. 1999b). Three are broad-absorption-line (BAL) quasars: F07599+6508, Mrk 231, and PG1700+518. All three are also low-ionization BALs (LoBALs), with broad absorption observed in Mg ii (Smith et al. 1995; Turnshek et al. 1997; Veilleux et al. 2016) and Na i D (Boroson & Meyers 1992; Rupke et al. 2005a).

The average value of ${L}_{H}/{L}_{H}^{* }$ in the sample is 2.8 (Veilleux et al. 2006, 2009b), which corresponds to a halo mass of log(${M}_{\mathrm{halo}}/{M}_{\odot }$) ∼12.8 based on abundance matching (Moster et al. 2013) and a baryon fraction f_b = 0.17. Thus, the sample probes well above the halo mass at which the stellar-to-halo mass fraction peaks; i.e., the region in which AGN feedback is thought to be necessary in quenching star formation (e.g., Somerville et al. 2008).

The 10 quasars were observed with the Integral Field Unit (IFU) in the Gemini Multi-Object Spectrograph (GMOS; Allington-Smith et al. 2002; Hook et al. 2004) on the Gemini North and South telescopes. Details of each observation are listed in Table 2. All observations used the one-slit mode, and 9/10 used the B600 grating; this setup yielded an approximately constant resolution element of full width at half-maximum (FWHM) 1.5 Å. All observations were dithered to increase the FOV. The resulting total FOV was centered on the quasar in each case except for PG1700+518. The wavelength range for the typical setup included most strong emission lines (Hα, [O iii] 4959, 5007 Å, [N ii] 6548, 6583 Å, [S ii] 6717, 6731 Å, [O i] 6300, 6364 Å). For two archival setups (PG1700+518 and F21219−1757) and a new Mrk 231 observation, the grating was tilted blueward, including only [O ii] 3727, 3729 Å through [O iii] 4959, 5007 Å. All cubes but one (PG1700+518) also included the Na i D 5890, 5896 Å absorption doublet. The image quality of the reduced data cubes, as determined from a PSF fit (Section 2.3.1), was 05–09. (Separate exposures for a given quasar were convolved to match the worst seeing in the data.)

Table 2.  Observations

Galaxy	PID	Dates	Grating	texp	PSF	Range	PA	FOV
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
I Zw 1	GN-2015A-DD-9	2015 Jul 10	B600	4 × 1800 s	06	4610–7490 Å	225°	51 × 78
F05189−2524	GS-2011B-Q-64	2011 Dec 03, 13, 31	B600	6 × 1800 s	06	4560–7430 Å	0°	52 × 50
F07599+6508	GN-2007A-Q-12	2007 Apr 08	R831	3 × 1440 s	05	6520–8640 Å	180°	30 × 45
Mrk 231	GN-2013A-Q-51	2013 Jun 16, 17	B600	5 × 1800 s	08	5620–6950 Å	215°	63 × 75
F13218+0552	GN-2012A-Q-15	2012 Apr 02, 08	B600	8 × 1800 s	08	5340–8230 Å	315°	39 × 45
F13342+3932	GN-2012A-Q-15	2012 May 17	B600	8 × 1800 s	06	5340–8230 Å	313°	42 × 51
PG1411+442	GN-2015A-DD-9	2015 Jul 22, 23	B600	4 × 1800 s	08	4610–7480 Å	335°	39 × 45
PG1613+658	GN-2012A-Q-15	2012 Apr 22, 24	B600	6 × 1800 s	09	5340–8220 Å	135°	51 × 78
PG1700+518	GN-2003B-C-5	2003 Sep 23	B600	2 × 1800 s	07	4320–7190 Å	14°	30 × 42
F21219−1757	GS-2005B-Q-50	2005 Sep 09	B600	2 × 1800 s	09	4060–6860 Å	0°	27 × 45

Note. Column 2: program ID. Column 3: UT dates of observations. Column 4: grating. Column 5: number of exposures × length of each exposure. Column 6: FWHM of quasar PSF in the data cube, except for F05189−2524, in which the number is an estimate of the seeing from observing logs. Column 7: observed wavelength range of combined data cube. Column 8: east of north position angle of IFU during observations. Column 9: size of combined FOV.

Download table as:  ASCIITypeset image

The Gemini IRAF package (version 1.12) was used to reduce the data, supplemented by revision 114 of the IFUDR GMOS package and IFSRED (Rupke 2014a). Data cubes were visualized with QFitsView (Ott 2012). In addition to the steps described in Rupke & Veilleux (2013b, 2015), scattered light was removed from the flat fields and the data prior to extraction. Because of the strong nuclear point source in most cubes, scattered-light removal was crucial for accurate continuum and emission-line fitting.

Scattered light could not be completely removed from two systems: I Zw 1 and PG1411+442. The dithering strategy adopted in these cases put the quasar near the edge of the FOV along the short axis of the IFU. Consequently, the quasar light was concentrated near either the edge or center of fiber blocks in the raw data. Because scattered-light subtraction relies on extrapolating from interblock regions, the off-centered quasar light in the blocks meant that this extrapolation was imperfect. The result is an undersubtraction of the scattered light and is visible as a vertical stripe emerging from the point source in the reconstructed host galaxy image (Section 2.3.1).

Each final data cube was rebinned to a spaxel size of 03 × 03, which balanced the sensitivity with sampling of the seeing disk (the sampling was Nyquist or better in almost every case).

The spectrum in each spaxel was modeled with IFSFIT (Rupke 2014b), which utilizes PPXF (Cappellari 2012) for starlight fitting and MPFIT (Markwardt 2012) for other continuum and emission-line fitting. In brief, IFSFIT masks emission-line regions, fits the continuum, and then simultaneously fits all emission lines in the continuum-subtracted spectrum. Subsequently, interstellar absorption lines are modeled.

2.3.1. Separation of Quasar and Host Galaxy Light

In most of this sample, the luminous AGN, or quasar, outshines the host galaxy. HST H-band data show that the quasar-light fraction is 0.45–0.89 (Veilleux et al. 2006, 2009a). Heavy nuclear obscuration in F05189−2524, and to a lesser degree in F13218+0552 and F13342+3932, lower this fraction in the optical. At HST resolution, the quasar contributes 45% of the H-band light of F05189−2524 but <1% of its I-band light (Kim et al. 2013), consistent with the appearance of a BLR only in the near-infrared (Veilleux et al. 1999b). On visual inspection of optical HST images (Figure 1), F13218+0552 and F13342+3932 show a small quasar-light fraction, while six other quasars dominate their hosts. The final quasar, F21219−1757, has no optical HST image.

Zoom In Zoom Out Reset image size

View extended figure (10 panels)

To visually compare HST and IFS data, the cubes were summed over wavelengths that overlap with the HST filter range (Figure 1). For the three systems with a small quasar-light fraction in the optical, the HST images were also Gaussian smoothed to match the seeing.

In 9/10 cases, the peak IFS spaxel in 01 sampling represents a pure quasar spectrum to very good approximation (Figure 2). For F05189−2524, the peak spaxel traces an unresolved point source at optical ground-based resolution (Veilleux et al. 2002; 47% contribution to total light), but it is a combination of quasar light and starlight. It has very strong, high-ionization emission lines from the NLR that are either spatially unresolved or dominated by an unresolved component (including [O iii] 4959, 5007 Å, He ii 4686 Å, [Fe vii] 5159 Å, [Fe vii] 5721 Å, [Fe vii] 6087 Å, and [Fe x] 6375 Å; Farrah et al. 2005).

Zoom In Zoom Out Reset image size

A multistep process is adopted to fit the quasar and host galaxy in each spaxel. This iterative approach is necessary for three reasons.

First, subtracting the starlight accurately is crucial for properly constraining the underlying interstellar component of the Na i D absorption feature. It also produces more accurate recombination-line fluxes by properly accounting for Balmer absorption. In dusty, infrared-luminous galaxies where Na i D is strong, the stellar continuum subtraction is much less important because interstellar absorption dominates the feature. However, in many of these (less dusty) quasars, the stellar and interstellar Na i D absorption are comparable in strength.

Second, simultaneously fitting the quasar, starlight, and emission lines was technically infeasible. The starlight must be fit in PPXF (to allow the velocity and dispersion of the stellar spectrum to be free parameters). In principle, the quasar spectrum could also be input into PPXF as a sky spectrum, since the sky spectrum remains unconvolved with the fitted stellar velocity dispersion. However, PPXF uses Legendre polynomials to correct the stellar fit for imperfectly calibrated data. Linear combinations of polynomials are too flexible because they have the freedom to produce negative continuum values. This is especially a problem in spaxels where the host galaxy contribution is uncertain due to domination by quasar light or the total flux is very low.

Third, there is a range of sensitivity in this data set. This is partly by design: galaxies were chosen over a range of redshifts. It is also by accident: galaxies that were targeted for this program were observed ∼2× longer than galaxies observed for different purposes and retrieved from the archive. There is also intrinsic variation in the strength of the host galaxy starlight compared to the quasar. Galaxies with a very high signal-to-noise ratio (S/N) in the stellar continuum (e.g., F05189−2524 and Mrk 231) can be fit using multiple stellar templates. However, the composition of the stellar population could not be constrained on a spaxel-by-spaxel basis in most of the sample. In one system (F07599+6508), the stellar population could not be constrained even in the spatially integrated spectrum.

A description of the fitting process follows, and it is illustrated in Figure 3. As discussed above, all fitting is performed with IFSFIT, which calls PPXF when fitting starlight spectra with stellar population synthesis (SPS) models. In the discussion below, "total" refers to a complete spectrum with all physical components.

• 1.  
    Fit total spectra with quasar + exponential starlight model + emission lines. Each spaxel n is fit with the sum of a scaled quasar spectrum, a starlight model, and emission lines:

  $$\begin{eqnarray}&&{I}_{\mathrm{total}}^{n}={I}_{\mathrm{quasar}}^{n}+{I}_{\mathrm{starlight},\,\exp .\ \mathrm{model}}^{n}+{I}_{\mathrm{emission}}^{n}.\end{eqnarray} \tag{ 1 }$$

  First, emission lines are masked; the quasar and starlight continuum components are then fit simultaneously using MPFIT; the continuum is subtracted from the data; and, finally, the emission lines are fit with MPFIT. The starlight model and quasar scaling in each spaxel are each the sum of four exponentials,

  $$\begin{eqnarray}&&{I}_{\mathrm{total}}^{n}=\displaystyle \sum _{{\rm{i}}=1}^{4}{I}_{{\rm{i}}}^{n}{I}_{\mathrm{quasar}}^{0}+\displaystyle \sum _{{\rm{j}}=1}^{4}{I}_{{\rm{j}}}^{n}+{I}_{\mathrm{emission}}^{n},\end{eqnarray} \tag{ 2 }$$

  where ${I}_{\mathrm{quasar}}^{0}$ is the nuclear quasar spectrum (Figure 2); the exponentials are

  $$\begin{eqnarray}&&{I}_{1}^{n}={a}_{1}^{n}{e}^{-{b}_{1}^{n}\bar{\lambda }},\end{eqnarray} \tag{ 3 }$$

  $$\begin{eqnarray}&&{I}_{2}^{n}={a}_{2}^{n}{e}^{-{b}_{2}^{n}(1-\bar{\lambda })},\end{eqnarray} \tag{ 4 }$$

  $$\begin{eqnarray}&&{I}_{3}^{n}={a}_{3}^{n}(1-{e}^{-{b}_{3}^{n}\bar{\lambda }}),\,\,\,\mathrm{and}\end{eqnarray} \tag{ 5 }$$

  $$\begin{eqnarray}&&{I}_{4}^{n}={a}_{4}^{n}(1-{e}^{-{b}_{4}^{n}(1-\bar{\lambda })});\end{eqnarray} \tag{ 6 }$$

  ${a}_{{\rm{i}}}^{n}\geqslant 0;$ ${b}_{{\rm{i}}}^{n}\geqslant 0;$ $\bar{\lambda }\equiv \tfrac{\lambda -{\lambda }_{\min }}{{\lambda }_{\max }-{\lambda }_{\min }};$ and $[{\lambda }_{\min }$,${\lambda }_{\max }]$ is the fit range. Exponentials are used because they are monotonic and can be made positive-definite. The four exponentials allow for all combinations of concave/convex and monotonically increasing/decreasing. Physically, the quasar multipliers correct for spatially varying errors due to imperfect scattered-light removal (or other calibration effects) and for the ${\lambda }^{-0.2}$ dependence of PSF size on atmospheric seeing. For a Gaussian PSF, the latter effect produces an exponential of the form ${e}^{-C{\lambda }^{0.8}}$, which is similar to I₁.
• 2.  
   Calculate starlight-only spectra. The quasar and emission-line components are subtracted from the total spectrum in each spaxel. The result is a data spectrum including only starlight:

  $$\begin{eqnarray}&&{I}_{\mathrm{starlight}}^{n}={I}_{\mathrm{total}}^{n}-{I}_{\mathrm{quasar}}^{n}-{I}_{\mathrm{emission}}^{n}.\end{eqnarray} \tag{ 7 }$$
• 3.  
    Spatially integrate spectra. The starlight and total spectra are spatially integrated:

  $$\begin{eqnarray}&&{I}_{\mathrm{starlight}}=\displaystyle \sum _{n}{I}_{\mathrm{starlight}}^{n}\,\mathrm{and}\end{eqnarray} \tag{ 8 }$$

  $$\begin{eqnarray}&&{I}_{\mathrm{total}}=\displaystyle \sum _{n}{I}_{\mathrm{total}}^{n}.\end{eqnarray} \tag{ 9 }$$
• 4.  
    Fit spatially integrated starlight spectrum. The spatially integrated starlight spectrum is fit with the González Delgado et al. (2005) high-resolution SPS models and Legendre polynomials P_l:

  $$\begin{eqnarray}{I}_{\mathrm{starlight}} & = & {I}_{\mathrm{starlight},\mathrm{SPS}\mathrm{model}}+{I}_{\mathrm{Leg}.\mathrm{poly}.}\\ & = & \displaystyle \sum _{t}{A}^{t}{I}_{\mathrm{starlight},\ \mathrm{SPS}\ \mathrm{model}}^{t}+\displaystyle \sum _{{\rm{l}}\,=\,0}^{{{\rm{l}}}_{\max }}{B}_{{\rm{l}}}{P}_{{\rm{l}}}.\end{eqnarray} \tag{ 10 }$$

  The models vary in age, t, but assume solar metallicity (Rupke et al. 2008). A Legendre polynomial of order ${l}_{\max }=20\mbox{--}50$ (small compared to ∼6200 pixels in wavelength) is used in order to account for imperfections in calibration (e.g., scattered-light subtraction) or the PSF subtraction. These imperfections appear as very shallow, broad ripples in the continuum or higher-frequency edge effects.
• 5.  
   Fit spatially integrated total spectrum. This fit proceeds as in steps 1, 2, and 4 but involves only a single invocation of IFSFIT. After masking emission lines, IFSFIT fits the continuum component of the spatially integrated total spectrum with a scaled quasar-plus-exponential model for starlight as in step 1:

  $$\begin{eqnarray}&&{I}_{\mathrm{total},\mathrm{cont}.}={I}_{\mathrm{quasar}}+{I}_{\mathrm{starlight},\exp .\mathrm{model}}.\end{eqnarray} \tag{ 11 }$$

  IFSFIT then calculates the starlight spectrum as in step 2,

  $$\begin{eqnarray}&&{I}_{\mathrm{starlight}}={I}_{\mathrm{total},\mathrm{cont}.}-{I}_{\mathrm{quasar}},\end{eqnarray} \tag{ 12 }$$

  and fits the result with the sum of the SPS and Legendre components as in step 4. Finally, the emission lines are unmasked and fit. The endpoint of this fit is represented by

  $$\begin{eqnarray}&&{I}_{\mathrm{total}}={I}_{\mathrm{quasar}}+{I}_{\mathrm{starlight},\mathrm{SPS}\mathrm{model}}+{I}_{\mathrm{emission}}.\end{eqnarray} \tag{ 13 }$$
• 6.  
    Choose best-fit SPS model from steps 4 and 5 and sum over ages. The fit (from either of the previous two steps) that shows the minimum contribution from a polynomial is typically chosen. The best-fit SPS model is summed over ages t, ignoring the polynomial component:

  $$\begin{eqnarray}&&{I}_{\mathrm{starlight},\,\mathrm{SPS}\,\mathrm{model}}^{\mathrm{fixed}}=\displaystyle \sum _{t}{A}^{t}{I}_{\mathrm{starlight},\,\mathrm{SPS}\,\mathrm{model}}^{t}.\end{eqnarray} \tag{ 14 }$$

  The "fixed" superscript indicates that this starlight model will not vary in subsequent steps.
• 7.  
   Fit total spectra with quasar + SPS starlight model + emission lines. Each spaxel is fit as in step 1 but with the (fixed) SPS starlight model from step 6 plus a sum of Legendre polynomials replacing the exponential model:

  $$\begin{eqnarray}&&\,{I}_{\mathrm{total}}^{n}={I}_{\mathrm{quasar}}^{n}+{C}^{n}{I}_{\mathrm{starlight},\,\mathrm{SPS}\,\mathrm{model}}^{\mathrm{fixed}}+\displaystyle \sum _{{\rm{l}}=0}^{{{\rm{l}}}_{\max }}{B}_{{\rm{l}}}{P}_{{\rm{l}}}+{I}_{\mathrm{emission}}^{n}.\end{eqnarray} \tag{ 15 }$$

  As in step 5, IFSFIT is invoked only once, but there are intermediate steps (within SPS) before the SPS model is fit. In particular, the entire spectrum in each spaxel is first fit with the quasar, an exponential starlight model, and emission lines. A starlight spectrum is then calculated by removing the best-fit quasar and emission components from the total spectrum, and the starlight-only spectrum is refit with the SPS model. The fixed SPS model, ${I}_{\mathrm{starlight},\,\mathrm{SPS}\,\mathrm{model}}^{\mathrm{fixed}}$, is allowed to vary from spaxel to spaxel in velocity, dispersion, and normalization (Cⁿ). The stellar population is thus implicitly assumed to be unvarying across the relatively small GMOS FOV. The Legendre polynomial is similar to the one used in step 4.
• 8.  
   Calculate starlight-only spectra. Repeat step 2, using the new fit from step 7. The new starlight-only spectrum can differ from the first one because the first fit (step 1) does not properly account for stellar absorption lines, which in turn affects how the emission lines are removed from the total spectrum. The new fit of each spaxel (step 7) accounts for stellar absorption lines.
• 9.  
   Spatially integrate starlight spectrum. Repeat step 3 using the new starlight spectrum from step 8.
• 10.  
   Fit spatially integrated starlight spectrum. Repeat step 4 using the new spatially integrated starlight spectrum from step 9.
• 11.  
   Compare SPS fits to starlight-only spectrum. Compare the results of steps 4 and 10 by calculating the fractional contribution to the total light of starlight in four age groups. If there is a ≥10% difference in any group, begin again with the new SPS model, starting from step 6. Iterate until the current SPS model differs from the previous by <10% in each of the four age groups.

Zoom In Zoom Out Reset image size

For most objects, the first iteration of steps 1–10 resulted in convergence. For three objects, a second iteration was required.

The final SPS fits to the spatially integrated starlight spectra are shown in Figure 4, and the percent contributions (by flux in the fitted wavelength range) of the stellar populations in four age bins are given in Table 3. The use of the González Delgado et al. (2005) templates imposes an upper rest wavelength limit of 7000 Å to the fits. Reddening of the stellar population, which is almost certainly present, was not included in the fits because of the additional degrees of freedom this requires. Inclusion of reddening would likely result in younger best-fit ages. The luminosity-weighted age mixtures from these fits vary from a combination of ages ≲10 Myr with ages ≳100 Myr to populations dominated by ages ≳1 Gyr. The precision of these ages is likely coarse, given the single-metallicity, unreddened populations and the lack of ages 10–100 Myr in the best-fit populations. However, the important outcome of these fits is not an accurate set of ages but rather an accurate fit to the underlying stellar continuum, particularly in the Balmer lines and Na i D, at modest spectral resolution, which can be produced by more than one combination of stellar populations.

Zoom In Zoom Out Reset image size

Table 3.  Stellar Population Synthesis Fits

	Ages (Gyr)
Galaxy	≤0.01	0.01–0.1	0.1–1	>1	Poly.
(1)	(2)	(3)	(4)	(5)	(6)
I Zw 1	38	0	0	23	39
F05189−2524	38	1	36	3	22
F07599+6508	46	0	0	19	35
Mrk 231	59	0	13	11	17
F13218+0552	22	0	34	13	31
F13342+3932	0	0	0	39	51
PG1411+442	26	0	1	26	47
PG1613+658	0	0	0	58	42
PG1700+518	37	0	0	4	59
F21219−1757	6	0	42	19	33

Notes. Columns 2–5: percent contribution to total flux (in the wavelength range of the fits to integrated host galaxy spectra) from solar metallicity stellar populations of different ages (Section 2.3.1). Column 6: percent contribution from linear combination of Legendre polynomials.

Download table as:  ASCIITypeset image

For the two sources where unsubtracted scattered light was most problematic (I Zw 1 and PG1411+442; Section 2.2), the starlight fits also included a separate very broad component in Hα and Hβ that was limited in wavelength but allowed to vary in flux and width. The scattered-light undersubtraction was pronounced in these two broad emission-line features because of their intrinsic strength in the quasar spectrum and the particular way in which light was scattered across the 2D GMOS raw data. Including these free components, which were discarded after fitting, significantly improved the fits.

The wavelength-integrated results of the quasar/galaxy decomposition with IFSFIT are shown in Figure 5. Moffat fits to the quasar PSF that illustrate the effectiveness of the spectral method in recovering a smooth PSF are also shown in this figure. The exponential model that multiplies the quasar spectrum in each spaxel is shaped primarily by the I₁ term, likely due to its similarity to the theoretical dependence of the seeing-only PSF on wavelength. In every case, the host galaxy (starlight plus emission lines) is visible as an excess over the PSF in the radial profile. In most cases, specific features in the host galaxy image are easily recognizable in the HST images (Figure 1).

Zoom In Zoom Out Reset image size

View extended figure (10 panels)

The Legendre polynomial components of the fits represent featureless continuum contaminants due to incomplete scattered-light removal, imperfect PSF removal, or unfit stellar continuum components. The percent contributions to the total fitted light are in the range 20%–50%, with a median of 37%. For the purposes of this analysis, the polynomial component is assumed to be part of the host galaxy, an assumption that is borne out by the qualitative comparisons of the host galaxy reconstructions with HST images. However, in one system with strong residual scattered light, I Zw 1, the polynomial component is identifiable as scattered light and is separated from the host emission (Figure 5).

The stellar velocity maps resulting from the starlight fits are shown in Figure 6.

Zoom In Zoom Out Reset image size

The fidelity of the host galaxy reconstructions is at a minimum within 1–2 spaxels of the nucleus (i.e., within 0.5–1 spatial resolution elements). However, even within this radius, some stellar and line properties can still be extracted accurately (see Figure 6 and the emission- and absorption-line maps presented below), with this accuracy depending on the relative strength of the PSF, the details of scattered-light contamination, and the strength of the relevant line features.

2.3.2. Emission-line Fitting

Resolved line emission was detected and fitted in each quasar, though in F21219−1757, the emission was too weak for the line properties to be robustly determined. Every data cube also included unresolved line emission from the NLR; this unresolved emission was present in the point-source spectrum. The NLR emission was not removed from the point-source spectrum prior to using it as a continuum template, because at the physical resolutions of the target and instrument combinations, the NLR emission was dominated by unresolved components. Tests showed that removing NLR emission from the point-source spectrum before fitting yielded very different fits. This was due to the emission-line fit in many spaxels being dominated by the unresolved point-source component.

Lines were modeled with 1–2 Gaussian velocity components in each spaxel. All strong emission lines were fit with the same velocity components. Only one component was fit in PG1411+442 and F21219−1757, and up to two components were fit in the rest of the sample. The initial guesses for the number of components per spaxel in a given data cube or region of a cube were set by hand beforehand and were tailored for each cube or region to produce better by-eye fits. Model line profiles were convolved with the spectral resolution before fitting. Each component was automatically checked for significance after the fit and kept only if it exceeded a 3σ flux threshold in at least one strong line. After fitting, each "total" emission line flux was checked and set to zero if the significance of that line was below 2σ.

If two components were fit, the components were sorted into two maps by wavelength or velocity dispersion. For a given galaxy, the first component (c1) had either a longer wavelength or a smaller velocity dispersion than component 2 (c2). The choice of sorting by wavelength or velocity dispersion was applied to the entire cube. The choice depended on which method by eye produced the smoothest rotation signal in c1. If an outflow was present, it was assumed to arise in c2. However, c2 was not assumed to be an outflow in every case. In I Zw 1, c2 has very low velocities and does not show velocity gradients consistent with a coherent outflow. This component is assumed to arise from other effects (e.g., tidal motions).

This method does not assign c1 or c2 to the galaxy-rotation model or "other" model a priori. Rather, c1 and c2 are assigned to galaxy rotation or "other" based on what makes the best galaxy-rotation model. This predisposes the analysis to a good disk model rather than having it predisposed toward the "other" model. When an outflow component is present in c2, there is thus higher confidence that such a feature exists.

Ionized outflow is represented by blueshifted emission lines in c2. It is possible that redshifted emission may in some cases represent the back side of the outflow; however, in most cases, c2 is dominated by high-velocity blueshifted motions. The only clear-cut cases of redshifted outflowing emission are PG1613+658 and PG1700+518; in these cases, the outflow is taken as the entirety of c2.

Maps of strong emission lines for each galaxy in which these lines are detected are shown in Figure 7. Here [N ii] 6583 Å (or [O iii] 5007 Å in the cases in which [N ii] 6583 Å is not present in the wavelength range) was chosen as a strong representative line. Each figure displays maps of the total line flux; the flux in each component; the velocity dispersion of each component, ${\sigma }^{{\rm{c}}1}$ and ${\sigma }^{{\rm{c}}2};$ and the central velocity of each component, ${v}_{50 \% }$^c1 and ${v}_{50 \% }$^c2.

Zoom In Zoom Out Reset image size

View extended figure (9 panels)

For galaxies in which more than one emission-line component was fitted, two other panels appear: a map of the flux density at a representative blueshifted velocity and a map of the velocity that encompasses 98% of the cumulative velocity distribution in each spaxel, ${v}_{98 \% }$^c2. This convention matches that used in some previous IFS studies (Rupke & Veilleux 2011, 2013a, 2013b). The cumulative percentages are derived in blueshifted spaxels by integrating from the red side of the velocity distribution, such that ${v}_{98 \% }$^c2 is representative of a "maximum" blueshifted velocity. In redshifted spaxels, ${v}_{98 \% }$^c2 is correspondingly calculated by integrating from the blue side of the line.

Because they are forbidden lines, [N ii] and [O iii] may reflect outflow or AGN-related physics more than the recombination-line maps. The Balmer lines are also more susceptible to errors in stellar or quasar subtraction, since Balmer lines feature strongly in these templates. Nonetheless, the maps of the other lines are consistent with the maps shown.

Besides the fluxes and kinematics, other quantities calculated from the emission lines on a spaxel-by-spaxel basis are reddening $E(B-V)$ (Figure 8), electron density n_e (Figure 9), and line ratios (Figure 10). These three quantities are calculated after fluxes have been summed over both components (c1+c2). The measurements for individual components are noisier than those for the flux or velocity measurements, and summing over both components produces more spatially coherent maps. However, the maps of $E(B-V)$, n_e, and line ratios for individual components (c1 or c2) are generally consistent with the c1+c2 maps.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View extended figure (9 panels)

$E(B-V)$ is calculated from the Balmer decrement of the c1+c2 line flux (Figure 8). Case B at 10⁴ K and the Cardelli et al. (1989) extinction curve are assumed. These reddening maps are not used to correct the line maps in Figures 7 and 10 because in many spaxels the reddening is too noisy to be useful (these spaxels are not shown in Figure 8). This is largely due to limited sensitivity to Hβ in some data cubes.

Electron densities are calculated from the c1+c2 fluxes and the [S ii] 6717 Å/[S ii] 6731 Å ratios (Figure 9). The fit of [S ii] 6717 Å/[S ii] 6731 Å versus n_e from Sanders et al. (2016) is used, but it is assumed that the ratio is not sensitive to n_e below 10 cm⁻² and above 10,000 cm⁻² as it asymptotes, and calculated densities outside of this range are set to the lower or upper limit. Spatially averaged values of n_e in each quasar range from 50 to 410 cm⁻², with a median of 150 cm⁻².

Flux differences among lines are displayed in standard line ratio plots in Figure 10. Resolved line ratio maps are typically used to distinguish the dominant excitation mechanism (stellar photoionization, shock ionization, or AGN photoionization; Veilleux & Osterbrock 1987; Kewley et al. 2006), determine the fraction of each mechanism that is ionizing the gas in each spaxel (Rich et al. 2011; Davies et al. 2016a), and measure the physical parameters of the interstellar medium and radiation field by comparison to theoretical models (Dopita & Sutherland 1995; Groves et al. 2004; Dopita et al. 2006). These plots are discussed further in Section 3.1.

2.3.3. Absorption-line Fitting

Half of the quasars in this sample show blueshifted resonant absorption from the Na i D doublet in the PSF template spectra (Figure 2); F07599+6508 and Mrk 231 show broad absorption (Boroson & Meyers 1992; Veilleux et al. 2013b, 2016), while F05189−2524 (Rupke & Veilleux 2015), F13342+3932, and PG1411+442 (Hamann et al. 2017, in preparation) show narrow absorption. Many also show possible resonant broad emission, which requires high column densities, though it is difficult to distinguish it from nearby He i 5876 Å emission without careful fitting (Thompson 1991).

In this work, the unresolved absorption is fit out as part of the PSF removal and discarded. The goal is instead to measure the resolved, extended absorption and emission that appears as part of the host galaxy spectrum.

The Na i D doublet and He i emission line were fit using the methods of Rupke et al. (2005b, for Na i D absorption) and Rupke & Veilleux (2015, if Na i D emission is present). These techniques robustly separate optical depth and covering factor in absorption, even for a blended doublet and multiple velocity components, by assuming a Gaussian in optical depth. They also handle cases of combined resonant absorption and emission. In all cases but one, the Na i D and He i lines were fit simultaneously, but the He i was constrained where needed by fits to other emission lines. In one case with very high S/N, F05189−2524, the He i line was fit with other emission lines rather than as part of the Na i D fitting.

In both the integrated and the spaxel-by-spaxel fits, one Na i D component was fit in all cases but PG1411+442 (see below), Mrk 231, and F05189−2524. The latter two are the nearest galaxies in this sample and have deep, high-S/N integrations. The F05189−2524 fits are similar to those previously published (Rupke et al. 2005a; Rupke & Veilleux 2015), though here the binning is coarser and the spatially unresolved continuum has been removed to match the treatment of the rest of the sample. The integrated spectrum requires three components in absorption, while the spaxel-by-spaxel fit uses two components. The Mrk 231 fits use deeper data than previously published (Rupke & Veilleux 2011), which enables two-component fitting.

Two cases of the strongest and highest-S/N absorption (F05189−2524 and Mrk 231) also show redshifted Na i D emission. The F05189−2524 emission is analyzed in Rupke & Veilleux (2015); the deep Mrk 231 data reveal the new Na i D emission, but it is not discussed further here.

For all quasar hosts but one, the host continuum for these fits is assumed to be a combination of the stellar and polynomial components. In one case with strong scattered light, I Zw 1, the polynomial component was successfully identified as scattered light and the continuum is assumed to be the purely stellar component. (The effect of removing the polynomial component is to increase the resulting values of the covering factor C_f and, consequently, to increase the global flow rates; Section 3.2.)

Fits to Na i D in the integrated quasar host spectra (Figure 11) reveal extended interstellar absorption in every case except PG1700+518 (for which the Na i D feature is not in the observed wavelength range).

Zoom In Zoom Out Reset image size

In the spaxel-by-spaxel fits, spaxels with fitted equivalent width S/N > 3 were considered significant detections. Of the nine quasar hosts with extended interstellar absorption, eight show this absorption on a spaxel-by-spaxel basis (Figure 12). The maximum extent of this absorption ranges from 1 to 10 kpc, and the structure of the absorption is varied. Wind structure is discussed further in the next section. Most equivalent widths are in the range 1–10 Å.

Zoom In Zoom Out Reset image size

View extended figure (9 panels)

Figure 12 also shows the velocity maps of the neutral outflows. In cases where two outflowing absorption-line components are fit, the lines in each spaxel are first combined into a cumulative velocity distribution weighted by optical depth; ${v}_{50 \% }$ and ${v}_{98 \% }$ are calculated with respect to this total distribution. The maps also show the velocity of the peak optical depth, which differs from ${v}_{50 \% }$ for two or more components. Finally, Figure 12 compares the fitted equivalent width with an empirical estimate from simple integration (Rupke & Veilleux 2015); this is a sanity check on the quality of the fits. Generally, the fits tend to recover more signal and are a better estimate of the line properties (Rupke & Veilleux 2015).

In PG1411+442 appears a weak but significant high-velocity extended component (−1300 $\mathrm{km}\,{{\rm{s}}}^{-1}$). Tests of integrated spectra of varying apertures show that this absorption has a maximum extent of 1–2 kpc from the quasar (in projection). The low-velocity absorption extends to larger radii. Thus, two components are fit within five spaxels of the nucleus, and one component outside of this.

For PG1613+658, the spaxel-by-spaxel accounting for the line does not reveal that it extends in every direction around the nucleus. The extended absorption that is detected on a spaxel-by-spaxel basis is a strong velocity component that arises in front of a small satellite galaxy that is ∼5 kpc from the quasar in projection (Figure 1) and is redshifted at a velocity similar to the ionized gas disk of the quasar host at that location.

The detected absorption is generally blueshifted; spaxels with ${v}_{50 \% }\lt 0$ $\mathrm{km}\,{{\rm{s}}}^{-1}$ are assumed to be outflowing. There is some absorption that is near systemic or slightly redshifted (up to ∼200 $\mathrm{km}\,{{\rm{s}}}^{-1}$). In three systems (F07599+6508, PG1411+442, and PG1613+658), these make up a substantial fraction of the spaxels with detected absorption. In PG1613+658, this is due to the nearby companion. The same may be true in PG1411+442. In F07599+6508, the origin of the redshifted absorption is unclear.

Of the radio-quiet Type 1 quasars in this sample, all contain large-scale (maximum sizes 3–12 kpc) outflows of ionized and/or neutral gas, as revealed by the presence of broad, blueshifted absorption lines and broad, blueshifted emission lines. If we break the detections down by phase, 70% have ionized outflows, 80% have neutral outflows, and 50% have both. This detection rate contrasts with the ∼10% detection rate of ionized outflows in a previous IFS study of Type 1 quasars (Husemann et al. 2013, hereafter H13), but is consistent with the very high detection rate in nearby Type 2s (Harrison et al. 2014; McElroy et al. 2015; Bae et al. 2017).

For the two objects in common between H13 and this work (I Zw 1 and PG1700+518), there is agreement on the ionized outflow properties. The significant disagreement on the detection rate could thus result from differences in sample selection, as H13 selected lower-luminosity AGNs. Part, though not all, of the discrepancy could also result from our choice to probe both neutral and ionized gas outflows, rather than just ionized outflows, which could increase the detection rate.

A third possibility is that the different detection rate between our study and H13 results from data quality differences. The data in H13 were taken with the Potsdam Multi-Aperture Spectrophotometer (PMAS) on the 3.5 m Calar Alto telescope. The quoted emission-line surface brightness limit of H13 (at lower spatial and spectral resolution) is a factor of 3–4 worse than the best values achieved in the current study (∼5 × 10⁻¹⁷ erg s⁻¹ cm⁻² arcsec⁻²; Figure 7). The 05 sampling of PMAS is larger than the 02 GMOS hexagonal spaxels, and the typical seeing of H13 is about twice that of the current study. The spectral resolution in H13 varies from 2× to 4× coarser than that of the current study. These factors make separating a second (typically broad, faint) outflowing emission-line component from a disk component more challenging. The exception is when the outflow dominates the line emission. In PG1700+518, the outflow completely dominates the line emission in [O iii], and both studies classify it as a large-scale outflow.

The approximate projected outflow boundaries and orientations are shown in Figure 13. The average outflow position in each phase is calculated in projected Cartesian coordinates; the position weighted by ${v}_{50 \% }$², ${v}_{98 \% }$², and W_eq or flux is also calculated. These different weightings give a sense of the uncertainty in computing the true orientation of the outflow.

Zoom In Zoom Out Reset image size

Loosely collimated outflows along their host galaxy's minor axes are common on kpc scales in many starburst and active galaxies observed with IFS (e.g., Rupke & Veilleux 2013b). The observed outflows in these quasars are also typically oriented near their hosts' minor axes. Clear exceptions to this are the outflows in F13218+0552, the neutral outflow in F07599+6508, and the ionized outflows in PG1613+658 and PG1700+518. One expects some exceptions due to projection effects, which will enhance misalignments for systems that are nearly face-on. PG1613+658 and PG1700+518 have outflows misaligned with the minor axis for the reasons described below. The outflows in F13342+3932 would at first glance appear to be exceptions as well. However, the highest velocity neutral gas lies right along the minor axis at 1–2 kpc scales, and the largest ionized fluxes are found along the minor axis on the same scales in the opposite direction. This behavior is not captured in the outflow positional averages, as the outflow is spatially dominated by an apparently less collimated flow at larger scales.

The opening angle of an outflow (between 0° and 180°) quantifies its collimation; loosely collimated outflows have large opening angles. In the bipolar cone model, the opening angle is simply the projected two-dimensional angle between the sides of each cone. The broad footprints of the outflows in this sample (Figure 13) suggest loosely collimated flows with large opening angles (perhaps approaching 180°).

In PG1700+518, the outflow lies along a known radio jet. This jet-driven outflow and its molecular counterpart are discussed in greater detail in Runnoe et al. 2018. This quasar is the most distant and radio-luminous in the sample with z = 0.2905 and log $({P}_{1.4\mathrm{GHz}})=24.8$ W Hz⁻¹. This places it at the knee of the radio luminosity function of a nearby AGN (Condon et al. 2002) but two orders of magnitude below the knee of the $z\sim 0.3$ quasar 1.4 GHz luminosity function (Condon et al. 2013). It is one of two quasars in the sample with a known jet (Kukula et al. 1998; Yang et al. 2012). (The other is Mrk 231; see Morganti et al. 2016 and references therein.)

In PG1613+658, there is an extended NLR or ionization cone, presumably aligned with and collimated by the AGN torus, that is roughly along the galaxy major axis and traced by velocity c2. It has very high [O iii]/Hα values (Figure 10) indicative of AGN photoionization. The flux and bulk velocity of c2 peak about 2 kpc from the nucleus. Though the projected bulk motions (as traced by ${v}_{50 \% }$^c2) of this component are similar to the disk motions in c1, the large values of ${\sigma }^{c2}$ and ${v}_{98 \% }$^c2 are indicative of an outflow, and it is classified as such. This outflow may arise from the torus being misaligned with the host minor axis. The bulk of the AGN radiation is pointed into the disk, accelerating disk material (as in nearby Seyferts; Fischer et al. 2013). The lack of high bulk motions could mean the outflow is in the plane of the sky, and the decrease in ${v}_{50 \% }$^c2 at large radii may mimic NLR deceleration in lower-luminosity systems (Fischer et al. 2013).

There is also an increased velocity dispersion and shock-like line ratios (high [N ii]/Hα, [S ii]/Hα, and [O i]/Hα) along the galaxy minor axis in PG1613+658, as seen in c1. These latter effects have been used as tracers of outflows in other resolved studies (Lehnert & Heckman 1996; Ho et al. 2016). However, the dense velocity contours along the minor axis of this system, coupled with the >2 kpc PSF, mean that the disk needs to be carefully modeled to account for beam smearing (Davies et al. 2011). In the absence of this detailed modeling, we do not include these spaxels as part of the outflow.

Strikingly, in the five systems where both neutral and ionized outflows are present, the neutral and ionized outflows are often oriented opposite each other. This ionized/neutral asymmetry is not obvious in other Type 2 quasars (Rupke & Veilleux 2013b). Spatially resolved observations show that Na i D absorption correlates with optical dust obscuration (Rupke & Veilleux 2013b, 2015). An anticorrelation in the outflow between line-emitting gas and dusty gas could reflect asymmetric obscuration in the outflow. On one side, the obscuration along the line of sight is small, and the ionized outflow is prominent. On the other, the obscuration is higher, and the neutral outflow dominates. One possible cause is the underlying gas disk, though the obscuration correlated with the wind is thought to be in the outflow rather than the disk (Rupke & Veilleux 2015).

The line ratio maps in Figure 10 (calculated on a spaxel-by-spaxel basis using fluxes summed over velocity components) are not used to determine whether or not an outflow is present, and a quantitative comparison of how the outflows relate to line excitation is reserved for future work. However, a qualitative understanding of the excitation of the ionized gas provides an overview of possible radiation physics in these outflows. The maps in Figure 7 and in the top panels of Figure 10 often visually reveal correlations between outflows and higher line ratios in many systems. However, this correlation is not one-to-one; for instance, only the redshifted side of the outflow in PG1700+518 has high excitation. And in some cases, there is no obvious correlation (Mrk 231).

Line ratio diagrams, seen in the bottom panels of Figure 10, typically show patterns suggestive of stellar and AGN photoionization in a sequence mixing pure stellar photoionization and pure AGN photoionization. This mixing sequence, which moves from the lower left to the upper middle of these diagrams, is due to the fractional contribution of stellar and AGN photoionization varying from one part of the galaxy to another (e.g., Davies et al. 2016a). Stellar photoionization dominates in the lower left and AGN photoionization in the upper middle of these diagrams. There is also evidence of a stellar photoionization sequence in metallicity in at least two hosts (I Zw 1 and PG1411+442); this sequence is visible as a curving locus of points from the middle left to the middle bottom (e.g., Kewley et al. 2001). Shock ionization, a common feature of galactic winds (Veilleux & Rupke 2002; Sharp & Bland-Hawthorn 2010; Rich et al. 2011), is intermixed with AGN and stellar ionization in other systems (F05189−2524, F07599+6508, F13218+0552, and PG1613+658) and is dominant in at least one system (F05189−2524). Shock ionization dominates the right-hand side of these diagnostics (Dopita & Sutherland 1995; Rich et al. 2011), and mixing with other ionization mechanisms is revealed as loci of points connecting regions of stellar, AGN, and shock ionization. Finally, evidence of AGN photoionization variations due to changing ionization parameters (Scharwächter et al. 2011; Davies et al. 2016b) are evident in at least one quasar (F13342+3932). These changes are visible as loci of points moving from the upper middle toward the lower right.

The key physical quantities that measure the potential impact of an outflow on its host are velocity, size (maximum projected radius R), mass (M), momentum (p), energy (E), and the outflow rates of M, p, and E in each gas phase. In Table 4, we list the mean and maximum values of FWHM^c2, ${v}_{50 \% }$^c2, and ${v}_{98 \% }$^c2. The averages are taken over spaxels in the outflowing component(s). The other quantities are listed in Table 5.

Table 4.  Outflow Velocity Statistics

		FWHMc2		${v}_{50 \% }^{{\rm{c}}2}$		${v}_{98 \% }^{{\rm{c}}2}$
Galaxy	Phase	avg	max	avg	max	avg	max
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
I Zw 1	neutral	202	654	−120	−326	−297	−635
F05189-2524	neutral	624	1728	−560	−1031	−1165	−2569
⋯	ionized	592	1591	−423	−989	−926	−1598
F07599+6508	neutral	275	682	−96	−315	−337	−793
⋯	ionized	310	641	−120	−839	−384	−1327
Mrk 231	neutral	395	1122	−416	−1112	−801	−1784
⋯	ionized	768	1397	−672	−1036	−1324	−2003
F13218+0552	neutral	189	386	−47	−185	−213	−524
⋯	ionized	641	1725	−169	−578	−714	−1897
F13342+3932	neutral	322	1027	−104	−432	−387	−1262
⋯	ionized	386	882	−144	−418	−473	−1167
PG1411+442	neutral	909	2148	−486	−1274	−1044	−1780
PG1613+658	ionized	409	609	−109	−578	−457	−704
PG1700+518	ionized	895	2179	−570	−854	−1331	−2627
F21219-1757	neutral	239	646	−230	−487	−441	−700

Notes. Column 2: gas phase. Columns 3−4: mean and maximum (over all spaxels) of FWHM of velocity component 2 (c2) in km s^–1. Columns 5−6: mean and maximum values (over all spaxels) of c2 central velocity in km s^–1. For the neutral phase, all components are combined in each spaxel into an optical depth–weighted cumulative velocity distribution. The averages are taken over spaxels with ${v}_{50 \% }\lt 0$ $\mathrm{km}\,{{\rm{s}}}^{-1}$. For the ionized phase, only c2 is used; with the exceptions of PG1613+658 and PG1700+518, the averages are taken over spaxels with ${v}_{50 \% }^{{\rm{c}}2}\lt 0$ $\mathrm{km}\,{{\rm{s}}}^{-1}$. For PG1613+658 and PG1700+518, spaxels with ${v}_{50 \% }^{{\rm{c}}2}\gt 0$ $\mathrm{km}\,{{\rm{s}}}^{-1}$ are included, but the velocities are multiplied by −1 so that the red and blue outflow components can be combined into a single average. Columns 7−8: mean and maximum of ${v}_{98 \% }\equiv {v}_{50 \% }-2\sigma$, where σ is the average velocity above and below ${v}_{50 \% }$ that encompasses 34% of the cumulative velocity distribution: the area between ${v}_{50 \% }\pm \sigma$ contains 68% of the distribution. As for ${v}_{50 \% }$, the redshifted ionized gas spaxels in PG1700+518 have their distribution flipped in sign and combined with the blueshifted components.

Download table as:  ASCIITypeset image

Table 5.  Mass, Momentum, and Energy

Galaxy	Phase	Robs	rwind	log[${f}_{{\rm{H}}\alpha }$/	log(M/	log[$({dM}/{dt})$/	log[p/	log[$(c\,{dp}/{dt})$/	log(E/	log[$({dE}/{dt})$/
		(kpc)	(kpc)	(erg s−1 cm−2)]	${M}_{\odot }$)	(${M}_{\odot }\,{\mathrm{yr}}^{-1}$)]	(dyne s)]	${L}_{\odot }$]	erg)	(erg s−1)]
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
I Zw 1	neutral	5.1	5.5	⋯	${8.77}_{-0.03}^{+0.06}$	${1.38}_{-0.06}^{+0.10}$	${49.41}_{-0.06}^{+0.10}$	${11.61}_{-0.09}^{+0.15}$	${56.79}_{-0.09}^{+0.12}$	${42.14}_{-0.12}^{+0.18}$
F05189-2524	neutral	2.7	3.0	⋯	${8.58}_{-0.02}^{+0.03}$	${1.98}_{-0.03}^{+0.05}$	${49.75}_{-0.03}^{+0.05}$	${12.66}_{-0.03}^{+0.06}$	${57.52}_{-0.03}^{+0.06}$	${43.58}_{-0.04}^{+0.07}$
⋯	ionized	2.8	3.0	−13.21	${7.36}_{-0.14}^{+0.06}$	${0.40}_{-0.14}^{+0.07}$	${48.16}_{-0.14}^{+0.07}$	${10.78}_{-0.14}^{+0.07}$	${55.55}_{-0.14}^{+0.07}$	${41.32}_{-0.13}^{+0.07}$
⋯	molecular	⋯	$\lt 1$	⋯	${8.28}_{-0.38}^{+0.38}$	${2.43}_{-0.29}^{+0.03}$	⋯	${12.72}_{-0.25}^{+0.08}$	⋯	${43.20}_{-0.25}^{+0.10}$
⋯	total	⋯	⋯	⋯	${8.77}_{-0.09}^{+0.16}$	${2.57}_{-0.19}^{+0.03}$	${49.76}_{-0.03}^{+0.05}$	${13.00}_{-0.11}^{+0.05}$	${57.53}_{-0.03}^{+0.06}$	${43.73}_{-0.07}^{+0.06}$
F07599+6508	neutral	7.5	8.1	⋯	${9.27}_{-0.09}^{+0.09}$	${1.27}_{-0.12}^{+0.39}$	${49.46}_{-0.12}^{+0.39}$	${11.20}_{-0.12}^{+0.19}$	${57.24}_{-0.24}^{+0.19}$	${41.81}_{-0.17}^{+0.44}$
⋯	ionized	7.0	8.1	−15.23	${7.50}_{-0.15}^{+0.06}$	$-{0.64}_{-0.19}^{+0.06}$	${47.56}_{-0.19}^{+0.06}$	${9.45}_{-0.33}^{+0.10}$	${54.65}_{-0.33}^{+0.10}$	${39.95}_{-0.32}^{+0.12}$
⋯	total	⋯	⋯	⋯	${9.27}_{-0.09}^{+0.09}$	${1.27}_{-0.12}^{+0.38}$	${49.47}_{-0.12}^{+0.38}$	${11.21}_{-0.12}^{+0.19}$	${57.24}_{-0.24}^{+0.19}$	${41.82}_{-0.17}^{+0.44}$
Mrk 231	neutral	3.4	3.6	⋯	${8.84}_{-0.02}^{+0.02}$	${2.00}_{-0.02}^{+0.03}$	${49.85}_{-0.02}^{+0.03}$	${12.54}_{-0.03}^{+0.05}$	${57.46}_{-0.03}^{+0.05}$	${43.35}_{-0.04}^{+0.06}$
⋯	ionized	2.9	1.3	−14.71	${5.88}_{-0.04}^{+0.05}$	$-{0.30}_{-0.05}^{+0.06}$	${47.10}_{-0.05}^{+0.06}$	${10.37}_{-0.05}^{+0.06}$	${54.78}_{-0.05}^{+0.06}$	${41.18}_{-0.05}^{+0.06}$
⋯	molecular	⋯	$\lt 1$	⋯	${9.04}_{-0.03}^{+0.07}$	${3.04}_{-0.06}^{+0.02}$	⋯	${13.34}_{-0.18}^{+0.06}$	⋯	${43.92}_{-0.16}^{+0.11}$
⋯	total	⋯	⋯	⋯	${9.25}_{-0.02}^{+0.05}$	${3.08}_{-0.06}^{+0.01}$	${49.85}_{-0.02}^{+0.03}$	${13.40}_{-0.15}^{+0.05}$	${57.47}_{-0.03}^{+0.05}$	${44.02}_{-0.12}^{+0.09}$
F13218+0552	neutral	11.3	12.0	⋯	${8.68}_{-0.39}^{+0.26}$	${0.93}_{-1.09}^{+0.35}$	${49.30}_{-1.09}^{+0.35}$	${11.08}_{-1.52}^{+0.39}$	${56.60}_{-1.01}^{+0.38}$	${41.48}_{-1.39}^{+0.41}$
⋯	ionized	11.5	12.0	−13.76	${9.21}_{-0.11}^{+0.04}$	${1.03}_{-0.09}^{+0.04}$	${49.40}_{-0.09}^{+0.04}$	${10.80}_{-0.07}^{+0.04}$	${56.18}_{-0.07}^{+0.04}$	${40.87}_{-0.06}^{+0.05}$
⋯	total	⋯	⋯	⋯	${9.32}_{-0.11}^{+0.08}$	${1.29}_{-0.24}^{+0.19}$	${49.65}_{-0.24}^{+0.19}$	${11.27}_{-0.44}^{+0.29}$	${56.74}_{-0.46}^{+0.30}$	${41.58}_{-0.64}^{+0.36}$
F13342+3932	neutral	9.2	10.7	⋯	${8.50}_{-0.07}^{+0.17}$	${0.71}_{-0.15}^{+0.22}$	${49.03}_{-0.15}^{+0.22}$	${10.91}_{-0.22}^{+0.24}$	${56.80}_{-0.32}^{+0.24}$	${41.83}_{-0.38}^{+0.27}$
⋯	ionized	7.9	10.7	−13.15	${9.17}_{-0.04}^{+0.03}$	${1.37}_{-0.04}^{+0.03}$	${49.69}_{-0.04}^{+0.03}$	${11.38}_{-0.04}^{+0.03}$	${56.71}_{-0.04}^{+0.03}$	${41.54}_{-0.05}^{+0.03}$
⋯	total	⋯	⋯	⋯	${9.25}_{-0.04}^{+0.04}$	${1.46}_{-0.04}^{+0.05}$	${49.78}_{-0.04}^{+0.05}$	${11.51}_{-0.06}^{+0.08}$	${57.06}_{-0.15}^{+0.15}$	${42.01}_{-0.21}^{+0.19}$
PG1411+442	neutral	2.8	3.2	⋯	${8.42}_{-0.12}^{+0.11}$	${1.68}_{-0.11}^{+0.15}$	${49.47}_{-0.11}^{+0.15}$	${12.57}_{-0.13}^{+0.20}$	${57.41}_{-0.13}^{+0.19}$	${43.65}_{-0.17}^{+0.24}$
PG1613+658	ionized	7.5	8.0	−14.11	${7.73}_{-0.30}^{+0.15}$	${0.04}_{-0.36}^{+0.16}$	${48.23}_{-0.36}^{+0.16}$	${9.98}_{-0.38}^{+0.17}$	${55.18}_{-0.38}^{+0.17}$	${40.11}_{-0.32}^{+0.15}$
PG1700+518	ionized	9.7	10.0	−12.89	${9.58}_{-0.09}^{+0.14}$	${2.36}_{-0.10}^{+0.14}$	${50.65}_{-0.10}^{+0.14}$	${12.87}_{-0.10}^{+0.14}$	${58.16}_{-0.10}^{+0.14}$	${43.53}_{-0.11}^{+0.15}$
F21219-1757	neutral	5.5	5.9	⋯	${8.79}_{-0.06}^{+0.09}$	${1.56}_{-0.12}^{+0.12}$	${49.62}_{-0.12}^{+0.12}$	${11.87}_{-0.19}^{+0.16}$	${57.04}_{-0.19}^{+0.16}$	${42.41}_{-0.25}^{+0.20}$

Note. Column 2: gas phase. Columns 3–4: maximum observed wind radius (in projection) and assumed wind radius, in kpc. Column 5: logarithm of the Hα flux in the wind, corrected for reddening using the Balmer decrement. For F07599+6508, the extincted flux is listed because the spectrum contains only one Balmer line. Columns 6–11: logarithms of the wind mass, momentum, energy, and their outflow rates, computed using a time-averaged thin-shell model that depends inversely on the assumed shell radius (Section 3.2; Rupke et al. 2005c; Shih & Rupke 2010; Rupke & Veilleux 2013b). The ionized gas values depend on electron density as ${n}_{e}^{-1}$. For quasars with [S ii] 6717, 6731 Å measurements, n_e is calculated on a spaxel-by-spaxel basis; for the others, n_e = 100 cm⁻² is assumed (Section 2.3.2). The neutral gas values depend on ionization fraction, Na abundance, metallicity, and Na dust depletion (Section 3.2 and Rupke et al. 2005b). Molecular gas values are from González-Alfonso et al. (2017).

Download table as:  ASCIITypeset image

Mass, momentum, and energy and their outflow rates are calculated using a time-averaged thin-shell model (the single radius free wind (SRFW) model of Rupke & Veilleux 2013b) in which these quantities depend inversely on the assumed shell radius (Rupke et al. 2005c; Shih & Rupke 2010; Rupke & Veilleux 2013b). The thin-shell model is a fiducial model for measuring such quantities. For spatially resolved data, it has the advantage of being able to associate each point in the plane of the sky with a unique three-dimensional radius, assuming that only the near side of the wind is observed. Real winds are almost certainly more complex. However, thick-shell models (e.g., Rupke et al. 2005c) require knowledge of the wind's density and/or radial distribution, which in turn requires more detailed spatial information than is available in the current data.

Some well-resolved winds are able to be fit with biconical or bipolar bubble models (e.g., Veilleux et al. 1994; Fischer et al. 2013; Rupke & Veilleux 2013b). Such models tend to better account for wind projection and often lead to higher values of mass, momentum, and energy (Rupke & Veilleux 2013b). There are indications of bipolar flows in at least two objects in our sample (PG1613+658 and PG1700+518), and there are further indications of minor-axis collimation in others (Section 3.1). However, the canonical sign of bipolar shells, line splitting, is noticeably absent. Detailed model fits of bipolar bubbles or cones would thus need to fully explore the parameter space of possible geometries to address parameter degeneracies. For simplicity and consistency with previous studies and within the sample, the thin shell is used exclusively here. Fits to bipolar models are reserved for future work, possibly using data of higher spatial resolution.

The measurements of M, p, and E may be lower limits to the true values due to the fact that any spatially unresolved outflow has been removed from the data. Thus, outflows at scales <0.5 kpc in the nearest systems, approaching 2–3 kpc in the most distant systems, are not considered here.

The maximum observed radius of each wind in two-dimensional projection (R_obs) is measured in each phase. The observed range of R_obs is 3–12 kpc. The (single) three-dimensional radius adopted for the SRFW model, r_wind, is a free parameter in the model. Adopting ${r}_{\mathrm{wind}}={R}_{\mathrm{obs}}$ exactly results in edge effects: spaxels with projected radii $R={R}_{\mathrm{obs}}$ yield formally infinite velocities when they are deprojected. The radius adopted is thus slightly larger (by ∼0.5 spaxel) than the observed radius. In most cases, ${r}_{\mathrm{wind}}^{\mathrm{neutral}}={r}_{\mathrm{wind}}^{\mathrm{ionized}}$ is adopted. One exception is Mrk 231; while there is blueshifted, AGN-photoionized emission at large radii, the line widths of this $R\sim 3\,\mathrm{kpc}$ emission are very narrow and lie atop the southern arc of star formation in this galaxy. It is thus assumed that these spaxels represent AGN-photoionized gas that is forming stars, rather than an outflow.

The ionized gas flux in each wind (except in F07599+6508, for which there is only one Balmer line in the present data) is corrected using the extinction computed from the Balmer decrement and listed in Table 5. As described in Section 2.3.2, Case B at 10⁴ K and the Cardelli et al. (1989) extinction curve are assumed. This correction (to the flux in c2) is done on a spaxel-by-spaxel basis where possible and uses the extinction computed from the combined flux in both components (c1+c2; Figure 8). Spaxels with uncertain extinction values are corrected using the median extinction value in the host.

The outflowing ionized gas mass in each spaxel depends inversely on electron density n_e. The electron density measured in each spaxel (calculated as described in Section 2.3.2 and shown in Figure 9) is used to calculate the gas mass. In the few spaxels where the [S ii] fluxes are too uncertain, a value of 100 cm⁻² is used (near the measured sample median; Section 2.3.2). This value is also assumed for the two systems without n_e measurements (Mrk 231 and PG1700+518).

The neutral gas masses are proportional to the measured column density of hydrogen (Rupke et al. 2005c). The hydrogen column density, in turn, is related to the measured Na column through the ionization fraction, the Na abundance relative to solar and the galaxy's metallicity, and the Na dust depletion (Rupke et al. 2005b). Solar metallicity is assumed (Rupke et al. 2008), and the values for ionization fraction, Na abundance, and Na dust depletion are taken from Rupke et al. (2005b).

Five quasar hosts have detections of both neutral and ionized outflows (though one, F07599+6508, does not have extinction-corrected outflow properties). Two hosts, F05189−2524 and Mrk 231, also have molecular phase outflow measurements of critical wind properties. The latest measurements from the modeling of multiple OH absorption lines observed with Herschel are adopted (González-Alfonso et al. 2017). These molecular outflows have predicted sizes under <1 kpc. Their properties are listed alongside the neutral and ionized properties in Tables 5 and 6.

The four hosts with extinction-corrected ionized outflow measurements and measurements in multiple gas phases show that the ionized gas phase dominates the mass of the wind in 2/4 cases but is negligible in the other two. However, the neutral and molecular phases typically contribute most to the energetics. This illustrates the importance of multiphase observations of winds, even in AGN outflows that are often presumed to be dominated by ionized gas.

Observational studies of outflow properties are useful for informing sub-grid prescriptions in cosmological simulations of galaxy formation and evolution if observed outflow properties correlate with host, AGN, or black hole properties. Previous studies have focused on how quasar-mode outflows depend on host or AGN properties (e.g., Sturm et al. 2011; Rupke & Veilleux 2013b; Spoon et al. 2013; Veilleux et al. 2013a; Cicone et al. 2014; Harrison et al. 2014; McElroy et al. 2015; Stone et al. 2016; Fiore et al. 2017; Sun et al. 2017). Since this sample of nearby luminous AGNs has substantial ancillary data, and thus useful constraints on the central black hole (Table 1), this paper focuses on the dependence of outflow properties on black hole properties.

To probe this dependence across a significant dynamic range, however, a larger range of black hole mass is required. Table 7 lists the outflow and host galaxy properties in 13 nearby Seyfert galaxies with spatially resolved outflows. The outflow properties are mostly based on ionized gas measurements of the NLR, except for those in NGC 1068, NGC 1266, and Circinus. The measurements of NGC 1266 and Circinus are of the molecular gas, which dominates the outflow mass and energy budgets. The measurements of NGC 1068 combine ionized and molecular gas. These NLR outflows have scales of 0.1 to a few kpc.

Table 7.  Large-scale Seyfert Galaxy Outflows

Galaxy	Type	dM/dt	log[$({dE}/{dt})/$	log(${M}_{\mathrm{BH}}/$	log[${L}_{\mathrm{AGN}}/$	σ	References
		(${M}_{\odot }\,{\mathrm{yr}}^{-1}$)	(erg s−1)]	${M}_{\odot }$)	(erg s−1)]	($\mathrm{km}\,{{\rm{s}}}^{-1}$)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
NGC 1068	2	76	42.87	${6.934}_{-0.015}^{+0.015}$	44.95 ± 0.50	176 ± 9	(b), (C), (D), (ff), (3), (7)
NGC 1266	2	110	41.04	${6.738}_{-0.441}^{+0.441}$	43.53 ± 0.30	94 ± 5	(a), (A), (dd), (1)
Mrk 79	1	10	40.68	${7.612}_{-0.565}^{+0.445}$	44.65 ± 0.18	129 ± 10	(d), (F), (gg), (2), (5)
NGC 2992	2	168	42.54	${7.703}_{-0.443}^{+0.442}$	44.48 ± 0.50	160 ± 17	(d), (D), (dd), (5), (8)
NGC 3783	1.5	$3.5\ \$	40.99	${7.367}_{-0.492}^{+0.442}$	43.70 ± 0.25	110 ± 19	(c), (D), (gg), (2), (5)
NGC 4151	1.5	3	41.63	${7.553}_{-0.456}^{+0.441}$	43.28 ± 0.28	96 ± 10	(b), (B), (aa), (2), (5)
NGC 4253	1	0.78	⋯	${6.822}_{-0.444}^{+0.443}$	43.66 ± 0.21	93 ± 32	(c), (G), (bb), (2), (5)
NGC 4388	2	0.02	37.30	${6.929}_{-0.010}^{+0.010}$	44.68 ± 0.52	99 ± 9	(b), (I), (ee), (3), (7)
Circinus	1	8.0	40.00	${6.230}_{-0.084}^{+0.092}$	43.30 ± 0.30	149 ± 18	(c), (J), (cc), (4)
NGC 5548	1.5	9.5	41.47	${7.718}_{-0.440}^{+0.440}$	44.29 ± 0.35	193 ± 11	(b), (H), (bb), (2), (5)
NGC 6814	1.5	10.5	41.05	${7.163}_{-0.493}^{+0.442}$	43.24 ± 0.34	106 ± 17	(c), (D), (bb), (2), (5)
NGC 7469	1	5.6	40.92	${6.979}_{-0.460}^{+0.441}$	44.62 ± 0.20	130 ± 6	(b), (D), (gg), (2), (5)
NGC 7582	1	3.4	⋯	${7.740}_{-0.097}^{+0.111}$	43.30 ± 0.30	148 ± 19	(b), (E), (hh), (6)

Notes. Column 2: Seyfert type. Columns 3–4: mass and energy outflow rates. Error bars are not consistently available for these measurements; we assume 0.3 dex. Most estimates are based on ionized gas measurements. The exceptions are NGC 1068 (molecular and ionized gas measurements have been summed), NGC 1266 (only a molecular gas measurement is available), and Circinus (the molecular gas value is much greater than the ionized gas value). In most cases, the ionized values are calculated geometrically and depend on the product of the electron density n_e and filling factor f, which is typically assumed to be n_e f = 5 cm⁻³. For consistency, we have normalized the geometric estimates to this value, except in the case of NGC 5548, for which f is calculated but not published. Some literature ionized gas estimates erroneously do not correct upward for mean particle mass; in these cases, we multiply the published rates by 1.4. Column 5: logarithm of the black hole mass. Errors from reverberation mapping (references aa, bb, and gg) have the uncertainty in the virial coefficient f added in quadrature ($\delta f=0.44;$ Woo et al. 2010). Column 6: logarithm of the bolometric AGN luminosity, as determined from the infrared, 5100 Å, or [O iii] luminosity, in order of preference. In one case, NGC 7582, the only available measurement is from X-rays. Errors for L_AGN determined from 5100 Å or [O iii] combine the observational error and the scatter in the bolometric correction. Errors for L_AGN determined from the infrared or X-rays are not listed; the errors are assumed to be 0.3 dex. Column 7: velocity dispersion from HyperLeda (Makarov et al. 2014) for all sources but NGC 4253 (Woo et al. 2010) and NGC 1266 (Cappellari et al. 2013). Column 8: references for data in columns 2–6.

References. Seyfert type: (a) Davis et al. (2012), (b) Osterbrock & Martel (1993), (c) Véron-Cetty & Véron (2006), (d) NED. Mass outflow rate: (A) Alatalo et al. (2015), (B) Crenshaw et al. (2015), (C) García-Burillo et al. (2014), (D) Müller-Sánchez et al. (2011), (E) Riffel et al. (2009), (F) Riffel et al. (2013), (G) Schönell et al. (2014), (H) Schönell et al. (2017), (I) Veilleux et al. (1999a), (J) Zschaechner et al. (2016). Black hole mass: (aa) Bentz et al. (2006), (bb) Bentz et al. (2010), (cc) Greenhill et al. (2003), (dd) Gültekin et al. (2009), (ee) Kuo et al. (2011), (ff) Lodato & Bertin (2003), (gg) Peterson et al. (2004), (hh) Wold et al. (2006). Bolometric AGN luminosity: (1) Alatalo et al. (2015), (2) Bentz et al. (2013), (3) Liu et al. (2009), (4) Moorwood et al. (1996), (5) Runnoe et al. (2012), (6) Rivers et al. (2015), (7) Schmitt et al. (2003), (8) Ward et al. (1980).

Download table as:  ASCIITypeset image

Figure 14 shows the total mass outflow rate (summed over all gas phases) versus black hole mass. Error bars on the IFS data reflect statistical uncertainties in the line fits to individual spaxels (as determined from Monte Carlo simulations in the case of absorption lines and least-squares errors in the case of emission lines), which have been propagated and summed in quadrature across all outflowing spaxels. However, unquantified systematic errors (due to, e.g., uncertain 3D spatial distribution) could be larger than the measurement uncertainties.

Zoom In Zoom Out Reset image size

Because of the significant intrinsic scatter in this plot that is not accounted for by the measurement errors, we estimate the strength and properties of a possible linear relationship using a Bayesian regression that fits the intrinsic scatter and uses a Gaussian mixture model for the intrinsic distribution of independent variables (LINMIX_ERR; Kelly 2007). The default mixture of three Gaussians and the Gibbs sampler are used. The resulting 95.5% (2σ) ranges for the correlation coefficient and slope are −0.01 to 0.88 and −0.02 to 1.37, respectively, with medians in the posterior distributions of 0.50 and 0.68. The median intrinsic scatter is 0.98, meaning that there are unaccounted-for variables that contribute significantly or that the error is underestimated. The best-fit linear regression line (using the medians and 1σ deviations from the median of intercept and slope) is then

$$\begin{eqnarray*}\mathrm{log}[({dM}/{dt})/({M}_{\odot }\,{\mathrm{yr}}^{-1})] & = & -({4.0}_{-2.5}^{+2.6})+({0.68}_{-0.34}^{+0.32})\\ & & \times \mathrm{log}({M}_{\mathrm{BH}}/{M}_{\odot }).\end{eqnarray*}$$

A robust fit that does not account for errors and bisects the fit of X on Y and Y on X gives a steeper slope of 1.2.

This statistical test points to a positive correlation between black hole mass and total mass outflow rate in this sample at 95.5% confidence, though the slope of this correlation is poorly constrained (95.5% range 0.0–1.4). It also suggests that unaccounted-for systematic errors or other variables play a significant role in their relationship. For instance, while the relationship among the different gas phases is not fully understood spatially, it may be that adding the phases together means that some phases are counted twice. The ionization correction to the Na i D measurement tries to account for all gas phases, but it is entirely possible that the gas clouds traced by Na i D are not spatially coincident with those detected by ionized or molecular gas tracers. In fact, they are often spatially disjoint (Figure 13). In this case, adding together the masses inferred from different phases is appropriate, since they trace different parts of the wind.

Selection effects may also play a role. Since most or all galaxies host nuclear black holes, they will appear on this plot somewhere. Studies of galaxy-scale outflows in complete, unbiased samples of galaxies will show whether the outflows in this sample are simply the tip of the iceberg (i.e., there is a broad distribution of mass outflow rates extending to zero) or whether the mass outflow rates measured here are representative of all AGNs with high Eddington ratio.

Because black hole mass correlates with velocity dispersion, the mass outflow rate should correlate with stellar velocity dispersion. Stellar velocity dispersion is known for only a subset of the current sample (Table 1). Figure 15 shows that there may be a correlation between these two quantities, but the correlation is not significant at the 2σ level. The 95.5% range of possible correlation coefficients and slopes is −0.08 to 0.97 and −0.7 to 10.6. The best-fit line is

$$\begin{eqnarray*}\mathrm{log}[({dM}/{dt})/({M}_{\odot }\,{\mathrm{yr}}^{-1})] & = & -({9.2}_{-5.5}^{+5.4})+({4.8}_{-2.5}^{+2.5})\\ & & \times \mathrm{log}[\sigma /(\mathrm{km}\,{{\rm{s}}}^{-1})].\end{eqnarray*}$$

Based on the correlation between log(dM/dt) and log(M_BH) and the M_BH–σ relation (Gültekin et al. 2009), the slope should be 2.9 ± 1.5, which is consistent with the fitted slope within the (large) errors.

Zoom In Zoom Out Reset image size

Energy outflow rates are also published for most of the Seyferts in Table 7. The correlation between dE/dt and M_BH or σ (Figure 16) is more significant than the correlations with mass outflow rate. The posterior distributions of correlation coefficients and slopes for dE/dt versus M_BH have a 95.5% range of 0.15–0.94 and 0.31–2.28, respectively. For dE/dt versus σ, these ranges are 0.08–0.98 and 1.1–17.3. The best-fit lines are

$$\begin{eqnarray*}\mathrm{log}[({dE}/{dt})/(\mathrm{erg}\,{{\rm{s}}}^{-1})] & = & -({32.0}_{-3.8}^{+3.7})+({1.27}_{-0.47}^{+0.48})\\ & & \times \mathrm{log}({M}_{\mathrm{BH}}/{M}_{\odot })\end{eqnarray*}$$

and

$$\begin{eqnarray*}\mathrm{log}[({dE}/{dt})/(\mathrm{erg}\,{{\rm{s}}}^{-1})] & = & -({23.0}_{-8.1}^{+7.8})+({8.7}_{-3.6}^{+3.7})\\ & & \times \mathrm{log}[\sigma /(\mathrm{km}\,{{\rm{s}}}^{-1})],\end{eqnarray*}$$

the best-fit correlation coefficients are 0.61 and 0.68, and the intrinsic scatter is 1.4 dex.

Zoom In Zoom Out Reset image size

The quasars in this sample host both a starburst and a luminous AGN, and it is important to determine which power source is driving the observed large-scale outflows. Previous work has argued that the Mrk 231 outflow is powered primarily by the AGN (e.g., Feruglio et al. 2010; Fischer et al. 2010; Rupke & Veilleux 2011; Feruglio et al. 2015; Morganti et al. 2016). In PG1700+518, the culprit is clearly the jet spawned by black hole accretion, while in PG1613+658, the outflow is driven by the AGN within the NLR.

In other systems (particularly F05189−2524 and PG1411+442 and, to a lesser extent, F13218+0552 and F13342+3932) the velocities ($\langle {v}_{98 \% }\rangle$) exceed those observed in local starbursts (e.g., Rupke et al. 2005c). Previous studies have found that the velocities of kpc-scale outflows in AGNs reach the highest velocities (∼1000 $\mathrm{km}\,{{\rm{s}}}^{-1}$ or higher) only at quasar-like luminosities (∼10⁴⁵ erg s⁻¹), and that there is a correlation between outflow velocity and AGN luminosity (Sturm et al. 2011; Rupke & Veilleux 2013b; Spoon et al. 2013; Veilleux et al. 2013a; Cicone et al. 2014; Harrison et al. 2014; Stone et al. 2016; Fiore et al. 2017; González-Alfonso et al. 2017; Sun et al. 2017). Figure 17 shows how velocity depends on AGN luminosity in IFS observations of ionized and neutral gas velocities in nearby Type 1 and 2 quasars and luminous starbursts (Rupke & Veilleux 2013b and this study). While there is still evidence that the highest outflow velocities only appear in quasars, there are many quasars with velocities comparable to those in starbursts. No correlation is evident. However, the dynamic range in luminosity of the current sample is relatively small.

Zoom In Zoom Out Reset image size

To further quantify how these outflows relate to their potential power sources, Table 8 lists the ratios of wind properties to properties of the star-forming host or AGN.

Table 8.  Power Sources and Feedback

Galaxy	Phase	log[$({dM}/{dt})/$	log[$({dM}/{dt})/$	log[$(c\,{dp}/{dt})/$	log[$({dE}/{dt})/$	log[$({dM}/{dt})/$	log[$(c\,{dp}/{dt})/$	log[$({dE}/{dt})/$
		SFR]	${({dM}/{dt})}_{\mathrm{SB}}$]	($3.5\times {L}_{\mathrm{SB}}$)]	${({dE}/{dt})}_{\mathrm{SB}}$]	${({dM}/{dt})}_{\mathrm{acc}}$]	LAGN]	LAGN]
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
I Zw 1	Neutral	0.40	1.38	0.12	−0.28	1.51	−0.43	−3.48
F05189-2524	Neutral	0.12	1.10	0.44	0.28	2.08	0.59	−2.08
⋯	Ionized	−1.47	−0.49	−1.45	−1.99	0.49	−1.30	−4.34
⋯	Molecular	0.57	1.55	0.50	−0.10	2.53	0.65	−2.45
⋯	Total	0.70	1.68	0.77	0.43	2.66	0.92	−1.92
F07599+6508	Neutral	−0.89	0.09	−1.32	−1.78	0.98	−1.25	−4.23
⋯	Ionized	−2.79	−1.82	−3.07	−3.64	−0.93	−3.01	−6.09
⋯	Total	−0.88	0.10	0.50	−1.78	1.03	−1.25	−4.22
Mrk 231	Neutral	−0.24	0.74	−0.06	−0.33	1.72	0.09	−2.68
⋯	Ionized	−2.54	−1.56	−2.23	−2.50	−0.58	−2.08	−4.85
⋯	Molecular	0.80	1.78	0.73	0.24	2.76	0.89	−2.12
⋯	Total	0.84	1.82	0.80	0.34	2.80	0.95	−2.01
F13218+0552	Neutral	−1.15	−0.17	−1.37	−2.04	0.50	−1.52	−4.70
⋯	Ionized	−1.05	−0.07	−1.64	−2.65	0.60	−1.80	−5.31
⋯	Total	−0.80	0.18	0.74	−1.94	0.91	−1.33	−4.61
F13342+3932	Neutral	−1.44	−0.46	−1.61	−1.77	0.55	−1.42	−4.09
⋯	Ionized	−0.78	0.20	−1.14	−2.05	1.21	−0.95	−4.37
⋯	Total	−0.70	0.28	0.74	−1.58	1.32	−0.82	−3.91
PG1411+442	Neutral	⋯	⋯	⋯	⋯	2.07	0.79	−1.71
PG1613+658	Ionized	−1.46	−0.48	−2.11	−2.83	0.00	−2.22	−5.68
PG1700+518	Ionized	0.23	1.21	−0.10	−0.03	1.50	−0.16	−3.07
F21219-1757	Neutral	−0.12	0.86	−0.18	−0.71	1.67	−0.20	−3.24

Notes. Column 2: gas phase. Column 3: logarithm of the mass outflow rate normalized to the star formation rate. Column 4: logarithm of the mass outflow rate divided by the hot gas mass production rate from a continuous starburst (Leitherer et al. 1999). Column 5: logarithm of the momentum outflow rate divided by the momentum injection rate from a continuous starburst (Leitherer et al. 1999; Veilleux et al. 2005; Heckman et al. 2015; González-Alfonso et al. 2017). Column 6: logarithm of the energy outflow rate divided by the mechanical energy production rate from a continuous starburst (Leitherer et al. 1999). Column 7: logarithm of the mass outflow rate divided by the accretion rate, calculated using a 10% efficiency of radiative energy released by accretion. Column 8: logarithm of the momentum outflow rate divided by the momentum input rate if the AGN bolometric luminosity accelerates the wind and each photon is intercepted once. Column 9: logarithm of the energy outflow rate divided by the AGN luminosity.

Download table as:  ASCIITypeset image

To compare to the star-forming host, the mass entrainment efficiency η (Rupke et al. 2005c) is computed as dM/dt normalized to the star formation rate, which comes, in turn, from the AGN fraction in Table 1 and the total infrared luminosity (Veilleux et al. 2009b). The mass, momentum, and energy outflow rates are also compared to the predicted rates from the starburst itself (Leitherer et al. 1999; Heckman et al. 2015; González-Alfonso et al. 2017). (The AGN fraction for PG1411+442 is formally 100%, so these quantities are not calculated.) The value of η ranges from 0.03 to 7, with a median value of 0.8. This median is a factor of several higher than in the purely neutral phase of starburst-driven winds in ULIRGs (Rupke et al. 2005c). The median increase in wind mass from entrainment, compared to the expected input from starbursts, is a factor of 7, which is again higher than in ULIRG starbursts (Rupke et al. 2005c). Seven out of nine systems require more momentum than can be input by the starburst (without multiple photon scatterings from the radiation field), and one of the remaining two is AGN-driven (PG1613+658). Three out of nine require more energy than produced by the starburst, and five out of nine require more than 20% of the starburst energy (meaning very little energy can be lost to radiation or dissipative shocks). Two systems, F07599+6508 and F13218+0552, have winds that could plausibly be explained by starbursts. Of the four remaining systems that require $\lt 10 \%$ of the starburst energy, one is clearly an AGN outflow (PG1613+658), and, in another, the ionized outflow rates are lower limits.

This conclusion that the large-scale outflows in this sample are largely powered by the AGN rather than star formation is strengthened by comparing to starbursts of similar mass. Though the errors on the slopes are large, the observed correlations between quasar + Seyfert outflow rates and σ may be steeper than the correlations between log(dM/dt) or log(dE/dt) and logarithm of circular velocity, log(v_c), determined for starbursts (slopes 2.9 ± 0.7 and 4.64 ± 0.8, respectively; Rupke et al. 2005c). Figures 15 and 16 show the relationship between log(dM/dt) and log(v_c) for starbursts after applying the definition ${v}_{{\rm{c}}}^{2}\equiv 2{\sigma }^{2}+{v}_{\mathrm{rot}}^{2}$. The result is insensitive to the choice of ${v}_{\mathrm{rot}}/2\sigma$, which is 0.42, on average, for ULIRGs (Dasyra et al. 2006b) and 0.58 for the Seyferts in Table 7 (based on HyperLeda data). The observed best-fit steep correlations with σ imply that while NLR outflows in nearby Seyferts are not too different from their starburst counterparts, galaxies hosting luminous black holes are ejecting mass and energy at rates many times those of starburst hosts of similar σ.

The conclusion is that most or all of the observed outflows are AGN-driven; however, star formation could play a role in driving the lowest-velocity and/or lowest-energy outflows.

Models for the powering of quasar-driven, large-scale outflows can be coarsely grouped into two categories: energy-conserving blast waves powered by small-scale accretion disk winds (e.g., King et al. 2011; Faucher-Giguère et al. 2012; Zubovas & King 2012; Wagner et al. 2013; Hopkins et al. 2016; Gaspari & Scadowski 2017) and momentum-conserving winds driven by multiple photon scatterings in dusty gas with high optical depths in the infrared (e.g., Ishibashi & Fabian 2015; Thompson et al. 2015; Costa et al. 2017). Besides the high observed velocities in quasar-mode winds, both types of models have striven to explain large "momentum boosts": the momentum outflow rate is a factor of a few to a few tens larger than the outflow rate expected from pure radiation pressure driving assuming one photon scattering (Sturm et al. 2011; Rupke & Veilleux 2013b; Spoon et al. 2013; Veilleux et al. 2013a; Cicone et al. 2014; Stone et al. 2016; González-Alfonso et al. 2017). However, these observational studies have contained few examples of Type 1 quasars.

To facilitate comparison to wind models and to assess feedback (Section 4.3), the outflow rate is compared to the accretion rate and the momentum and energy outflow rate are compared to the AGN luminosity in Table 8.

The momentum boosts in the neutral and ionized phases of nearby quasar winds (Table 8 and Figure 18) show a wide dynamic range, from $(c\,{dp}/{dt})/{L}_{\mathrm{AGN}}=0.001$ to 6. For comparison, molecular gas measurements for this sample are taken to be total momentum outflow rates and maximum velocities from González-Alfonso et al. (2017). Momentum boosts tend to be higher but velocities lower in the molecular phase of individual systems. When all phases are summed, half of the IFS quasar sample shows $(c\,{dp}/{dt})/{L}_{\mathrm{AGN}}\gt 1$, with a range of 0.01–20.

Zoom In Zoom Out Reset image size

Roughly half the sample lies 2–3 orders of magnitude below predictions of the momentum boost from energy-conserving blast-wave models (Faucher-Giguère et al. 2012). This implies either that the blast wave is not energy-conserving in these systems, that some other mechanism drives the wind (such as radiation pressure), or that another phase dominates the momentum budget. Measurements of one of these systems, I Zw 1, have not detected a molecular outflow, though the neutral outflow velocities in this system are comparable to the disk velocities, and the outflow could be hidden in the core of a molecular-line profile (Veilleux et al. 2013a; Cicone et al. 2014).

Comparison of the observed velocities at a given L_AGN (Figure 17) with energy-conserving models (Zubovas & King 2012) shows that the models can reproduce the spatially averaged highest velocities in most, though not all, quasars. These models also match the best-fit slope in the dM/dt versus M_BH plane (Figure 14): 0.63 (theory) vs. 0.68 (observed). However, these models overpredict the average velocity in most systems and the highest velocities in some Type 1 quasars. They also overpredict the observed mass outflow rates by 1–2 orders of magnitude (Figure 14). However, the observed mass outflow boost in AGN-driven winds compared to starburst-only winds (0.01–10; Section 4.1) is of the same order of magnitude as predicted by some simulations (Hopkins et al. 2016).

The efficiencies with which the AGNs power these winds are high in the sense that the mechanical energy required is typically a very small fraction of the AGN luminosity. The range of $({dE}/{dt})/{L}_{\mathrm{AGN}}$ for both Seyferts and quasars is large, with most sources in the range 0.002%–2%; such values are low enough for models to successfully drive AGN outflows (e.g., Hopkins & Elvis 2010; Zubovas & King 2012) and are similar to efficiencies assumed in numerical simulations (e.g., Hopkins et al. 2016). There is no correlation of this efficiency measure with M_BH (Figure 19), L_AGN, or Eddington ratio.

Zoom In Zoom Out Reset image size

The conclusion is that it is difficult to distinguish among competing models or constrain individual models. However, none of these models make predictions about how outflow properties scale with black hole mass—such predictions would optimally leverage the results of the current study.

A key remaining question is how these large-scale, quasar-driven outflows will affect the star formation and AGN activity in their hosts. This question is difficult to answer conclusively with the current data.

It is clear that in an absolute sense, the most massive black holes produce the most massive and energetic outflows. High-velocity, high-momentum boost feedback, as observed in part of this sample, may be the most effective at shaping the M_BH–σ relation at high black hole masses in simulations as well (Anglés-Alcázar et al. 2017). In essence, these outflow rates illustrate how well the black hole, through its current accretion rate, is moving away gas that could at some future time feed the black hole. This future potential depends on (1) the uncertain timescales for gas at large radii to inflow to the accretion disk and then to the black hole and (2) what fraction of gas at large radii will actually fall into the black hole at some future time.

A possible relative measure of feedback on the AGN itself is the mass outflow rate divided by the present black hole accretion rate. In the quasar + Seyfert sample, this ratio is mostly in the range 1–1000, which is larger than at radii <100 pc in the simulations of Hopkins et al. (2016). Figure 20 shows how $\tfrac{{({dM}/{dt})}_{\mathrm{out}}}{{({dM}/{dt})}_{\mathrm{acc}}}$ depends on black hole mass and Eddington ratio. There is no significant correlation with M_BH. However, the Seyfert subsample, on average, has higher $\tfrac{{({dM}/{dt})}_{\mathrm{out}}}{{({dM}/{dt})}_{\mathrm{acc}}}$ and lower M_BH than the quasar subsample. The correlation with Eddington ratio is almost 2σ (using Bayesian fitting as in Section 3.3), with a median slope of −1.0 and a 95.5% range of −2.4 to 0.1. However, the Bayesian test does not account for the fact that the variables are not completely independent, since both depend on L_AGN (${dM}/{{dt}}_{\mathrm{acc}}\sim {L}_{\mathrm{AGN}}$ and Eddington ratio is proportional to L_AGN).

Zoom In Zoom Out Reset image size

We conclude that there is weak evidence of inverse correlations of $\tfrac{{({dM}/{dt})}_{\mathrm{out}}}{{({dM}/{dt})}_{\mathrm{acc}}}$ with M_BH and Eddington ratio. If such correlations exist, then one consequence could be that relative feedback (gas out on large scales over gas in at small scales) is higher in systems with lower Eddington ratio and black hole mass. Alternatively, it could be that the assumption of a constant radiative efficiency from accretion (how well the black hole converts accretion energy into luminosity) is incorrect. Instead, the efficiency may be higher in higher-mass black holes (Davis & Laor 2011). Higher radiative efficiency means that a lower accretion rate ${({dM}/{dt})}_{\mathrm{acc}}$ produces the same luminosity L_AGN, which could remove any correlation of $\tfrac{{({dM}/{dt})}_{\mathrm{out}}}{{({dM}/{dt})}_{\mathrm{acc}}}$ with M_BH. A higher efficiency could result from higher black hole spin, since in standard thin-disk accretion (Shakura & Sunyaev 1973), the efficiency depends on the spin (Reynolds 2013). However, current measurements of black hole spin do not show higher spin in higher-mass black holes (Reynolds 2013).

The issue of timescale is also an important consideration. Quasar outflows are powered by current and past AGN activity. The dynamical timescales of these flows are of order 1–100 Myr, while AGN accretion may have a smaller duty cycle of ∼0.1 Myr (e.g., King & Nixon 2015; Schawinski et al. 2015). Thus, the accretion rate (and thus the AGN luminosity and Eddington ratio) may change significantly on shorter timescales than the outflow. Some of the scatter (and the correlation itself) in the right panel of Figure 20 could result from short-timescale variations of L_AGN if ${dM}/{{dt}}_{\mathrm{out}}$ and M_BH indeed vary on longer timescales.

As with the other correlations found in this work, selection biases likely exist, and larger samples are needed. For instance, the lower regions of Figure 20 could be occupied by galaxies with lower mass outflow rates but similar accretion rates and black hole masses (e.g., the radio-mode feedback regime in radiatively inefficient accretion modes), which could also fill in other parts of Figures 14–16. An open question is whether the distributions in outflow properties, taken over all galaxies with accreting black holes, are continuous or significantly bimodal. For instance, in a hypothetical on-or-off model, outflows are either present and look like the ones present here, or they are absent and any outflow rates are exactly zero. The result would be a bimodal distribution in outflow rates. In this on-or-off model, the correlations presented here may accurately describe active feedback modes.