© The Authors 2024

1. Introduction

The evolution of massive galaxies is related to the accretion of matter onto super massive black holes (SMBHs). SMBHs are thought to self-regulate gas accretion by releasing energy, which has an impact on the host galaxy (e.g., Fabian 2012), as gas cannot cool down efficiently to form stars (Dubois et al. 2016). This feedback from active galactic nuclei (AGN) is required in cosmological simulations to regulate the growth of galaxies and match the observations (e.g., Schaye et al. 2015; Nelson et al. 2018). Outflows of gas are a manifestation of feedback, and they are capable of modifying the distribution of gas in the central hundreds of parsecs of galaxies (e.g., García-Burillo et al. 2021; Ramos Almeida et al. 2022), and altering the chemical enrichment of the circumgalactic medium up to hundreds of kiloparsecs from the SMBHs (e.g., Martín-Navarro et al. 2021).

The main feedback channels are the so-called quasar or radiative mode and the radio or kinetic mode. The first is supposed to act in AGN with high accretion rates (Fabian 2012), while the second one is generally associated with powerful radio galaxies, with a lower accretion rate, where the gas is mechanically compressed and accelerated by jets (McNamara & Nulsen 2007). However, this distinction of the two main feedback mechanisms is too simplistic, and both modes can act simultaneously. In fact, in radio-quiet type-2 quasars (QSO2s, L_[O III] > 10^8.5 L_⊙; Reyes et al. 2008) low-power jets have been proposed to accelerate outflows (e.g., Jarvis et al. 2019; Girdhar et al. 2022, 2024; Speranza et al. 2022; Venturi et al. 2023), even though the quasar mode should dominate due to their high accretion rates. In this context, QSO2s are good laboratories to test whether compact low-power jets can efficiently drive gas outflows (e.g., Audibert et al. 2023). This can be explored by inspecting the spatial distribution and geometry of the radio emission with respect to the outflowing gas and comparing their energetics.

To detect signatures of ionized outflows, it is common to characterize the [O III]λ5007 Å emission line which, being a forbidden transition, never has a broad component from the broad line region (BLR) and therefore is a good tracer of the kinematics of the narrow line region (NLR), from the central parsecs to kiloparsec scales (e.g., Osterbrock 1989). Any broad [O III]λ5007 Å profiles with prominent and asymmetric wings can be associated with turbulent and fast outflowing gas that is not ordinarily rotating, as claimed in early (e.g., Heckman et al. 1981; Feldman et al. 1982; Whittle 1985) and recent studies (e.g., Cresci et al. 2015; Brusa et al. 2015; Carniani et al. 2015; Venturi et al. 2018; Speranza et al. 2021; Peralta de Arriba et al. 2023). These ionized outflows are commonly observed in nearby QSO2s (e.g., Liu et al. 2013; Mullaney et al. 2013; Harrison et al. 2014; Villar Martín et al. 2014; Fischer et al. 2018) where the contamination from the continuum and the BLR is naturally masked by nuclear dust (Ramos Almeida & Ricci 2017), although large-scale obscuration might be important as well.

However, the ionized emission represents only a single phase of the gas, and to fully characterize the impact of outflows a multiphase characterization is needed (Cicone et al. 2018). From observational scaling relations, the cold molecular outflows are expected to drive more massive outflows than the ionized ones (e.g., Fiore et al. 2017) and the outflow properties, including densities and velocities, have been observed to be different depending on the gas phase (e.g., Finlez et al. 2018; Herrera-Camus et al. 2019; Fluetsch et al. 2021; García-Bernete et al. 2021). However, the majority of works are based on a single outflow phase, providing a partial view of AGN feedback in which the interplay between the different phases of the gas is still to be unveiled.

Here we present new integral field unit (IFU) data of five QSO2s observed with the MEGARA¹ instrument of the 10.4 m Gran Telescopio Canarias (GTC). This instrument provides optical spectra in which, at the redshift of the targets (z ∼ 0.1), the [O III]λ5007 Å emission line is detected at a high spectral resolution (R ∼ 5850 using the low-resolution grating). These QSO2s have already been observed with the Atacama Large Millimeter/submillimeter Array (ALMA) at 0.2″ resolution and their cold molecular gas kinematics have been characterized by Ramos Almeida et al. (2022) using the CO(2–1) emission line. These authors reported the presence of cold molecular outflows for the five QSO2s, four of them spatially resolved with projected radii of ∼0.2–0.7″ (∼0.4–1.3 kpc). Here we aim to study the ionized kinematics of the QSO2s and compare them with their cold molecular gas counterparts. The paper is organized as follows: in Sect. 2 we present the main properties of the sample. Section 3 includes the details on the observations and data reduction. In Sect. 4 we show the methodology applied to analyze the nuclear and extended emission of the QSO2s and to measure their outflow properties, and we describe the results. In Sect. 5 we discuss the implication of our findings, which we summarize in Sect. 6. Throughout this work we assume the following cosmology: H₀ = 70.0 km s⁻¹ Mpc⁻¹, Ω_M = 0.3, and Λ = 0.7.

2. Sample selection and characteristics

