Ionized Outflows in Nearby Quasars are Poorly Coupled to their Host Galaxies

We analyze Multi Unit Spectroscopic Explorer observations of nine low-redshift (z<0.1) Palomar-Green quasar host galaxies to investigate the spatial distribution and kinematics of the warm, ionized interstellar medium, with the goal of searching for and constraining the efficiency of active galactic nucleus (AGN) feedback. After separating the bright AGN from the starlight and nebular emission, we use pixel-wise, kpc-scale diagnostics to determine the underlying excitation mechanism of the line emission, and we measure the kinematics of the narrow-line region (NLR) to estimate the physical properties of the ionized outflows. The radial size of the NLR correlates with the AGN luminosity, reaching scales of $\sim 5\,$kpc and beyond. The geometry of the NLR is well-represented by a projected biconical structure, suggesting that the AGN radiation preferably escapes through the ionization cone. We find enhanced velocity dispersions ($\sim 100\,$km$\,$s$^{-1}$) traced by the H$\alpha$ emission line in localized zones within the ionization cones. Interpreting these kinematic features as signatures of interaction between an AGN-driven ionized gas outflow and the host galaxy interstellar medium, we derive mass outflow rates of $\sim 0.008-1.6\, M_\odot \,$yr$^{-1}$ and kinetic injection rates of $\sim 10^{39}-10^{42} \,$erg$\,$s$^{-1}$, which yield extremely low coupling efficiencies of $\lesssim 10^{-3}$. These findings add to the growing body of recent observational evidence that AGN feedback is highly ineffective in the host galaxies of nearby AGNs.