The five QSO2s studied here are part of the Quasar Feedback (http://research.iac.es/galeria/cra/qsofeed/) sample (Ramos Almeida et al. 2022; Pierce et al. 2023), selected from the Reyes et al. (2008) compilation of narrow-line AGN to have L_{[O III]} > 10^8.5 L_⊙ and redshift of z < 0.14. These selection criteria lead to a final number of 48 QSO2s with an average stellar mass and standard deviation of M_* = 10^{11.1 ± 0.2} M_⊙, obtained from the 2MASS Extended Source Catalogue (XSC) K-band magnitudes (Ramos Almeida et al. 2022). Our subset of five QSO2s have redshifts 0.08 ≤ z ≤ 0.12 and their main properties are listed in Table 1. We selected these five QSO2s because they were detected in the CO(2–1) ALMA observations reported in Ramos Almeida et al. (2022). Another two QSO2s, hosted in red early-type galaxies, were observed with ALMA at the same angular resolution, but they were not detected in CO(2–1) and therefore have molecular gas masses < 7 × 10⁸ M_⊙. Thus, the five QSO2s studied here are representative of the gas-rich QSO2s in the QSOFEED sample, as they have total molecular gas masses of 4–18 × 10⁹ M_⊙. Similar gas masses, of ∼1 × 10¹⁰ M_⊙, were reported by Jarvis et al. (2020), based on APEX CO(2–1) observations a sample of QSO2s at z < 0.2 hosting kpc-scale ionized outflows and jets.

In Fig. 1 we show how our subset of five QSO2s compares with the whole QSOFEED sample in terms of luminosity distance (D_L), bolometric luminosity, optical morphology, 1.4 GHz luminosity, and [O III] kinematics (W₈₀ and vel₀₅; i.e., velocity dispersion and outflow velocity).

Fig. 1. Orange and red histograms correspond to all QSO2s in the QSOFEED sample and to the five QSO2s studied here. From left to right and from top to bottom: luminosity distance, AGN luminosity, optical morphology (L=lenticular, E=elliptical, S=spiral, M=merger, I=interacting), integrated radio luminosity from FIRST 1.4 GHz data, W80, and vel05 of the [O III] line profile measured from SDSS optical spectra of the targets (Bessiere et al., in prep.). SDSS J134733.36+121724.3 was excluded from the W80 and vel05 histograms because of its extreme kinematics (W80 = 2500 km s−1 and vel05 = −2500 km s−1).

The luminosity distances of the five QSO2s are intermediate compared with those of the QSOFEED sample, and the bolometric luminosities, derived from the extinction-corrected [O III] luminosities from Kong & Ho (2018) by applying the correction factor of 474 from Lamastra et al. (2009), are among the highest (log L_bol = 45.5 − 46.0 erg s⁻¹; see top panels of Fig. 1). The five QSO2s have diverse optical morphologies, with three of them showing signatures of galaxy interactions and mergers. In this context, they are representative of the QSOFEED sample, which includes at least 65% of interacting and merging galaxies (Pierce et al. 2023). As can be seen for the middle left panel of Fig. 1, we lack undisturbed lenticular and elliptical galaxies.

Our subset of QSO2s have radio luminosities that are intermediate (log L_1.4 GHz = [23.7,24.4] W Hz⁻¹) within the range probed by the QSOFEED sample, although the majority of the QSO2s have radio luminosities in the range [22.5,23.5] W Hz⁻¹. It is noteworthy that four of the QSO2s in the QSOFEED sample qualify as radio-loud objects, having log L_1.4 GHz > 24.5 W Hz⁻¹ (see the middle right panel of Fig. 1) and L_1.4 GHz/L_{[O III]} ratios above the Xu et al. (1999) division. The other 44 QSO2s in the QSOFEED sample are radio-quiet, including the five QSO2s studied here. The relatively high radio luminosities of the latter put them above the radio-FIR correlation of star-forming galaxies (Bell 2003). According to Jarvis et al. (2019), in the case of the four QSO2s that we have in common (J1010, J1100, J1356, and J1430), less than 7% of their radio emission is associated with star formation, with nonthermal AGN emission being responsible for the rest. Ramos Almeida et al. (2022) estimated the SFRs of the QSO2s from their FIR luminosities, following Kennicutt (1998), corrected to a Chabrier initial mass function (i.e., dividing by a factor of 1.59; Chabrier 2003). These SFRs are supposed to be less affected by the AGN contribution than estimations from mid-infrared or radio, and therefore we adopt these values here. The SFRs range between 12 and 69 M_⊙ yr⁻¹, with uncertainties of ∼0.3 dex, and they are comparable with the values reported by Jarvis et al. (2019) for J1010, J1100, J1356, and J1430 (8–84 M_⊙ yr⁻¹), measured from the IR luminosity due to star formation, derived from spectral energy distribution (SED) fitting.

Finally, the bottom panels of Fig. 1 provide a comparison of the [O III] kinematics of the QSOFEED and MEGARA samples. To ensure a fair comparison, we fit the [O III] lines detected in SDSS spectra of the 48 objects in the QSOFEED sample, including the five MEGARA QSO2s, using the same methodology presented in Sect. 4.3, where we define W₈₀ and vel₀₅. The MEGARA QSO2s are at the high-velocity dispersion and high-velocity end of the two histograms, particularly in the case of W₈₀. Therefore, the five QSO2s studied here are representative of the QSO2s having the highest AGN luminosities and the most extreme ionized gas kinematics in the QSOFEED sample.

The ionized outflows of the five QSO2s have been studied using data from different facilities operating in the optical or near-infrared (NIR): the Gemini-South GMOS IFU (J1010, J1100, J1356, and J1430; Harrison et al. 2014), the Low-Dispersion Survey Spectrograph (J1356; Greene et al. 2012), HST/STIS (J1100; Fischer et al. 2018), VLT/SINFONI (J1430; Ramos Almeida et al. 2017), and GTC/EMIR (J1509; Ramos Almeida et al. 2019). In the case of the last two datasets, the NIR spectra enable the simultaneous characterization of the ionized gas emission from the Paα and [Si VI]λ1.963 μm lines and the warm molecular gas from the H₂ lines. However, only J1509 shows a clearly detected warm molecular outflow. Finally, Ramos Almeida et al. (2022) reported signatures of cold molecular outflows in all the QSO2s, four of them spatially resolved (J1010, J1100, J1356, and J1430).

3. MEGARA observations and data reduction

The sample was observed during four nights on 2021, between March and April (see Table 2; program ID: GTC62-20B, PI: C. Ramos Almeida), using the optical IFU mode of MEGARA that is installed on the Nasmyth A focus of the GTC. The IFU Field-of-view (FOV) of MEGARA is 12.5″ × 11.3″ (∼23.0 × 20.8 kpc² at z = 0.1) and it is composed by 567 hexagonal fibers, each one circumscribed in a circle of 0.62″ in diameter. Additional 56 fibres, divided in 8 minibundles, are located at 1.75–2 arcmin from the center of the FOV to simultaneously measure the sky background during observations. We made use of one of the low resolution volume phase holographic gratings (VPH570-LR; hereafter LR-V), which has a wavelength coverage of 5144–6168 Å with a nominal spectral resolution of R = 5850 (∼5 km s⁻¹) at the central wavelength of λ_c = 5687.63 Å. Using our own dataset, we measured a median FWHM = 0.9 Å from the ArXe lamps that corresponds to a spectral resolution of ∼45 km s⁻¹. We use the corresponding values measured for each source from the lamps as our instrumental broadening.

The integration times for each target are reported in Table 2. We took three exposures of a maximum of 590 s of integration to reduce the impact of cosmic rays. During the entire observation run the conditions were clear (i.e., less than 10% of the sky covered by cirrus), except for a photometric night when J1010 was observed. We calculated the seeing full width at half maximum (FWHM) using the photometric standard stars observed just before each target. The seeing was FWHM = 1.1″ (∼2 kpc) for the whole sample, except for J1100, which has a slightly lower value of 1.0″ (see Table 2).

The data reduction was performed using the MEGARA pipeline (Pascual et al. 2020, 2023). This pipeline performs bias subtraction, bad-pixels masking, tracing fibres, wavelength calibration (including versions corrected and uncorrected of barycentric velocity), flatfield correction, flux calibration, sky subtraction, and cube reconstruction. As recommended by the pipeline developers, the final data cubes were produced with spaxels of 0.3″ in size.

4. Methodology and results

4.1. The nuclear [O III] emission

In order to study the ionized gas kinematics we firstly extracted the nuclear spectra of the five QSO2s in a circular aperture of 1.2″ (∼2.2 kpc) in diameter. We did this to match the angular resolution set by the seeing size (see Table 2). Thanks to the LR-V wavelength coverage of MEGARA, Hβ and the [O III]λλ4959,5007 Å doublet are detected at the redshift of our sample. Considering that the [O III]λ5007 Å line is the most prominent, we extracted the nuclear spectra centering them on its peak, that always coincides with the continuum peak except in J1356, where it is off-setted by one spaxel.

We characterized the [O III] doublet by simultaneously fitting multiple Gaussians after subtracting the contribution of the continuum. We selected the continuum in the blue and red side adjacent to the [O III] doublet, so as to exclude any flux contribution from the emission lines, and we performed a linear fit of the selected regions to then subtract it. The width of each of these bands and the distance from the line peak depends on the broadness of the doublet and on the proximity of the Hβ emission line. The multiple Gaussians were fit using the same in-house Python program described in Sect. 3.1 of Speranza et al. (2022), that is based on the modeling module of Astropy (Astropy Collaboration 2013, 2018). For this work, we select the best fit as the model with the minimum number of Gaussian components for which the reduced-χ² improves more than 10% by adding the last component, as in Bessiere & Ramos Almeida (2022). We fit two Gaussians in the case of J1509, four in J1010, J1100, and J1430, and five in J1356. Data points having a standard deviation higher than 10 times the standard deviation of the residuals were removed (see the fits of J1356 and J1509 in Fig. 2). This threshold was selected by performing a preliminary fit of the [O III] doublet from which we evaluated the points clearly detached from the line profiles, which most likely correspond to hot pixels.

Fig. 2. Fits of the [O III] doublets observed in the nuclear spectra extracted with an aperture of 1.2″ in diameter. The observed spectra are shown in blue, the model including all the components in black, and the residuals from the fit in brown. Same colors correspond to the same physical components: narrow (n; orange dashed), intermediate (i; green dot-dashed), and broad (b; red solid). Noisy pixels not considered in the fits are marked with red crosses.

In all the five QSO2s the best fits include at least two Gaussians: a narrow component with FWHM between 150 and 460 km s⁻¹, and a broad (FWHM ≥ 1300 km s⁻¹) and component blueshifted from the systemic velocity set by the narrow component(s) (velocity shift, v_s, between −40 and −430 km s⁻¹). In four of the QSO2s we also detect intermediate components with FWHM = 340–890 km s⁻¹, most of them blueshifted (between −20 and −370 km s⁻¹) but also redshifted (26 and 280 km s⁻¹ in J1430 and J1356, respectively). We note that J1356 has an intermediate component that is narrower than the narrow components of J1100 and J1509 (see Table 3). Here we define the intermediate components as those being at least 100 km s⁻¹ broader in FWHM than the narrowest component fit to the same line, but still significantly narrower than the broadest component. All the fits are shown in Fig. 2 and the parameters are listed in Table 3.

In general, we associate the narrow components with gas in the NLR. In all cases but J1509 we resolve two narrow components that might correspond to the approaching and receding sides of a rotating gas distribution (see Sect. 4.3), although some of the QSO2s show complex kinematics (J1356 and J1430), suggesting that ongoing mergers and/or the presence of jet-ISM interactions might be altering or preventing gas rotation. On the other hand, the broad components correspond to turbulent outflowing gas, since they have FHWM ≥ 1300 km s⁻¹ and they are blueshifted with respect to the systemic velocity. We note that these blueshifts are small in the case of the nuclear spectra because we are looking at the integrated spectrum of the central 1.2 arcsec of the galaxies, but taking advantage of the spatial information provided by the IFU, we measure velocities of ≥800 km s⁻¹ for all the QSO2s (see Appendix A). Therefore, we detect the previously reported ionized outflows in the five QSO2s using this high spectral resolution dataset. The case of the intermediate components is more ambiguous since in some cases their FWHMs and v_s are consistent with those of the NLR (e.g., J1356), and in others with outflowing gas (e.g., J1100). For this reason in Sect. 4.3 we use a nonparametric method to measure the outflow properties, rather a parametric one. By doing this we do not have to ascribe any physical meaning to each kinematic components.

In Table 3 we report the results from the [O III] fits performed using Gaussians. In particular, we list the FWHM, the velocity shift from the narrow component (v_s), and the integrated flux. In the case of double peaked profiles (i.e., two narrow components detected with comparable fluxes), the systemic velocity corresponds to the flux weighted value between the two peaks (J1010 and J1430). Instead, when a dominant narrow component is fit, it is used to obtain the systemic velocity (J1100, J1356, and J1509). In Table 3 and throughout the text, we refer only to the parameters of [O III]λ5007 Å, since it shares the same number of components, each with identical FWHM and v_s, and an intensity ratio of 2.98 with [O III]λ4959 Å (Storey & Zeippen 2000), as they are part of the same doublet. The uncertainties of the parameters were measured by means of a Monte Carlo simulation. We simulate mock spectra by varying the flux of each spectral element of the emission line profile by adding random values extracted from a normal distribution with an amplitude given by the noise of the pixels. The uncertainties of the parameters are then computed as 1σ of each parameter distribution obtained from 100 mock spectra.

4.2. A spatially resolved analysis

From the analysis of the nuclear spectra presented in Sect. 4.1, we find signatures of ionized outflows in the five QSO2s. However, before exploring the kinematics of the gas across the MEGARA FOV, it is important to evaluate whether the [O III] emission is spatially resolved or not. In seeing-limited observations such as ours, beam smearing can replicate the nuclear spectrum throughout all the spaxels included within the seeing disk.

To do so we used the PSF-subtraction method proposed in Carniani et al. (2013) and applied in several studies (e.g., Carniani et al. 2015; Kakkad et al. 2020, 2023). This method consists of comparing the observed spectra across the MEGARA FOV with the spectrum extracted at the AGN position. The emission line profiles are expected to change with respect to the central one if the [O III] emission is spatially resolved. Therefore, we model the spectrum of the central five spaxels and we subtract it spaxel by spaxel, after applying a normalization factor to scale the amplitude and to minimize the residuals between the nuclear fit and the observed spectra of the other spaxels. The normalization factor is measured in the entire wavelength range covered by the model (see the right panels of Fig. 3). Before applying this procedure, we masked spaxels having less than three points in the spectrum detected at < 3σ over the continuum leftward and redward of the [O III] peak. Finally, we generate the cube of the residuals collapsing the entire wavelength range shown in the right panels of Fig. 3 and we inspected them to see whether they are close to zero or show a significant excess.

Fig. 3. Residuals obtained from the PSF-subtraction analysis performed across the MEGARA FOV. The left panels show the residual maps generated by subtracting the normalized model of the central spectrum, extracted in the five spaxels around the peak of the [O III] emission (black or white crosses), from the spectra of all the spaxels. North is up and east to the left. The black squares indicate the excess regions (≥3σ and in 0.6″ × 0.6″ apertures) from which we extract the corresponding spectra shown in the right panels. The dot-dashed and the dashed lines correspond to the [O III]λ4959 and 5007 Å centroids of the central model. For J1509, the grey band in the bottom right panel indicates bad pixels that were excluded when collapsing the spectral range shown on each panel to build the residual maps.

In Fig. 3 we show the residual maps (left panels) and the corresponding spectra (right panels) extracted from the excess regions indicated with black boxes in the corresponding maps. The residual maps are in ±3σ units, where the sigma corresponds to the standard deviation of the background measured in the edges of each source’s FOV, using an aperture of 1.5″ × 1.5″. We consider the [O III] emission resolved when residual structures are observed at ≥3σ. Instead, noisy maps with corresponding low-signal spectra are indicative of unresolved emission. This is the case of J1010, whose [O III] emission appears spatially unresolved in our data. Instead, we find it to be spatially resolved in J1100, J1356, J1430, and tentatively resolved in J1509. Therefore, in the following, we only consider outflow measurements for the QSO2s where the [O III] emission is resolved, although in Sect. 4.3 we describe the kinematics maps of J1010 as well.

4.3. The [O III] kinematics of the extended emission

In order to characterize the [OIII] kinematics of the QSO2s, we performed a nonparametric analysis of the line profiles by measuring the velocity percentiles, that is, velocities associated with a percentage of the flux of the emission line (e.g., Whittle 1985; Harrison et al. 2014; Speranza et al. 2021; Bessiere & Ramos Almeida 2022). We adopt this approach to automatize the fitting procedure across the MEGARA FOV and to overcome the complexity of the line profiles. In fact, some of the QSO2s are part of merging/interacting systems (e.g., J1356) showing a wide diversity of [O III] profiles at different distances from the AGN, that makes it challenging to obtain an acceptable physical description using a multicomponent model. To accomplish the nonparametric approach we fit multiple Gaussians to describe the emission line profiles, as we did in Sect. 4.1, but without ascribing any physical meaning to the single components. To perform the fits we use the nonlinear least-squares minimization and curve fitting for Python (LMFIT, Newville et al. 2014, 2016). We fit from a minimum of one to a maximum of six Gaussian components to each line of the [O III] doublet, including a new Gaussian to the fit only when the reduced χ² improves more than 10% with respect to the previous fit. As for the nuclear analysis presented in Sect. 4.1, we fix the kinematics and flux ratios (i.e., [O III]λ5007/[O III]λ4959 = 2.98) between the lines of the doublet. We set the width of the Gaussians to be larger than the instrumental broadening and to have amplitudes ≥3σ, in order to avoid bad fits. Once we obtain the best fit for the [O III] doublet in a given spaxel, we remove the model of the [O III]λ4959 line to obtain an isolated model of the [O III]λ5007 emission line. An example of the fitting procedure is shown in the left panel of Fig. 4 (for a random spaxel of J1430), with the corresponding [O III] model shown in the right panel. In order to avoid spaxels with poor or non detection of the [O III] line during the fitting process, we excluded those having less than three points in the spectrum detected at < 3σ over the continuum leftward and redward of the emission line peak. This is the same criterion used for the PSF-subtraction analysis described in Sect. 4.2.

Fig. 4. Example of nonparametric analysis of a random spaxel of J1430. Left panel: [O III]λλ5459,5007 Å continuum subtracted emission lines, with the observed spectrum shown in blue, the best fit model in black, and the corresponding residuals in grey. The same number of Gaussians (four in this case) with the same kinematics and a flux ratio of 2.98 are fit for each line of the doublet. Right panel: corresponding model of the [O III]λ5007 Å line, obtained from the fit shown in the left panel, with W80 indicated as the pink shaded region, vel05 and vel95 indicated by the pink vertical lines, and the upper (vel02 and vel98) and lower bounds (vel10 and vel90) marked by the dot-dashed vertical gray lines. The dashed vertical black line represents the velocity at the [O III] peak (velp).

Applying this process spaxel-by-spaxel, we can extract the parameters used to build the kinematics maps of [O III]λ5007. These parameters are shown in the right panel of Fig. 4 and they are the following.

• The velocity peak (vel_p), which is the velocity shift of the peak of the line from the systemic velocity measured in the nuclear spectrum (i.e., from λ_c in Table 3). vel_p mainly traces the narrow component, as can be seen from the top right panel of Fig. 5.

• The 5th and 95th velocity percentiles (vel₀₅ and vel₉₅), which are the velocities including the 5% and 95% of the emission line flux. They are representative of the blueshifted and redshifted outflow velocities, respectively (see middle panels of Fig. 5).

• which is the line region including 80% of the flux (W₈₀ = vel₉₀−vel₁₀). This parameter is indicative of how turbulent the gas is and, in the case of a line profile that can be described with a single Gaussian, it approximately corresponds to 1.09 × FWHM (see bottom left panel of Fig. 5). We note that the values of W₈₀ have not been corrected from instrumental broadening.

In Fig. 5 we present the kinematic maps of J1430, and the rest of the maps are presented and described in Appendix A (see Figs. A.1–A.4). The top left panel of Fig. 5 shows the [O III]λ5007 integrated flux. Due to the peculiar morphology of the ionized gas, J1430 has been nicknamed the Teacup (Keel et al. 2012). This morphology has been also detected in the radio (Harrison et al. 2015) and in the X-rays (Lansbury et al. 2018). The ionized bubble reaching 12 kpc northeast from the nucleus is detected in our map, as well as part of the southwest structure. As shown in the top right panel, negative velocities are seen to the southwest and positive velocities to the northeast, up to ±300 km s⁻¹. This is in agreement with previous results (Gagne et al. 2014; Harrison et al. 2014; Ramos Almeida et al. 2017; Villar-Martín et al. 2018; Moiseev & Ikhsanova 2023; Venturi et al. 2023). In the case of this QSO2, strong jet-ISM interaction is inducing enhanced turbulence in both the molecular and ionized gas, and the ionized and molecular major axes do not coincide (Audibert et al. 2023; Venturi et al. 2023), possibly due to the past merger even that J1430 experienced. Recently, Moiseev & Ikhsanova 2023 performed kinematic modeling of the ionized gas and claimed that most of it is rotating. However, another possibility is that the gas that we are tracing with vel_p is part of the outflow itself (see below).

Fig. 5. Kinematic maps obtained from the nonparametric analysis of the [O III] emission line profiles of J1430. Top panels: integrated flux on the left, with the ∼1″ resolution 6 GHz VLA contours at 3,5,10,20, and 60σ from Jarvis et al. (2019) superimposed (white solid lines), and velp on the right, which traces the narrow component. Middle panels: vel05 (left) and vel95 (right), which trace blueshifted and redshifted outflowing gas. The same 6 GHz contours shown in the flux map are shown in cyan, and the ∼0.25″ resolution 6 GHz VLA contours at 3, 5, 15, 30, and 60σ from Jarvis et al. (2019) are shown in solid white. The black contours correspond to the brightness temperature ratio (T32/T21) from ALMA CO observations at 0.5″ resolution (Audibert et al. 2023). Bottom panels: W80 on the left, with the same contours as in the vel05 map, and the number of Gaussians fit in each spaxel on the right. White and black crosses indicate the AGN location, defined as the maximum of the [O III] emission. North is up and east to the left.

The middle panels of Fig. 5 show high values of blueshifted and redshifted outflow velocities in the northeast and southwest direction respectively, with vel₀₅ of up to ∼–1000 km s⁻¹ and vel₉₅ of up to ∼800 km s⁻¹. These high velocities are detected within the same region where W₈₀ increases over ∼600 km s⁻¹ (see the left bottom panel of Fig. 5). The maximum values of W₈₀ are co-spatial with the most blueshifted velocities (i.e., to the northeast). This is in agreement with Harrison et al. (2014, 2015) and Ramos Almeida et al. (2017). Keel et al. (2017) reported blue [O III] wings eastward of the nucleus and red wings westward, with velocities of up to ±1000 km s⁻¹, using GMOS IFS observations. Thanks to the higher spectral resolution of our data, to the nonparametric analysis we performed, and to the larger FOV of MEGARA, which is more than twice as big as the FOV of the GMOS IFU data used by Harrison et al. (2014) and Keel et al. (2017), we clearly resolve the approaching and receding sides of the outflow in J1430, in the form of expanding shells/bubble. It is also possible that the reverse velocity pattern that we see in the vel_p map of J1430 corresponds to the inner and outer sides of northeast and southwest hollow cones, respectively, but here we do not perform any kinematic modeling that might allow us to confirm this. The last panel of Fig. 5 shows the number of Gaussians fit in each spaxel, which varies from five with increasing distance from the nucleus.

From the other maps presented in Appendix A, which are described there, we detect gas rotation in in J1010, J1100, and J1509, based on the similarity with the CO gas distribution (Ramos Almeida et al. 2022). In J1356 and J1430, the kinematics are clearly affected by the mergers, which increase the turbulence. J1100 shows disturbed rotation, which might be related to the presence of the stellar bar (Fischer et al. 2018; Ramos Almeida et al. 2022). J1430 and J1356 are the only objects in the sample in which we observe the approaching and receding sides of the outflows, with both QSO2s having similar outflow velocities of ±800 km s⁻¹. However, in J1356, the distribution of the outflowing gas is irregular, showing the highest blueshifted and redshifted velocities and W₈₀ ≈ 1000 km s⁻¹ in the southwest region (PA ∼ 200°). For the rest of the targets, we only detect negative outflow velocities (vel₀₅) coinciding with high velocity dispersion (W₈₀) regions. In J1100 the most turbulent gas is observed in the northeast (PA ≈ 60°), where W₈₀ ∼ 1600–1800 km s⁻¹ and vel₀₅ ≥ − 1500 km s⁻¹. In J1509, the highest values of W₈₀, of 1600–1700 km s⁻¹ and vel₀₅, of up to −1500 km s⁻¹, are measured in the northwest. Finally, for J1010 we detect high values of vel₀₅ and vel₉₅, but similar values of W₈₀ are measured across almost the entire extent of the [O III] emission-line region, consistent with this being unresolved (see Sect. 4.2). Despite this, we detect a rotation pattern in the vel_p map of J1010 because approximately half of the instrument fibers are probing redshifted velocities and the other half blueshifted velocities. In cases of severe beam smearing, the amplitude of rotation is usually smaller than the real one.

4.4. Outflow properties

From the analysis of the [O III] emission and kinematics presented in Sects. 4.2 and 4.3, we concluded that they are spatially resolved in J1100, J1356, J1430, and J1509. Therefore, in the following we describe and measure the outflow properties of these QSO2s, because our data do not allow us to constrain the outflow properties of J1010.

The velocity maps shown in Figs. 5, A.2, A.3, and A.4 reveal different patterns: J1430 shows an expanding shell/bubble with negative and positive velocities, J1100 and J1509 present high blueshifted velocities coincident with the highest W80 values, and J1356 shows an irregular and complex gas morphology having positive and negative outflow velocities. Therefore, we need to analyze each target separately to estimate the outflow properties and their impact on the surrounding environment, as well as the mechanism(s) that might be driving them.

We first need to measure the electron densities (n_e) and outflow fluxes to estimate the outflow masses (M_out), and the outflow radii (R_out) and velocities (v_out) to quantify the mass outflow rates (Ṁ_out) and kinetic powers (Ė_kin) of the winds. The largest uncertainty is given by n_e, which can change the results by several orders of magnitude (e.g., Rose et al. 2018; Davies et al. 2020; Holden & Tadhunter 2023). To evaluate how Ṁ_out varies with n_e, we measure the outflow density using two different methods involving different emission lines detected in the SDSS spectra available for the QSO2s. To measure the line fluxes we fit a first-order polynomial to the continuum and a single Gaussian profile to each line. A multiple Gaussian fit including broad and/or intermediate components was not done because of the low signal-to-noise of some of the emission lines involved. Therefore, the two sets of densities correspond to total densities, and they are most likely lower limits to the outflow density (e.g., Mingozzi et al. 2019; Fluetsch et al. 2021).

The first method is based on the density-sensitive ratio of [S II]λλ6716,6731 Å (Osterbrock & Ferland 2006) and the temperature-sensitive ratio of [O III]λλ4363,5007 Å. Using these lines and the Pyneb tool (Luridiana et al. 2015) we obtain densities between 2.5 ≲ log n_e≲ 3 cm⁻³ (see Table 4). However, this classical technique suffers for saturation at high densities (log n_e ≥ 3.5 cm⁻³; Rose et al. 2018) and it could lead to an overestimation of the outflow masses and, consequently, of Ṁ_out. To overcome this problem the trans-auroral lines can be used instead (Holt et al. 2011; Rose et al. 2018; Ramos Almeida et al. 2019; Holden & Tadhunter 2023; Holden et al. 2023). This second method makes use of the [S II]λλ6716,6731 Å and [O II]λλ3726,3729 Å doublets as well as of the trans-auroral [O II]λλ7319,7331 Å and [S II]λλ4068,4076 Å doublets. The electron densities are derived by comparing the [O II] and [S II] ratios TR([O II]) = F(3726+3729)/F(7319+7331) and TR([S II]) = F(4068+4076)/F(6717+6731) with a grid of photoionization models computed with Cloudy (Ferland et al. 2013), using a spectral index of α = 1.5 and a ionization parameter of log U=−2.3. The resulting densities obtained from the trans-auroral lines are 3.2 ≲ log n_e≲ 4 cm⁻³. We list the individual values obtained with both methods in Table 4, and for further details on the estimations of the electron densities we refer the reader to Sect. 3.2 of Speranza et al. (2022), where the two methods are deeply explained.

To measure the outflow flux we integrate the [O III] outflow emission at ≥3σ (black contours in Fig. 6), as in Speranza et al. (2022). To do so, we selected the high-velocity [O III] gas, which varies from case to case: −200 and −2000 km s⁻¹ for J1100 and J1509, and −400 and −2000 km s⁻¹ and 400 and 1800 km s⁻¹ for J1356 and J1430, for which the redshifted counterpart of the outflow is also detected. The interval of velocities are relative to vel_p of each spaxel. Since here we are using a nonparametric analysis, the flux is computed using the modeled [O III]λ5007 Å profile (see Fig. 4 for reference). This has the advantage that we do not have to worry about other emission lines close to [O III]λ5007 Å or noisy continua, making it possible to integrate the flux at large velocities. Once we have the outflow flux maps, we fit an ellipse to the 3σ contour (red dashed lines in Fig. 6) and we integrate the flux within the ellipse. The ellipse is fit using a Python script presented by Hill (2016) in the online supplementary material. We calculate the errors on the flux using ellipses fit to the 4σ and 2σ contours, shown as orange lines in Fig. 6. The [O III] fluxes are then corrected from extinction using the A_v values reported by Kong & Ho (2018) for our QSO2s and the extinction law of Cardelli et al. (1989). Using the previously described outflow densities and fluxes, we can measure the corresponding outflow masses using Eq. (5) in Carniani et al. (2015) (see also Osterbrock & Ferland 2006):

$$\begin{array}{c}\hfill M_{\mathrm{out}}=4\times {10}^7{M}_{\odot }\left(\frac{\text{C}}{{10}^{\mathrm{O}/\mathrm{H}}}\right)\left(\frac{L_{\mathrm{out}}}{{10}^{44}}\right){\left(\frac{⟨n_{\mathrm{e}}⟩}{{10}^3}\right)}^{-1}\end{array}$$(1)

Fig. 6. [O III] flux maps corresponding to high-velocity gas only, with velocities between −2000 < vout <  − 200 km s−1, −2000 < vout < − 400 km s−1, and 400 < vout < 1800 km s−1 depending on the source. The black solid line corresponds to the flux contour at 3σ and the red dashed line to the ellipse fit to it. The white cross indicates the [O III] emission peak, and the red cross corresponds to the most distant point from the white cross within the fit ellipse. The distance between the two crosses is the maximum outflow radius (Rout) reported in Table 4. The two orange lines are the flux contours at the 2σ and 4σ levels used to estimate the errors of the outflow radius and flux.

where L_out = 3 × L_{[O III]} to consider the total mass of ionized outflow as in Fiore et al. (2017), n_e is the electron density, $C=\langle n_{\text{e}}^2\rangle -{\langle n_{\text{e}}\rangle }^2=1$ is the condensation factor, and O/H = [O/H]–[O/H]_⊙ = 0 (assuming that the gas clouds have the same electron density and solar metallicity as in Bischetti et al. 2017).

To measure the outflow radius (R_out) we use the same outflow flux maps described above. We obtain R_out as the maximum distance between the [O III] emission peak and the ellipse that we fit to the 3σ contours (see Fig. 6). By doing this we simultaneously obtain the outflow radius and PA, measured from the north. This is the same method proposed in Speranza et al. (2022), which allows us to include spaxels with low surface brightness but high values of velocity and velocity dispersion that would be missed by measuring the FWHM of the high-velocity flux maps (as in Ramos Almeida et al. 2017). Here, we also estimate the errors by measuring the corresponding extensions at 4σ and 2σ (see Fig. 6). Following Kang & Woo (2018), we then subtract in quadrature the half width at half maximum (HWHM) measured from the standard stars observed after each target, which is half of the seeing FWHM measured for each QSO2 (shown in Table 2), from the outflow radius. With this method, we measure a range of R_out values for the QSO2s, from 2″ (∼3.1 kpc) to 5.7″ (∼12.6 kpc). For comparison, Ramos Almeida et al. (2017, 2019) reported radial sizes of 1.1 ± 0.1 kpc and 1.3 ± 0.2 kpc for the ionized outflows detected in J1430 and J1509, respectively, whereas here we measure $R_{\mathrm{out}}=3.7{\pm }_{0.2}^{0.4}$ kpc for the blueshifted side of the outflow in J1430 and $R_{\mathrm{out}}=3.8{\pm }_{0.2}^{0.3}$ for J1509.

Finally, to measure the outflow velocity within the outflow radius, we averaged the values of vel₀₅ and vel₉₅ (if the redshifted outflow counterpart is detected) from the outflow velocity maps (see middle panels of Fig. 5). We did this by selecting only those spaxels within the 3σ ellipse that we use to measure the corresponding R_out. The errors of vel₀₅ were obtained by averaging the vel₀₂ and vel₁₀ maps and calculating the difference with vel₀₅, and for vel₉₅, we did the same but using vel₉₈ and vel₉₀. These percentiles are representative of the least and most conservative outflow velocities, respectively. All the resulting velocities are listed in Table 4.

Then, using the outflow mass, radius, and velocity, and assuming a constant average volume density, as in Fiore et al. (2017), the outflow mass rate can be calculated as:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{out}}=3\times v_{\mathrm{out}}\times \frac{M_{\mathrm{out}}}{R_{\mathrm{out}}}\end{array}$$(2)

and the kinetic power as:

$$\begin{array}{c}\hfill {\dot{E}}_{\mathrm{kin}}=\frac{{\dot{M}}_{\mathrm{out}}}{2}(v_{\mathrm{out}}^2+3\sigma )\end{array}$$(3)

where σ∼ FWHM/2.355, with FWHM = W80/1.09, and W₈₀ is averaged within the red ellipses shown in Fig. 6. In the cases of J1356 and J1430, for which we detect both sides of the outflow, Ṁ_out and Ė_kin correspond to the sum of the individual values measured for the blueshifted and redshifted counterparts.

In Table 4 we show all the outflow properties described in this Section (i.e., R_out, PA, v_out, M_out, Ṁ_out, and Ė_kin), based on the two electron density measurements presented above. Using the [S II]-based n_e we obtain Ṁ_out = 3.3 − 6.5 M_⊙ yr⁻¹ and Ė_kin = 10^41.1−42.0 erg s⁻¹, whereas using the transauroral-based densities, they are a factor of 3–4 lower: Ṁ_out = 0.7–1.6 M_⊙ yr⁻¹ and Ė_kin = 10^40.6−41.0 erg s⁻¹. We also list the corresponding coupling efficiencies in Table 4, ξ = Ė_kin/L_bol, which are 0.001–0.01% if we use the [S II]-based densities, and 0.0003–0.001% if we use the trans-auroral-based densities.