INTRODUCTION
The accretion of matter onto supermassive black holes (BHs) residing within the nuclei of galaxies is widely accepted to be the source of the prodigious energy released in active galactic nuclei (AGNs; Rees 1984). The net amount of energy emitted during the growth of the BH greatly exceeds the binding energy of its host galaxy (Fabian 2012), positioning the AGN-host galaxy inter-on galaxy-wide scales that are powerful enough to heat or expel the interstellar medium (ISM) from the host galaxy (e.g., Silk & Rees 1998;Somerville et al. 2008;Schaye et al. 2015;Sijacki et al. 2015;Lacey et al. 2016), thereby suppressing star formation (Dubois et al. 2016) and keeping the host galaxy quiescent (Fabian 2012).
Despite the ubiquitous evidence of multi-phase gas outflows in AGN hosts, the role of AGN feedback in establishing the host galaxy properties and its ability to quench ongoing star formation activity are still intensively debated. Local low-luminosity AGNs exhibit global star formation rates (SFRs) lower than those measured from the inactive, star-forming galaxy population (Ellison et al. 2016;Leslie et al. 2016;Jackson et al. 2020), and some systems are effectively passive (Ho et al. 2003). On the other hand, the more luminous nearby AGN host galaxies are similar to the inactive galaxy population in several key properties, such as SFR, atomic hydrogen gas content, molecular gas mass, and dust emission and/or dust attenuation (Maiolino et al. 1997;Evans et al. 2001;Scoville et al. 2003;Evans et al. 2006;Bertram et al. 2007;Ho et al. 2008;Fabello 1 While some recent studies use the term "extended narrow-line region" (ENLR), we prefer to use the more traditional terminology of "NLR" and reserve "ENLR" to refer to gaseous regions outside of the host galaxy that are ionized by the AGN radiation (e.g., Fu & Stockton 2009  The open circles show the sample at z < 0.5 . We focus on the nearest (z < 0.1) host galaxies with available MUSE data. These targets mostly sample the L5100 10 45 erg s −1 and MBH 10 8 M range, with PG 1426+015 being the only exception (Table 1). et al. 2011;Harrison et al. 2012;Geréb et al. 2015;Zhu & Wu 2015;Husemann et al. 2017;Bernhard et al. 2019;Ellison et al. 2019;Shangguan et al. 2018;Shangguan & Ho 2019;Shangguan et al. 2020a;Grimmett et al. 2020;Jarvis et al. 2020;Yesuf & Ho 2020;Zhuang et al. 2021). Far from being quenched, stars are actively forming within the host galaxies of luminous AGNs , and in some the star formation efficiencies are comparable to those seen in starburst systems (Kirkpatrick et al. 2020;Shangguan et al. 2020b;Zhuang & Ho 2020;Xie et al. 2021). These observations suggest that negative AGN feedback, if present, likely produces only a localized effect on the host galaxies.
The most luminous AGNs, the quasars, are the ideal laboratories to scrutinize the interplay between AGN feedback and the host galaxy ISM. Within the context of the merger-induced evolution scenario of quasars (Sanders et al. 1988), the birth of a largely unobscured quasar is produced after the overwhelming release of energy by the AGN expels the enshrouding gas and dust near the nucleus Treister et al. 2010). However, studying active galaxies is challenging because the radiation from the AGN contaminates the emission from the underlying host, a problem that is most severe for type 1 (unobscured, broad-line) sources, although it is not entirely negligible even in type 2 (obscured, narrow-line) systems. Integral-field unit (IFU) observations offer the possibility to overcome such a limitation by spatially resolving the host galaxies (e.g., Jahnke et al. 2004;Lipari et al. 2009;Husemann et al. 2014;Harrison et al. 2016;Husemann et al. 2017;Ilha et al. 2019;Feruglio et al. 2020;Kakkad et al. 2020;Lacerda et al. 2020;Riffel et al. 2021;Scholtz et al. 2021). The AGN emission, as a point-like source that follows the observation point-spread function (PSF), in principle can be deblended from the extended emission of the host using sophisticated computational procedures (e.g., Husemann et al. 2013;Rupke et al. 2017;Husemann et al. 2022) to produce cleaner measurements of the host galaxy properties. This work reports on seeing-limited and ground-layer adaptive optics-aided (FWHM seeing ∼ 0. 6 − 1 , which corresponds to physical scales 1 kpc) Very Large telescope (VLT) Multi Unit Spectroscopic Explorer (MUSE) observations targeting the optical light toward the host galaxies of nine Palomar-Green (PG) quasars (Schmidt & Green 1983). The low-redshift (z < 0.5) subset of the PG survey contains 87 type 1 quasars (Boroson & Green 1992). This sample is one of the most studied quasar surveys to date, enjoying a rich repository of multi-wavelength data for the AGN and host galaxy, encompassing X-ray observations (Reeves & Turner 2000;Piconcelli et al. 2005;Bianchi et al. 2009), optical long-slit spectra (Boroson & Green 1992;Ho & Kim 2009), Hubble Space Telescope (HST) highresolution imaging of the host galaxy (Kim et al. 2008(Kim et al. , 2017Zhang et al. 2016;Zhao et al. 2021), dust properties for both the torus and host galaxy (Shi et al. 2014;Petric et al. 2015;Shangguan et al. 2018;Zhuang et al. 2018), global SFRs (Xie et al. 2021), molecular gas masses measured from CO (Shangguan et al. 2020a,b), and radio properties (Kellermann et al. 1989(Kellermann et al. , 1994. The main goal of this study is to characterize the effects of AGN feedback on the host galaxy traced by optical nebular emission lines. We focus on a small subset of the nearest (z < 0.1) PG quasars with L 5100 ≡ λL λ (5100Å) 10 45 erg s −1 and M BH 10 8 M (Figure 1). The upper limit on redshift is chosen to avoid mixing observations with large differences in physical scale resolution, but trading off the quasar luminosity and BH mass range covered by the sample. We employ no further constraints to select our targets. A companion publication will examine the star formation activity and star formation efficiency of the host galaxies.

SAMPLE AND OBSERVATIONS
We analyze optical light VLT-MUSE wide-fieldmode IFU observations of nine low-redshift (z < 0.1) PG quasar host galaxies (Table 1). We observed three targets during our European Southern Observatory (ESO) program 0103.B−0496(B; PI: F. Bauer). We complement those observations with archival VLT-MUSE data of PG quasar host galaxies taken from the ESO programs 094.B−0345(A), 095.B−0015(A), 097.B−0080(A), 0101.B−0368(B), and 0104.B−0151(A). All observations were carried out from January 2015 to October 2019. The field-of-view of a single MUSE observation is nearly 1 ×1 with a pixel size sampling of 0. 2 × 0. 2, which generates ∼ 90, 000 spectra per pointing. The MUSE spectral coverage ranges from ∼4700 to ∼9350Å, with a wavelength sampling of 1.25Å channel −1 and at a mean full width at half maximum (FWHM) resolution of R ≈ 3000 (2.65Å). We analyze the MUSE archival data cubes, 2 which correspond to a combination of multiple observing blocks (OBs) typically ranging from 2 to 4 for single program observations. In the particular case of PG 0050+124, we used the "MUSE-DEEP" 3 reduced data cube, which corresponds to a combination of 15 single OBs focused on maximizing the output signal contrast. The benchmark is to obtain a better signal-to-noise (S/N) ratio than any individual observing run data cube. MUSE and MUSE-DEEP data cubes are reduced following the same pipeline framework outlined in Weilbacher et al. (2020). We note that all the observations were carried out in good weather conditions, with most of the OBs graded "A." A few OBs were graded "B" because the sky conditions did not meet the requested observation seeing (e.g., PG 1126−041). Only one OB was graded "C" because the data reduction pipeline did not converge when performing the master sky line fit 4 . This means that only the sky lines were not properly modeled and subtracted due to an inaccurate characterization during the run of the data reduction pipeline. This single OB (from a total of four) pertains to observations of PG 1426+015, and we included it when analyzing this host galaxy data cube. The observations were taken under average seeing conditions of ∼ 1 (Table 2). Five targets were observed under natural seeing conditions, while the remaining four MUSE observations were assisted by the ground-layer adaptive optics (GLAO) module.
We processed the combined data cubes by applying the Zurich Atmosphere Purge (ZAP) sky subtraction tool (version '2.1.dev'; Soto et al. 2016) with cfwidthSP = 5 (the other parameters were not modified) to remove residual sky subtraction features. We masked all the pixels that were affected by spectra flux saturation along the full spectral range. We also masked the spectra at wavelengths where strong sky-line residuals are present in the original data cubes (∼5578.5, 5894.6, 6301.7, 6362.5, and 7640Å) due to the uncer- (4) Image quality of the observations. In the case of the AO-aided observations we report the spatial resolution values measured from the whitelight image (data cube header keyword "SKY RES"), while for the "seeing-limited" observations we report the seeing value. (5) MUSE instrument mode. All the observations were taken in wide-field-mode (WFM), but a few (WFM-AO) were aided by GLAO. Observations taken in naturalseeing conditions (WFM-noAO) were carried out before the GLAO system became operational for MUSE.
tain ZAP correction. We note that in the AO-aided observations the MUSE spectra have been previously masked by the data reduction pipeline at the wavelength range surrounding NaD emission (5840-5940Å). The data cube for each source was corrected for Galactic reddening assuming a Cardelli et al. (1989) extinction law and extinction values from Green et al. (2019). We further checked the flux calibration of the MUSE observations by comparing the host galaxy MUSE-based Sloan Digital Sky Survey (SDSS; Abolfathi et al. 2018) r-and i-band total magnitudes with available SDSS Petrosian magnitudes, finding good agreement between both quantities, with average magnitude difference of ∼ 0.2 (MUSE magnitudes are lower) and 0.2 dex scatter. In this work we adopt the MUSE line-spread function (LSF) parametrization of Guérou et al. (2017) to correct the line widths for instrumental resolution.

METHODS
Our main goal is to study AGN feedback effects on their host galaxies by analyzing the ionized gas properties. We characterize the host galaxies using two different binning approaches, separately optimized for the stellar component and nebular emission lines. When  Figure 2. MUSE white-light images of the nine PG quasar host galaxies. The yellow-shaded regions represent the circular masks that we manually applied to the individual data cubes to avoid additional sources from contaminating the host galaxy continuum across the line-of-sight. modeling the host galaxy stellar component, we focus on maximizing the recovery of the S/N of the continuum, even at the expense of degrading the spatial resolution. On the contrary, when characterizing the emission lines we only slightly downgrade the observation spatial resolution to boost the S/N of the data. Preserving the highest possible spatial resolution is critical to deblending the nuclear emission from that of the host, as the nucleus is, by far, the brightest component in the MUSE data cube, and any subtraction error leads to large systematics propagated throughout data analysis.
As with imaging studies of AGN host galaxies (e.g., Zhao et al. 2021), in the IFU data the unresolved emission coming from the nuclear region is seen as a blurred, point-like source that follows the shape of the PSF. This unresolved component needs to be rigorously characterized to not misinterpret the extended emission component. The main difference with respect to the imaging studies is that we have to account for a nuclear spectrum instead of a single flux value per pixel, implying that a correct deblending procedure must include an accurate characterization of this spectrum, which is to be subtracted from the data cube. Two different strategies can be adopted to deblend the AGN emission from that of the host galaxy. The first option is to follow an imaging-based strategy, where the PSF and its variation with wavelength are characterized by modeling pseudo broad-band images derived from the data cube. The model PSF is scaled by the nuclear spectra flux density at each wavelength, and then subtracted from the data cube (e.g., QDeblend 3D Husemann et al. 2022). The alternative strategy is to follow a spectra fittingbased approach, where the nuclear spectrum is modeled to then be used as a template when pixel-wise fitting the data cube. This strategy is similar to that of using stellar spectral templates when modeling the galaxy stellar emission. In this work we adopt the latter option as this procedure does not rely on modeling the PSF. Indeed, the PSF is obtained as a by-product because the nuclear spectrum template characterizes the unresolved AGN emission spread over the data cube by the PSF. The nuclear spectrum comprises three components: (1) a featureless continuum from the hot accretion disk, which arises from scales of less than a few light days (e.g., Homayouni et al. 2019;Guo et al. 2022;Jha et al. 2022); (2) broad emission lines from the BLR, which has a size that scales with the luminosity of the source, but even in the most luminous sources studied to date rarely exceeds a light-year (e.g., Kaspi et al. 2000;Li et al. 2021); and (3) a compact component associated with the photoionized NLR near the vicinity of the nucleus, on scales much smaller than the spatial resolution of our observations (∼ kpc scale). We emphasize that subtracting the inner NLR component is critical. The [O iii] λλ4959, 5007 emission produced near the nucleus can achieve a S/N even higher than that for broad Hβ, meaning that in every pixel in which we observe the blurred broad Hβ emission we also observe nuclear [O iii] emission due to the observation finite spatial resolution (see also Husemann et al. 2016). Not considering this component would result in an overestimation of the [O iii] emission and consequently an inaccurate classification of the underlying source powering the nebular emission lines. It would also lead to the confusion of this [O iii] component as an additional kinematic signature throughout the host galaxy. Keeping in mind this observational effect, we now proceed to explain in detail our data analysis and quasar deblending technique.

Stellar Component Modeling
We first characterize the spatial extent of the host galaxies in the MUSE data cubes by collapsing them across the spectral axis and deriving "white-light images." We use the background2D task of the photutils Python package (Bradley et al. 2020) with an 1 width smoothing kernel to determine the background level of the white-light images and to mask any emission with S/N < 5. We then use the photutils task detect sources to build segmentation maps and to isolate the host galaxies from other sources within the MUSE field-of-view. At this step, we manually mask any other contaminating object covered by the host galaxy segmentation map sources (Figure 2), including, for instance, small galaxy companions (e.g., PG 1426+015), field stars (e.g., PG 2130+099), or both (e.g., PG 0050+124).
We further apply Voronoi binning to the host galaxy segmentation map after imposing a target criterion of S/N = 50 (Cappellari & Copin 2003) at 9010-9280Å in the observer frame avoiding sky-line residuals. We employ the penalized pixel-fitting method (pPXF; Cappellari & Emsellem 2004) and model the spectra of the Voronoi cells using a stellar template library generated from the extended MILES stellar population synthesis models (MIUSCAT; Vazdekis et al. 2012). The MIUSCAT composite templates are based on the MILES and calcium triplet (CaT) empirical stellar spectral libraries for single stellar populations (SSPs) of a given age and metallicity, covering the wavelength range 3465-9469Å at an uniform spectral resolution of FWHM = 2.51Å (for more details, see Vazdekis et al. 2012). The MIUSCAT templates are broadened in velocity space by taking into account their spectral resolution difference with respect to the MUSE LSF width. We mask wide 4800-5100Å and 5400-6740Å spectral windows in rest-frame and mainly covering the Hβ, [O iii], Hα, [N ii] λλ6548, 6583, and [S ii] λλ6716, 6731 emission lines among others (e.g., He i λ5875). We also mask the wavelength ranges 7230-8070Å and 8270-8510Å, where additional sky lines are observed. We note that in our observations the stellar kinematics are constrained mainly by CaT λλ8498, 8542, 8662 and Mg i λλ5167, 5173, 5184. Owing to the modest S/N and spectral resolution of the MUSE spectra, we only model the absorption lines as Gaussian functions, avoiding more sophisticated fits using a Gauss-Hermite series (Cappellari & Emsellem 2004).
The pPXF procedure allows the inclusion of additive and multiplicative polynomials during the fit, aiming to correct imperfect sky subtraction or scattered light with the former, and inaccuracies in spectral calibration or mismatches in dust reddening correction with the latter (Cappellari 2017). We experimented with both. Our results are insensitive to the inclusion of a multiplicative polynomial of order 3 or less. Conversely, the inclusion of an additive polynomial preferentially downweights the young (age 1 Gyr) SSP template, fails to produce sufficiently strong hydrogen absorption lines, and hence introduces systematic biases (overestimates) in the resulting Hα and Hβ emission-line fluxes. In light of these considerations, we do not include any additive or multiplicative polynomials during the stellar continuum modeling.

Template for the Nuclear Spectra
We extract the nuclear spectrum for each target by spatially collapsing the data cube within a circular region with radius equal to the observation spatial resolution (Table 2), centered on the peak of the quasar emission. Following standard practice (e.g., Boroson & Green 1992;Ho & Kim 2009), we decompose the continuum of the nuclear spectrum into three components: a featureless continuum to represent accretion disk emission, a "pseudo-continuum" produced by many blended broad iron emission lines, and starlight from the host galaxy. The featureless continuum is modeled as a power law characterized by an amplitude, spectral index, and pivot wavelength. We constrain the spectral index to be negative, and the pivot wavelength is restricted to lie within 4600-6000Å in the rest frame. We use the Fe ii template of PG 0050+124 (I Zw I) provided by Boroson & Green (1992) to model the iron pseudo-continuum, allowing it to adjust in amplitude, velocity width, and velocity shift (Hu et al. 2008).
The central stellar component is very poorly constrained because it is outshone by the quasar emission, whose blue, featureless continuum mimics and is degenerate with a young SSP. To circumvent this difficulty, we adopt an empirical template of the stellar continuum extracted from regions of the host galaxy that are relatively uncontaminated by the AGN. The critical assumption, of course, is that the stellar population at the off-nuclear position is similar to that of the nuclear position. We quantify the AGN contamination as a function of radius by studying the radial variation of the flux density at ∼ 4650Å. Assuming that the MUSE PSF can be described by a Moffat (1969) function with index β = 2.2 (e.g., Bacon et al. 2017;Guérou et al. 2017) and using the image quality values recorded during the observations (Table 2), we estimate that the AGN continuum on average drops to ∼ 4 % (range ∼ 3 % − 7 %) of its central peak value 5 at a projected distance of 5 times the seeing radius. Hence, for each host galaxy, we extract an AGN-free spectrum at an annulus centered at 5 times the seeing radius from the nucleus. The best-fit pPXF model then serves as the starlight template for the nuclear spectrum.   . Nuclear spectra for the nine PG quasar host galaxies. The spectra are extracted from the MUSE data cubes using a circular region centered on the quasar location and with aperture equal to the observation spatial resolution (Table 2). We only show the wavelength range covered by our spectral modeling routine. The flux densities are normalized with respect to the peak value of the Hα emission. The iron pseudo-continuum component is shifted above the power-law component to improve visualization. The broad-line components correspond to the sum of all the Gaussian sub-components (typically four in each case). The residuals (data − model) are below the ∼2.5 % limit in all cases.  The emission-line spectrum contains blends of many features, which need to be decomposed carefully. The permitted lines (e.g., Hβ and Hα) are especially complicated because they comprise heavily blended contributions from the BLR and NLR, each of which can have complex kinematics. Following common practice (e.g., Greene & Ho 2005), we first fit the well-isolated, strong narrow forbidden lines of the [O iii] λλ4959, 5007 doublet, fixing their known separation and relative intensity (Osterbrock & Ferland 2006), and assuming that both components share the same profile. The lines, though narrow, often have non-trivial profiles. We model them using as many Gaussians as necessary to achieve a satisfactory fit, although in practice two suffice: a narrow core plus a broader, often blueshifted wing. We then use the best-fitting model for [O iii] λ5007 as an empirical template to extract the flux of the narrow component of Hβ and Hα, as well as [N ii] λλ6548, 6584, which normally is also severely blended with Hα. Note that in principle [S ii] λλ6716, 6731 gives a more suitable template for treating the Hα+[N ii] complex (Ho et al. 1997), but in our sample [S ii] has much lower S/N than [O iii], and we use the [O iii] template throughout our analysis. The profile of [O iii] may be a poor representation of the narrow Hα and Hβ emission lines if the NLR is highly stratified in density (Greene & Ho 2005), although, this is usually not the case (Whittle 1985;Veilleux 1991). Nevertheless, the nuclear spectra fitting procedure is not particularly sensitive to this election. We tested the use of the [O iii] narrow core as narrow line shape template, finding that accurate nuclear spectra models can also be obtained in this case. A single Gaussian is sufficient for most of the weaker forbidden lines ([Fe vii] λ5720, [Fe vii] λ6087, , and [Ar iv] λ7170), while two components are usually required for the helium lines (He ii λ4686, He i λ5875, He i λ6678, and He i λ7065). With the NLR component thus constrained, we fit broad Hβ and Hα with multiple Gaussians, but assign no physical meaning to any of the sub-components; typically four sub-components suffice. We restrict the fit to 4450-7600Å in the rest-frame. The lower limit is set to avoid Hγ λ4340, which is not included in the MUSE spectral coverage for all sources, while the upper range is limited by the lack of an iron template at longer wavelengths.
Although the stellar continuum accounts for only a minor fraction of the total flux in the nucleus, we confirm, with the aid of the Bayesian information criterion (BIC; Schwarz 1978), that including it in the fit yields statistically improved models for six out of the nine cases (∆BIC > 10). In one case there is not enough evidence to decide, while in the two cases where the model without the stellar component is favored by the BIC, the differences in the residuals are negligible. Therefore, for the sake of consistency, all the final fits for the nuclear spectra include the stellar continuum component. The nuclear spectra and their best-fit models are shown in Figure 3.

Characterization of the Host Galaxy Ionized Gas Emission
We implement a pixel-wise fitting approach to characterize the principal bright diagnostic lines, Hα, Hβ, We preferentially preserve the spatial resolution of the observations regardless of the S/N of the lines, in contrast to the Voronoi binningbased strategy adopted when characterizing the stellar component (Section 3.1).
For each pixel, we average the spectrum by considering a box size comparable to the seeing (e.g., Swinbank et al. 2012). The average spectrum is modeled by the sum of a scaled nuclear template, a stellar continuum, plus nebular emission lines. The nuclear template is multiplied by an amplitude free parameter. We assign to each pixel the pPXF best-fit stellar continuum model of the corresponding Voronoi cell spectrum (see Section 3.1). Note that due to the strong quasar emission in the center, we do not have a stellar continuum model for the central pixels; for these central pixels, we follow the approach used when deriving the nuclear template, by adopting the best-fit pPXF solution of the AGN-free spectra as the stellar continuum template (see Section 3.2 for details). The stellar continuum component is simply scaled in amplitude. The emission lines are modeled by using Gaussian functions; in almost all cases only a single component is required. At this step we consider the MUSE instrumental resolution by adding the LSF line width in quadrature within the line model. We use the BIC to determine whether the spectra are well-characterized by our model by comparing the best-fit solution with a simple straight line fit (i.e., no emission present).
After obtaining the initial set of spectral models, we repeat our procedure using a new set of initial parameters to mitigate the known sensitivity of the leastsquares minimization technique to the initial guess. For a given pixel, we take a new set of initial parameters from the best-fit models obtained from the neighboring pixels in the previous run, and from these new solutions we select the model associated with the lowest BIC value. Finally, we use Monte Carlo re-sampling to derive the parameter uncertainties. We measure the average level of the root-mean-square of the spectrum for each pixel from the residuals, add simulated noise to the  observed spectrum assuming a normal distribution, and then repeat the fit. We iterate 100 times to obtain a probability distribution for each parameter, from which we estimate the 1 σ uncertainties from the 16th and 84th percentiles. Figure 4 presents example fits of individual pixels of PG 0923+129, and partly show the models quality that we achieve. Naturally, different correlations among the spectra model parameters arise depending on the target spectra shape. If the blurred nuclear spectra dominates the pixel emission (1st row in Figure 4), then its amplitude scale parameter tends to be anti-correlated with the Hα and Hβ line intensities and line widths. As expected, the nuclear spectrum and stellar continuum scale parameters are anti-correlated when both components significantly contribute to the target spectra emission (2nd row in Figure 4). When the stellar continuum is significant, its scale amplitude parameter tends to be correlated with the Hβ line intensity because of the presence of the Hβ stellar absorption line. The [O iii] line intensity tends to be correlated with the [O iii] line width, but only in regions where the target spectra continuum is faint (3rd row in Figure 4). Our quasar deblending approach suggests that the blueshifted wing seen in [O iii] arises from an unresolved region. However, we note that a more sophisticated analysis may be needed to fully constrain the size of the outflow traced by the blueshift seen in [O iii] (e.g., Singha et al. 2022). We emphasize that given the diversity of nuclear spectra shapes observed across our PG quasar sample (Figure 3), plus the plethora of combinations between these multiple spectra components digitized in each data cube, the parameter correlations associated with pixel-wise modeling the data must be considered only as a rough description.   (Baldwin et al. 1981; hereinafter BPT). The two diagrams provide different levels of differentiation in terms of the ionization mechanism. The [O iii]/Hβ-[N ii]/Hα diagram designates three broad classes ( Figure 5): objects below the dashed line are powered by star formation (Kauffmann et al. 2003); objects above the "maximum starburst" solid line are powered by AGNs (Kewley et al. 2001); and those in between these two boundaries have a composite source of ionization (mixture of young stars and AGNs). Classifications based on the [O iii]/Hβ-[S ii]/Hα diagram use the demarcations given by Kewley et al. (2001Kewley et al. ( , 2006, which distinguish star-forming objects from Seyferts and low-ionization nuclear emissionline regions (LINERs). The properties of low-luminosity AGNs, including the physical distinction between highionization (Seyfert) and low-ionization (LINER) sources have been extensively reviewed by Ho (2008Ho ( , 2009) and will not be repeated here. We assign a classification to all pixels whose emission lines were detected at S/N > 3, and to those whose ionization source can be classified unambiguously even with just emission-line flux limits. For the former pixels, all emission-line fluxes are corrected for internal reddening based on the observed Balmer decrement of the narrow lines, assuming an intrinsic value of Hα/Hβ = 2.86 for the star-forming pixels and Hα/Hβ = 3.1 for the others (Osterbrock & Ferland 2006), assuming the extinction curve of Cardelli et al. (1989).

Narrow-line Region
Often characterized by its extension and geometry, the NLR reflects how the ionizing radiation from the central AGN permeates through the ISM of the galaxy host. It also traces the deposition of kinetic energy and momentum when outflows are present (e.g., Greene et al. 2012;Husemann et al. 2013Husemann et al. , 2016Husemann et al. , 2019bStorchi-Bergmann et al. 2018;Sun et al. 2018) Figures 5 and A1 reveal that the PG quasar hosts have vast ISM zones in which the emission is consistent with being powered by the AGN, allowing us to characterize in detail their NLR. 6

NLR Size-Luminosity Correlation
The correlation between the NLR size and the AGN [O iii] luminosity (L [O III] ) is usually employed to constrain the properties of the gas that is photoionized by the central source (e.g., Bennert et al. 2002;Liu et al. 2013Liu et al. , 2014Hainline et al. 2014;Dempsey & Zakamska 2018 Chen et al. (2019) and fit the NLR radial profiles using Sérsic functions (Sérsic 1963) to characterize the NLR radius by the isophote at an intrinsic [O iii] surface brightness of 10 −16 erg s −1 cm −2 arcsec −2 (1+z) −4 and measured from the Seyfert-classified pixels (R NLR 16,fit ); and (2) we estimate the maximum NLR radial size (R NLR 16,max ) by computing the deprojected location of the most distant AGN-ionized region emitting above the 10 −16 erg s −1 cm −2 arcsec −2 (1+z) −4 [O iii] surface brightness threshold (e.g., Husemann et al. 2022). In both cases we account for the effect of cosmic dimming, enabling the comparison of AGN hosts observed at different redshifts (Liu et al. 2013). We account for inclination correction by measuring the minor-to-major axis ratio of the host galaxy from the white-light image using photutils. We report those estimates in Table 3. Figure 6 compares our R NLR 16,fit measurements of the PG quasars hosts with the large sample of 152 Seyfert galaxies (both type 1 and type 2, but mostly type 2) studied by Chen et al. (2019) using the Mapping Nearby Galaxies at Apache Point Observatory (MaNGA) survey (Bundy et al. 2015). Also included are luminous type 1 and type 2 quasars at z ≈ 0.5 from Liu et al. (2013Liu et al. ( , 2014. With R NLR 16,fit ≈ 1.6 − 4.7 kpc (Table 3), the PG quasars follow the best-fit NLR size-L [O III] relation previously reported by Chen et al. (2019). We only give an upper limit for the NLR size of PG 0050+124  (Chen et al. 2019). The sky-blue and blue squares correspond to the z ≈ 0.5 type 2 and type 1 quasar IFU data presented in Liu et al. (2013Liu et al. ( , 2014) and re-analyzed by Chen et al. (2019). The dashed line corresponds to the best linear fit (slope ∼ 0.42) computed by Chen et al. (2019). The arrow shows our upper limit estimate for the NLR isophote size in PG 0050+124. The PG quasar host galaxies follow the NLR luminosity-size relation.
because we cannot detect extended [O iii] emission from its central R 2.5 kpc zone (see Appendix A).
Interestingly, among all the PG quasar hosts, the zones where the AGN mainly contributes to the ionization state of the ISM coincide with the regions where Hα emission is faint and the ionized gas follows the global rotation pattern of the galaxy (Figures 5 and  A1), suggesting that this gas was uplifted from the disk due to local secular processes before being ionized by the AGN (Husemann et al. 2019a), or simply that the AGN-ionized gas emission is outshined by the H ii region ionizing photons where the star formation activity takes place. Geometrical effects may play a substantial role when determining the sizes of the NLR from simple profile fits. To highlight this, in Figure 6 we show in open circles the R NLR 16,max values for the PG quasars, and we compare with the R NLR 16,fit values (orange circles). Both estimates are in agreement within the uncertainties for most of the cases, but three outliers can be clearly seen. The PG quasar host galaxies with consistent R NLR 16,fit and R NLR 16,max estimates tend to show smooth [O iii] surface brightness maps, and with the brighter zones being detected closer to the quasar location. In one of the outliers, PG 0923+129, R NLR 16,max ≈ 0.7 kpc is significantly smaller than R NLR 16,fit ≈ 1.8 kpc. Careful scrutiny reveals that R NLR 16,fit is given by an extrapolation of the radial profile in a location where this host galaxy contains an inner star-forming ring seen in Hα ( Figure 5) and in CO (Molina et al. 2021), suggesting that the smaller R NLR 16,max value measured for this system may be produced by the AGN-ionized gas emission being outshone by the ionizing photons coming from the underlying star formation activity in that local ISM sub-structure. Another possibility could be that the dust on the star-forming ring might be attenuating the number of AGN ionizing photons that can permeate to larger radii. In the other two cases we have the opposite evidence, where R NLR 16,max values are significantly higher than the R NLR 16,fit estimates. One case corresponds to PG 1426+015, a host galaxy that has two tidal features that are being ionized by the AGN, with the most distant feature ∼ 26 kpc away from the nucleus presenting [O iii] emission above the surface brightness threshold adopted to define the NLR size. Similarly, PG 0050+124 also displays a set of distant (∼ 18 kpc) ionized clouds emitting [O iii] above the surface brightness threshold, which is surprising given that we only found a R NLR 16,fit upper limit value for this target. With log [N ii]/Hα ≈ −0.5 and log [O iii]/Hβ ≈ 1.0, the gas composition appears to be metal-poor, not unlike the situation found in HE 1353−1917 (Husemann et al. 2019a). We note that some of the systems analyzed by Chen et al. (2019) also have NLR sizes significantly above their best-fit relation ( Figure 6). They rejected merger pairs when designing their sample. Nearly half of their systems with R NLR 16,fit 5 kpc tend to be z ≈ 0.10 − 0.15 galaxies suggesting that the MaNGA observation's spatial resolution (∼ 2. 5) may be a concern.
Overall, our findings suggest that, as expected, the surface brightness of the NLR depends of the central AGN radiation power, but the NLR extension is mainly determined by geometrical factors. The exemplary cases are PG 0050+124 and PG 1426+015 where the NLR emission is detected abnormally bright at the outskirts of the host galaxies ( Figure 5). In the PG quasars, the NLR extension is probably limited by the availability of gas that can be ionized by the AGN radiation (matterbounded case), and not by the lack of AGN ionizing photons at large galactocentric radius (ionization-bounded case). This finding is in line with the reports of other examples of distant AGN-photoionized gas debris and filaments, as well as large-scale gas outflows located beyond the main disk of the host galaxy disk (e.g., Greene et al. 2012  . Azimuthal angular difference between the ionization cone side position angles for the PG quasar host galaxies with two-sided ionization cones. The dashed line indicates the ideal 180 • angle difference between both sides. The shaded zone represents the azimuthal angular size bin (22.5 • ) used to compute the C f values from the two-dimensional maps; for each ionization cone PA value we assume an uncertainty equal to half the bin size. The data are clustered around the expected azimuthal angle difference, consistent with the bipolar geometry of an ionization cone.

NLR Geometry
To explore the NLR geometry we analyze the ionization energy source maps in polar projection (Figure 7). Aiming to identify possible bipolar ionization cones, we first bin the maps in 16 azimuthal bins, and for each individual bin and we compute the covering fraction of the AGN-classified pixels as, where N AGN pix is the number of pixels classified as AGNlike according to the [N ii]/Hα map and as Seyfert-like according to the [S ii]/Hα map, and N tot pix is the number of pixels in which we detect any emission that can be classified in the BPT maps (Section 3.4). Figure 7 presents the trends of C f with respect to the azimuthal angle for both BPT classification schemes for example object PG 0923+129, while the rest of our sample is presented in Appendix B. For most of the systems we find that C f peaks at different azimuthal angle values in both BPT maps, but we caution that the presence of star-forming sub-structures (e.g., star-forming ring in PG 0923+129 or spiral arms in PG 0050+124) may outshine the AGN-induced emission and/or shape the NLR geometry. We focus on the C f values estimated from the [S ii]/Hα map, as this classification scheme better separates between emission produced by AGN photoionization and shocks. 7 Six of the objects show an increase of C f at particular bins of azimuthal angle, suggesting that the AGNs ionize gas along preferred position angles. For the remaining sources we do not find any significant trend with respect to azimuthal angle. One 7 As reviewed in Ho (2008), LINER emission far from the central region of galaxy can be powered by a variety of sources, including shocks. This should not be conflated with the proposition that LINERs, especially in low-luminosity AGNs, bear no relation to BH accretion. As this study concerns highly accreting, luminous quasars, it is reasonable to attribute LINER-like emission to nonnuclear sources of excitation, such as shocks, and only attribute to AGN photoionization those regions classified as Seyfert-like.
of these is PG 1426+015, a merger for which we cannot determine any NLR geometry given its complex morphology. For the two remaining systems PG 1011−040 and PG 1244+026, we speculate that the observed round geometry of its NLR may be consistent with the source being viewed through a LOS close to the direction of the ionization cone (pole-on view; e.g., Storchi-Bergmann et al. 2018). Observing larger samples of type 1 and type 2 AGNs would help to understand whether this possibility is consistent with the expectations implied by the AGN unification model from a statistical point of view. We characterize the sources presenting preferred AGN ionization direction by fitting the C f values using up to two Gaussians. During this modeling we account for the azimuthal binning implemented to estimate the C f values, and we also add a constant free parameter value that helps to deal with cases where few Seyfert-classified pixels are detected in all directions (e.g., PG 1126−041). We assume that the centroids of the Gaussian trace the position angles of the ionization cone sides, while the width of the Gaussian represents the opening angle of the cone (likely to be lower limits due to surface brightness bias). Hence, each Gaussian is assumed to represent one side of the ionization cone, and a one-sided ionization cone is well-represented by only one Gaussian component. The inferred location of the projected ionization cones is overlaid on the [S ii]/Hα map in Figures 7 and B1. We find two-sided ionization cones in PG 0050+124, PG 0934+013, PG 1126−041, and PG 1211+143, and more subtle evidence is suggested in PG 0923+129 and PG 2130+099. To determine if these two-sided ionization cones are consistent with a biconical geometry, in Figure 8 we show the difference between the ionization cone sides position angles (∆PA cone ) for these six quasar hosts. We find a mean position angle difference of ∼ 187 • , in close agreement with the expected difference of 180 • for a bipolar cone geometry. We measure a scatter of ∼ 23 • , which is comparable to the azimuthal bin size (22.5 • ) employed to tessellate the host galaxy maps. Admittedly, our NLR geometry analysis is limited to the systems where the ionization cones are misaligned with respect to the observation line-of-sight. While C f is conveniently defined to be immune to the host galaxy inclination angle, it may not be useful to study the NLR geometry in highly inclined objects.

Ionized Gas Outflows
The association of disturbed ionized gas kinematics with ionization cones has been interpreted as evidence for ionized gas outflows in active galaxies (e.g., Sun et al. 2017;Kang & Woo 2018). While ionization cones can ex-tend up to several kpc within the host galaxy or beyond, the ionized gas outflows tend to be primarily detected on smaller scales (e.g., Karouzos et al. 2016;Fischer et al. 2018;Husemann et al. 2019a). Our MUSE maps of PG quasar host galaxies offer us the opportunity to explore the possible contribution of outflows to AGN feedback. Outflows can prevent or delay star formation if they succeed in evacuating or relocating a significant fraction of the ISM (Fabian 2012). Conversely, outflows can promote star formation through gas compression by shocks (Cresci et al. 2015;Maiolino et al. 2017).

Local Observational Signatures
Visual inspection of the host galaxy maps (Figures 5  and A1) suggests that the lower Hα velocity dispersions (σ V,Hα ) tend to be measured in zones where the star-forming regions mainly contribute to the ionization state of the ISM (pixels enclosed by green contours in the σ V,Hα maps). Interestingly, in some sources (e.g., PG 0050+124 and PG 0934+013) the pixels with higher σ V,Hα are located in zones that overlap with the direction of the ionization cone. This is not the case when observing the [O iii] velocity dispersions (σ V, [OIII] ), where higher values do not spatially coincide with the high σ V,Hα estimates in some cases (e.g., PG 1126−041). However, the [O iii] emission line is observed in the bluer region of the MUSE spectra coverage implying that the low σ V, [OIII] estimates are less reliable when compared to the σ V,Hα values due to the MUSE LSF variation with wavelength, i.e., the [O iii] lines are observed at coarser spectral resolution and subject to a more uncertain line width deconvolution correction. Hence, in this study we focus on studying the apparent connection between the Hα velocity dispersion with respect to the ionization state of the gas. We also inspect the gas electron density (n e ) variation across the ionization cone direction, as this quantity may reflect the compression of the gas by shocks (Osterbrock & Ferland 2006). To this end, we take advantage of the respective sensitivity of both BPT maps, using the [S ii]/Hα map to isolate the effect of AGN ionization and the [N ii]/Hα map to highlight the star-forming regions. Figure 9 examines the deprojected radial variation of σ V,Hα and [S ii] λ6716/[S ii] λ6731 across the six detected ionization cones. We convert the ratio of the [S ii] doublet to n e adopting an electron temperature of 10 4 K following Eq. 3 of Proxauf et al. (2014). Consequently, we assume that the [S ii] emission only provides a reliable estimate of n e within the 40 − 10 4 cm −3 range. However, we caution that the n e values correspond to the line-of-sight integrated densities instead of local ISM estimates. Across all six ionization cones Pixel-wise Hα velocity dispersion (left) and [S ii] emission-line intensity ratio (right) across the axis of the detected ionization cones. In both panels the data are color-coded following the [S ii]/Hα maps of Figure 5, and the open black circles represent azimuthal average values computed from the pixels labeled as star-forming (according to the [N ii]/Hα classification scheme), if any. In the right panels, the unshaded region denotes the range for which the [S ii] λ6716/[S ii] λ6731 ratio is sensitive to ne (∼ 40 − 10 4 cm −3 ), with the corresponding ne scale given in the right axis, following the conversion of Proxauf et al. (2014). All objects exhibit regions with enhanced σV,Hα, although ne increases only in PG 0934+013 and PG 1126−041.
we find an increase in gas velocity dispersion up to σ V,Hα ≈ 100 − 200 km s −1 at specific radii. The only exception is PG 1126−041, whose σ V,Hα profile continuously increases as a function of radius. In all the cases, the regions associated with an increase of σ V,Hα tend to be ionized by the AGN, according to the [S ii]/Hα diagram classification. In contrast, the pixels classified as star-forming from the [N ii]/Hα map have relatively smooth radial trends in velocity dispersion in the range σ V,Hα ≈ 20 − 45 km s −1 (we omit PG 1126−041, which has too few star-forming pixels), ∼ 5 times lower than those in the AGN-classified locations. We find no clear trends in terms of [S ii] λ6716/[S ii] λ6731. PG 0934+013 is the only system in which an increase in σ V,Hα spatially coincides with a decrease in [S ii] λ6716/[S ii] λ6731 (increase of n e ); the most turbulent, AGN-dominated regions have densities 200 times higher than those in the star-forming locations (n e 40 cm −3 ). That σ V,Hα only increases at specific radii within the host suggests that the axis of the ionization cone is not pointed along the gas disk, but instead is slightly offset, such that it affects the ISM at large vertical scale heights where diffuse ionized gas lies. Such a disk-outflow geometric configuration may imply that the outflows are unable to significantly affect the gas located in the disk midplane where the star-forming complexes lie. The outflows do not dominate the ionization, so that the fluxweighted emission coming from co-spatial H ii regions mainly reflects the star formation process, and hence the Seyfert-classified regions are mainly detected where the Hα emission generated by the star formation activity is faint (see Figures 5 and A1). Instead of a homogeneously filled cone aligned with respect to the disk midplane, the outflow geometry may be more akin to that of an inclined expanding shell, which is capable of producing only a local enhancement in gas density and turbulence at the shock front interface with the gas clouds (e.g., expanding hot gas cocoon in Figure 5 of Husemann et al. 2019a). We caution that we cannot rule out the possibility that we are missing some inner faint emission and kinematic fingerprint induced by the outflow because the weak signal is swamped by brighter emission from star-forming complexes across the same LOS (e.g., PG 0923+129 and PG 0934+013). Our ability to detect kinematically subtle signatures toward the inner regions of the host galaxies is limited by the accuracy of our deblending technique.

Mass Outflow Rates and Energetics
We now proceed to compute the mass outflow rates (Ṁ out ) and kinetic injection rates (Ė kin ), following Fiore et al. (2017). Assuming that the kinematic features For the latter case, the vertical bar plotted in the bottom-right corner represents the systematic uncertainty associated with our assumption for the shell thickness. The dashed line shows the best fit reported by Fiore et al. (2017). We also show the mass outflow rates given by Baron & Netzer (2019) for a sample of less powerful AGNs, using their "log U " method. (Right) Same as the left panel but using ne estimated from the observed [S ii] doublet (Section 4.2.1). For PG 2130+099, we only present upper limits (assuming ne = 40 cm −3 ), which are highlighted by upward-pointing arrows attached to the corresponding data points. We do not show PG 1211+143 because it has an unconstrained value of ne for the outflow. In all cases the mass outflow rates are low,Ṁout 1.6 M yr −1 .
noted in Section 4.2.1 correspond to genuine outflow signatures, we isolate the emission associated with those pixels in order to quantify their outflow properties (e.g., Feruglio et al. 2020). We select the Seyfert-classified pixels (according to the [S ii]/Hα map) along the ionization cones that have σ V,Hα ≥ 50 km s −1 , which lies comfortably above the velocity dispersion threshold of the star-forming regions (open circles in Figure 9). For each individual pixel, we remove the parent Voronoi cell stellar continuum component ( § 3.1) from its spectrum (recall that we have already deblended the nuclear spectrum), and spatially collapse the spectra of all pixels to obtain a pure emission-line spectrum for each side of the ionization cone. To characterize the outflow spectra, we fit double-Gaussian profiles to the emission lines to account for line asymmetries. Following several authors in the literature, we compute the outflow velocity V out = V cone +2 σ cone , where V cone and σ cone are given by the best-fit Hα line centroid and line width (e.g., Rupke & Veilleux 2013;Fiore et al. 2017;Feruglio et al. 2020). The associated ionized gas mass is computed as M out = 3.4 × 10 6 100 cm −3 n e L Hα 10 41 erg s −1 M , (2) with L Hα the extinction-corrected Hα luminosity. Two sets of values for M out are calculated, one using n e derived from [S ii] λ6716/[S ii] λ6731, and another assuming n e = 200 cm −3 to make a consistent comparison with the quantities reported by Fiore et al. (2017). 8 We derive n e values in the range of ∼ 100 − 360 cm −3 with the exception of PG 1211+143 and PG 2130+099. For the former system we were unable to derive a reliable n e estimate, while for the later case we report n e = 40 cm −3 upper limits. For a homogeneously filled cone geometry, the mass outflow rate (Fiore et al. 2017) is given by, where R out is the radius at which the mass outflow rate is computed. We set this value equal to the mean radius at which we detect the enhancement of the Hα velocity dispersion. However, as we suggested in Section 4.2.1, a homogeneously filled cone geometry may not be the most suitable model to describe our data. If an expanding shell is a better approximation, where ∆R shell is a shell thickness that must be assumed (Husemann et al. 2019b). We use ∆R shell = 200 pc and adopt a range of 20-500 pc (Husemann et al. 2019b) to bracket the systematic error budget. Then, the kinetic energy injection rate is simply estimated asĖ kin = 0.5V 2 outṀout . Table 3 lists the estimates of M out ,Ṁ out , andĖ kin for the ionized gas outflows. Figure 10 presents the mass outflow rates as a function of the AGN bolometric luminosity, L bol = 10 L 5100 (Richards et al. 2006). We consider two separate cases: (1) adopting n e = 200 cm −3 , as assumed by Fiore et al. (2017) and (2) using the actually measured values or upper limits of n e . In both cases, we showṀ out computed assuming an outflow geometry of a homogeneously filled cone and a shell-like cone. We find that all six of the PG quasars have very low mass outflow rates (Ṁ out ≈ 0.008 − 0.5 M yr −1 ) when compared with the sample collected by Fiore et al. (2017) at similar bolometric luminosities (L bol ≈ 10 45 − 10 46 erg s −1 ), for mass outflow rates computed assuming identical n e and outflow geometry. Our values ofṀ out are similar to those reported for less powerful AGNs by Baron & Netzer (2019), who derived mass outflow rates based on the ionization parameter and AGN luminosity (their "log U " method). We also find similarṀ out values compared to those reported by Rojas et al. (2020), who used the nuclear [O iii] line shape and luminosity to estimatė M out for a sample of nearby X-ray-selected AGNs with bolometric luminosities comparable to those of our PG quasars.
The mass outflow rates derived here are extraordinarily sensitive to the assumptions of the outflow geometry, which is essentially unknown. For example, the mass outflow rates increase systematically by a factor of ∼ 10 if we adopt a shell-like outflow cone geometry, and, in view of the linear scaling with shell-thickness, another order of magnitude boost inṀ out can be achieved if ∆R shell ≈ 20 pc. In the face of the other sources of systematic uncertainty, by lowering the Hα velocity dispersion threshold to ∼30 km s −1 when selecting the Seyfertclassified pixels we find thatṀ out increases only by a factor of ∼ 4 on average, where both the increase of Hα luminosities and the decrease of the electron densities contribute evenly to this boost factor budget. On the other hand, the two assumptions of n e that we explored do not make a big difference. However, we should recognize that density does pose a major source of uncertainty. The choice of the methodology and tracer to calculate n e affect the values of mass outflow rate and quantities that depend on it (e.g., Baron & Netzer 2019;Riffel 2021). It is highly unlikely that a single density characterizes the line-emitting gas, nor is [S ii], which has a low critical density for collisional de-excitation, the most appropriate tracer for all the gas if a large range of densities exists. Baron & Netzer (2019) suggest that the [S ii] emission mainly arises close behind the ionization front where most of the gas is neutral and n e is effectively lower (see also Figure 10 of Davies et al. 2020). The bulk of the ionized gas is poorly traced by [S ii], and using this line doublet as an electron density tracer may overestimateṀ out (Revalski et al. 2022). This systematic uncertainty would also affect the mass outflow rates presented by Fiore et al. (2017) and showed in Figure 10. Indeed, Davies et al. (2020) suggest that theṀ out −L bol correlation for the more luminous AGNs weaken when correcting the data by n e (we refer to Davies et al. 2020, for more details).
Plotting the kinetic injection rate of the outflow versus the bolometric luminosity of the AGN (Figure 11), we findĖ kin /L bol 10 −3 , values that are systematically lower than those reported by Fiore et al. (2017) for AGNs of similar power and studied under the same model assumptions. Our values, however, are in better agreement with the estimates reported by Baron & Netzer (2019) and Rojas et al. (2020) for the less powerful AGNs and X-ray-selected AGNs with similar bolometric luminosities, respectively. Our findings suggest that ionized gas outflows withĖ kin /L bol 10 −3 are inefficient in terms of impacting the global ISM of the host galaxy. Note that our estimates ofṀ out andĖ kin /L bol should be regarded as upper limits if sources other than AGNdriven outflows contribute to the local enhancement in gas dispersion (σ V,Hα ). The same applies regarding the use of [S ii] as an electron density estimator (Baron & Netzer 2019;Davies et al. 2020;Riffel 2021;Revalski et al. 2022).

SUMMARY AND DISCUSSION
The role of AGN feedback on the evolution of galaxies remains a topic of intense debate. Theoretical arguments abound concerning the need for negative AGN feedback to solve a suite of otherwise intractable problems in galaxy evolution. Theoretical and numerical works argue that outflows with kinetic power of a few percent the energy radiated by the AGN (Ė kin /L bol ≈ 0.05 − 0.5) can significantly impact the evolution of the host galaxy (Di Matteo et al. 2005;Springel et al. 2005;Hopkins & Quataert 2010;Zubovas & King 2012;Hopkins et al. 2016). Indeed, literature estimates based on observations of ionized, molecular, and Xray tracers collected by Fiore et al. (2017) suggest that the kinetic power of outflows lie mainly in the rangė E kin /L bol 0.1, with only a few reaching values as high asĖ kin /L bol ≈ 0.1 − 1 (see also Shimizu et al. 2019;Kakkad et al. 2022). A major source of uncertainty pertains to the degree to which outflows actually couple to the cold gas reservoir of the host galaxy. If the outflow expands in a direction perpendicular to the disk plane of the galaxy, then the cold gas, which resides largely in the disk, would be hardly touched (Gabor & Bournaud 2014). But the accretion disk is not typically aligned with the host galaxy (Wilson & Tsvetanov 1994).
While there has been much tantalizing evidence of multi-phase outflows in AGNs reported in the literature (e.g., Cicone et al. 2014;Feruglio et al. 2015;Karouzos et al. 2016;Morganti 2017;Cicone et al. 2018;Fluetsch et al. 2019;Veilleux et al. 2020), we have a poor understanding of how efficiently the outflows actually affect the host galaxy (e.g., Harrison et al. 2018;Kakkad et al. 2022). The naive expectation that AGN feedback expels and depletes the cold gas reservoir is certainly not borne out by much of the evidence. A variety of observations have demonstrated that active galaxies, covering a wide gamut of accretion power and evolutionary states, are categorically not gas-deficient compared to inactive, star-forming of similar stellar mass (e.g., Bertram et al. 2007;Krips et al. 2012;Xia et al. 2012;Husemann et al. 2017;Shangguan et al. 2018Shangguan et al. , 2020bShangguan & Ho 2019;Jarvis et al. 2020;Yesuf & Ho 2020;Koss et al. 2021;Salvestrini et al. 2022). The cold gas shows no signs of being substantially disturbed in terms of its global kinematics (e.g., Ho et al. 2008;Shangguan et al. 2020b;Molina et al. 2021), nor does it have any trouble forming stars (e.g., Bernhard et al. 2019;Grimmett et al. 2020;Jarvis et al. 2020;Kirkpatrick et al. 2020;Xie et al. 2021;Zhuang et al. 2021).
mate the physical properties of the ionized gas outflows. Our main results are: • The NLR is spatially resolved for all the host galaxies. Down to an [O iii] surface brightness level of 10 −16 erg s −1 cm −2 arcsec −2 , the NLR has radial extensions in the range 1.6 − 4.7 kpc, which correlate with the total extinction-corrected [O iii] luminosity. In some cases, the [O iii] emission can reach distances beyond the physical extent of the host galaxy.
• The radiation field of the AGN ionizes the host galaxy ISM in preferred directions, producing an NLR morphology that resembles a bipolar ionization cone.
• The velocity dispersion of the ionized gas, as traced by Hα, can reach values 100 km s −1 at specific galactocentric radii (mostly ∼ 5 − 10 kpc) across the AGN ionization cones. These velocity dispersions are significantly higher than the average value ( 40 km s −1 ) measured from starforming regions located at similar distances from the galaxy center.
• Assuming that the local velocity dispersion enhancements reflect the interaction between the AGN-driven outflow and the host galaxy ISM, we compute ionized gas mass outflow ratesṀ out ≈ 0.008 − 1.6 M yr −1 and kinetic injection powerṡ E kin /L bol 10 −3 . These values are substantially lower than those expected from the empirical correlations of Fiore et al. (2017), but they agree better with the estimates of Baron & Netzer (2019) and Rojas et al. (2020) for AGNs with similar and lower bolometric luminosities.
Our IFU analysis of the present sample of PG quasars offers some clues on the conundrum of why the host galaxies of luminous AGNs in the local Universe seem to show so little evidence of negative AGN feedback. Although the NLRs follow the expected luminosity-size relation, implying that the AGN radiation can permeate across a solid angle not significantly misaligned with respect to the disk plane of the host, the emission detected across the host galaxy is faint and seems to be easily swamped by the Hα emission generated by the star formation activity (Figure 7). More importantly, the extremely low kinetic powers ofĖ kin /L bol 10 −3 suggest that the ionized gas outflows are in general very poorly coupled to the ISM of the host galaxy. This casts serious doubt on the efficiency of negative AGN feedback and its global impact on the host galaxies of nearby quasars. Although the present results only pertain to the ionized gas in a small number of objects, we note that Shangguan et al. (2020b) also placed strong limits on the absence of cold, molecular outflows for a subset of the low-redshift PG quasars. Koss et al. (2021) found similar results when analyzing the molecular gas content in a large sample (∼ 200) of nearby X-ray-selected AGNs.
Aside from the low effectiveness of negative AGN feedback, an intriguing characteristic revealed by the IFU observations is the prevalence of enhanced velocity dispersion for the ionized gas at certain zones located within the ionization cone. The increase in velocity dispersion is associated, at least in one instance (PG 0934+013; Figure 9), with an accompanying increase in electron density. A local increase in the velocity widths of the nebular emission lines has been observed also in the inner kpc-scale regions of nearby AGNs (Ruschel-Dutra et al. 2021), whose kinematic disturbance has been regarded as an indicator of outflows within the NLR (Kang et al. 2017;Sun et al. 2018;Husemann et al. 2019b). Jet-like radio emission has also been reported (Husemann et al. 2013;Villar-Martín et al. 2017;Santoro et al. 2020;Venturi et al. 2021). The available radio data (∼ 0.685 − 5 GHz) for the PG quasars suggest unresolved ( 0. 5) low-powered ratio jet emission in PG 1011−040, PG0934+013 and PG1426+015, while non-conclusive reports can be made for the other host galaxies as the radio fluxes could also be produced by the optically thin synchrotron emission from AGN-and/or starburst-driven winds (Silpa et al. 2020). We interpret the localized velocity dispersion enhancements as a signature of the shock front interface between the ISM and an expanding, shell-like outflow propagating in a direction offset from the disk midplane of the host galaxy.

ACKNOWLEDGMENTS
We thank to the anonymous referee for her/his helpful comments and suggestions.
We acknowledge support from the National Science Foundation of China (11721303, 11991052, 12192220, 12192222, 12011540375) Figure A1. Two-dimensional moment maps and ionization diagnostic diagrams for the PG quasar host galaxies. The panels are organized and color-coded following Figure 5. The host galaxy has two extended and clumpy spiral arms that extend out to the outskirts of the stellar disk. The emission coming from the spirals can be unambiguously associated with underlying star formation activity according to the BPT maps. On the other hand, the emission coming from the inter-arm regions are much fainter and likely produced by the AGN ionization and/or shocks. The LOS velocity map clearly shows a disk-like rotational pattern, while the velocity dispersion map suggests more complex ISM dynamics. The spiral arms can be associated with regions that show lower values of σ V,Hα within the host galaxy, while the inter-arm regions tend to show larger σ V,Hα .
PG 0923+129: The Hα emission exhibits an inner disk structure from which two clumpy spiral arms extend out to the outskirts of the host galaxy, which shows ubiquitous star formation activity. The inner Hα emission seems to trace a ring-like structure. This is interesting, as the galaxy appears to be an S0 system when seen in HST or MUSE broad-band images. A main region with Seyfert classification seen toward the south-east suggests the presence of an ionization cone. The velocity map shows the characteristic rotation pattern of a disk galaxy, while the velocity dispersion map suggests that star formation occurs in zones with lower σ V,Hα .
PG 0934+013: Two outer, clumpy, spiral armlike regions can be clearly seen in the Hα map. The central zone, where a bar is present, also emits Hα, which, according to the BPT maps is powered by star formation activity. The zones ionized by the AGN are mainly distributed from the north-east to the south-west in the form of a biconical ionization cone. The [S ii]/Hα BPT diagram, however, suggests that LINER-like emission is ubiquitous across the host galaxy. Similar to PG 0050+124 and PG 0923+129, regions with LINER characteristics have higher σ V,Hα than those dominated by star formation. The LOS velocity map shows a clear disk-like rotation pattern.
PG 1011−040: The Hα emission traces two main zones, a nuclear region surrounded by ring-like clumpy structure. Five large star-forming clumps are clearly seen inside the ring. The LOS velocity field is complex, but the velocity gradient varies smoothly across the field-of-view, sug-gesting that both structures are spatially connected. This is also supported by the rest-frame I-band HST/WF3 image and MUSE white-light image, which reveals a stellar bar crossing from the ring-like structure from one extreme to the other through the nuclear region. The central zone has larger σ V,Hα .
PG 1126−041: This system clearly shows a disklike morphology and rotational pattern traced by the Hα emission. A clumpy morphology is seen toward the galaxy outskirts, while in the central zone the Hα emission seems to be distributed in a bar-like or outflow structure. The velocity dispersion map is complex, with an enhancement of σ V,Hα toward the north-west and south-east, zones that are ionized by the AGN, although no clear substructures are visible in the Hα and [O iii] intensity maps. Ongoing star formation activity is necessary to explain the BPT maps.
PG 1211+143: The MUSE observation shows a compact system in which some evidence of rotation can be appreciated from the LOS velocity map. The data suggest two major clumps toward the north of the galaxy; they are consistent with star formation once the [N ii] upper limits are considered. The lack of [N ii] but the detection of Hα, Hβ, and [O iii] suggests that these two zones are metal-poor.
PG 1244+026: The Hα emission traces a compact ionized gas distribution and a star-forming clump present toward the south. Interestingly, while the central emission of the clump is consistently classified as "composite," it is surrounded by pixels with "Seyfert" classification according to the [N ii]/Hα map. This suggests that PSF blurring is inducing emission-line flux mixing. A rotation pattern is evident in the LOS velocity map. No clear trend is seen in the velocity dispersion map.
PG 1426+015: The host galaxy is an ongoing merger (Kim et al. 2017). Although the Hα emission map looks patchy, we can see clearly three different structures: a large tidal feature toward the east, a second tidal feature toward the west, and a central, rotating disk-like component according to the Hα LOS velocity map. The BPT maps suggest that star formation activity is probably taking place within the tidal features and in some small clumps located in the central disk-like component. The higher σ V,Hα observed toward the eastern zone of the inner disk-like component may suggest the presence of an outflow component and/or shocks affecting the ionized gas, judging by the [S ii]/Hα diagram.
PG 2130+099: The Hα emission traces an outer spiral-like clumpy structure that extends out to the outskirts of the galaxy. The clumpy regions are consistently classified as star-forming zones. A central clumpy disk-like region is seen toward the center. The [N ii]/Hα map suggests ongoing star formation, although shocks cannot be completely discarded, according to the [S ii]/Hα map. The LOS velocity maps clearly show a rotating disk-like pattern. However, the velocity dispersion map suggests more complex dynamics. While σ V,Hα increases overall toward the central disk-like region, higher velocity dispersions are seen in the zones between the spiral and the central disk-like region.