5. Discussion

In Sect. 4 we presented the results from our analysis of the kinematics of the nuclear (i.e., the central ∼2.2 kpc of the QSO2s) and extended ionized gas emission, traced by [O III], of the five QSO2s observed with GTC/MEGARA. From this analysis we concluded that four of the five QSO2s show spatially resolved [O III] kinematics and we characterized their outflow properties. Hereafter we discuss the mechanisms that are either responsible or contribute to driving them, their impact on the host galaxies, and the interplay between the ionized and cold molecular gas phases.

5.1. The role of radio jets in driving the outflows

As described in Sect. 2, the five QSO2s are radio-quiet and their AGN luminosities are relatively high (10^{45.5 − 46} erg s⁻¹). Considering this, the ionized gas outflows that we studied here should in principle be driven by radiation pressure, generally referred to as “quasar-mode” feedback (Fabian 2012). However, it has been extensively claimed in the literature that low-power compact jets in radio-quiet AGN can compress and shock the surrounding gas as they expand, contributing to drive both ionized and molecular outflows (see e.g., Venturi et al. 2018; Jarvis et al. 2019; Speranza et al. 2022; Audibert et al. 2023). Alternatively, it has been suggested that the extended radio structures often seen in radio-quiet AGN might be a consequence of AGN winds impacting a higher-density ISM and producing shocks (Fischer et al. 2019, 2023).

In order to explore possible jet-gas interactions in our sample, we compared our outflow flux maps with VLA data at 6 GHz of similar angular resolution. All the targets, except J1509, have observations available at high (∼0.25″ beam) and low-angular resolution (∼1″ beam) from Jarvis et al. (2019). These authors reported that in the cases of J1010, J1100, J1356, and J1430, only between 3% and 6% of their radio emission can be accounted for by star formation. The VLA data of J1010 and J1100, at both high and low angular resolutions, do not reveal any extended component. For J1010, Jarvis et al. (2019) reported a deconvolved size of 0.11″ (∼200 pc), but for J1100 they did not do the same because they relied on the detection of an extended structure only seen in e-MERLIN data. In Figs. A.1 and A.2 we show the low-angular resolution VLA contours of J1010 and J1100. On the other hand, J1356 and J1430 show extended radio emission in the VLA data that Jarvis et al. (2019) identified with a jet or lobe in the case of J1430, and inconclusive in the case of J1356. For these two galaxies, they reported that the extended radio structures were associated with morphologically and kinematically distinct features detected in [O III], indicative of interaction between them. Here we take advantage of the larger FOV covered by GTC/MEGARA to further investigate this.

J1356 shows irregular and messy gas kinematics, as it can be seen from Fig. A.3. This is because this QSO2 is hosted in an ongoing merger between two galaxies, the north and south components (see Sect. 5.1.3 in Ramos Almeida et al. 2022). We detect a blueshifted component of the outflow expanding toward the south up to ∼12.6 kpc from the nucleus. Indeed, this is the maximum radius that we can measure from our data, as we are limited by the FOV (see Fig. 6). This large-scale blueshifted high-velocity [O III] gas was first reported by Greene et al. (2012) and referred to as “the bubble”. These authors reported a total size of ∼20 kpc for this outflowing gas. Here we measure a redshifted counterpart of the outflow expanding towards the southwest (PA ∼ 200°) with a maximum extension of ∼6.8 kpc (see Fig. 6 and Appendix A.3 for more details). This redshifted emission, with velocities between 400 and 1800 km s⁻¹ relative to the corresponding velocity peak, is perfectly aligned with the 6 GHz radio contours at 1″ resolution, as can be seen from the left panel of Fig. 7. The radio contours, which have the same angular resolution as our GTC/MEGARA data (∼1″) extend to the southwest to then bend towards the south, up to −2.5″ (∼5.5 kpc) from the [O III] emission peak (the black cross in Fig. 7). At the same distance from the center we observe a sharp decrease of the outflowing flux, of one order of magnitude (from 10⁻¹⁵ to 10⁻¹⁶ erg s⁻¹ cm⁻² Å ⁻¹). The same contours are superimposed to the vel₀₅, vel₉₅, and W₈₀ maps in Fig. A.3. The highest outflow velocities and W₈₀ are measured within the radio contours and/or immediately next to them, particularly in the case of the redshifted counterpart of the outflow. We argue that this might be observational evidence for a jet-like structure accelerating the ionized gas as it progresses through the galaxy’s ISM (see Audibert et al. 2023 and references therein). However, it is worth noting that such extended radio structure lacks a high resolution counterpart (i.e., in the ∼0.2″ VLA data). This is the reason behind the inconclusive classification of this feature reported by Jarvis et al. (2019), although they favored the jet or lobe scenario. Another possibility would be a shocked-wind origin for this extended radio structure. Radiatively driven winds could be generating synchrotron emitting shocks in the ISM (e.g., Zakamska et al. 2016; Hwang et al. 2018; Fischer et al. 2023). This is plausible for radio-quiet AGN with bolometric luminosities comparable to J1356 (i.e., L_bol ≈ 10⁴⁵ erg s⁻¹; Faucher-Giguère & Quataert 2012).

Fig. 7. Continuum subtracted [O III] flux maps of J1356 and J1430 extracted from the cubes in different velocity intervals, indicated in the top right corner of each panel. The black cross corresponds to the [O III] peak, the black contours to the 6 GHz VLA contours at ∼1″ resolution (at 3, 10, 15, 30, and 60σ for J1356 and at 3, 5, 10, 20, and 60σ for J1430), and the white contours to the 6 GHz VLA contours at ∼0.25″ resolution (at 3, 5, 15, 30, and 60σ).

Focusing on J1430, our kinematic maps reveal an extended [O III] outflow, with blueshifted velocities increasing to the northeast (R_out = 3.7 kpc and PA = 65°) and redshifted velocities towards the southwest (R_out = 3.1 kpc and PA = 198°; see middle panels of Fig. 5). The W₈₀ map shown in the bottom left panel of the same figure shows maximum values of ∼900 km s⁻¹ in the same region where the maximum blueshifted velocities are detected. The high-angular resolution VLA observations at 6 GHz reveal a compact jet-like structure with a total extent of ∼0.8 kpc (PA ∼ 60°; Harrison et al. 2015; Jarvis et al. 2019; see right panel of Fig. 7). Here we find that the PA of the blueshifted outflow emission coincides with the jet direction (see white contours in the middle left panel of Fig. 5), in agreement with the results from VLT/MUSE data presented in Venturi et al. (2023). In a recent work, Audibert et al. (2023) showed that this compact radio jet, which subtends a small angle with the disk of cold molecular gas, is compressing and accelerating the molecular gas, and driving a lateral outflow that shows enhanced velocity dispersion and higher gas excitation. Our data also shows enhanced velocity dispersion (W₈₀) in the direction perpendicular to the jet, reported by Venturi et al. (2023), with maximum values at the jet’s forefront (see the bottom left panel of Fig. 5), in agreement with the results obtained for the molecular gas in this object (Audibert et al. 2023), and those found for the warm ionized outflows in some nearby Seyfert galaxies (Venturi et al. 2021). This type of jet-gas interaction is typical of the so-called frustrated-jets, that is, low-power jets that make slow progress through the ISM, disrupting and causing greater damage than powerful collimated jets (Wagner & Bicknell 2011; Nyland et al. 2018). Hydrodynamical simulations show that clouds of material with densities of n_e = 300 cm⁻³ can efficiently confine jets with P_jet ≤ 10⁴⁴ erg s⁻¹ (Mukherjee et al. 2016; Bicknell et al. 2018; Talbot et al. 2022). In the case of J1430, we measure outflow densities > 1700 cm⁻³ using the trans-auroral lines, and Audibert et al. (2023) reported a jet power of P_jet ∼ 10⁴³ erg s⁻¹. For J1356 we measure outflow densities > 1600 cm⁻³ using the same method, and if we consider that the extended radio structure detected in the low-resolution VLA data is a jet, from these data we measure P_jet∼10⁴⁴ erg s⁻¹, following the same prescription as in Audibert et al. (2023). Therefore, J1430, and possibly J1356 as well if the extended VLA feature is confirmed as a jet, could be examples of radio-quiet quasars with a low-power jet slowly advancing through the galaxy ISM, shocking and accelerating the surrounding gas, hence contributing to drive the ionized and molecular outflows. A jet-driven outflow is consistent with the kinetic powers that we measure for both QSO2, of Ė_kin ∼ 10⁴¹x erg s⁻¹, two-three orders of magnitude lower than the jet power.

In the middle and bottom left panels of Fig. 5 we show the contours of the brightness temperature ratio (T₃₂/T₂₁) derived from the CO(3–2) and CO(2–1) lines by Audibert et al. (2023) for J1430. The highest values are distributed in the perpendicular direction to the radio jet, slightly offset from the maximum velocities of the CO velocity dispersion (see Fig. A.1 in Audibert et al. 2023). We note that the T₃₂/T₂₁ contours extend up to where the ionized gas shows the highest outflow velocities and W₈₀ values, that is, at ∼ − 0.5″ in the x-axis. This might suggest that the ionized outflow is being accelerated by the combined action of the compact radio jet and the molecular gas that it is entraining.

On larger scales, we report an enhancement of vel₉₅ and W₈₀ that overlaps with the extended radio emission traced by the low-resolution 6 GHz data to the southwest (see middle right and bottom left panels of Fig. 5, where the 1″ resolution radio contours are shown in cyan). In this region the [O III] emission reaches the highest positive outflow velocities, between 600 and 800 km s⁻¹, with corresponding W₈₀ values of 600–800 km s⁻¹.

The results that we report here for J1430, and tentatively for J1356 as well, contribute to strengthen the idea that, even in radio-quiet sources, low-power jets can successfully transfer mechanical energy to the ISM, contributing to drive and/or accelerate ionized gas outflows (e.g., Cresci et al. 2023). Therefore, radiation-driven winds might not be the only channel contributing to launch kpc-scale outflows in luminous AGN such as our QSO2s, but also low-power compact jets.

5.2. The impact of electron density on scaling relations

Harrison et al. (2014) performed a nonparametric analysis of the [O III] kinematics for a sample of 16 radio-quiet QSO2s that included J1010, J1100, J1356, and J1430. They used GMOS IFU data, which has a 5″ × 3.5″ FOV, smaller than MEGARA’s (12.5″ × 11.3″), and assumed the same spherical or multicone outflow geometry that we are using here to calculate the mass outflow rates (i.e., Eq. (2)). They measured the electron densities from the [S II] doublet detected in the SDSS spectra of the QSO2s, as we also do here, and they found values between 200 and 1000 cm⁻³, fully consistent with ours. However, Harrison et al. (2014) assumed a value of n_e = 100 cm⁻³ for all the targets, and the total, extinction-uncorrected Hβ luminosities to work out the outflow masses. For this reason they reported values of M_out = (2 − 40)×10⁷ M_⊙, which are one order of magnitude larger than ours, of (4–35)×10⁶ M_⊙, when we use the [S II]-based n_e. Assuming σ = W₈₀/2.355 and R_out=R_[OIII], they reported outflow mass rates between 3 and 70 M_⊙ yr⁻¹ for their sample. For our QSO2s, we measured outflow mass rates between 3.3 and 6.5 M_⊙ yr⁻¹ using the [S II]-based n_e (see Table 4). If we use the trans-auroral-based n_e values (1600–9800 cm⁻³), reported in the same table, we measure outflow mass rates between 0.7 and 1.6 M_⊙ yr⁻¹.

A mass outflow rate of ≤0.46 M_⊙ yr⁻¹ was reported by Ramos Almeida et al. (2019) for J1509 based on near-infrared data from GTC/EMIR. In fact, the outflow rate is an upper limit because the electron density was derived from the total fluxes of the trans-auroral lines, as we are doing here. For J1509 here we measure $1.1{\pm }_{0.4}^{0.3}{M}_{\odot }{\mathrm{yr}}^{-1}$. For J1430, Venturi et al. (2023) reported an ionized mass outflow rate of 40–130 M_⊙ yr⁻¹, measured from VLT/MUSE data on the second radial bin they considered (i.e., between 1 and 2 kpc). The discrepancy between our ionized mass outflow rate and that reported by Venturi et al. (2023) mainly comes from 1) the lower values of the electron density employed in the latter work (average value of ∼500 cm⁻³ in the second radial bin), which were used to calculate the outflow mass on a spaxel-by-spaxel basis; and specially 2) the different methods adopted to estimate the outflow rate. In Venturi et al. (2023), the authors integrated the flux of the broad components fit to the line profiles in each spaxel, which is significantly higher than the outflow flux estimated from the nonparametric method we used here (see Hervella Seoane et al. 2023), and they computed the outflow rate on each spaxel by using ΔR_out equal to the size of each spaxel (0.2″). Then, they produced radial profiles of Ṁ_out by summing, in each radial bin, the values measured on each spaxel and then rescaling using ΔR_out/ΔR_bin, with ΔR_bin = 1 kpc. In addition, we note that Venturi et al. (2023) reported a decrease in Ṁ_out by a factor of 3–20 if they use v_out = v_s instead of v_out = v_s + FWHM/2, where v_s is the velocity shift of the broad components from the systemic velocity.

To put our results in a broader context, we compare them with those reported by Fiore et al. (2017), who compiled and put together several ionized outflow measurements from the literature, establishing a scaling relation between Ṁ_out and the AGN luminosity. This comparison is shown in Fig. 8, where we also represent how our results vary by using different values of n_e. When we adopt n_e = 200 cm⁻³, as done by Fiore et al. (2017), we obtain the highest mass outflow rates, of 10–32 M_⊙ yr⁻¹ (green circles in Fig. 8). These values are consistent with the observational scaling relation, defined by the fit of AGN at z < 0.5 (i.e., the purple line in the same figure). However, when the electron densities are calculated individually for each source, the mass outflow rates decrease significantly. Using the n_e measured from the [S II] doublet, our QSO2s (the cyan squares in Fig. 8) fall well below the scaling relation. This offset becomes of one order of magnitude when we use the densities derived from the trans-auroral method (see green triangles in Fig. 8), a trend also shown in the study conducted by Santoro et al. (2020). These results highlight the importance of having reliable estimates of the outflow densities to derive mass outflow rates (Harrison et al. 2018; Rose et al. 2018; Davies et al. 2020; Revalski et al. 2022; Holden & Tadhunter 2023), and also, of revisiting the widely used empirical relations. These are likely an upper limit to the real, and much more scattered relation between the mass outflow rate and AGN luminosity (Ramos Almeida et al. 2022; Lamperti et al. 2022).

Fig. 8. Mass outflow rate versus AGN luminosity. The purple line corresponds to the linear fit of the ionized outflow mass rates calculated by Fiore et al. (2017) for a sample of AGN at z < 0.5 (purple squares) assuming a value of the density of 200 cm−3. The black dot-dashed line corresponds to the linear fit from Davies et al. (2020), that includes the mass rates of Fiore et al. (2017), but corrected by using individual electron densities. The other symbols correspond to the four QSO2s with spatially resolved [O III] outflows, with different shapes and colors indicating the values of ne adopted to perform the outflow mass rate calculations: green circles correspond to densities of 200 cm−3, and cyan squares and light-green triangles to those measured from the [S II] doublet and the trans-auroral lines, respectively. The outflow mass rates of the QSO2s vary up to 2 orders of magnitudes depending on how the ne is measured.

In line with the latter, Baron & Netzer (2019) did not find a linear correlation between the outflow mass rate and AGN luminosity using a sample of 234 local type-2 AGN with L_bol = 10^{43.3 − 45.8} erg s⁻¹. They estimated individual electron densities using optical line ratios at the location of the wind (see also Revalski et al. 2018a,b), and they reported Ṁ_out of 1–2 orders of magnitude lower than previous estimates. Davies et al. (2020) also reported lower mass outflow rates than those predicted from the Fiore et al. (2017) relation for a sample of 11 Seyfert 2 galaxies using individual n_e values determined as in Baron & Netzer (2019). By putting together their outflow mass rates, those from Baron & Netzer (2019), and those from Fiore et al. (2017) after calculating individual n_e values for the latter, they proposed a modified version of the Fiore et al. (2017) scaling relation, shown as a dot-dashed line in Fig. 8. They concluded that, even with a larger scatter, a linear correlation is still observed across 5 orders of magnitude in L_bol, but with outflow mass rates that are lower by about a factor 3, due to the higher n_e values. Our mass rates, both computed using [S II] and trans-auroral densities, are in agreement with this corrected scaling relation from Davies et al. (2020), which is also shown in Fig. 8.

As mentioned in Sect. 4.4, here we used two different methods to estimate the outflow densities, but they correspond to total values, as one single Gaussian component was fit to the corresponding emission lines measured within the SDSS fiber (1.5″ radius). In Speranza et al. (2022) we were able to fit the [S II] and the trans-auroral lines with a narrow and a broad component for the QSO2 J0945, and we measured a significantly larger n_e value for the latter, being the total value more similar to that of the narrow component (see Table 2 and Fig. 5 in Speranza et al. 2022). Therefore, it is likely that our QSO2s outflow densities will be even higher than the total n_e values employed here, which implies that the outflow masses and mass rates will be even lower.

Finally, it is worth mentioning that here we are focusing on the influence of the outflow density on the mass outflow rates, but there are other parameters, such as the outflow radius, velocity, and assumed geometry, that strongly impact the outflow rates. Furthermore, the bolometric luminosities also have large uncertainties, of up to ∼1 dex, because of the systematic uncertainty on the bolometric correction (Jarvis et al. 2019). Nevertheless, Ramos Almeida et al. (2023) reported good agreement between the AGN luminosities derived from extinction-corrected [O III] luminosities and hard X-rays for the QSO2 Mrk 477 (∼0.2 dex).

5.3. A multiphase view of outflows

The ionized gas emission corresponds to the warm phase of gas, having typical temperatures of T ∼ 10⁴ K. As we have shown in Sect. 5.2, QSO2s typically have total ionized gas masses of (2–40)×10⁷ M_⊙ (Harrison et al. 2014). These gas masses are much lower than the molecular gas masses measured from CO lines detected in the millimiter range, which range between ∼10⁸ and 10¹⁰ M_⊙ for QSO2s at z < 0.2 (Jarvis et al. 2020; Ramos Almeida et al. 2022; Molyneux et al. 2024). It has been observationally demonstrated that molecular outflows carry more mass than their ionized counterparts, at least in the local universe, and at AGN luminosities below 10⁴⁷ erg s⁻¹ (Fiore et al. 2017; Fluetsch et al. 2021). However, as it happens with the ionized outflow mass rates, the values used to build the molecular gas counterpart of the scaling relations (Cicone et al. 2014; Fiore et al. 2017) are higher than those found for the general population of quasars (Ramos Almeida et al. 2022), ULIRGs (Lamperti et al. 2022), and Seyfert galaxies (Alonso-Herrero et al. 2019; Domínguez-Fernández et al. 2020; Zanchettin et al. 2021).

In this work, we aim to provide a more comprehensive view of quasar-driven outflows by putting together the cold molecular and ionized outflow measurements available for our QSO2s. Using the CO(2–1) transition observed with ALMA at 0.2″ resolution, Ramos Almeida et al. (2022) detected spatially resolved cold molecular outflows for the same four QSO2s where we detect spatially resolved ionized outflows (i.e., J1100, J1356, J1430, and J1509). They report higher mass outflow rates, of ${\dot{\mathrm{M}}}_{\mathrm{out}}^{{\mathrm{H}}_2}=8-16M_{\odot }{\mathrm{yr}}^{-1}$, than our ionized outflow mass rates, which are Ṁ_out = 3.3 − 6.5 M_⊙ yr⁻¹ if we use the [S II]-based n_e, and 0.7–1.6 M_⊙ yr⁻¹ if use the higher densities estimated from the trans-auroral lines. It is worth mentioning that this comparison is not straightforward because of (1) the different angular resolution of the observations used (0.2″ in the case of the ALMA observations and 1.2″ for the MEGARA observations) and (2) the different methodologies employed to derive the outflow properties. In the case of the ALMA CO(2–1) observations, Ramos Almeida et al. (2022) performed kinematic modeling of the data using ^3DBAROLO (Di Teodoro & Fraternali 2015) to disentangle rotation from noncircular motions, and integrated the flux of the high-velocity gas detected along the kinematic minor axis, using an slit-simulated aperture of 0.2″ width (see Ramos Almeida et al. 2022 for further details). Instead, here we integrated all the flux of the high-velocity ionized gas within ellipses defined by the 3σ contours. However, despite these differences in angular resolution and methodology, we note that the molecular outflow measurements are rather conservative, and different methodologies will most likely result in even higher molecular outflow masses and/or outflow rates. As an example, Audibert et al. (2023) considered four different scenarios, from most to least conservative, to calculate the outflow mass rates of J1430 using the same CO(2–1) observations as Ramos Almeida et al. (2022), and reported values from 6.7 to 44 M_⊙ yr⁻¹.

Higher molecular outflow rates are found because, even though the molecular outflow velocities are lower (≤300 km s⁻¹) than the ionized outflow velocities that we measure here (500 ≤ |v_out| ≤ 1300 km s⁻¹), the cold molecular outflows carry more mass and are more compact. For the latter Ramos Almeida et al. (2022) reported masses of (1–10) × 10⁷ M_⊙ and radii ≤3 kpc, whereas here we measure (1–35) × 10⁶ M_⊙ and radii of 3.1–12.6 kpc for the ionized outflows (see Table 4). Thus, for our sample, the molecular mass outflow rates are on average 2–3 times larger than the ionized mass outflow rates. However, from the comparison between the ionized and cold molecular scaling relations in Fiore et al. (2017), the gap between the two gas phases at the bolometric luminosities of the QSO2s is expected to be almost of two orders of magnitude, that is a factor of ∼50. This is shown in Fig. 9, where the observational scaling relations for the molecular and ionized outflows are the pink and purple lines, respectively. The molecular outflow mass rates measured by Ramos Almeida et al. (2022) are the yellow circles, and the corresponding ionized mass outflow rates, derived using the [S II]-based densities, are the cyan squares. If we consider the ionized outflow rates obtained using the trans-auroral-based densities, they are ∼12 times lower than the molecular outflow rates. Using 0.2″ resolution CO(2–1) observations of a sample of nearby ULIRGs, Lamperti et al. (2022) measured molecular outflow rates of 16–145 M_⊙ yr⁻¹ for six AGN-dominated ULIRGs with similar bolometric luminosities as our QSO2s, and upper limits of ≤3 − 32 M_⊙ yr⁻¹ for another five AGN-dominated ULIRGs (see gray circles in Fig. 9).

Fig. 9. Ionized and cold molecular gas mass outflow rates versus AGN luminosity. The purple and cyan squares are the same as in Fig. 8, and the pink circles are the cold molecular outflow measurements from Fiore et al. (2017). The corresponding linear fits are indicated with the same colors. The yellow circles are the cold molecular mass outflow rates reported by Ramos Almeida et al. (2022) for the same QSO2s analyzed in this work, and the grey circles AGN-dominated ULIRGs from Lamperti et al. (2022).

These results highlight that the Ṁ_out versus L_bol empirical relation might be just the upper envelope of a much more scattered relation between these two quantities. Based on our multiphase analysis of this small sample of QSO2s, we find that the cold molecular outflows are indeed carrying more mass and are more compact than their ionized counterparts. This introduces a significant offset in Ṁ_out between the two gas phases. A larger sample with multiphase outflow measurements in the same targets, which is the ultimate goal of the QSOFEED project (Ramos Almeida et al. 2022), is needed to derive statistically significant conclusions. If the linear trend is still observed in a larger sample, it might just need to be rescaled (Davies et al. 2020). Otherwise, as discussed in Ramos Almeida et al. (2022), outflow mass rates might not depend only on AGN luminosity, but in other factors such as the coupling between the jets and/or the winds and the ionized and molecular gas, and the amount and distribution of gas in the central regions. A combination of these two scenarios is also possible, in which case the scaling relation would persist, albeit with larger scatter.

Using the SFRs of the QSO2s, derived from their far-infrared luminosities (see Table 1), we can calculate the mass loading factors (η = Ṁ_out/SFR) of the ionized outflows, which are η = 0.1–0.2. These values cover a smaller range than those measured for the cold molecular phase (η_H2 = 0.1–1.3; Ramos Almeida et al. 2022). If we sum the loading factors of the two phases, η_tot = η + η_H2, we get η_tot = 0.3 + 1.3 = 1.6 for J1430, and η_tot = 0.2–0.5 for the other QSO2s. Values of η_tot ≤ 1 imply, in principle, that outflows are less effective removing gas than star formation, and also, that star-formation could be potentially driving them. However, a few caveats must be mentioned in this respect. First, we argue that the outflows, both molecular and ionized, are AGN-driven. The high outflow velocities that we measure cannot be accounted for by star formation (Speranza et al. 2022). Second, the SFR measurements are calculated from the total far-infrared luminosity, and therefore are galaxy-wide SFRs. The outflows, on the other hand, are more compact, specially the molecular outflows (r_out = 0.5–3 kpc; Ramos Almeida et al. 2022). Indeed, despite the relatively low η_H2 values of the QSO2s, Ramos Almeida et al. (2022) and Audibert et al. (2023) reported that these outflows, in some cases in combination with compact jets, are modifying the molecular gas distribution in the central kiloparsec of the galaxies. High angular resolution IFU observations are needed to compare the spatial distribution of recent star formation and AGN outflows to really evaluate the impact of the latter.

Another caveat is the different timescale of the star formation probed by the far-infrared luminosity and the outflows. Considering that the QSO2 outflows have dynamical timescales of ∼1–10 Myr (Bessiere & Ramos Almeida 2022; Ramos Almeida et al. 2022), their properties should be compared with recent star formation. This can be probed using resolved stellar population analyses (Bessiere & Ramos Almeida 2022) or the polycyclic aromatic hydrocarbon (PAH) features detected in the mid-infrared, which are good probes of young and massive stars, but whose properties vary with distance from the AGN (García-Bernete et al. 2022; Ramos Almeida et al. 2023). These types of studies can reveal whether recent star formation is prevented or triggered by the outflows, or if both can simultaneously happen (Cresci et al. 2015; Carniani et al. 2016; Bessiere & Ramos Almeida 2022). Indeed, for the five QSO2s studied here, we will study the spatial distribution of the PAH features using recently awarded Cycle 2 James Webb Space Telescope (JWST) observations (PI: Ramos Almeida, proposal 3655), and we will compare it with the spatially resolved outflows presented here. From our measurements of the outflow radii (R_out = 3.3–12.7 kpc) and velocities (v_out = 500–1300 km s⁻¹) we obtain dynamical timescales for the ionized outflows of t_dyn = R_out/v_out = 3–20 Myr, consistent with those of the molecular outflows and the young stellar populations that we will probe with JWST.

6. Summary and conclusions

We have investigated the [O III]λ5007 Å emission of five nearby QSO2s (z ∼ 0.1) from the QSOFEED sample using seeing-limited, high spectral resolution observations obtained with GTC/MEGARA. The IFU capabilities of this instrument allowed us to explore the nuclear emission from the central 1.2″ (∼2.2 kpc) of the galaxies, and the ionized extended emission across a FOV of 12.5″ × 11.3″ (∼23 × 21 kpc²). We characterize the [O III] kinematics of the QSO2s and, thanks to a recent high angular resolution study of their CO(2–1) emission (Ramos Almeida et al. 2022), we provide a comprehensive multiwavelength view of their gas content and kinematics. The main conclusions from our work are listed as it follows.

• The nuclear spectra of the five QSO2s present signatures of outflowing gas. Broad components with 1300 ≤ FWHM ≤ 2200 km s⁻¹ are detected, all with blueshifted velocities ranging from −40 to −430 km s⁻¹. In all the QSO2s but in J1509, we also fit intermediate components to reproduce the line profiles, with 340  ≤  FWHM  ≤  900 km s⁻¹, some of them blueshifted (with velocities of up to −365 km s⁻¹) and others redshifted (up to 280 km s⁻¹) with respect to the systemic velocity.

• We detect kinematically disturbed gas in the five QSO2s (high velocities and turbulence) that we attribute to outflowing gas. These outflows are spatially resolved in all the QSO2s except in J1010. For the four QSO2s in which the outflows are resolved, we measure radii between 3.1 and 12.6 kpc. In J1356 and J1430 we detect both sides of the outflow, approaching and receding, and in J1100 and J1509, only the approaching side.

• The receding side of the outflow in J1356, and the approaching side of the one in J1430 are well aligned with the extended radio structures detected in VLA observations. This co-spatiality might be indicative of compact, low-power jets driving or contributing to accelerate the ionized outflows in these radio-quiet QSO2s. Alternatively, radiativelly driven winds could be inducing shocks in the ISM and producing the extended structures, specially in the case of J1356, where the nature of the extended radio structure is unclear.

• The electron densities measured from the [S II]λλ6716, 6713 doublet range between 300 and 1000 cm⁻³, and the corresponding mass outflow rates are 3.3–6.5 M_⊙ yr⁻¹. If instead we use the higher electron densities measured from the trans-auroral lines (1600–9800 cm⁻³), the mass outflow rates decrease by a factor of 3–4 (0.7–1.6 M_⊙ yr⁻¹).

• The ionized and cold molecular mass outflow rates measured for the four QSO2s with spatially resolved outflows are well below the empirical scaling relation from Fiore et al. (2017). The cold molecular outflows carry more mass and are more compact than their ionized counterparts, but the expected gap between the two phases at the bolometric luminosities of the QSO2s decreases from a factor of ∼50 to 2–12 depending of the method used to calculate the electron densities. We argue that empirical scaling relations might need to be rescaled by using precise multiphase outflow measurements of the same targets.

• Using our ionized and cold molecular outflow mass rates and the total SFRs measured from the far-infrared luminosities we find a total mass loading factor ${\eta }_{\mathrm{tot}}=({\dot{M}}_{\mathrm{out}}+{\dot{M}}_{\mathrm{out}}^{{\mathrm{H}}_2})/$SFR = 1.6 in J1430 (the Teacup), and η_tot = 0.2 − 0.5 in the other three QSO2s. Despite these low values, we argue that the outflows are AGN-driven, based on different considerations (e.g., the high outflow velocities).

We do not find a significant impact of quasar-driven outflows on the global star formation rates when considering the energy budget of the molecular and ionized outflows together. However, considering the dynamical timescales of the outflows, of 3–20 Myr in the case of the ionized gas and 1–10 Myr for the molecular gas, spatially resolved measurements of the recent star formation in these targets are needed in order to evaluate this impact fairly. Our forthcoming JWST observations of the five QSO2s studied here will permit us to do so and thus advance in our understanding of AGN feedback in the local universe.

Acknowledgments

Based on observations made with the GTC, installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, in the island of La Palma. The authors are extremely grateful to the GTC staff for their constant and enthusiastic support. This work has been supported by the EU H2020-MSCA-ITN-2019 Project 860744 “BiD4BESt: Big Data applications for black hole Evolution STudies”. CRA and GS acknowledge the project “Feeding and feedback in active galaxies”, with reference PID2019-106027GB-C42, funded by MICINN-AEI/10.13039/501100011033. CRA also acknowledges the project “Quantifying the impact of quasar feedback on galaxy evolution”, with reference EUR2020-112266, funded by MICINN-AEI/10.13039/501100011033 and the European Union NextGenerationEU/PRTR. OGM acknowledge support from the PAPIIT/DGAPA project IN109123 from UNAM and the “Frontera de la Ciencia” project CF-2023-G-100 from CONACyT. AL is partly supported by the PRIN MIUR 2017 prot. 20173ML3WW 002 ‘Opening the ALMA window on the cosmic evolution of gas, stars, and massive black holes’. CRA thanks the Kavli Institute for Cosmology of the University of Cambridge for their hospitality while working on this paper, and the IAC Severo Ochoa Program for the corresponding financial support. GS and CRA thank Begoña García Lorenzo, Giacomo Venturi, and Chris Harrison for helpful discussions and Patricia Bessiere for sharing the SDSS outflow measurements of the QSOFEED sample. AL thanks Maria Vittoria Zanchettin for helpful discussions. We finally thank the anonymous referee for useful and constructive suggestions.

References

Appendix A: Kinematics maps

In Section 4.3 we present and describe the kinematic maps of J1430. Here we do the same for the rest of the sample, comparing our results with previous studies. The outflow orientations and maximum extensions are measured as described in Section 4.4. For a more detailed description of the five QSO2s, including their molecular gas content, distribution, and kinematics, we refer the reader to Ramos Almeida et al. (2022).

A.1. J1010

This QSO2 is hosted in an early-type galaxy which is interacting with a smaller companion located at 13 kpc to the southwest. The CO(2-1) morphology of the QSO2, studied in Ramos Almeida et al. (2022) at 0.8″ resolution, shows two peaks that are 1.25 kpc apart. This peculiar CO morphology is most likely produced by gas removal and/or excitation in the northsouth direction. We show the kinematic maps of J1010 in Fig. A.1. Its [O III] emission is round and centrally peaked, having a diameter of ∼4″ (7.2 kpc), directly estimated from the kinematic maps. The vel_p map shows gas rotation, with negative velocities in the southeast and positive velocities to the northwest. This rotation pattern was first reported by Harrison et al. (2014), who used Gemini/GMOS observations of the central 5×3.5 arcsec² of the QSO2. The vel₀₅ map shows the highest velocities towards the northwest, and the vel₉₅ map in the opposite direction. However, for this QSO2 we find that the [O III] emission is not spatially resolved (see Section 4.2), explaining the detection of high values of W₈₀ (≥1200 km s⁻¹) across almost the entire FOV (see bottom left panel of Fig. A.1). Therefore, based on our analysis we find that our seeing-limited data does not allow us to resolve the ionized outflow or disentangle its orientation. However, Ramos Almeida et al. (2022) inspected archival VLT/MUSE data (PI: G. Venturi, 0104.B-0476) and reported higher [O III] FWHMs in the northwest region of the galaxy.

Fig. A.1. Same as in Fig. 5, but for J1010. The ∼1″ resolution 6 GHz VLA contours at 3,10,15,20, and 25σ from Jarvis et al. (2019) are superimposed to the vel05, vel95, and W80 maps.

A.2. J1100

J1100 is a barred spiral galaxy, moderately inclined (i = 38°), the south being the near side (Ramos Almeida et al. 2022). Fischer et al. (2018) presented HST/ACS imaging observations of J1100 that show rounded [O III] emission, of 1.1″ (∼2 kpc) in radius, with a small tail towards the southeast. Our GTC/MEGARA kinematic maps, shown in Fig. A.2, also show a round [O III] morphology, with a radius of ∼3″ (∼5.5 kpc).

Fig. A.2. Same as in Fig. 5, but for J1100. The ∼1″ resolution 6 GHz VLA contours at 3,10,15,30, and 40σ from Jarvis et al. (2019) are superimposed to the vel05, vel95, and W80 maps.

The velocity peak map shows some deviations from regular rotation, also reported from GMOS/IFU data by Harrison et al. (2014), which are most likely due to the stellar bar present in this galaxy. Despite that, we can distinguish positive velocities to be dominant in the northeast region and negative in the southwest. This rotation pattern coincides with the kinematics of the cold molecular gas studied with ALMA at 0.2″ resolution (Ramos Almeida et al. 2022).

Harrison et al. (2014) reported broad [O III] profiles across the entire FOV and velocity blueshifts increasing towards the southeast, but with large uncertainty. According to our analysis, presented in Section 4.2, the [O III] emission is spatially resolved, and thus we can clearly identify high blueshifted velocities (≥1200 km s⁻¹) on the northeast side of the galaxy that correspond to broad line profiles (W₈₀ ≥ 1600 km s⁻¹). Using HST/STIS data of J1100 and a slit orientation of −19°, Fischer et al. (2018) reported the presence of blueshifted velocities to the southeast, with FWHM of 1780 km s⁻¹. This is in agreement with our results if we look at the values of the kinematic maps at PA=−19°, but thanks to the IFU capabilities of MEGARA, we can see from the maps that the preferential orientation of the outflow (i.e., the maximum velocities and FWHMs) is to the northeast, with PA = 63° and R_out = 2.8″ (∼5.1 kpc). Considering this new information, for the approaching side of the outflow to be detected in the north, which is the far side of the galaxy, it must subtend a large angle relative to the galaxy disk.

A.3. J1356

J1356 is a merging system with two nuclei, north and south, observed in the optical (Greene et al. 2012; Harrison et al. 2014) and in molecular gas (Sun et al. 2014; Ramos Almeida et al. 2022). From our GTC/MEGARA data we find the north nucleus to be the brightest, as it is also the case for CO. According to Ramos Almeida et al. (2022), it contains ∼55% of the total molecular mass in the system. The kinematic maps shown in Fig. A.3 are then centered at the [O III] peak of the north nucleus.

Fig. A.3. Same as in Fig. 5, but for J1356. The red cross in the top left panel marks the AGN position, that is offset by one pixel from the [O III] peak (indicated by the white cross). The ∼1″ resolution 6 GHz VLA contours at 3,10,15,30, and 60 σ from from Jarvis et al. (2019) are superimposed to the vel05, vel95, and W80 maps.

The ionized gas shows positive velocities northeast of the north galaxy, and negative to the southwest, but the major axis is not aligned with the CO major axis (PA = 110°) reported by Ramos Almeida et al. (2022). Overall, the kinematics are complex and disturbed because of the ongoing merger. In the case of this QSO2 we detect both the approaching and receding sides of the outflow, with higher values of vel₀₅ and vel₉₅ corresponding to an enhancement in W₈₀. We detect the [O III] outflowing “bubble” proposed by Greene et al. (2012), which is expanding towards the south up to ∼12 kpc from the north nucleus (maximum R_out∼ 12.6 kpc), with a PA∼197°. In fact, this outflow radius is limited by the size of the GTC/MEGARA FOV, but according to Greene et al. (2012), the size of the bubble is ∼20 kpc. The widths measured in this south region are W₈₀ ≈ 800 − 1000 km s⁻¹. A receding outflow component is also detected to the south with a similar PA (∼200°) but less extended, up to 6.8 kpc, and with the bulk of the flux confined within 2.5″ (∼5.5 kpc). We find this redshifted outflow component to be most likely accelerated by extended radio-structure detected in VLA at 1″ resolution (see left panel of Fig. 7), as discussed in Section 5.1. The blueshifted bubble and the redshifted outflow might correspond to different outflow episodes, considering their different spatial and dynamical timescales (we measure ∼ 20 Myr for the blueshifted bubble and ∼15 Myr for the redshifted outflow). Ramos Almeida et al. (2022) also reported redshifted outflow velocities for the molecular outflow that they detected in this QSO2, albeit much more compact (R_out = 0.4 kpc; t_dyn = 1.4 Myr).

A.4. J1430

The maps corresponding to J1430 are shown in Fig. 5. This QSO2 has been extensively studied in different wavelength ranges, and it is known as the Teacup. As mentioned in Section 4.3, our kinematic maps are similar to those presented in Harrison et al. (2014), where Gemini/GMOS data were used, in Harrison et al. 2015, based on VLT/VIMOS data, in Ramos Almeida et al. (2017), based on VLT/SINFONI data, and most recently in Venturi et al. (2023), using VLT/MUSE observations. In the vel_p map, we see positive velocities to the northeast and negative to the southwest, having a total velocity gradient of 600 km s⁻¹ (see the top right panel of Fig. 5). Ramos Almeida et al. (2017) reported a spatially resolved ionized outflow with a PA∼72 − 74°, which is comparable to the PA that we measure in this study (≈65°) for the approaching side of the outflow, which exhibits velocities reaching up to -1000 km s⁻¹. This is consistent with the blueshifted outflow reported by Harrison et al. (2015), oriented as the central small-scale radio jet. Venturi et al. (2023) reported an enhancement of the velocity dispersion perpendicular to this radio jet, also found for the molecular gas (Ramos Almeida et al. 2022; Audibert et al. 2023). Here we confirm these results and, furthermore, our observations reveal the receding side of the outflow, extending in the opposite direction, with a PA of approximately 200°. This was first reported by Keel et al. (2017), who measured velocities of up to ±1000 km s⁻¹ in the same area where we detect blue and red wings. According to Ramos Almeida et al. (2022), considering that the east is the far side of the galaxy, for the ionized outflow to be blueshifted to the east, it must subtend a large angle relative to the galaxy disk.

The CO(2-1) morphology of J1430 exhibits a double peaked structure perpendicular to the radio jet, as reported by Ramos Almeida et al. (2022). Differently, the optical emission displays a single peak, as observed in previous studies by Keel et al. (2012), Harrison et al. (2015), Ramos Almeida et al. (2017), Venturi et al. (2023) and our own work. Ramos Almeida et al. (2022) suggested that this difference in the distribution of molecular gas compared to ionized gas may be attributed to the combined action of the molecular outflow and the jet. They found the cold molecular outflow to be redshifted to the east and blueshifted to the west, that is the opposite trend compared to the ionized outflow, where blueshifted velocities are detected to the nort-east and redshifted to the southwest. However, it noteworthy that Ramos Almeida et al. (2022) studied the molecular outflow emission along the minor axis only. In the more detailed study of Audibert et al. (2023), the blue side of the outflow would be oriented to the southeast, and the red to the northwest. Therefore, the ionized and cold molecular outflow in J1430 do not share either orientation or PA.

A.5. J1509

The [O III] morphology of J1509, a barred galaxy whose near side is the north, shows a compact disk of ∼1.5″ (3 kpc) in radius, as directly measured from the MEGARA maps shown in Fig. A.4, and a fish tail-like feature to the southeast. This is the only QSO2 in the sample not included in the kinematic study of Harrison et al. (2014) and lacking VLA radio observations. Indeed, to the best of our knowledge, these are the first optical IFU observations of this QSO2. The kinematic maps shown in Fig. A.4 reveal gas rotation, with negative velocities to the east and positive to the west. This coincides with the kinematics of the molecular gas studied with ALMA by Ramos Almeida et al. (2022).

Fig. A.4. Same as in Fig. 5, but for J1509.

A multiphase outflow was reported in this QSO2 by Ramos Almeida et al. (2019) using near-infrared GTC/EMIR long-slit observations. They detected blueshifted ionized and warm molecular gas, with FWHM ∼ 1500–1700 km s⁻¹. The slit was placed roughly along the minor axis of the galaxy, PA = −16°. Our vel₀₅ map, shown in the middle-left panel of Fig. A.4, indeed shows maximum blueshifted velocities to the northwest, which is also the region with the highest values of W₈₀ = 1500–1600 km s⁻¹. Therefore, we confirm the detection of the approaching side of the outflow to the northwest, with a PA ∼ 40° and R_out = 3.8 kpc. The ALMA CO(2–1) observations also show a very compact and coplanar molecular outflow, blueshifted to the north, and redshifted to the south. We do not observe the redshifted side of the ionized outflow, since the vel₉₅ map does not show any structure ascribable to an outflow. If present, it must be hidden by the galaxy disk. Thus, for the blueshifted side of the ionized outflow to the detected in the north, which is the near side, the outflow must subtend a small angle with the galaxy disk.

All Tables

All Figures

Fig. 1. Orange and red histograms correspond to all QSO2s in the QSOFEED sample and to the five QSO2s studied here. From left to right and from top to bottom: luminosity distance, AGN luminosity, optical morphology (L=lenticular, E=elliptical, S=spiral, M=merger, I=interacting), integrated radio luminosity from FIRST 1.4 GHz data, W80, and vel05 of the [O III] line profile measured from SDSS optical spectra of the targets (Bessiere et al., in prep.). SDSS J134733.36+121724.3 was excluded from the W80 and vel05 histograms because of its extreme kinematics (W80 = 2500 km s−1 and vel05 = −2500 km s−1).

In the text

Fig. 2. Fits of the [O III] doublets observed in the nuclear spectra extracted with an aperture of 1.2″ in diameter. The observed spectra are shown in blue, the model including all the components in black, and the residuals from the fit in brown. Same colors correspond to the same physical components: narrow (n; orange dashed), intermediate (i; green dot-dashed), and broad (b; red solid). Noisy pixels not considered in the fits are marked with red crosses.

In the text

Fig. 3. Residuals obtained from the PSF-subtraction analysis performed across the MEGARA FOV. The left panels show the residual maps generated by subtracting the normalized model of the central spectrum, extracted in the five spaxels around the peak of the [O III] emission (black or white crosses), from the spectra of all the spaxels. North is up and east to the left. The black squares indicate the excess regions (≥3σ and in 0.6″ × 0.6″ apertures) from which we extract the corresponding spectra shown in the right panels. The dot-dashed and the dashed lines correspond to the [O III]λ4959 and 5007 Å centroids of the central model. For J1509, the grey band in the bottom right panel indicates bad pixels that were excluded when collapsing the spectral range shown on each panel to build the residual maps.

In the text

Fig. 4. Example of nonparametric analysis of a random spaxel of J1430. Left panel: [O III]λλ5459,5007 Å continuum subtracted emission lines, with the observed spectrum shown in blue, the best fit model in black, and the corresponding residuals in grey. The same number of Gaussians (four in this case) with the same kinematics and a flux ratio of 2.98 are fit for each line of the doublet. Right panel: corresponding model of the [O III]λ5007 Å line, obtained from the fit shown in the left panel, with W80 indicated as the pink shaded region, vel05 and vel95 indicated by the pink vertical lines, and the upper (vel02 and vel98) and lower bounds (vel10 and vel90) marked by the dot-dashed vertical gray lines. The dashed vertical black line represents the velocity at the [O III] peak (velp).

In the text

Fig. 5. Kinematic maps obtained from the nonparametric analysis of the [O III] emission line profiles of J1430. Top panels: integrated flux on the left, with the ∼1″ resolution 6 GHz VLA contours at 3,5,10,20, and 60σ from Jarvis et al. (2019) superimposed (white solid lines), and velp on the right, which traces the narrow component. Middle panels: vel05 (left) and vel95 (right), which trace blueshifted and redshifted outflowing gas. The same 6 GHz contours shown in the flux map are shown in cyan, and the ∼0.25″ resolution 6 GHz VLA contours at 3, 5, 15, 30, and 60σ from Jarvis et al. (2019) are shown in solid white. The black contours correspond to the brightness temperature ratio (T32/T21) from ALMA CO observations at 0.5″ resolution (Audibert et al. 2023). Bottom panels: W80 on the left, with the same contours as in the vel05 map, and the number of Gaussians fit in each spaxel on the right. White and black crosses indicate the AGN location, defined as the maximum of the [O III] emission. North is up and east to the left.

In the text

Fig. 6. [O III] flux maps corresponding to high-velocity gas only, with velocities between −2000 < vout <  − 200 km s−1, −2000 < vout < − 400 km s−1, and 400 < vout < 1800 km s−1 depending on the source. The black solid line corresponds to the flux contour at 3σ and the red dashed line to the ellipse fit to it. The white cross indicates the [O III] emission peak, and the red cross corresponds to the most distant point from the white cross within the fit ellipse. The distance between the two crosses is the maximum outflow radius (Rout) reported in Table 4. The two orange lines are the flux contours at the 2σ and 4σ levels used to estimate the errors of the outflow radius and flux.

In the text

Fig. 7. Continuum subtracted [O III] flux maps of J1356 and J1430 extracted from the cubes in different velocity intervals, indicated in the top right corner of each panel. The black cross corresponds to the [O III] peak, the black contours to the 6 GHz VLA contours at ∼1″ resolution (at 3, 10, 15, 30, and 60σ for J1356 and at 3, 5, 10, 20, and 60σ for J1430), and the white contours to the 6 GHz VLA contours at ∼0.25″ resolution (at 3, 5, 15, 30, and 60σ).

In the text

Fig. 8. Mass outflow rate versus AGN luminosity. The purple line corresponds to the linear fit of the ionized outflow mass rates calculated by Fiore et al. (2017) for a sample of AGN at z < 0.5 (purple squares) assuming a value of the density of 200 cm−3. The black dot-dashed line corresponds to the linear fit from Davies et al. (2020), that includes the mass rates of Fiore et al. (2017), but corrected by using individual electron densities. The other symbols correspond to the four QSO2s with spatially resolved [O III] outflows, with different shapes and colors indicating the values of ne adopted to perform the outflow mass rate calculations: green circles correspond to densities of 200 cm−3, and cyan squares and light-green triangles to those measured from the [S II] doublet and the trans-auroral lines, respectively. The outflow mass rates of the QSO2s vary up to 2 orders of magnitudes depending on how the ne is measured.

In the text

Fig. 9. Ionized and cold molecular gas mass outflow rates versus AGN luminosity. The purple and cyan squares are the same as in Fig. 8, and the pink circles are the cold molecular outflow measurements from Fiore et al. (2017). The corresponding linear fits are indicated with the same colors. The yellow circles are the cold molecular mass outflow rates reported by Ramos Almeida et al. (2022) for the same QSO2s analyzed in this work, and the grey circles AGN-dominated ULIRGs from Lamperti et al. (2022).

In the text

	Fig. A.1. Same as in Fig. 5, but for J1010. The ∼1″ resolution 6 GHz VLA contours at 3,10,15,20, and 25σ from Jarvis et al. (2019) are superimposed to the vel05, vel95, and W80 maps.
In the text

	Fig. A.2. Same as in Fig. 5, but for J1100. The ∼1″ resolution 6 GHz VLA contours at 3,10,15,30, and 40σ from Jarvis et al. (2019) are superimposed to the vel05, vel95, and W80 maps.
In the text

	Fig. A.3. Same as in Fig. 5, but for J1356. The red cross in the top left panel marks the AGN position, that is offset by one pixel from the [O III] peak (indicated by the white cross). The ∼1″ resolution 6 GHz VLA contours at 3,10,15,30, and 60 σ from from Jarvis et al. (2019) are superimposed to the vel05, vel95, and W80 maps.
In the text

	Fig. A.4. Same as in Fig. 5, but for J1509.
In the text