© ESO 2022

1. Introduction

Different modes of active galactic nucleus (AGN) feedback are considered to be important processes driving the evolution of massive galaxies, by regulating black hole and galaxy growth (Di Matteo et al. 2005; Harrison 2017). If AGN feedback is not considered in cosmological simulations, gas can cool efficiently and galaxies keep forming stars, growing too big (e.g., Dubois et al. 2016). Additionally, the predicted galaxy-halo mass relations do not match observations either (Silk & Rees 1998; Croton et al. 2006; Moster et al. 2010). Currently there is plenty of observational evidence showing that AGN or black hole feedback has an impact on very different scales, which go from the central tens to hundreds of parsecs (e.g., García-Burillo et al. 2021) to hundreds of kiloparsecs (e.g., Rupke et al. 2019; Martín-Navarro et al. 2021). We now need to understand how AGN feedback works in relation to AGN and host galaxy properties and the coupling between the two while considering the short timescales associated with nuclear activity (Martini et al. 2004; Hickox et al. 2014).

The most widely studied manifestation of AGN feedback are gas outflows. Radio jets and/or AGN winds can drive multiphase outflows (Mukherjee et al. 2018; Wylezalek & Morganti 2018; Jarvis et al. 2019) that have an impact on the star formation efficiency of galaxies, as they can remove, heat up, and/or disrupt the gas available to form stars. An important source of uncertainty affecting outflow studies is that the contribution from the different gas phases entrained in the winds has not been determined in representative AGN samples (Cicone et al. 2018). If the neutral atomic, ionized, and molecular gas phases are accounted for, AGN-driven outflows might be massive and energetic enough to regulate star formation, at least in the central region of galaxies (e.g., Sánchez et al. 2018). In the case of low-to-intermediate luminosity AGN (log L_bol ∼ 42–46 erg s⁻¹), the molecular gas phase appears to be dominant in terms of mass over the other interstellar medium (ISM) phases (e.g., Fiore et al. 2017; Fluetsch et al. 2019, 2021). However, at log L_bol > 46 erg s⁻¹, ionized outflows might be as massive as their molecular counterparts (Fiore et al. 2017; Bischetti et al. 2019a; Fluetsch et al. 2021), although this has to be explored in larger AGN samples with multiphase outflow measurements. In Ramos Almeida et al. (2017, 2019) we used near-infrared (NIR) spectroscopy of nearby (z ∼ 0.1) obscured quasars with log L_bol = 45.5–46.5 erg s⁻¹ to study quasar-driven outflows using ionized (Paα and [Si VI]) and warm molecular emission lines (H₂). We found that the outflow properties were different in the two gas phases, as previously claimed by Rupke & Veilleux (2013) for the same type of source. Ionized outflows are faster, but the molecular outflows appear to carry the bulk of the mass (e.g., Fluetsch et al. 2021). However, estimating total molecular outflow masses from warm molecular gas depends on assuming warm-to-cold gas mass ratios, which are affected by large uncertainties (Ramos Almeida et al. 2019).

Cold molecular outflows in the nearby universe have been studied in low-to-intermediate luminosity AGN (i.e., Seyfert galaxies with L_IR < 10¹² L_⊙) and ultra-luminous infrared galaxies (ULIRGs), mainly, but not only, using hydroxyl (OH) and carbon monoxide (CO) transitions (see Veilleux et al. 2020 for a recent review). Focusing on CO-based studies, it has been found that the two types of objects show very different outflow properties, going from maximum velocities of ∼100–200 km s⁻¹ and outflow rates of ≲10 M_⊙ yr⁻¹ in Seyfert galaxies (e.g., Alonso-Herrero et al. 2019; Domínguez-Fernández et al. 2020; García-Bernete et al. 2021), to > 300 km s⁻¹ and hundreds of M_⊙ yr⁻¹ in ULIRGs (e.g., Cicone et al. 2014; Pereira-Santaella et al. 2018; Fluetsch et al. 2019, 2021; Herrera-Camus et al. 2020). It is noteworthy that jetted Seyfert galaxies like NGC 1068 and IC 5063 (García-Burillo et al. 2014, 2019; Morganti et al. 2015) show much higher outflow velocities and masses than those without, indicating the strong influence that jets might have in launching/accelerating molecular outflows (Wylezalek & Morganti 2018; Jarvis et al. 2019). Indeed, even in lower radio-power Seyfert galaxies, small-scale jets can drive nuclear molecular outflows (e.g., Aalto et al. 2016; Audibert et al. 2019; Fernández-Ontiveros et al. 2020; García-Bernete et al. 2021). Studying molecular outflows in nearby quasars, which have higher AGN luminosities than ULIRGs, and are hosted in galaxies of different morphologies (Bessiere et al. 2012; Pierce et al. 2022), might help us to identify which factors, including AGN luminosity, are relevant for producing more or less massive molecular outflows.

In the most powerful quasars and ULIRGs in the local universe, AGN-driven outflows have the potential to deplete the host galaxies of molecular gas on timescales of ∼10–50 Myr (Feruglio et al. 2010; Cicone et al. 2014; Bischetti et al. 2019a). However, as constrained from simulations and observations, only a small percentage of the outflowing gas is capable of escaping the galaxy (Negri & Volonteri 2017; Fluetsch et al. 2019), with the rest raining back down onto it. This result remains valid at high redshifts (z ∼ 6) and at high AGN luminosities (e.g., Carniani et al. 2019; Bischetti et al. 2019b). Regardless, it has been shown that AGN feedback is capable of influencing the distribution of molecular gas in the central kpc of nearby galaxies (Rosario et al. 2019; García-Burillo et al. 2021) and/or injecting energy in the haloes, preventing hot gas from cooling (Trussler et al. 2020). This gentler form of AGN feedback, sometimes referred to as “maintenance mode”, could be enough to regulate galaxy growth.

Here we explore for the first time the cold (T_k < 100 K) molecular gas content traced by the 2–1 line of CO, and adjacent continuum emission, of a sample of nearby (z ∼ 0.1) type-2 quasars (QSO2s) at an angular resolution of ∼0.2″ (370 pc). QSO2s (L_[OIII] > 10^8.5 L_⊙; Zakamska et al. 2003; Reyes et al. 2008) are optically obscured quasars (i.e., buried AGN) showing narrow emission lines (FWHM < 2000 km s⁻¹) of high equivalent widths. They have been extensively studied in the optical (e.g., Zakamska & Greene 2014; Villar-Martín et al. 2014; Woo et al. 2016) and to a much lesser extent in the NIR (Rupke & Veilleux 2013; Villar Martín et al. 2015; Ramos Almeida et al. 2017, 2019). QSO2s are excellent laboratories to search for outflows and study their influence in their host galaxies in the optical and NIR, because the strong AGN continuum and the broad permitted lines produced in the broad-line region are naturally obscured by dust. This makes it easier to detect (1) broad components associated with outflowing gas and (2) stellar absorption features, required to study the stellar populations of the host galaxies. In addition, QSO2s might constitute a crucial phase in the coevolution of AGN and their host galaxies in which the AGN is clearing up gas and dust to eventually shine as a type-1 QSO (e.g., Sanders et al. 1988; Hopkins et al. 2008; Hickox et al. 2009). However, in contrast with what would be expected from this scenario, Shangguan & Ho (2019) showed that the gas content of nearby QSO2s appears indistinguishable from that of type-1 quasars of the same luminosity.

From unresolved millimeter observations carried out with the Plateau de Bure Interferometer, IRAM 30 m and APEX, gas masses of ∼7–25×10⁹ M_⊙¹ have been reported for small samples of QSO2s at z < 0.4 (Krips et al. 2012; Villar-Martín et al. 2013; Jarvis et al. 2020). Due to the low sensitivities available in the past, a limited number of AGN have been studied in the high-luminosity regime. Most of them are ULIRGs, and hence, CO-bright (Downes & Solomon 1998; Wilson et al. 2008; Cicone et al. 2014). The Atacama Large Millimeter/sub-millimeter Array (ALMA) now enables spatially resolved CO maps of galaxy samples to be obtained with reasonable integration times. High angular resolution observations are particularly important in the case of luminous quasars and QSO2s in particular. They are commonly found in galaxy groups (Ramos Almeida et al. 2013) and are often seen to be interacting with other galaxies (Bessiere et al. 2012; Pierce et al. 2022). These companion galaxies might be included in the large apertures of single-dish observations, resulting in larger molecular gas masses/fractions. Furthermore, the high angular resolution is key for identifying molecular outflows, since even in nearby ULIRGs they appear generally compact (r ∼ 1–2 kpc; e.g., Feruglio et al. 2010; Cicone et al. 2014), although these radii can extend to a few kpc if low surface brightness components are considered (Feruglio et al. 2013; Herrera-Camus et al. 2019; Cicone et al. 2020).

In Sect. 2 we describe the sample selection and properties of the QSO2s. Section 3 describes the ALMA observations and data reduction. In Sect. 4 we explain the methodology that we followed to interpret the molecular gas kinematics. Section 5 includes the results found for the individual galaxies (Sect. 5.1), the millimeter continuum emission of the QSO2s (Sect. 5.2) and their molecular gas content (Sect. 5.3). In Sect. 6 we discuss the results on the molecular gas reservoirs and molecular outflows, and in Sect. 7 we summarize the findings of this work. We assume a cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7. The measurements from other works discussed here have been converted to this cosmology.

2. Sample selection and properties

Our QSO2 sample was drawn from Reyes et al. (2008), one of the largest compilations of narrow emission line AGN. We selected all the QSO2s with L_[OIII] > 10^8.5 L_⊙ (L_bol > 10^45.6 erg s⁻¹ using the bolometric correction of 3500 from Heckman et al. 2004) and redshifts z < 0.14. These constraints leave us with a sample of 48 QSO2s, hereinafter referred to as the Quasar Feedback (QSOFEED) sample, with 45.6 < log L_bol < 46.5 erg s⁻¹ (average value of 45.9 ± 0.2 erg s⁻¹). The optical selection of the targets prevents biases in dust and gas content, unlike in the case of infrared-selected samples. QSO2s might have, in principle, higher gas and dust masses than type-1 quasars, but recent results do not support this hypothesis (Shangguan & Ho 2019). These luminous QSO2s have stellar masses ranging from 10^10.7 to 10^11.6 M_⊙ (average M_* = 10^{11.1 ± 0.2} M_⊙), calculated from 2MASS Extended Source Catalogue (XSC) K-band magnitudes. We followed the procedure applied in Pierce et al. (2022), which uses the Bell (2003) equations, extinction and k-corrected K-band magnitudes, but considering a Chabrier IMF and (B − V) = 0.95. Considering that these are type-2 AGN, we do not expect the AGN contamination of the stellar mass estimates to be high. Indeed, the recent study of QSO2 hosts by Shangguan & Ho (2019), in which they performed SED fitting for a sample of 86 optically selected QSO2s at z < 0.5, shows that only five objects show strong hot dust emission in the NIR, one of them being one of our targets (J1100; see Tables 1 and 2). Indeed, the average stellar mass of the QSOFEED sample is similar to the value of M_* = 10^{10.9 ± 0.2} M_⊙ reported by Shangguan & Ho (2019).

From the QSOFEED sample we selected a subset of 7 QSO2s with redshifts 0.07 ≤ z ≤ 0.12. These QSO2s are representative of the whole sample in terms of AGN and radio luminosity, stellar mass, galaxy morphology and ionized outflow properties. These and other properties are listed in Tables 1–3. The 7 QSO2s show a range of different galaxy morphologies in the optical SDSS images (see Fig. 1) including barred spirals, early-type galaxies (ETGs) and interacting, merging and post-merger systems (see Table 2). The two seemingly undisturbed ETGs in the sample, J0232 and J1152, show rest-frame colors M_u − M_g ∼ 1.4. These are much redder and closer to typical red-sequence galaxies (M_u − M_g > 1.5; Blanton 2006) than the other five QSO2s, which show (M_u − M_g)∼0.95–1.18. In the following, we refer to J0232 and J1152 as red ETGs. J1152 has been revealed as a post-merger system from deep optical imaging (see Pierce et al. 2022). Thus, four of the seven QSO2s (57%) show disturbed morphologies that are indicative of a past interaction. This is consistent with the percentage of disturbance found for the QSOFEED sample based on deep optical imaging taken with the Isaac Newton Telescope (INT), in La Palma, which is 65% (Pierce et al., in prep.).

Fig. 1. SDSS gri color composite images of the seven QSO2s. Their optical morphologies include barred spiral galaxies (J1100 and J1509), red ETGs (J0232 and J1152) and interacting, merging and post-merger galaxies (J1010, J1356, and J1430). North is up and east to the left. The green horizontal bars at the top left of each panel correspond to 5″, which, at the average redshift of the targets (z = 0.1), corresponds to ∼9 kpc. The images are 40″ × 40″ (74 × 74 kpc2).

All the QSO2s are luminous infrared galaxies (LIRGs; see Table 2 and Fig. 2) except the red ETGs. From their FIR luminosities (L_FIR), calculated from 60 and 100 μm fluxes (Helou et al. 1985), we estimated the IR luminosities (L_IR; 8–1000 μm) by multiplying L_FIR by a factor 1.82 ± 0.17. We determined this value using the sample of LIRGs at z < 0.1 presented in Greve et al. (2014), for which both L_FIR and L_IR are available. From the total sample of 68 targets, we selected the 43 LIRGs with 10^10.2 L_⊙ ≤ L_FIR ≤ 10^11.5 L_⊙, which is the FIR luminosity range covered by the 7 QSO2s studied here. We then measured the star formation rates (SFRs) of the QSO2s using Eq. (4) in Kennicutt (1998), corrected to a Chabrier IMF (i.e., dividing by a factor of 1.59; Chabrier 2003). SFRs estimated from FIR fluxes are more appropriate for massive galaxies with high dust contents, as is the case for our QSO2s (Bell 2003; Shangguan et al. 2020a). In addition, the FIR emission of quasars is less affected by the AGN contribution than the mid-infrared or the radio, and it is not attenuated by dust, as is the case for the ultraviolet and Hα fluxes. All the QSO2s but the red ETGs have SFRs = ∼12–69 M_⊙ yr⁻¹, which place them between 0.75 and 1.5 dex above the main sequence (MS) of local SDSS DR7 star-forming galaxies (Saintonge et al. 2016). For J1010, J1100, J1356, and J1430, Jarvis et al. (2020) reported SFRs of 35, 34, 84 and 8 M_⊙ yr⁻¹, calculated from the IR luminosity due to star formation, excluding the AGN contribution. These values are in good agreement with ours (30, 34, 69, and 12 M_⊙ yr⁻¹; see Table 2), and thus we are confident that the AGN contribution to the FIR emission of these QSO2s is small. The case of the red ETGs might be different because dust heating from older stars or AGN might be responsible for a non-negligible fraction of the FIR emission (Kennicutt 1998), which is much lower than in the case of the LIRG QSO2s. For J0232 and J1152 we measure SFRs of 3 and 3.7 M_⊙yr⁻¹, which place them 0.2 and 0.3 dex above the MS. These are high SFRs for red ETGs, and indeed, the fit of the IR SED of J1152 reported in Shangguan & Ho (2019) reveals an important contribution from AGN-heated dust to the FIR, unlike in the cases of J1100, J1356, and J1430.

Fig. 2. Radio (1.4 GHz) versus IR luminosities (8–1000 μm) for the QSO2s in our sample. The radio-FIR correlation of star-forming galaxies from Bell (2003) is shown as a solid line (slope of 1.05 ± 0.04), with a correlation scatter of 0.26 dex. All the QSO2s in our sample show a radio excess. Different colors indicate whether the QSO2s are red ETGs, spirals, or interacting galaxies (see Table 2).

The black hole masses and Eddington ratios reported in Table 1 are from Kong & Ho (2018) and were estimated using stellar velocity dispersions measured from SDSS spectra and the M_BH − σ_* relation. The black hole masses range from 10^7.5 to 10^8.6 M_⊙ and the Eddington ratios (f_Edd) between 0.1 and 1, which are around the median values (10^8.2 M_⊙ and 0.2) reported by Kong & Ho (2018) from 669 of the QSO2s in Reyes et al. (2008). Thus, our seven QSO2s are near-Eddington to Eddington-limit obscured AGN in the local universe, unlike Seyfert 2 galaxies, which have typical Eddington ratios of f_Edd ∼ 0.001–0.1.

The influence of radio jets on the gas properties of radio-quiet AGN is another open question that requires further investigation (e.g., Villar-Martín et al. 2014, 2021; Jarvis et al. 2019, 2021). Our targets are radio-quiet in terms of their L_1.4GHz/L_[OIII] values (Xu et al. 1999), but they span a wide range of radio luminosities (log L_1.4GHz = 22.7–24.4 W Hz⁻¹; see Fig. 2 and Table 1). These values are representative of the QSOFEED sample, for which we measure an average luminosity of log L_1.4 GHz = 23.4  ±  0.7 W Hz⁻¹. These radio luminosities are intermediate between those of Seyfert galaxies (log L_1.4GHz < 22.5 W Hz⁻¹) and luminous high-excitation radio galaxies (HERGs; log L_1.4 GHz > 25.0 W Hz⁻¹). We note that the red ETGs have the lowest radio luminosities in our sample (log L_1.4 GHz = 22.7–23.0 W Hz⁻¹), whilst the other QSO2s have log L_1.4 GHz = 23.7–24.4 W Hz⁻¹ (see Fig. 2). The 7 QSO2s are well above the radio-FIR correlation of star-forming galaxies (Bell 2003), indicating an excess of radio emission unrelated to star-formation (see Fig. 2). In fact, according to Jarvis et al. (2019), in J1010, J1100, J1356, and J1430, star formation accounts for ≤10% of the radio emission. This, together with the steep radio spectra and the radio morphologies measured from 1–7 GHz VLA data, led them to conclude that jets are responsible for the high radio luminosities derived for these QSO2s.

The 7 QSO2s also cover a wide range of ionized outflow properties, as can be seen from Table 3. The red ETGs, J0232 and J1152, have slower and less turbulent outflows of ionized gas than the other QSO2s. This is not surprising considering the well-known connection between radio power and ionized gas kinematics (Mullaney et al. 2013; Zakamska & Greene 2014). AGN with log L_1.4 GHz > 23 W Hz⁻¹ are five times more likely to have [O III] line profiles with FWHM > 1000 km s⁻¹ than AGN with lower L_1.4 GHz, with the FWHM peaking at log L_1.4 GHz ∼ 24 W Hz⁻¹ (Mullaney et al. 2013). Using the SDSS spectra that are available for all the QSO2s, we measure an average FWHM = 1100 ± 450 km s⁻¹ for the broadest component of the [OIII] lines. We find FWHMs ranging from 370 to 2300 km s⁻¹. By performing a non-parametric analysis, we measure an average W80 (i.e., the width that contains 80% of the line flux) of 780 ± 310 km s⁻¹, ranging from 345 to 1800 km s⁻¹. Thus, the emission line kinematics of the 7 QSO2s studied here are representative of the QSOFEED sample.

3. ALMA observations

We targeted the ¹²CO(2−1) emission line (rest frequency 230.538 GHz) and its underlying continuum emission (rest frequency 220–240 GHz, which corresponds to λ ∼ 1.2–1.3 mm) in the 7 QSO2s with ALMA during Cycle 6 (project 2018.1.00870.S; PI: C. Ramos Almeida). Observations were performed between October 2018 and August 2019 using the C43-3 and C46-6 antenna configurations, with maximum baseline lengths of 0.5 km and 2.5 km. The typical on-source times were ∼0.3 h and ∼0.6 h per source, respectively. Band 5 or Band 6 receivers were selected depending on the redshift of each galaxy, providing us with four spectral windows of 1.875 GHz width and a spectral resolution of ∼10 km s⁻¹. One spectral window was centered on the CO(2−1) expected frequency, while the others accounted for continuum emission. We used a single pointing with a field-of-view (FOV) of 26.3″ in the case of the Band 6 observations, and 28.7–29.3″ for the Band 5 observations.

We created datacubes without spectral averaging, i.e., keeping the native spectral resolution of ∼10 km s⁻¹. The visibilities from the two configurations, which were calibrated using the CASA 5.4.0 software (McMullin et al. 2007) in the pipeline mode, were combined for all sources except for J1010, which was observed in the compact configuration only. Moreover, visibilities from all spectral windows were averaged to measure the millimeter continuum emission, by excluding the spectral region associated with the CO(2−1) emission line (|v − v_sys|< 500 km s⁻¹ for J1100, J1010 and J1509, and |v − v_sys|< 600 km s⁻¹ for the other QSO2s). To model the continuum emission next to the line, we fitted a first order polynomial in the uv plane to channels with velocities 500(600) < |v − v_sys|< 2000 km s⁻¹. The subtraction of this fit provided us with continuum-subtracted visibilities. We thus produced continuum-subtracted data-cubes by using the tclean CASA task in non-interactive mode, using natural weighting. The Hogbom cleaning algorithm was used in combination with a threshold of three times the rms sensitivity. Similarly, imaging of continuum maps was performed. The resulting rms sensitivity and beamsize of our QSO2 observations are listed in Table 4. In addition to the rms noise, there is a 10% flux calibration error. We note that the combination of the C43-3 and C43-6 configurations allowed us to reach angular resolutions of ∼0.18–0.25″ (300–500 pc at z ∼ 0.1), while recovering extended emission on scales of up to ∼6.5–8.6″ (12–16 kpc). Such a combination was fundamental to recover the total CO emission of the QSO2s. As shown in Sect. 5.3, four of the QSO2s have APEX CO(2–1) luminosities measured in an aperture of ∼28″ (Jarvis et al. 2020), which we have compared with our total CO fluxes (see Sect. 5.3.2).

In the case of J0232, the observations in the extended configuration were classified as semi-pass, and the data reduction was manually optimized by the Observatory. Specifically, the spectral index of the phase calibrator was identified from closely paired (within 1–2 days) Band 3 and Band 7 measurements as close as possible in time to the actual observing date (∼20 days). The flux of the calibrator was then extrapolated in time from the two Band 3 measurements made before and after the observations. Using this value, of 0.475 Jy at 97.475 GHz, along with the spectral index (−0.45) we established the flux scaling. These values were used as input for the gfluxscale calibration step.

4. Methodology

In order to interpret the molecular gas kinematics, which we describe in detail in Sect. 5.1, we attempted to model them with a simple disk model using ^3DBAROLO (Di Teodoro & Fraternali 2015). The aim of this modeling is to help us identify noncircular motions that could be associated with inflows or outflows (see Fig. 3 and e.g., Sirressi et al. 2019; Domínguez-Fernández et al. 2020). The fits were done following the methodology described in Alonso-Herrero et al. (2018) and Domínguez-Fernández et al. (2020). We first ran the model for each galaxy by fixing the kinematic center to the coordinates of the 1.3 mm continuum and the disk scale height to the default thin disk approximation describing a CO disk whose vertical structure is mostly unresolved at our spatial resolution. We allowed the systemic and rotation velocities, velocity dispersion, disk inclination and PA to vary and used uniform weighting. We then ran the model again fixing the systemic velocity, inclination and PA to the average values derived from the first run, so only rotation velocity and velocity dispersion are allowed to vary. By subtracting the model velocity map from the observed mean velocity field, we obtained mean-velocity residual maps that we use to investigate deviations from circular motions. We did the same with the velocity dispersion maps.

Fig. 3. Flowchart of the methodology used here for interpreting the molecular gas kinematics of the QSO2s. We note that the coplanar outflow/vertical inflow and coplanar inflow/vertical outflow degeneracy only holds if the vertical motions are spatially extended. Otherwise, line splitting would be detected.

Finally, we produced position-velocity (PV) diagrams along the kinematic minor and major axis of the CO distribution using ^3DBAROLO. These diagrams were extracted using a slit width of approximately the beam size. The PV diagrams along the minor axis are better suited for studying noncircular motions because along this axis the only gas motions with non-null projection in the plane of the galaxy are inward or outward radial motions (i.e., inflows or outflows). On the other hand, PV diagrams along the major axis serve to further characterize gas rotation. In the case of pure rotation, the PV diagram along the minor axis would just show emission around the systemic velocity with a width of tens of km s⁻¹. This width is a combination of beam smearing and cloud-cloud velocity dispersion (García-Burillo et al. 2014).

In case of detecting radial motions along the minor axis, they can either be within the virial range (i.e., on the order of the observed circular motions) or exceed it. In the latter case we would be most likely witnessing outflowing molecular gas, as the velocities would be too large for an inflow. On the other hand, if the radial motions are on the order of the rotational velocities, we cannot assign them unambiguously to purely radial outflows without a further careful scrutiny of the gas kinematics. In that case, to interpret the CO(2−1) velocity residuals and PV diagrams we need to know the PA and orientation of the galaxies relative to the plane of the sky (i.e., to determine the near and far sides using optical and/or NIR data). This can be inferred from analysis of the morphology of the gas response (in terms of trailing spirals or leading edges of the bar), and/or observed colors/dust extinction. Fortunately, the majority of the QSO2s have been observed with the Hubble Space Telescope (HST) and/or have optical/NIR integral field observations that we used for gathering this information.

In the case of barred spirals such as J1100 and J1509, we also need to know the bar sizes and derive the best guesses for the position of their corotation resonances. The primary bar regulates the distribution and kinematics of the molecular gas on scales going from hundreds of pc (the gas follows the so-called x₂ orbits) to several kpc (x₁ orbits). We need to estimate the extent of the corotation region of the bar to disentangle whether the molecular gas motions being studied are inside/outside of this corotation region. The corotation radius can be estimated from the bar radius as R_cr ∼ 1.2 ± 0.2 × R_bar (Athanassoula 1992).

Once we determine which are the far and near sides of the disks, we can easily interpret our molecular gas observation in terms of radial motions. These radial motions will correspond to inflows or outflows depending on whether they are coplanar or vertical with the CO disks. In Fig. 3 we summarize the methodology described above.

In the case of quasars, in principle we expect any molecular outflow to be mostly coplanar with the CO disk (see Fig. 4). This is supported by previous observational evidence in the case of spatially extended AGN-driven molecular outflows, for example NGC 1068 (García-Burillo et al. 2014, 2019), NGC 3227 (Alonso-Herrero et al. 2019), IC 5063 (Morganti et al. 2015; Mukherjee et al. 2018), NGC 5643 (Alonso-Herrero et al. 2018; García-Bernete et al. 2021), NGC 4388, NGC 5506, and NGC 7582 (García-Burillo et al. 2021). As shown in Fig. 4, the AGN wind, jet and/or ionized outflow produces a radial expansion in the CO disk (i.e., a molecular outflow), which will be more or less massive depending on its orientation relative to the CO disk (strong and weak coupling, scenarios A and B in Fig. 4). As we know from detailed studies of nearby AGN such as NGC 1068, some of the molecular gas is forced to leave the plane of the galaxy and adopt a rather shell-like geometry (see García-Burillo et al. 2019), but in order to interpret the CO kinematics here we consider that the bulk of the molecular outflows are either coplanar or vertical with the CO disk.

Fig. 4. AGN wind/jet/ionized outflow geometries considered in this work and the coplanar molecular outflows that they produce. Left panel: example of strong coupling (scenario A), in which the wind/jet/ionized outflow are coplanar with the disk, launching a more massive molecular outflow. Right panel: example of weak coupling (scenario B), in which the wind/jet/ionized outflow subtends a certain angle from the CO disk, having a small intersection region (shown in dark gray) and launching a more modest molecular outflow.

To derive outflow mass rates (Ṁ_out) we need to estimate outflow masses (M_out) and assume a certain outflow geometry. The outflow masses were calculated by integrating the CO(2−1) emission along the minor axis and within the regions and velocities indicated in the PV diagrams shown in Sect. 5.1. These outflow regions and velocities were selected from inspection of the PV diagrams along the minor axis and the velocity residual maps (see Sect. 5.1 for details). We assumed that the CO emission is thermalized and optically thick, so R₂₁ = CO(2−1)/CO(1−0) = 1 (Braine & Combes 1992; Solomon & Vanden Bout 2005), and α_CO = 0.8 ± 0.5 M_⊙(K km s⁻¹ pc²) ⁻¹ (Downes & Solomon 1998). The latter is more appropriate for outflowing molecular gas than the commonly assumed Galactic factor (see e.g., Morganti et al. 2015).

We then assume a time-averaged thin expelled shell geometry for deriving outflow mass rates:

(1)

This geometry is commonly adopted for molecular outflow calculations in the local universe (e.g., Audibert et al. 2019; Fluetsch et al. 2019; Lutz et al. 2020), and it is more conservative than the multi-conical outflow geometry uniformly filled by outflowing clouds. This assumption for the outflow geometry implies a constant mass rate since the outflow started (Lutz et al. 2020). In Eq. (1), M_out, v_out and r_out are the outflow mass, velocity and radius, and α is the angle between the molecular outflow and the line of sight. In the case of a coplanar molecular outflow, tan(α) = tan(90–i) = 1/tan(i), where i is the inclination angle of the CO disk. Thus, Eq. (1) corresponds to the deprojected mass outflow rate, and the tan(α) factor results from deprojecting the outflow velocity, v_out/cos(α), and the outflow radius, r_out/sin(α).

5. Results

5.1. Individual galaxies

We detect continuum emission (rest-frame λ∼1.2–1.3 mm) at > 3σ in the seven QSO2s observed with ALMA, and CO(2−1) emission in five of them. The CO moment maps were generated by integrating over the spectral range where line emission is detected above 3σ in the continuum-subtracted cubes. The continuum emission is shown as images for the QSO2s without CO detection (J0232 and J1152, for which the results are presented in Appendices A and B). For the other five QSO2s, it is shown as contours starting from 3σ overlaid on the corresponding CO moment 0 maps. To interpret the molecular gas kinematics, we followed the procedure described in Sect. 4. The regions with outflow motions have been identified in the corresponding PV diagrams.

5.1.1. SDSS J101043.36+061201.4 (J1010)

The SDSS optical image of this QSO2 resembles an ETG morphology (S0a; Willett et al. 2013) clearly interacting with a small galaxy at ∼7″ (13 kpc) SW that it is also detected in CO (see Appendix C). The SFR estimated from the IR luminosity, of 30 M_⊙ yr⁻¹, is comparable to those of the spiral galaxies J1100 and J1509, placing J1010 1.15 dex above the MS (see Sect. 2). This is likely a result of the interaction with the small companion. Despite the high SFR, J1010 is well above the radio-IR correlation of star-forming galaxies (see Fig. 2). Jarvis et al. (2019) estimated that only 2.6% of the FIRST 1.4 GHz luminosity of J1010 can be accounted for by star formation. High angular resolution VLA data (∼0.25″) at 6 GHz shows a compact and rounded morphology with a deconvolved major axis of ∼200 pc and PA ∼ 180° (Jarvis et al. 2019). The ALMA 1.3 mm continuum image shows the same morphology (see pink contours in Fig. 5), with a deconvolved size of 0.25″ × 0.24″ (450 × 434 pc²) and PA = 41 ± 69°. The position of the AGN is derived from the peak of this continuum emission.

Fig. 5. CO(2−1) moment 0, 1, and 2 maps (i.e., flux, velocity, and velocity dispersion) of J1010. The 3DBAROLO models of the moment 1 and 2 maps and corresponding residuals are also shown. Continuum contours starting at 3σ are shown in pink in the moment 0 map (σ = 0.037 mJy beam−1), and the beam size (0.80″ × 0.69″) is shown in the bottom left corner of each panel. The major axis of the CO kinematics is shown as a black dashed line in the moment 1 maps. North is up and east to the left.

The CO(2−1) moment maps of J1010 are shown in Fig. 5. The maximum radius of the CO emission is R_CO = 1.3″ (2.4 kpc). The moment 0 map reveals a double-peaked structure, with a separation of ∼0.7″ (1.25 kpc) between peaks in the E–W direction. The maximum of the 1.3 mm continuum lies in the middle of the two CO peaks. This morphology could be explained by a highly inclined disk with a ring-like in-plane distribution of molecular gas, but our ^3DBAROLO modeling of the data is not compatible with such disk orientation (see below). Alternatively, this peculiar CO morphology could be the result of AGN feedback, either produced by molecular gas removal in the N–S direction, or by CO being excited to higher levels (Rosario et al. 2019). There is no optical counterpart to this double-peaked morphology. An optical image generated from a publicly available VLT/MUSE data cube with the same angular resolution as our ALMA data (FWHM ∼ 0.79″; PI: G. Venturi) shows a single nucleus coincident with the peak of the 1.3 mm continuum.

Kinematics. The kinematics of the ionized gas were studied by Harrison et al. (2014) using Gemini/GMOS in integral field mode. They reported the presence of a compact rotating [OIII] disk with major axis PA ∼ 299°= −61° in addition to broad [OIII] wings of FWHM ∼ 1000–1200 km s⁻¹ blueshifted by −100 km s⁻¹ (see also Villar-Martín et al. 2014). Using the MUSE datacube mentioned before, we detect broader and more intense [O III] emission lines to the NW, as well as maximum blueshifts. We also generated two continuum images centered at 6200 and 8500 Å to obtain a color map (see Appendix C). In both images we detect a dust lane SW from the nucleus, which produces redder colors. According to this, the N would be the far side of the galaxy. Thus, for the approaching side of the ionized outflow to be detected in the NW, it has to subtend a large angle relative to the galaxy/CO disk (scenario B in Fig. 4).

The CO velocity field shown in Fig. 5 resembles rotation, blueshifted to the E and redshifted to the W. Indeed, this velocity map can be reproduced with a rotating disk of major axis PA ∼ 290° = −70°, in good agreement with the [OIII] kinematic major axis reported by Harrison et al. (2014) and with the galaxy major axis (−67°; see Table 2). This is indicative of rotation-dominated gas kinematics. The inclination of the rotating disk model fitted to the CO data is i ∼ 36°, also consistent with the galaxy inclination (35°; see Table 2). Thus, despite the tidal interaction evidenced by the optical image, the optical and CO disks share similar kinematics. The velocity dispersion map peaks at the position of the nucleus, showing values of up to 120 km s⁻¹. The ^3DBAROLO model accounts to some extent for the large line width feature identified close to the minor axis region resulting from the beam smearing of the motions attributable to circular rotation and turbulence. This is better visualized in the modeled minor axis PV diagram shown in Fig. 6.

Fig. 6. PV diagrams along the CO kinematic major axis (PA = 290° = −70°) and minor axis (PA = 200° = 20°) of J1010 extracted with 3DBAROLO. We used a slit width of ∼0.8″. Blue contours correspond to the observed CO(2−1) emission above 2σ and red contours to the 3DBAROLO rotating disk model (inclination of i = 36°) shown in Fig. 5. Yellow dots are the average model velocities at different radii. Blue contours outside the boundaries of the model correspond to noncircular motions. Orange lines indicate the positive/negative pattern that, in the case of this QSO2, might be produced by either a coplanar inflow or a spatially extended vertical outflow.

In addition to ordered rotation along the major axis (v_rot ∼ 200 km s⁻¹; see top panel of Fig. 6), the PV diagram along the kinematic minor axis (PA ∼ 200° = 20°) reveals noncircular motions in the central 1.6″ (2.9 kpc) of J1010, mostly blueshifted to the N and redshifted to the S (see bottom panel of Fig. 6). This can be also seen in the residual velocity map shown in Fig. 5.

Thus, in the case of this QSO2, for which the N is the far side, we are witnessing either a coplanar inflow or a vertical outflow (see Fig. 3). In the latter case, the outflow could be star-formation driven in the case of this interacting galaxy with a high SFR. An outflow could explain the double-peaked CO morphology shown in the moment 0 map. However, if this is the case, the vertical outflow must be spatially extended, as otherwise we would expect significant line splitting, something that we do not see in the PV diagrams. Due to the limited angular resolution of the ALMA observations of this QSO2, we cannot fully resolve these radial motions. Higher angular resolution data are required to further investigate this interesting system.

5.1.2. SDSS J110012.39+084616.3 (J1100)

J1100 has a SFR of 34 M_⊙ yr⁻¹ and it shows a clear radio-excess in Fig. 2, with only ∼4% of the 1.4 GHz luminosity being due to star formation as estimated by Jarvis et al. (2019). VLA 6 GHz images of this QSO2 at the same angular resolution as our ALMA data (∼0.2″) reveal a compact and rounded morphology, very similar to J1010, with a linear size of 800 pc and PA ∼ 170° (Jarvis et al. 2019). Lower resolution VLA data do not show evidence for extended emission either. The ALMA 1.3 mm continuum emission of J1100 is shown as pink contours in the moment 0 map in Fig. 7. It appears very compact at 3σ, with the emission peak indicating the AGN position. We measure a deconvolved size of 0.30″ × 0.25″ (555 × 463 pc²) with PA = 72 ± 76°.

Fig. 7. Same as in Fig. 5 but for J1100. Continuum contours starting at 3σ are shown in pink in the moment 0 map (σ = 0.017 mJy beam−1), and beam size is 0.24″ × 0.20″.

The QSO2 host is a large and moderately inclined (i = 38°) barred galaxy (SBb; Willett et al. 2013). The optical continuum SDSS and HST/ACS images available for this QSO2 show strong and wound spiral arms, and an undisturbed morphology (Fischer et al. 2018). The [OIII] emission appears very compact in the HST/ACS images, nearly circular (1.1″ radius ≃2 kpc), and showing a small tail to the SE. The ALMA CO(2−1) moment maps shown in Fig. 7 reveal a large disk of molecular gas with R_CO = 2.7″ (∼5 kpc). The CO morphology follows the stellar bar and the inner part of the spiral arms seen in the HST/ACS optical continuum image. From this image we estimated a bar radius of R_bar ∼ 2.5″ (4.5 kpc) in the E–W direction (PA ∼ 85°). In the CO maps, from the west edge of the bar we see the beginning of one of the spiral arms, which shows negative velocities. From the eastern edge of the bar we see another spiral arm showing positive velocities.

Kinematics. The [OIII] kinematics are dominated by rotation, with a kinematic major axis PA ∼ 74° (Harrison et al. 2014). Based on HST/STIS spectroscopy and using a slit orientation of −19° (almost coincident with the galaxy minor axis; see Table 2), Fischer et al. (2018) reported the presence of an [OIII] outflow of FWHM ∼ 1780 km s⁻¹, velocities of up to −350 km s⁻¹ and projected radius of ∼460 pc. The bulk of this outflowing gas is found toward the SE, in agreement with the Gemini/GMOS analysis of the [OIII] kinematics presented in Harrison et al. (2014).

The ^3DBAROLO model of the moment 1 and 2 maps and corresponding residuals are shown in the middle and right panels of Fig. 7. In the case of this QSO2 we fixed the inclination to match that of the galaxy (i ∼ 38°) to prevent the fitted model from having an unrealistic orientation of 70°. This happens because ^3DBAROLO tries to fit the most conspicuous emission, which resembles an apparently more inclined disk. The kinematic major axis of the model (PA ∼ 69°) is very similar to the major axis of the galaxy (PA ∼ 67°; see Table 2) and of the [OIII] gas (PA ∼ 74°), consistent with rotation-dominated gas kinematics.

As explained in Sect. 4, the bar regulates the CO distribution and kinematics within the corotation radius. In the case of J1100, R_CR ∼ (1.2 ± 0.2)×R_bar ∼ 3″ (5.5 kpc). The morphology of the molecular gas within this radius is in agreement with the canonical gas response to a stellar bar with a radially extended inner Lindblad resonance (ILR) region. The gas appears to concentrate along the two leading edges of the stellar bar. This is the expected pattern when gas transits from the outer x₁ orbits to the inner x₂ orbits in the presence of an ILR region. This can also be seen in the PV diagrams shown in Fig. 8. Along the major axis (top panel), part of the gas follows the x₁ orbits of the bar, that is the smooth rotation pattern modeled with ^3DBAROLO. The corotation region would correspond to offsets of ±3″ (5.5 kpc), but we see gas rotation extending even further. In the PV diagram along the major axis we also see high-velocity gas (±200 km s⁻¹) distributed in a nuclear disk (the inner 1″ ∼ 1.8 kpc), which follows the x₂ orbits of the bar.

Fig. 8. Same as in Fig. 6 but for J1100. The CO kinematic major axis is PA = 69°, the minor axis PA ∼ 159° = −21° and the inclination i = 38°. We used a slit width of ∼0.2″.

The PV diagram along the minor axis (PA ∼ 159° = −21°; bottom panel of Fig. 8) shows noncircular gas motions with maximum velocities of ±200 km s⁻¹. The bulk of these fast noncircular motions is found at ∼0.5″ (0.90 kpc) to the S (i.e., negative offsets), although forbidden velocities of less than 100 km s⁻¹ are also seen to the N and S up to a radius of ∼2″ (3.7 kpc), following the classical positive/negative pattern. All these noncircular gas motions happen well inside the corotation radius of the bar and therefore could correspond to gas inflows. In order to assess whether or not this is the case, we first need to assume that the spiral arms trail galaxy rotation (Pasha & Smirnov 1982). In this case the galaxy rotation has to be clockwise. Considering this and the CO velocity field shown in Fig. 7, the N has to be the far side. Thus, for the approaching side of the ionized outflow (Harrison et al. 2014; Fischer et al. 2018) to be detected in the S, it must be almost coplanar with the galaxy disk.

Most of the molecular gas along the minor axis appears blueshifted in the S, although there is also redshifted gas in the N. This can be also seen from the residual velocity map shown in the middle right panel of Fig. 7, and it is the opposite behavior of what it is expected from a coplanar inflow of gas induced by a stellar bar within the corotation radius. Instead, our analysis of the kinematics suggests the presence of a coplanar outflow with an average outflow velocity of ±115 km s⁻¹ and radius r_out ∼ 0.7″ (1.3 kpc). These are the average values of the regions where the molecular gas shows velocities larger than those of the rotating disk model at a given radius. We discard the vertical inflow scenario (see Fig. 3) because this is highly unlikely in an undisturbed disk galaxy such as J1100. Besides, we do not see any evidence of line splitting.

By integrating the CO emission along the minor axis within these regions (see Fig. 8), we estimate an outflow gas mass of 10.5 × 10⁷ M_⊙ and an outflow rate of 9.5 M_⊙ yr⁻¹. Considering the deprojected outflow rate, 12.2 M_⊙ yr⁻¹ and the SFR, of 34 M_⊙ yr⁻¹, we estimate a mass loading factor of η = Ṁ_out/SFR ∼ 0.3. Thus, the CO kinematics in the central kpc of J1100 are peculiar. Overall the molecular gas follows the canonical response to the bar (i.e., falling inward), but the quasar-driven wind and ionized outflow perturb the molecular gas in the disk and drive it outward.

5.1.3. SDSS J135646.10+102609.0 (J1356)

This QSO2 is hosted in a spectacular merger system showing a completely distorted optical morphology. There are two nuclei (N and S), separated by only 1.1″ (2.4 kpc). The N nucleus is the QSO2, whose host galaxy is the dominant member of the merger in both the optical and molecular gas (Sun et al. 2014). It is a massive ETG with PA = 156° and i = 55° (see Table 2). Although the underlying stellar population is dominated by old stars (Greene et al. 2009), the SFR derived from its IR luminosity is the largest in our sample (69 M_⊙ yr⁻¹; see Table 2), likely a consequence of the ongoing merger. Other components of the merging system, as identified from the ALMA and HST morphologies by Sun et al. (2014), are the S nucleus and the western arm (W arm hereafter). The latter is a huge stellar feature containing a large percentage of the molecular gas in the galaxy.

We detect the N nucleus in CO(2−1), which is shown in the top panel of Fig. 9, as well as the W arm. Part of the latter is seen in the moment 0 map as blueshifted emission westward of the N nucleus. This feature extends up to a maximum distance of ∼2.3″ (5 kpc) from the quasar. By integrating the CO emission at velocities < 0 km s⁻¹ we detect the S nucleus and also the W arm (see Appendix D and Sun et al. 2014). The N galaxy represents ∼55% of the total molecular gas mass in the system, the W arm ∼32%, and the S nucleus ∼13%.

Fig. 9. Same as in Fig. 5 but for J1356 (N nucleus). Part of the W arm is also shown and labeled in orange in the moment 0 map. Continuum contours starting at 3σ are shown in pink in the moment 0 map (σ = 0.015 mJy beam−1), and beam size is 0.25″ × 0.24″.

The 1.3 mm continuum contours at 3σ follow the CO emission, but they also show narrower structures extending up to ∼0.75″ (1.7 kpc) NE and ∼0.5″ (1.1 kpc) NW. These structures have different orientations than the major axis of the 6 GHz continuum emission (PA ∼ 20°, coincident with the kinematic minor axis of the CO disk; see below) measured from low-resolution VLA data (∼1″; Jarvis et al. 2019). This 6 GHz emission has the appearance of a bent jet of 5.6 kpc to the SW (see Fig. 5 in the previously mentioned work). The VLA data at ∼0.25″ resolution shows a fairly compact morphology of ∼300 pc. From our ALMA continuum image we measure a deconvolved size of 0.27″ × 0.20″ (600 × 440 pc²) with PA = 61 ± 17°. Using the VLA fluxes reported in Jarvis et al. (2019) and the ALMA 1.3 mm flux we measure a spectral index α = −0.79.

Kinematics. Based on long-slit optical spectroscopy, Greene et al. (2012) reported the presence of a ∼20 kpc ionized outflow in the form of an expanding bubble. The base of this outflow would be located at ∼3.5″ (7.7 kpc) S of the QSO2 (see Fig. 1 in Sun et al. 2014). Our analysis of the CO(2−1) data does not reveal a molecular counterpart to the ionized bubble, as it was also the case for the CO(1−0) and CO(3−2) transitions reported by Sun et al. (2014).

The CO(2−1) moment 1 and 2 maps are shown in Fig. 9. The N nucleus shows high average velocity dispersion values, of up to 160 km s⁻¹. We fitted the CO kinematics with a rotating disk model of PA = 110° and i = 52°. This indicates that the galaxy and CO disk would be coplanar, but the PAs do not coincide (see Table 2), which is not surprising in an ongoing major merger. Harrison et al. (2014) analyzed the [OIII] kinematics of the central ∼3.5″ × 5″ (7.7 × 11 kpc²) of the galaxy and reported a major axis with PA = 225° = 45°. These authors found an average blueshift of −215 km s⁻¹ for the galaxy-integrated [OIII] profile, with a maximum blueshift of −600 km s⁻¹ at ∼2″ (4.4 kpc) SW of the nucleus and a projected outflow radius of ≤3.1 kpc. The outflow orientation coincides with the PA of the radio jet inferred from VLA data (20°), but it is possible that the blueshifted gas detected to the SW is associated with the W arm and the S nucleus, which are blueshifted relative to the systemic velocity (see Appendix D). We note that the analysis of the molecular gas kinematics that we are describing here focuses on the N nucleus.

Using the HST/WFC3 images from Comerford et al. (2015) we constructed a B–I color map of the galaxy. A dust lane crosses the N nucleus with a PA coincident with the CO major axis. On the small spatial scales considered in Fig. 9, we detect redder optical colors NE of the nucleus (see Appendix D). Based on this, the N would be the near side. Since the ionized outflow is detected to the SW, it must subtend a large angle relative to the CO disk (scenario B in Fig. 4, i.e., weak coupling).

The PV diagram along the minor axis (PA = 20°; see bottom panel of Fig. 10), shows tentative evidence of noncircular motions within the central ∼0.6″ (1.3 kpc) of the QSO2. We do not see the positive/negative pattern characteristic of outflowing/inflowing gas as in J1100, but this could be due to the complex gas kinematics of ongoing merger systems like this one. If the N is the near side, molecular gas outflowing in an almost coplanar geometry should be blueshifted to the N and redshifted to the S. The high-velocity gas with v_max ∼ 350 km s⁻¹ detected within a radius of 0.2″ (0.4 kpc) to the S would then correspond to outflowing gas (see also the redshifted residuals along the minor axis in Fig. 9). This is supported by the order of magnitude of the radial velocities. A vertical inflow scenario is, in principle, also compatible with having redshifted gas detected in the far side (see Fig. 3). In a major merger, a vertical inflow driven by, for example, a tidal tail would be plausible, but in that case we would detect line splitting in the PV diagrams (see e.g., García-Burillo et al. 2019), something that does not happen.

Fig. 10. Same as in Fig. 6 but for J1356. The CO kinematic major axis is PA = 110°, the minor axis PA = 200° = 20° and the inclination angle i = 52°. We used a slit width of ∼0.25″. The orange box indicates the region and velocities used to estimate the mass of the outflowing gas, which in the case of this QSO2 does not show the characteristic positive/negative pattern.

A high-velocity component of ∼400 km s⁻¹, detected in CO(3−2) in the N nucleus, was reported by Sun et al. (2014). They identified this with an outflow of M_out ∼ 7 × 10⁷ M_⊙ by assuming = . This outflow mass is larger than our measurement from CO(2−1), of (1.4 ± 1.2)×10⁷ M_⊙. We note, however, that this value has been obtained from integrating the high-velocity gas to the S and along the minor axis only (i.e., within the orange box in Fig. 10). If we consider all the redshifted high-velocity gas, we measure a mass of 7.1 × 10⁷ M_⊙, as in Sun et al. (2014). We also detect some high-velocity blueshifted gas in the PV diagram along the minor axis (see Fig. 10), but not in the residual map shown in Fig. 9. Thus, we prefer to be conservative and only use the high-velocity redshifted gas to the S along the minor axis to work out the outflow mass. Using this mass, the radius indicated above, and an average outflow velocity of 310 km s⁻¹, we estimate an outflow rate of 10.0 M_⊙ yr⁻¹. Considering the high SFR estimated for this QSO2, of 69 M_⊙ yr⁻¹, and the deprojected Ṁ_out = 7.8 M_⊙ yr⁻¹, the mass loading factor is η ∼ 0.1.

5.1.4. SDSS J143029.88+133912.0 (J1430; The Teacup)

Often referred to as the Teacup galaxy, J1430 has been widely studied in different wavelength ranges including the X-rays (Lansbury et al. 2018), optical (Keel et al. 2012; Harrison et al. 2015; Villar-Martín et al. 2018), NIR (Ramos Almeida et al. 2017) and radio (Jarvis et al. 2019). The host galaxy is bulge-dominated, showing clear signatures of a past interaction. These include a spectacular system of concentric shells and nuclear dust lanes (see Fig. 1 and Keel et al. 2015).

Using VLA data at 0.3″ resolution, Harrison et al. (2015) reported the presence of two compact steep radio sources: one corresponding to the QSO2 nucleus (HR-A) and another (HR-B) at ∼0.5″ (0.8 kpc) NE (PA ∼ 60°). These authors proposed that this could be a compact radio jet that might be accelerating the ionized gas in the central kpc and possibly driving the large-scale radio bubbles. Our ALMA data reveal continuum emission at 1.3 mm peaking in the middle of two CO(2−1) blobs, as can be seen from the top panel of Fig. 11. This continuum emission looks more jet-like than in the case of the other QSO2s in our sample and it has an almost E–W orientation, perpendicular to the CO emission. From the analysis of the continuum image we derive a deconvolved size of 0.37″ × 0.23″ (590 × 370 pc²) and PA = 80 ± 12°, almost coincident with the major axis of the 1.4 GHz continuum emission measured from FIRST (77°; Harrison et al. 2014).

Fig. 11. Same as in Fig. 5 but for J1430. Continuum contours starting at 3σ are shown in pink in the moment 0 map (σ = 0.015 mJy beam−1), and beam size is 0.21″ × 0.18″.

The morphology of the cold molecular gas at 0.2″ resolution looks very different to that of the warm molecular gas observed with VLT/SINFONI at 0.5″ resolution (Ramos Almeida et al. 2017). Instead of the single-peaked disk-like structure observed in H₂, the CO moment 0 map shown in Fig. 11 shows a double-peaked morphology, with the two peaks separated by ∼0.8″ (1.3 kpc) with PA ∼ −10°. This galaxy has been observed in the optical with HST (angular resolution of 0.1″) and it shows a single nucleus in both continuum and [OIII] (Keel et al. 2012; Harrison et al. 2015).

Kinematics. The moment 1 map shows a distorted rotation pattern, redshifted and blueshifted to the N and S, respectively (see Fig. 11). This is also the case for the H₂ (Ramos Almeida et al. 2017) and [OIII] velocity fields (Harrison et al. 2014). The moment 2 map reveals higher values of the velocity dispersion (100–120 km s⁻¹) across ∼1″ (1.6 kpc) in the direction perpendicular to the jet. This enhancement of the velocity dispersion has been observed in ionized gas in four nearby AGN from the MAGNUM survey (Venturi et al. 2021), and interpreted as due to the action of the jets perturbing the gas in the galaxy disk.

The ^3DBAROLO model of the moment 1 and 2 maps and corresponding residuals are shown in the middle and bottom rows of Fig. 11. Since it is difficult to define a clear major axis of the CO distribution in the case of this QSO2, we first let the PA vary freely and fixed the inclination to 38° (i.e., coincident with the galaxy inclination measured from optical images; see Table 2). This is a post-merger system and the inclination of the CO and stellar disk might be different, but if we try to fix the PA to the average value of the previous step (PA = 4°) and let the inclination vary, we get values between i = 35°and 41°. These values produce lower residuals than higher/lower inclinations, so we can safely assume i = 38° for the CO disk. The kinematic major axis of the CO disk (PA = 4°) is different from the [OIII] major axis derived from the analysis of Gemini/GMOS data (PA = 27°; Harrison et al. 2014) and from the galaxy’s PA of −19°. This is likely due to the past merger and/or to the action of the jet on the ionized and molecular gas.

Ionized outflows detected in the form of blueshifted broad components of optical and NIR emission lines have been reported for the Teacup (Harrison et al. 2015; Ramos Almeida et al. 2017). A spatially resolved outflow with radial extent of ∼1 kpc was found by Ramos Almeida et al. (2017) using the Paschen α and [Si VI] lines in the NIR. They reported maximum outflow velocities of −1100 and −900 km s⁻¹ and PAs of 75°and 72° respectively. The approaching side of the ionized outflow is on the NE side, within ∼1″ (1.6 kpc), as also found for [OIII] in the optical (Harrison et al. 2014, 2015). Considering this, if the ionized outflow subtends a large angle relative to the CO disk (scenario B in Fig. 4), the E would be the far side. This is confirmed by the presence of strong dust lanes westward of the nucleus (Keel et al. 2015). These dust lanes produce redder colors on that side of the galaxy, as shown by the r–i color map shown in Appendix E.

The PV diagram along the minor axis (see bottom panel of Fig. 12) shows strong noncircular motions within r_out ∼ 0.3″ (0.5 kpc) of the galaxy. These motions have maximum velocities of 250 and −180 km s⁻¹, mainly redshifted to the E and blueshifted to the W. This can be seen also from the residual map shown in Fig. 11. Since the E is the far side, this corresponds to either outflowing gas in the CO disk plane or vertical inflowing gas. As in the case of J1356, we favor the first scenario based on the high radial velocities, of up to 250 km s⁻¹, and the lack of line splitting in the PV diagrams.

Fig. 12. Same as in Fig. 6 but for J1430. The CO kinematic major axis is PA = 4°, the minor axis PA = 94° and the inclination angle i = 38°. We used a slit width of ∼0.2″. The orange boxes indicate the regions used to estimate the mass of the outflowing gas.

By integrating the high-velocity CO(2−1) emission along the minor axis and within the regions labeled in Fig. 12, we measure an outflow mass of M_out = 3.12 × 10⁷ M_⊙ and an outflow rate of 12.3 M_⊙ yr⁻¹. Considering the deprojected outflow mass rate, of Ṁ_out = 15.8 M_⊙ yr⁻¹ and the SFR, of 12 M_⊙ yr⁻¹, the mass loading factor is η ∼ 1.3, the largest in our sample.

We did not detect a warm molecular counterpart of this outflow using NIR data from VLT/SINFONI (Ramos Almeida et al. 2017). Unlike the nuclear Paschen α and [Si VI] lines, for which both narrow and broad Gaussian components were necessary to reproduce the line profiles, the H₂ could be fitted with a single Gaussian component, although slightly blueshifted (−50 km s⁻¹) relative to the narrow component of Paα. We note that the cold molecular outflow is barely resolved even at the high angular resolution of the ALMA data, making it challenging to detect its warm molecular gas counterpart with the 0.5″ resolution of the SINFONI observations. Furthermore, cold molecular gas is much more abundant than the warm molecular gas traced with the NIR H₂ lines, which might be another reason for detecting the molecular outflow in CO but not in H₂. A detailed comparison of the warm and cold molecular gas in this galaxy will be the subject of forthcoming work (Audibert et al., in prep.).

The combined action of the cold molecular outflow and the jet could explain the double-peaked morphology shown in the flux map of Fig. 11. This morphology is similar to that of the Seyfert 2 galaxy NGC 2110, which also has a radio jet oriented perpendicularly to an area depleted of CO(2–1). The molecular gas might have being pushed outward, resulting in a depletion of the central region. Alternatively, the molecular gas in the region more affected by the jet might have been excited to higher-J transitions (Rosario et al. 2019). Thus, the Teacup could be another example of AGN feedback shaping the molecular gas reservoir in the central kpc of an AGN (see also Rosario et al. 2019; García-Bernete et al. 2021; García-Burillo et al. 2021).

5.1.5. SDSS J150904.22+043441.8 (J1509)

This galaxy is classified as a barred spiral (SBa) in Galaxy Zoo (Willett et al. 2013), with PA = 94° and inclination i = 44° (see Table 2). From the r-band SDSS image we measure a bar orientation of ∼ − 30° and R_bar ∼ 2.25″ (5.6 kpc).

Unlike the majority of QSO2s in our sample, there are no high angular resolution optical and radio data of J1509 to be compared with the ALMA continuum and molecular gas distribution. The 1.3 mm continuum emission of J1509 is quite peculiar (see pink contours in Fig. 13). It follows the inner oblong structure of molecular gas, but from this central region an arc-like structure directed to the SE that does not follow the CO distribution is exhibited. By fitting a 2D Gaussian to the continuum emission at ≥3σ we measure a deconvolved size of 0.74″ × 0.55″ (1.5 × 1.1 kpc²) with PA = 135 ± 22°. This structure cannot be associated either with dust or free-free emission, since it would otherwise follow the CO morphology. It should then correspond to synchrotron emission, and the most likely explanation for the peculiar morphology is a radio jet that in its interaction with the surrounding environment has been forced to bend. Alternatively, if the jet precesses about a defined axis it can result in the jet being curved as observed in the plane of the sky (Begelman et al. 1984). The presence of a jet is further supported by the spectral index measured from the 1.4 GHz and 200 GHz fluxes, which corresponds to a steep spectrum (α = −0.64; see Sect. 5.2).

Fig. 13. Same as in Fig. 5 but for J1509. Continuum contours starting at 3σ are shown in pink in the moment 0 map (σ = 0.013 mJy beam−1), and beam size is 0.26″ × 0.24″.

The moment 0 map in Fig. 13 shows an oblong structure in CO, almost in the E–W direction and with a major axis of ∼1.7″ (3.4 kpc). Outside this compact elongated gas disk there are two spiral arms that develop inside the stellar bar out to r ∼ 2″ (4 kpc) to the NW and SE (PA ∼ −30°). This feature mimics the leading edges signature typical of the gas response to the bar in the presence of an extended ILR region. The leading edges are joint at the northern end of the stellar bar by a gas arc at r ∼ 2″. This arc lacks a southern counterpart, which gives it the appearance of an asymmetric ring.

Kinematics. In Ramos Almeida et al. (2019) we analyzed a NIR long-slit spectrum (slit PA = −16°, almost coincident with the minor axis of the CO disk) of this QSO2 obtained with the instrument EMIR on the 10.4 m Gran Telescopio CANARIAS (GTC). We detected blueshifted ionized and warm molecular gas within radial sizes of 1.3 ± 0.2 and 1.5 ± 0.2 kpc. The maximum velocity that we measured for the warm molecular outflow is −750 km s⁻¹, with a FWHM  ∼ 1500 km s⁻¹. For the ionized gas, using the Paα and Brδ lines, we measured maximum velocities of −1200 km s⁻¹ and FHWMs ∼ 1800 km s⁻¹. We also detected a blueshifted broad component in the coronal line of [Si VI]λ1.963 μm, with v_max = −850 km s⁻¹ and FWHM = 1500 km s⁻¹. This emission line can only be produced by AGN photoionization or shocks (Rodríguez-Ardila & Fonseca-Faria 2020). Thus, the NIR data clearly shows that we are witnessing a multiphase AGN-driven outflow in J1509.

The CO velocity field (moment 1 map in Fig. 13) shows a rotating distribution, blueshifted to the E and redshifted to the W. Considering this, and assuming that the spiral arms trail galaxy rotation, rotation has to be clockwise and thus, the N is the near side. Our ^3DBAROLO models of the moment 1 and 2 maps and corresponding residuals are shown in the middle and right panels of Fig. 13. The orientation of the kinematic major axis (PA = 82°) and inclination of the CO disk (i = 43°) are very similar to those measured from the optical image (see above), indicating rotation-dominated CO kinematics. The stellar bar regulates the CO distribution and kinematics within the corotation radius, which we estimate as R_CR ∼ (1.2 ± 0.2)×R_bar ∼ 3.3″. The PV diagram along the major axis displayed in Fig. 14 shows that the bulk of the gas is rotating within the inner 2″ (4 kpc). Beyond this region, we also see rotation up to distances of ∼2.7″ to the E and ∼2″ to the W. This corresponds to gas transiting from the x₁ to the x₂ orbits of the bar.

Fig. 14. Same as in Fig. 6 but for J1509. The CO kinematic major axis is PA = 262°= 82°, the minor axis PA = 352° = −8° and the inclination i = 43°. We used a slit width of ∼0.25″. The orange boxes indicate the regions and velocities used to estimate a lower limit to the mass of the outflowing gas. The gas shows the characteristic positive/negative pattern but we do not resolve the high-velocity gas within the outflow.

The gas beyond the central oblong structure is concentrated along the two leading edges of the stellar bar and it follows the canonical response. If we look at the residual velocities in Fig. 13, once we subtract our rotating disk model we see mostly blueshifted velocities to the N and redshifted to the S along the minor axis. These residuals correspond to outflowing gas in a mostly coplanar geometry, considering that the N is the near side. We also see gas falling onto the galaxy center (i.e., the redshifted blob in the center of the velocity residual map). The PV diagram along the minor axis (PA = −8°; bottom panel of Fig. 14) shows evidence for outflowing gas at low velocities (up to ±70 km s⁻¹), showing the characteristic positive/negative velocity pattern to the S and N. This low-velocity gas is clearly resolved up to r ∼ 2″ (4 kpc), but we do not resolve the high-velocity gas in spite of the good angular resolution (see the ^3DBAROLO contours in the bottom panel of Fig. 14).

Thus, the CO kinematics of this QSO2 suggest a competition between the gas that follows the canonical response to the bar and consequently falls inward, and the gas that is being pushed away by the jet, wind and/or ionized outflow in a coplanar geometry. We discard the vertical inflow scenario to explain the blueshifted/redshifted gas to the N/S (see Fig. 3) because (1) it is highly unlikely in an undisturbed disk galaxy, and (2) we do not observe any signature of line splitting. The case of J1509 is different from the Teacup because we see evidence for both a warm and cold molecular outflows, as found in the case of the Seyfert 2 galaxy IC 5063 (Tadhunter et al. 2014; Morganti et al. 2015), where the jet and the multiphase outflow are coplanar with the galaxy and CO disks (Mukherjee et al. 2018). This could also be the case for J1509, and in this case we should see the blueshifted side of the ionized outflow in the N side, being almost coplanar with the H₂/CO disk (i.e., scenario A in Fig. 4).

In Ramos Almeida et al. (2019) we estimated a warm molecular outflow mass of (1.0 ± 0.2) × 10⁴ M_⊙, assuming local thermal equilibrium. Using the conversion factor of 6 × 10⁻⁵ reported by Emonts et al. (2014) for two nearby LIRGs observed in the NIR and millimeter ranges we estimated a total molecular gas mass in the outflow of (1.7 ± 0.4) × 10⁸ M_⊙. This is consistent with the lower limit of M_out≥6.8 × 10⁷ M_⊙ that we measure by integrating the CO(2−1) emission along the minor axis and within the regions indicated in Fig. 14. These regions are the ones where the molecular gas shows higher velocities than the rotating disk model at these radii. The cold molecular outflow is therefore slower and less turbulent than its warm molecular and ionized counterparts, but dominant in terms of mass. From this lower limit of the outflow mass, we estimate an outflow rate ≥1.03 M_⊙ yr⁻¹. Considering the deprojected value, of ≥1.11 M_⊙ yr⁻¹, and the SFR derived from the IR luminosity, of 34 M_⊙ yr⁻¹, we measure a mass loading factor of η ≥ 0.03.

5.2. Continuum emission of the QSO2s

As can be seen from Sect. 5.1 and Appendices A and B, most of the continuum maps are either compact (J0232, J1010, and J1100; having deconvolved sizes of 0.2, 0.4, and 0.5 kpc, respectively) or slightly asymmetric (J1152 and J1356; 0.3 × 0.2 and 0.6 × 0.4 kpc²), with the peak indicating the AGN position. This is not the case for the Teacup (J1430) and J1509, which show jet-like morphologies with sizes and PAs of 0.6 × 0.4 kpc² and ∼80°, and 1.5 × 1.1 kpc² and ∼135°, respectively (see Table 4). These millimeter morphologies are similar to the centimeter morphologies depicted from 6 GHz VLA data at 0.25″ resolution for the four QSO2s in common with Jarvis et al. (2019). J1010 and J1100 appear compact in the VLA images, whilst J1356 and J1430 are extended. According to these authors, the extended radio structures most likely correspond to jets (see Table 4). In the case of J1356 the structures that we detect at 3σ in our ALMA continuum image have different PA than the radio axis derived from VLA data (see Sect. 5.1.3 and Table 4). For J1430, the PA from the ALMA continuum map is very similar to the PAs measured from the VLA data at 6 GHz (60°) and also from FIRST (77°; Harrison et al. 2014). All the continuum sizes, PAs, flux densities and corresponding rms values measured from our ALMA data are reported in Table 4, together with the PAs and morphologies measured from the VLA observations of the QSO2s (Jarvis et al. 2019).

For J1010, J1100, J1356, and J1430 we can calculate spectral indices using the ALMA 200–220 GHz observed fluxes and the 1.5, 5.2, and 7.2 GHz VLA fluxes reported in Jarvis et al. (2019). For the other three sources we used the 1.4 GHz FIRST integrated fluxes instead. We note that the difference in angular resolution between the FIRST and ALMA data could have an impact on the determination of the spectral index of the QSO2s with extended radio morphologies (i.e., J1509 and, to a lesser extent, J1152). We obtain spectral indices (α, with S_ν ∝ ν^α) between −0.55 and −0.90, which correspond to steep spectrum radio sources (see Table 4). Thus, even the QSO2s with compact radio morphologies in the sample (J0232, J1010 and J1100) have steep spectra, similar to compact HERGs of intermediate radio power at z < 0.1 (Pierce et al. 2020). This suggests that the observed millimeter continuum emission of the QSO2s most likely corresponds to synchrotron radiation from particles accelerated by shocks and/or small-scale jets.

5.3. Molecular gas content of the QSO2s

5.3.1. CO morphologies

We detect CO(2−1) emission in the five QSO2s included in Sect. 5.1, and the CO morphologies are very diverse. The two spiral galaxies, J1100 and J1509, are the most extended, showing CO emission at ≥3σ up to distances of 5 and 4 kpc from the AGN position, respectively. Their CO emission follows the spiral arms and bars, as is normally the case for this type of galaxy (see e.g., Fig. 1 in Bolatto et al. 2017). In the case of J1356, if we only consider the N nucleus of J1356 we measure R_CO ∼ 1.7 kpc, but we detect CO at ≥3σ in the W arm up to a distance of 5 kpc from the QSO2 nucleus (see Table 5). Finally, J1010 and J1430 show double-peaked morphologies that do not have optical counterpart (see top panels of Figs. 5 and 11), with peak separations of 1.25 and 1.3 kpc, respectively, and almost aligned with the CO kinematic major axis. In the case of J1430, the jet is perpendicular to the two peaks, whereas in the case of J1010 the continuum morphology appears rather compact and round (0.4 × 0.4 kpc²). This CO morphologies could have been produced by the action of AGN feedback in the central kpc of these QSO2s (see Sects. 5.1.1 and 5.1.4).

The integrated CO spectra are shown in Fig. 15. We extracted them from the line emitting regions at ≥3σ in the combined data. The two spiral galaxies and J1010 show double-peaked CO profiles, typical of rotating disks. J1430 shows a broad single-peaked CO line profile with a full width at half maximum (FWHM) of 485 km s⁻¹. Finally, J1356 shows an asymmetric CO profile, with a prominent red wing and a blueshifted peak. This profile corresponds to emission from the N nucleus (broad CO component) and the W arm (blueshifted peak). Since J1356 is an ongoing merger system, in Table 5 we also report the CO(2−1) flux measured for the N nucleus only, which is the one shown in Figs. 9 and 10. For the red ETGs, J0232 and J1152, we estimated 3σ upper limits for the fluxes by assuming FWHM = 430 km s⁻¹ and considering a circular aperture of R_CO = 2.9 kpc, which are the average values of the three ETGs with CO detections (J1010, J1356, and J1430). All the CO(2−1) fluxes are reported in Table 5. Corresponding errors include the uncertainty associated with the measurement of the spectra and the 10% of flux calibration error.

Fig. 15. Integrated spectra of the five QSO2s with CO(2−1) detection. J1010, J1100 and J1509 show double-peaked CO profiles characteristic of rotating disks. J1356 and J1430 show very broad CO profiles, indicating the presence of different emission-line components. For J1356 we also include the spectrum of the N nucleus only (J1356N). In all panels, v = 0 km s−1 corresponds to the CO frequencies listed in Table 5.

Finally, we measured the flux of the central kpc of the galaxies (r = 0.5 kpc; except for J1010, for which we used r = 0.72 kpc instead because of its lower angular resolution) to estimate the percentage of the total emission that it represents. We find that the central kpc of the spiral galaxies contains ∼5–12% of the total molecular gas, whereas in the interacting, merging and post-merger systems it represents between 18 and 25% (in the case of J1356 it is 32% if we consider the N nucleus only; i.e., J1356N in Table 5). Merging systems are thus more centrally concentrated than the spiral galaxies in our sample. Large nuclear molecular gas concentrations of between 38 and 75% were also reported for a sample of four PG-quasars at z ∼ 0.06 and L_bol ∼ 10⁴⁵ erg s⁻¹ observed in CO(2−1) with ALMA (Izumi et al. 2020). These concentrations correspond to the gas in the central 700 pc as compared with the inner 2 kpc of the PG-quasars. Larger AGN and control samples need to be observed with high angular resolution millimeter data to investigate any relation between the distribution of molecular gas and galaxy morphology, AGN luminosity, Eddington ratio and/or outflow properties (García-Burillo et al. 2021).

5.3.2. Integrated CO luminosities

Using the fluxes reported in Table 5 and following Eq. (3) in Solomon & Vanden Bout (2005) we calculated the CO luminosities reported in Table 5.

In Fig. 16 we show L_IR (8–1000 μm) versus for our QSO2s, with different symbols and colors indicating whether the galaxies are spirals, red ETGs or interacting, merging and post-merger systems. Quasars of roughly the same bolometric luminosity (log L_bol = 45.7–46.3 erg s⁻¹) and hosted in galaxies of similar stellar masses (log M_* = 10.9–11.3 M_⊙) have CO luminosities spanning over more than one order of magnitude. This is in contradiction with the idea that all low-redshift, optically-selected quasars reside in gas-rich host galaxies and not in ellipticals (e.g., Scoville et al. 2003; Shangguan et al. 2020b).

Fig. 16. IR (8–1000 μm) versus CO(2−1) luminosities. The QSO2s in our sample are shown with same colors and symbols as in Fig. 2. The solid line is the relation (log LIR = 1.05 log ) found for LIRGs and ULIRGs at z ≤ 0.1 and dusty star forming galaxies at z > 1 from Greve et al. (2014). The dashed lines indicate the correlation scatter, of 0.27 dex. The dot-dashed line is the fit of star-forming galaxies with z = 0–2 from Genzel et al. (2010). This relation (log LFIR = 1.15 log ) has a scatter of 0.31 dex and it corresponds to CO(1−0) instead of CO(2−1), so we assume r21 = 1 and use LIR = 1.3 × LFIR to estimate total IR luminosity as indicated in Genzel et al. (2010). Pink plus symbols are the PG quasars at z < 0.1 from Shangguan et al. (2020b), open diamonds the PUMA ULIRGs at z < 0.16 from Pereira-Santaella et al. (2021), and blue asterisks the QSO2s at z ∼ 0.1 from Jarvis et al. (2020). The APEX of J1010 is higher than the corresponding ALMA value because it includes two companion galaxies, it is lower in the case of J1100, equal for J1430, and for J1356 is an upper limit.

For comparison, in Fig. 16 we show the IR-CO(2−1) relation of Greve et al. (2014), derived from a sample of nearby LIRGs and ULIRGs at z ≤ 0.1 and submillimeter galaxies at z > 1. We also include ALMA Compact Array (ACA) measurements for PG quasars at z < 0.1 with ⟨L_bol⟩ = 44.27 ± 0.36 erg s⁻¹ from Shangguan et al. (2020b). These bolometric luminosities are closer to the values measured for Seyfert galaxies than for quasars, and indeed, they have lower IR and CO luminosities than the five QSO2s with CO(2−1) detections. The beamsize of the ACA CO(2−1) observations is 7.4″ × 4.8″ (13.7 × 8.9 kpc² at z = 0.1), and the IR luminosities were estimated from SED fitting. We also include in Fig. 16 the nine QSO2s at z ∼ 0.1 from Jarvis et al. (2020). These targets have [OIII] luminosities above 10^8.4 L_⊙, 1.4 GHz luminosities ≥10^23.5 W Hz⁻¹ and ionized outflows with FWHM ≳ 700 km s⁻¹. These criteria exclude QSO2s in red ETGs like J0232 and J1152, which show slower ionized outflows and lower radio luminosities. Four of our five radio luminous QSO2s are in Jarvis et al. (2020) and thus they have APEX CO(2−1) fluxes measured in an aperture of ∼28″ (52 kpc at z = 0.1). With ALMA we measure exactly the same flux for J1430, and a higher flux for J1100 (see Fig. 16). For J1356 an APEX upper limit is reported in Jarvis et al. (2020). Finally, for J1010 the APEX CO(2−1) flux is overestimated by a factor of 2 because the large aperture includes two companion galaxies that also emit in CO (see Appendix C and Table 5). We inspected the optical morphologies of the other five QSO2s studied in Jarvis et al. (2020) using color-combined SDSS images and all of them are also disks and/or interacting galaxies. This is why they all show large CO and IR luminosities in Fig. 16, occupying the same region as our spirals, interacting, and merging QSO2s. As we discuss in Sect. 6.1, selecting luminous QSO2s that have strong ionized outflows and high radio powers biases the samples to have large molecular gas reservoirs and high SFRs.

These gas-rich QSO2s lie close to the Greve et al. (2014) relation, which is reasonable considering that all of them are LIRGs. Indeed, most of the galaxies at z ≤ 0.1 in Greve et al. (2014) are LIRGs. For further comparison, in Fig. 16 we also include the PUMA sample of ULIRGs at z < 0.16 from Pereira-Santaella et al. (2021). These 23 ULIRGs, observed with ALMA in CO(2−1) at an angular resolution of ∼400 pc, are above the Greve relation. Thus, according to their IR and CO luminosities, QSO2s show intermediate values between those of MS galaxies (dot-dashed line in Fig. 16) and nearby ULIRGs.

5.3.3. Molecular gas masses

Assuming that the CO emission is thermalized and optically thick, the CO luminosity is independent of J and of rest frequency, and thus the brightness temperature ratio, R₂₁ = CO(2−1)/CO(1−0) = 1 (Braine & Combes 1992; Solomon & Vanden Bout 2005). Since we do not have CO(1−0) fluxes measured for our QSO2s, we assume R₂₁ = 1 instead of using average values from the literature, which usually range between 0.5 and 1.2 (see e.g., Shangguan et al. 2020b and references therein) and are missing for the particular case of QSO2s. The only exception is J1356, for which a total CO(1−0) luminosity of (1.03 ± 0.32)×10⁹ K km s⁻¹ pc² is reported in Sun et al. (2014), measured from ALMA Cycle 0 data. This value is considerably lower than the CO(2−1) luminosity reported in Table 5, probably due to the best sensitivity of our combined ALMA Cycle 6 data.

In order to estimate molecular gas masses, we also need to assume a CO-to-H₂ conversion factor (α_CO). For easier comparison with the literature, we chose the Milky Way value α_CO = 4.35 ± 1.30 M_⊙(K km s⁻¹ pc²) ⁻¹ from Bolatto et al. (2013). The same or similar values have been used to estimate the gas masses of the COLD GASS (Saintonge et al. 2012, 2016) and ATLAS^3D surveys (Young et al. 2011), the type-1 quasars in Husemann et al. (2017) and the QSO2s in Jarvis et al. (2020). We discuss results from these surveys in Sect. 6.1. Moreover, the Galactic factor is very close to the peak value of the distribution of conversion factors estimated for the xCOLD GASS high-mass sample (log M_*/M_⊙ > 10.0; see Fig. 7 in Accurso et al. 2017). Furthermore, this α_CO distribution is narrow, indicating that molecular gas scaling relations should not change substantially for massive galaxies at low-redshift. Corresponding gas masses (M_H2), depletion timescales due to star formation (t_dep = M_H2/SFR) and molecular gas fractions (f_H2 = M_H2/M_*) for the QSO2s in our sample are given in Table 5.

The two spiral galaxies have molecular gas masses of ∼2 × 10¹⁰ M_⊙ and the interacting, merging and post-merger galaxies ∼4–11 × 10⁹ M_⊙. These total gas masses are in good agreement with the values reported in the literature for small samples of QSO2s at z < 0.3, of ∼7–25 × 10⁹ M_⊙ (Krips et al. 2012; Villar-Martín et al. 2013; Jarvis et al. 2020), and larger than the gas masses measured for COLD GASS (purple and blue squares in Fig. 17). On the other hand, the red ETGs have gas masses < 4–7 × 10⁸ M_⊙, consistent with the values reported for bulge-dominated type-1 quasars in Husemann et al. (2017) and also with those of the ETGs in the ATLAS^3D survey (red diamonds in Fig. 17).

Fig. 17. SFR versus molecular gas masses for the QSO2s, shown as circles of the same colors as in Fig. 2. Average errors are indicated in the bottom right corner. Dotted, dashed and dot-dashed lines correspond to constant depletion times of 108, 109, and 1010 yr. Data from different samples are shown as follows: ATLAS3D are red diamonds, COLD GASS are blue/purple squares (disk/bulge-dominated galaxies), type-1 QSOs from Husemann et al. (2017) are triangles of the same colors as our QSO2s (orange=bulge-dominated, yellow=disk-dominated, light turquoise=major mergers), and QSO2s from Jarvis et al. (2020) are open circles.

6. Discussion

6.1. The molecular gas reservoirs of nearby QSO2s

Husemann et al. (2017) studied the CO(1−0) and CO(2−1) emission of 14 nearby type-1 AGN at z < 0.2 with log L_bol = 44–46 erg s⁻¹ using NOEMA data². Although some of the bolometric luminosities are more Seyfert- than QSO-like, in the following we refer to them as type-1 QSOs. Husemann et al. (2017) claimed that galaxy morphology has an influence on the CO content of quasars. Whilst disk-dominated and merging quasars show gas masses typical of star-forming galaxies of the same stellar mass (M_H2 ∼ 10^9 − 10 M_⊙), bulge-dominated quasars have M_H2 ≤ 10⁹ M_⊙ (see Fig. 17). This dependence of gas mass on galaxy morphology is also apparent when comparing the disk-dominated massive galaxies in the COLD GASS survey and the bulge-dominated galaxies in the COLD GASS and ATLAS^3D surveys, also shown in Fig. 17.

As can be seen from the same figure, we also find that our interacting, merging and post-merger QSO2s and the spirals have the largest gas masses (4–18 × 10⁹ M_⊙), occupying the upper right corner of Fig. 17. This is also the case of the 5 QSO2s in Jarvis et al. (2020) not included in our sample, shown as open circles. On the other hand, the red ETGs have smaller CO masses (< 6.7 × 10⁸ M_⊙). We note, however, that J1152 is a post-merger ETG, showing bright tidal features similar to these observed in J1430. The difference between both J0232 and J1152 and the other ETGs in our sample are their redder optical colors, lower SFRs (3 and 3.7 M_⊙yr⁻¹), lower radio luminosities (log L_1.4 GHz ≤ 23 W Hz⁻¹) and slower ionized outflows (see Table 4).

The red ETGs, interacting and merging QSO2s show the shortest depletion times (< 220 Myr), as also found by Husemann et al. (2017) for their sample of type-1 quasars using Hα-based SFRs (see Fig. 17). In the bulge-dominated quasars, the SFR is enhanced relative to non-active ETGs of the same stellar mass (e.g., COLD GASS bulge-dominated galaxies, shown as purple squares in Fig. 17), and the molecular gas masses are small (gas fractions f_H2< 0.01), leading to short t_dep. However, as mentioned in Sect. 2, the IR-based SFRs that we measure for the red ETGs are likely overestimated, which would translate into longer depletion timescales. In these objects we might be witnessing either later-stage or less intense AGN feedback, in which outflows are weaker, the molecular gas reservoirs smaller, and the last bursts of star formation are taking place.

For the five QSO2s with CO detections, we measure M_H2 = 4–18 × 10⁹ M_⊙, t_dep = 130–520 Myr and f_H2 = 0.04–0.20. These values are similar to those reported for the QSO2s observed with APEX by Jarvis et al. (2020) if we use R₂₁ = 1 to convert their CO(2−1) luminosities into CO(1−0) values: M_H2 = 6–25 × 10⁹ M_⊙, t_dep = 100–800 Myr and f_H2 = 0.08–1.0. These gas fractions are higher than ours, possibly due to the different methods employed to estimate the stellar masses (SED fitting versus NIR magnitudes). Indeed, for the four QSO2s that we have in common, the stellar masses reported by Jarvis et al. (2020) are lower. Potentially, our stellar masses could be contaminated with AGN emission, but they are consistent with those reported by Shangguan & Ho (2019), also derived from SED fitting (see Sect. 2). These values of M_H2, t_dep and f_H2 are representative of disk-dominated, interacting and merging QSO2s, but not of QSO2s in red ETGs such as J0232 and J1152. QSO2s with fast ionized outflows and high radio luminosities (log L_1.4 GHz > 23.5 W Hz⁻¹) have large reservoirs of molecular gas and high SFRs, as also discussed in Jarvis et al. (2020). Thus, even though AGN feedback is taking place in the form of ionized and molecular gas outflows (see Sect. 6.2), there is still plenty of molecular gas and star formation in the spirals and merging QSO2s.

The gas fractions measured for the five QSO2s with CO(2−1) detections (f_H2 = 0.04–0.20) are in between those reported in Saintonge et al. (2016) for massive galaxies (log M_*/M_⊙ = 10.6–11.2) with specific-SFRs (sSFRs)³ of log sSFR = −9.95 ± 0.01 yr⁻¹, f_H2 = 0.13 ± 0.01, and those of galaxies in the MS, f_H2 = 0.03–0.04. J1100 and J1509 (i.e., the spirals) show gas fractions of 0.15 and 0.20, whilst the merging QSO2s have f_H2 = 0.04–0.06. The spirals and J1430 have t_dep ∼ 500 Myr, and J1010 and J1356 have t_dep ∼ 130–170 Myr. These depletion times are shorter than those of massive galaxies (log M_*/M_⊙ > 10.8) in the MS (1.2–1.5 Gyr; Saintonge et al. 2016) and also of massive galaxies with −10.4 < logsSFR < −9.6 (i.e., similar to the QSO2s), of ∼1 Gyr.

The interacting and merging QSO2s show the shortest depletion times. This is also the case for the two quasars in major mergers in Husemann et al. (2017), shown as turquoise triangles in Fig. 17, and for the interacting QSO2s in Jarvis et al. (2020), as e.g. J0945+1737, J1000+1242, and J1316+1753, all of them having t_dep < 300 Myr. Even shorter depletion timescales due to SFR only, of tens of Myr, are estimated for ULIRGs (e.g., Cicone et al. 2014). As mentioned in Sect. 5.3.2, our QSO2s show intermediate properties between ULIRGs and MS galaxies in terms of molecular gas content and SFR.

Saintonge et al. (2017) compared the molecular gas content of BPT-selected active and non-active galaxies in the xCOLD GASS sample with log M_*/M_⊙ > 10 and matched in sSFR. They found slightly lower molecular gas fractions in AGN (f_H2 ∼ 0.014) than in the matched non-AGN sample (f_H2 ∼ 0.021). Similar molecular gas contents measured from single-dish radio telescopes were reported by Rosario et al. (2018) for the LLAMA sample of low-to-intermediate luminosity AGN and matched control sample galaxies, and also by Shangguan et al. (2020a) for a representative sample of 40 z < 0.3 PG quasars. Nevertheless, for the few luminous AGN in xCOLD GASS, Saintonge et al. (2017) reported higher gas fractions than those of weaker AGN and similar to sSFR-matched inactive galaxies. Recently, Koss et al. (2021) found that X-ray selected AGN in nearby massive galaxies have higher gas fractions than inactive galaxies matched in stellar mass. Moreover, they find that the more luminous the AGN and the higher its Eddington ratio, the larger the molecular gas mass. This is consistent with luminous AGN being more frequently found in more massive star-forming and merging/interacting galaxies (e.g., Shimizu et al. 2015), although here we showed that this is not the case for all luminous quasars.

These and other works searched for differences between the gas content of AGN and non-active galaxies in the hope of finding evidence for quenching induced by AGN feedback. However, most of these comparisons show that the AGN gas fractions are equal or larger than those of non-active galaxies matched in stellar mass, indicating that when the AGN is on, generally there is gas and star formation (although this depends on galaxy morphology; see for example the QSO2s in red ETGs). This does not mean that AGN feedback is not having an effect on the molecular gas reservoirs of galaxies, but considering the short AGN lifetimes, and the possibility that every massive galaxy goes through one or several active phases, spotting the action of feedback on the total gas fractions becomes challenging. The region we should focus on if we aim to evaluate the effect of the current episode of nuclear activity is the central kpc. This is the region having a dynamical timescale comparable to the periods of nuclear activity and the outflows (∼1–100 Myr; see e.g., Martini et al. 2004). Furthermore, reduced gas fractions relative to star-forming regions in the same host galaxies have been recently measured in the central kpc of AGN using CO observations (Ellison et al. 2021). In our QSO2s, we find smaller concentrations of molecular gas in the central kpc of the spiral galaxies, which are the ones showing more massive molecular outflows and higher Eddington ratios (see Sect. 6.2). In the merging QSO2s, we see tentative evidence of AGN feedback modifying the distribution of CO in the central 1–2 kpc of the galaxies (i.e., double-peaked CO morphologies in J1010 and J1430). This could be consistent with an inside-out quenching scenario (e.g., Sánchez et al. 2018).

In summary, our ALMA data of nearby QSO2s, along with the comparison with the QSO2s in Jarvis et al. (2020) and the QSO1s in Husemann et al. (2017), reveal the existence of at least two different scenarios regarding the molecular gas content. First, the disk-dominated and merging quasars have large molecular gas masses and extended CO morphologies. The gas fractions (f_H2 = 0.04–0.20) are consistent with those of massive non-active galaxies with similar sSFRs, but the depletion timescales of the QSO2s (t_dep = 100–500 Myr) are shorter. The merging QSO2s show more concentrated CO morphologies than the spirals. Second, the quasars in red ETGs show smaller molecular gas reservoirs, slower ionized outflows and lower radio luminosities. These red ETGs might be undergoing the last episodes of star formation, leading to short depletion timescales of < 200 Myr. High resolution millimeter observations of a larger sample of quasars are required to confirm/discard these trends.

6.2. Cold molecular outflows

From our analysis of the CO kinematics, we find that the bulk of the molecular gas is rotating in the five QSO2s with CO(2−1) detection. In addition to this, we detect noncircular gas motions that are on the order of the circular motions. Therefore, after close inspection of the PV diagrams and residual maps shown in Sect. 5.1, we find that four of the QSO2s show evidence of coplanar outflow motions (J1100, J1356, J1430, and J1509). The alternative interpretation are vertical inflows, a scenario that is not justified in the central regions of molecular gas disks, where the inflow motions are known to take place inside the plane (driven by angular momentum transport by bars and/or spirals). Only in the case of J1356, which is a completely distorted merging system, might it be possible to have a vertical inflow driven by, for example, a tidal tail, but in that case we would observe line splitting in the PV diagrams, something that does not happen. In the case of J1010 the PV diagram along the minor axis indicates that we are witnessing either a coplanar inflow or vertical outflow. Vertical molecular outflows are usually starburst-driven, as in the case of M82 (Walter et al. 2002; Leroy et al. 2015). J1010 has strong star formation, and it is an interacting galaxy, making it possible to have a vertical outflow. However, such a vertical outflow would produce line splitting in the corresponding PV diagrams unless spatially extended beyond the projected radius of the disk (García-Burillo et al. 2019). In the case of J1010 the inflow/outflow is barely resolved, so we cannot rule out the vertical outflow scenario.

For the coplanar molecular outflows we measure projected radii r_out ∼ 0.4–1.3 kpc, velocities v_out ∼ 115–310 km s⁻¹, and masses M_out ∼ 1.4–10.5 × 10⁷ M_⊙ (see Table 6). They are compact, we barely resolve some of them despite the high angular resolution of the data. For example, in the case of J1509 we do not resolve the high-velocity gas within the outflow, but we can estimate a lower limit of its mass from the low-velocity gas (see Fig. 14). The outflow mass rates range between 8 and 16 M_⊙ yr⁻¹, and ≥1.1 M_⊙ yr⁻¹ for J1509 (see Table 6). Considering the SFRs reported in Table 2, we measure mass loading factors η = Ṁ_out/SFR ∼ 0.1–1.3. Except in the case of J1430, the SFRs are larger than the outflow rates, indicating that these molecular outflows might not be efficient in removing molecular gas (η ≤ 1) despite the high AGN luminosity of the QSO2s (Cicone et al. 2014). Using the galaxies total molecular gas masses we estimate depletion timescales associated with the outflows () that range between 400 Myr and 1.3 Gyr (see Table 6)⁴. These are longer depletion timescales than those associated with SFR (100–500 Myr).

Taken at face value, the low mass loading factors and long depletion timescales that we find for the QSO2s could be indicating that the molecular outflows are not AGN-driven but star formation-driven. However, the fact that precisely in these four QSO2s the outflows are coplanar with the CO disks is indicative of AGN-driven outflows (García-Burillo et al. 2021). Moreover, we note that our molecular outflow masses are rather conservative, as we only consider noncircular gas motions consistent with outflowing gas along the minor axis (see Sects. 4 and 5.1). For example, in the case J1356, Sun et al. (2014) reported an outflow mass rate of 115 M_⊙ yr⁻¹ (using Eq. (1)), whereas here we are measuring 7.8 M_⊙ yr⁻¹. As indicated in Sect. 5.1.3, if we consider all the high-velocity gas (v_out ≥ 300 km s⁻¹) that we detect in CO(2−1), we measure the same outflow mass as Sun et al. (2014), of ∼7 × 10⁷ M_⊙. However, according to our analysis of the kinematics, this mass would include rotation and/or tangential noncircular motions. Apart from the larger outflow mass, Sun et al. (2014) used the highest outflow velocity (v_out = 500 km s⁻¹) and r_out = 0.3 kpc to compute the outflow mass rate.

6.2.1. QSO2 molecular outflows in a broader context

The outflow radii, velocities and mass rates of the QSO2s are intermediate between those of the cold molecular outflows reported for Seyfert galaxies and those of AGN in ULIRGs (see Table 7). In Seyfert galaxies, the outflow masses represent 0.1–1% of the total gas mass in the galaxies, very similar to the QSO2s (0.2–0.7%; see Table 7). It is noteworthy that the jetted Seyfert galaxies NGC 1068 and IC 5063 (García-Burillo et al. 2014; Morganti et al. 2015) show faster and more massive molecular outflows. The spatial coincidence between the ionized/molecular gas outflow and the radio jet in these two Seyfert galaxies indicates that the molecular gas is being pushed by the combined action of the AGN-driven wind and the jet. These are examples of strong coupling driving massive molecular outflows, which represent between 3–5% of the total gas mass in the galaxies.

The other extreme in this comparison of outflow properties are AGN in ULIRGs (e.g., Feruglio et al. 2010; Cicone et al. 2014; Fluetsch et al. 2019) and extremely powerful quasars such as PDS 456 (Bischetti et al. 2019a; Herrera-Camus et al. 2020). These objects show the most powerful AGN-driven molecular outflows in the local universe, having radii of 0.5–1.2 kpc, maximum velocities of 600–1200 km s⁻¹ and accounting for 2–27% of the total molecular gas mass (see Table 7). The AGN bolometric luminosities of these ULIRGs are slightly lower than those of our QSO2s (10^{44.7 − 45.7} erg s⁻¹), and their morphologies are consistent with major mergers, having SFRs = 6.5–140 M_⊙ yr⁻¹. Thus, as discussed in Sect. 6.1, galaxy morphology and offset from the star-formation MS are linked to the molecular gas content of galaxies, and apparently, also to the strength of the molecular outflows. The molecular outflow radii of the ULIRGs can extend to a few kpc when low surface brightness components are considered (Feruglio et al. 2013; Herrera-Camus et al. 2019; Cicone et al. 2020). Indeed, it could be possible that having observed our QSO2s in CO(2−1) prevents us from detecting a more diffuse neutral or molecular gas phase of the outflows. The neutral component might be detected in HI (n_H ∼ 100 cm⁻³), and the diffuse molecular component in CI(1−0) and CI(2−1), which have critical densities of 500 and 1000 cm⁻³ (Papadopoulos et al. 2004), or CO(1−0), with n_H = 440 cm⁻³ for a t_kin = 20 K (Yang et al. 2010).

If we plot the mass outflow rates measured for the QSO2s in a Ṁ_out versus L_bol diagram (see Fig. 18), the QSO2s are significantly below the linear fit of the CO outflow mass rates compiled by Fiore et al. (2017). According to this observational scaling relation, at the average luminosity of the QSO2s (log L_bol ∼ 46 erg s⁻¹), we would expect outflow rates of ∼500 M_⊙ yr⁻¹. Even if we use the more conservative values of L_bol of J1100, J1356, and J1430, derived from SED fitting by Jarvis et al. (2019), we would expect Ṁ_out values ranging from ∼100 to 500 M_⊙ yr⁻¹. This suggests that is not only AGN luminosity what defines how powerful molecular outflows are. Indeed, the powerful quasar PDS 456 (log L_bol = 47 erg s⁻¹ and Ṁ_out ∼ 290 M_⊙ yr⁻¹; Bischetti et al. 2019a) also has a much lower mass outflow rate than what it would be expected from extrapolating this scaling relation (∼2800 M_⊙ yr⁻¹). The Seyfert galaxies NGC 3227 and Mrk 1066, included in Table 7, also fall well below the blue dashed line. In fact, this scaling relation might be the upper boundary of the Ṁ_out versus L_bol relation, according to Fig. 18, mostly driven by the powerful outflows of the ULIRGs and the jetted-Seyferts. Other factors different from AGN luminosity, including jet power, coupling between wind/jet/ionized outflows and CO disks, and amount/geometry of dense gas in the nuclear regions might be also important for efficiently driving massive molecular outflows.

Fig. 18. Molecular mass outflow rate as a function of AGN luminosity. Circles correspond to the QSO2s. Horizontal dotted lines indicate the position that they would occupy if we use the more conservative Lbol values derived from SED fitting from Jarvis et al. (2019). Blue squares are the CO outflows of Seyferts and ULIRGs compiled by Fiore et al. (2017), but dividing the outflow rates by a factor of 3 to match our outflow geometry (see Eq. (1)). Blue dashed line is the linear fit of the blue squares. Pink diamonds are the galaxies included in Table 7 that are not in Fiore et al. (2017) except NGC 1068 and IC 5063, for which they report different AGN luminosities and outflow rates.

Focusing on the individual outflow and galaxy properties of the QSO2s, we find that the spiral galaxies, J1100 and J1509, have more massive outflows, of M_out = 10.5 and ≥6.8 × 10⁷ M_⊙, respectively (0.7% and ≥0.4% of the total gas mass), than J1356 and J1430. For them we measure M_out = 1.4 and 3.1 × 10⁷ M_⊙ (0.2% and 0.5% of the total gas mass). From the analysis of the kinematics, described in Sect. 5.1, we found that the ionized gas outflows reported for J1100 and J1509 must be almost coplanar with the galaxy disks (i.e., more favorable orientation for dragging molecular gas outward; scenario A in Fig. 4), whereas those of J1010, J1430, and J1356 would subtend a relatively large angle (scenario B in Fig. 4). This could explain the more massive molecular outflows measured for the spiral galaxies, whose ionized outflows, in combination with the jets and winds, might be contributing to drag a larger amount of molecular gas outward in an almost coplanar geometry. We have to bear in mind, however, that the merging QSO2s show more chaotic kinematics than the spirals, and thus, their molecular outflows could have a more 3D geometry than those of the spirals. In that case, the outflow masses of J1356 and J1430 would be underestimated from the methodology used here. A more detailed study of the molecular gas kinematics of these QSO2s will be the subject of a forthcoming work (Audibert et al., in prep.).

The more massive outflows measured for the spirals, with r_out ∼ 1.3 and ≤3 kpc, could be responsible for the smaller molecular gas concentrations that we measure in the central kpc (r = 0.5 kpc) of the galaxies (5–12%), in comparison with the merging systems (18–25%). The radii and dynamical timescales of the QSO2 outflows (r_out/v_out; see Table 6), between 1.4 and 11 Myr (and ≤66 Myr for J1509), are consistent with them being driven by the current AGN episode and having an effect on the central kpc of the galaxies. On smaller spatial scales (the inner ∼200 pc), a radial redistribution of the molecular gas induced by AGN feedback has been reported by García-Burillo et al. (2021) from a comparison between low and intermediate-luminosity Seyferts (NUGA and GATOS samples). The most luminous Seyferts, which also have the highest Eddington ratios and show evidence for molecular outflows, have smaller molecular gas concentration in the inner 50 pc of the galaxies relative to the inner 200 pc. On this respect, it is noteworthy that the Eddington ratios of the two QSO2s in spiral galaxies are higher (f_Edd = 1.1 and 0.6; see Table 1) than for the merging systems (f_Edd = 0.15, 0.09, and 0.4). Larger samples of AGN are required to continue exploring these tantalizing trends.

7. Conclusions

We have explored for the first time the cold molecular gas content traced by the 2–1 line of CO, and adjacent continuum emission, of a sample of nearby (z ∼ 0.1) and luminous L_[OIII] > 10^8.5 L_⊙ QSO2s at an angular resolution of ∼0.2″ (370 pc). The ALMA observations permit us to study the molecular gas content, morphology and kinematics of these QSO2s. The main results are summarized as follows.

1. Quasars of roughly the same bolometric luminosity (log L_bol = 45.7–46.3 erg s⁻¹) and hosted in galaxies of similar stellar masses (log M_* = 10.9–11.3 M_⊙) have molecular gas masses ranging from 4–18 × 10⁹ M_⊙ to < 4–7 × 10⁸ M_⊙. This is in contradiction with the idea that all low-redshift, optically-selected quasars reside in gas-rich host galaxies and not in ellipticals.

2. Galaxy morphology and color, radio luminosity, and outflow properties are related with the total molecular gas content of quasars. QSO2s in disks and/or merging systems with high radio luminosities and fast ionized outflows have larger gas fractions than red ETGs with lower radio luminosities and slower ionized outflows.

3. QSO2s show intermediate properties between those of local ULIRGs and MS galaxies regarding molecular gas content and SFRs. They have depletion times associated with star formation ranging from ∼100 to 500 Myr.

4. The CO(2−1) morphologies of QSO2s in the merging systems are more centrally concentrated than those in the spiral galaxies. The central kpc (r = 0.5 kpc) of the galaxies contains ∼5–12% of the total molecular gas in the case of the spirals, and between 18 and 25% in the interacting, merging and post-merger systems.

5. We find more massive and extended molecular outflows in the spiral galaxies than in the merging systems. The spirals have ionized outflows that are mostly coplanar with the CO disks/galaxies, whereas in the case of the merging systems the ionized outflows would subtend a larger angle relative to the CO disks/galaxies. Furthermore, the QSO2s in spirals have larger Eddington ratios than the merging QSO2s. This could explain the smaller molecular gas concentrations in the central kpc of the spirals.

6. The outflow mass rates that we measure (8–16 M_⊙ yr⁻¹) are much lower than expected from observational scaling relations with AGN luminosity (e.g., Fiore et al. 2017). This implies that is not only the AGN luminosity that defines how powerful molecular outflows are. Other factors including jet power, coupling between wind, jet, and/or ionized outflows and the CO disks, and amount/geometry of dense gas in the nuclear regions might be also relevant.

7. We use the total molecular gas masses and deprojected outflow mass rates to estimate depletion timescales associated with the outflows. We measure ranging between 400 Myr and 1.3 Gyr. These are larger depletion timescales than those associated with SFR, leading to low mass loading factors of η ∼ 0.1–1.3. Except in the case of the Teacup, star formation is more effectively consuming molecular gas than the outflows.

8. Taken at face value, the low mass loading factors and long depletion timescales that we find for the QSO2s could be indicating that the molecular outflows are not AGN-driven but star formation-driven. However, the fact that in four of the five QSO2s these outflows are coplanar with the CO disks is indicative of AGN-driven outflows.

We do not find evidence for a significant impact of the QSO2s’ molecular outflows, or more generally, of quasar feedback, on the total molecular gas reservoirs and star formation rates. However, they appear to be modifying the distribution of cold molecular gas in the central kiloparsec of the galaxies. Similar high-resolution mm observations as those presented here of a larger sample of QSO2s are required to confirm this.

Acknowledgments

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.00870.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. CRA and MB thank Anna López (IRAM) and Luke Maud (ESO) for their support with the data processing. CRA thanks R. Morganti, J. Falcón-Barroso, and J. A. Acosta-Pulido for useful discussions. CRA also thanks B. Husemann for kindly providing the data to produce Fig. 17 and G. Pérez Díaz (SMM IAC) for designing Figure 4. CRA acknowledges financial support from the European Union’s Horizon 2020 research and innovation programme under Marie Skłodowska-Curie grant agreement No 860744 (BiD4BESt). CRA and PSB acknowledge support from the project “Feeding and feedback in active galaxies”, with reference PID2019-106027GB-C42, funded by MICINN-AEI/10.13039/501100011033. CRA and AA acknowledge the projects “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, and from the Consejería de Economía, Conocimiento y Empleo del Gobierno de Canarias and the European Regional Development Fund (ERDF) under grant “Quasar feedback and molecular gas reservoirs”, with reference ProID2020010105, ACCISI/FEDER, UE. MB and CF acknowledge support from PRIN MIUR project “Black Hole winds and the Baryon Life Cycle of Galaxies: the stone-guest at the galaxy evolution supper”, contract 2017PH3WAT. AAH and SGB acknowledge financial support from grant PGC2018-094671-B-I00 (MCIU/AEI/FEDER,UE). SGB acknowledges support from the research project PID2019-106027GA-C44 of the Spanish MICINN. AAH work was done under project No. MDM-2017-0737 Unidad de Excelencia “María de Maeztu”– Centro de Astrobiología (INTA-CSIC). MPS acknowledges support from the Comunidad de Madrid through the Atracción de Talento Investigador Grant 2018-T1/TIC-11035 and PID2019-105423GA-I00 (MCIU/AEI/FEDER,UE). JP acknowledges support from STFC [ST/V000624/1]. We finally thank the anonymous referee for her/his constructive report.

References

Appendix A: SDSS J023224.24-081140.2 (J0232)

This QSO2 is hosted in the most compact galaxy of our sample, with a major axis of 6″ (11 kpc) measured from the r-band SDSS image (see Table 2). Its optical morphology resembles an undisturbed elliptical or lenticular galaxy (ETG; see Figure 1). From the SDSS magnitudes we measure M_u-M_g = 1.45, which is closer to the typical colors of red-sequence galaxies. Together with J1152 it is one of the least powerful radio and FIR emitters in the sample (see Fig. 2). The SFR that we estimate from the extrapolated FIR luminosity (see Table 2 for details), of 1.7 M_⊙ yr⁻¹, places it in the MS of local SDSS DR7 star-forming galaxies (see Sect. 2). Regarding the ionized gas kinematics, the [OIII] lines in the SDSS spectrum show a blue wing that can be reproduced with a Gaussian of FWHM∼770 km s⁻¹, blueshifted by -370 km s⁻¹ relative to the narrow component (Villar-Martín et al. 2014).

Our ALMA observations reveal compact 1.3 mm continuum emission at ≥3σ, with a deconvolved size of 0.13″ × 0.09″ (240 × 166 pc²) and PA = 145 ± 50° (see Fig. A.1 and Table 4). Using our continuum flux and the FIRST 1.4 GHz flux we measure a spectral index α = −0.67, which corresponds to a steep spectrum. We do not detect CO(2−1) emission above 3σ, but from this non-detection we can estimate an upper limit on the gas mass by assuming R₂₁ = 1 and α_CO = 4.35 M_⊙(K km s⁻¹ pc²) ⁻¹, which yields M_H2< 6.7×10⁸ M_⊙ (see Sect. 5.3 for details).

Fig. A.1. 1.48 mm (203 GHz) continuum map of J0232 with contours at 2, 3, 4, 5, and 6σ overlaid in black (σ = 0.013 mJy beam−1). The beam size (0.18″ × 0.15″) is shown in the bottom left corner. North is up and east to the left.

Appendix B: SDSS J115245.66+101623.8 (J1152)

The host galaxy of this QSO2 is classified as an elliptical galaxy in the Galaxy Zoo project (Willett et al. 2013), based on the SDSS images (see Fig. 1). Fischer et al. (2018) reported no signs of galaxy interactions from the analysis of the HST/ACS continuum image. However, deeper optical imaging observations obtained in January 2020 with the Wide Field Camera on the 2.5 m Isaac Newton Telescope (La Palma, Spain), reveal spectacular shells and a tidal tail of several kpc (Pierce et al. 2022). These features are the result of a past merger or interaction with another galaxy, and they coincide with the faint structures toward the NE (shell) and SW (fan and tail) shown in the SDSS image (see Fig. 1). Finally, from inspection of the HST/ACS image of J1152 we detect a dust lane crossing the galaxy W of the nucleus. The HST/ACS [OIII] image reveals a spectacular biconical morphology (PA = 10°) with bubbles and ripples showing a remarkable resemblance to the Teacup galaxy (Keel et al. 2012). HST/STIS spectroscopy with a slit PA = 10° reveals mainly rotation-dominated kinematics (positive velocities to the N and negative to the S) with a high-amplitude rotation curve (peak velocities ± 200 km s⁻¹). noncircular motions are restricted to a radius of 130 pc NE of the nucleus, with FWHM = 360 km s⁻¹ and velocities of up to -300 km s⁻¹ (Fischer et al. 2018). This is the slowest ionized outflow in our sample (see Table 3).

This QSO2 was observed in band 6 with ALMA because of its redshift (z = 0.07). We detect continuum emission at > 3σ, although weaker than in the other QSO2s (except J0232). The continuum in the case of J1152 is centered at 1.2 mm (246.7 GHz) rest-frame because of the spectral setup chosen. At ≥3σ the continuum appears compact, although there is a fainter elongated structure toward the NW detected at 2.5σ (see Fig. B.1), but similar structures at the same level are detected across the whole field-of-view. We measure a deconvolved size of 0.20″ × 0.12″ (267 × 160 pc²) with PA = 142 ± 80°. In the cm regime, J1152 has the lowest 1.4 GHz luminosity in our sample, but it still shows a radio excess in the radio-FIR diagram shown in Fig. 2. Using our ALMA continuum flux and the FIRST 1.4 GHz flux we measure a spectral index α = −0.65, very similar to J0232.

Fig. B.1. 1.35 mm (222 GHz) continuum map of J1152 with contours at 2, 2.5, 3, and 3.5σ in black (σ = 0.023 mJy/beam). The beam size (0.23″ × 0.15″) is shown at the bottom left corner. North is up and east to the left.

The host galaxy appears depleted of CO(2−1), which is somehow surprising considering the bright morphological features detected in the optical image, which are the product of a past gas-rich interaction with another galaxy (see Ramos Almeida et al. 2011 and references therein). As we did for J0232, from this non-detection we can estimate an upper limit to the gas mass of M_H2< 3.7×10⁸ M_⊙, by assuming R₂₁ = 1 and α_CO = 4.35 (see Sect. 5.3 for details).

Appendix C: SDSS J101043.36+061201.4 (J1010)

The QSO2 host galaxy is in clear interaction with a small companion at ∼7″ (13 kpc) SW (see Fig. 1). A long tidal tail of several kpc in size linking them is also evident. As can be seen from Fig. C.1, the small companion and another faint source at ∼9.5″ (17 kpc) SW are detected in CO(2−1). The faintest companion lacks an optical counterpart in either the SDSS or PanSTARRS images. None of the CO-emitting companions show continuum emission at 1.3 mm. Figure C.2 corresponds to the R-I color map constructed from the VLT/MUSE datacube (see Table 3), with the CO(2−1) contours overlaid. The optical images have been shifted to match the ALMA astrometry. Redder colors are seen SW from the nucleus, produced by the dust lane detected in the two optical images.

Fig. C.1. Large-scale CO(2−1) map of J1010 and two fainter companion sources at the same redshift. Contours correspond to 3, 5, 10, and 20σ, with σ = 0.09 Jy/beam. The beam size (0.80″ × 0.69″) is shown at the top left corner of the figure. North is up and east to the left.

Fig. C.2. CO(2−1) contours (0.2, 0.35, 0.5, 0.65, 0.8, 0.95 times the peak of CO emission) of J1010 overlaid on the R-I color map, shown in color scale, derived from continuum images obtained from VLT/MUSE data (see Table 3), and centered at 6200 and 8500 Å. The star corresponds to the peak of the ALMA continuum image and the color map corresponds to the central 6″ × 6″. North is up and east to the left.

Appendix D: SDSS J135646.10+102609.0 (J1356)

J1356 is the only QSO2 in our sample with published spatially resolved CO and adjacent continuum observations to date. Sun et al. (2014) observed the CO(1−0) and CO(3−2) emission lines with ALMA in Cycle 0 (beam size 1.9″ × 1.3″) and Cycle 1 (0.35″ × 0.29″), respectively. The peak of the 3 mm continuum as measured from the Cycle 0 data coincides with the position of the QSO2 (the N nucleus). This N nucleus is detected in CO(1−0) and CO(3−2), as well as the “western arm” (W arm) and the S nucleus.

When we integrate over the broad CO(2−1) line (see Fig. 15) we only detect the N nucleus and the W arm. In Fig. 9 the latter appears as blueshifted emission at ∼1″ westward of the N nucleus, extending to a maximum distance of 2.3″ (5 kpc) from the N nucleus (see left panel of Fig. D.1). It coincides with the location of the base of an extended stellar feature in the HST continuum images (Sun et al. 2014). If we integrate over the narrower and blueshifted CO emission detected in this region we see three blobs, including the N nucleus, along PA = −90°, all of them surrounded by diffuse emission (see left panel of Fig. D.1). In Table D.1 we report the integrated flux, FWHM, and velocity of CO(2−1) extracted from the 4σ contours at the positions of the two western blobs shown in Fig. D.1. The molecular gas at these positions is blueshifted by 200 km s⁻¹, with a FWHM of 150–260 km s⁻¹, in agreement with the CO(1−0) measurements reported by Sun et al. (2014). These authors suggested that the W arm could be a tidal feature resulting from the merger, slightly offset from its stellar counterpart. Offsets between the gas and the stars have been observed in tidal features (e.g., Hibbard et al. 2000).

Fig. D.1. CO(2−1) moment 0 maps of J1356 integrating over the emission detected in the “W arm,” westward of the N nucleus (left panel; black contours go from 2 to 7σ with σ = 0.0196 mJy beam−1), and in the S nucleus, ∼1.1″ south of the N nucleus (right panel; white contours go from 2 to 5σ with σ = 0.0156 mJy beam−1). The N and S nuclei are connected by a CO-emitting bridge or spiral arm. The beamsize is indicated in the bottom left corner of each panel.

We also detect CO(2−1) emission at ≥3σ at the position of the S nucleus when we integrate over the narrower and slightly blueshifted CO profile there detected at ≥3σ (see right panel of Fig. D.1). Our data also reveal a bridge/arm connecting the S and N nuclei. Corresponding CO line measurements for the S nucleus and S arm are reported in Table D.1. We measure a blueshift of -45 km s⁻¹, in good agreement with the blueshift estimated from the 2σ CO(1−0) detection at the S nucleus (Sun et al. 2014). No CO(3–2) emission was found at this position.

Sun et al. (2014) measured a total molecular gas mass in the merging system of 4.9 × 10⁹ M_⊙ if we assume an α_CO = 4.35 M_⊙(K km s⁻¹ pc²) ⁻¹. This mass is smaller than our value of 11.5×10⁹ M_⊙, but consistent within the large uncertainties. According to Sun et al. (2014), one third of this molecular gas would be distributed in a compact rotating disk at the position of the QSO2, with a radius of 300 pc, and half of the molecular gas in the W arm. The S nucleus would contain ∼10% of the molecular gas mass in the system. From our CO(2−1) data, we find that the N nucleus represents ∼55% of the total mass (6.4×10⁹ M_⊙), the W arm ∼32% (after including the diffuse emission at ≥3σ), and the S nucleus and S arm ∼13%.

Figure D.2 shows the B-I color map constructed from the HST/WFC3 F438W and F814W images of J1356 described in Comerford et al. (2015), with the CO(2−1) contours overlaid. The HST images have been shifted to match the ALMA astrometry. A dust lane crosses the N nucleus, and slightly redder colors are observed toward the N of the CO disk.

Fig. D.2. Same as in Figure C.2 but with the CO(2−1) contours of J1356 overlaid on the B-I color map derived from HST/WFC3 F438W and F814W images. The color map corresponds to the central 4″ × 4″.

Appendix E: SDSS J143029.88+133912.0 (J1430)

In Figure E.1 we show the r-i color map constructed from the HST/WFC3 F621M and F763M images of J1430 described in Keel et al. (2015), with the CO(2−1) contours overlaid. The HST images have been shifted to match the ALMA astrometry. An intricate system of shells can be clearly seen in the two images westward of the nucleus, where red colors are observed.

Fig. E.1. Same as in Figure C.2 but with the CO(2−1) contours of J1430 (0.15, 0.2, 0.35, 0.5, 0.65, 0.8, and 0.95 times the peak of CO emission), overlaid on the r-i color map derived from HST/WFC3 F621M and F763M images. The color map corresponds to the central 4″ × 4″.

All Tables

All Figures

Fig. 1. SDSS gri color composite images of the seven QSO2s. Their optical morphologies include barred spiral galaxies (J1100 and J1509), red ETGs (J0232 and J1152) and interacting, merging and post-merger galaxies (J1010, J1356, and J1430). North is up and east to the left. The green horizontal bars at the top left of each panel correspond to 5″, which, at the average redshift of the targets (z = 0.1), corresponds to ∼9 kpc. The images are 40″ × 40″ (74 × 74 kpc2).

In the text

	Fig. 2. Radio (1.4 GHz) versus IR luminosities (8–1000 μm) for the QSO2s in our sample. The radio-FIR correlation of star-forming galaxies from Bell (2003) is shown as a solid line (slope of 1.05 ± 0.04), with a correlation scatter of 0.26 dex. All the QSO2s in our sample show a radio excess. Different colors indicate whether the QSO2s are red ETGs, spirals, or interacting galaxies (see Table 2).
In the text

	Fig. 3. Flowchart of the methodology used here for interpreting the molecular gas kinematics of the QSO2s. We note that the coplanar outflow/vertical inflow and coplanar inflow/vertical outflow degeneracy only holds if the vertical motions are spatially extended. Otherwise, line splitting would be detected.
In the text

Fig. 4. AGN wind/jet/ionized outflow geometries considered in this work and the coplanar molecular outflows that they produce. Left panel: example of strong coupling (scenario A), in which the wind/jet/ionized outflow are coplanar with the disk, launching a more massive molecular outflow. Right panel: example of weak coupling (scenario B), in which the wind/jet/ionized outflow subtends a certain angle from the CO disk, having a small intersection region (shown in dark gray) and launching a more modest molecular outflow.

In the text

Fig. 5. CO(2−1) moment 0, 1, and 2 maps (i.e., flux, velocity, and velocity dispersion) of J1010. The 3DBAROLO models of the moment 1 and 2 maps and corresponding residuals are also shown. Continuum contours starting at 3σ are shown in pink in the moment 0 map (σ = 0.037 mJy beam−1), and the beam size (0.80″ × 0.69″) is shown in the bottom left corner of each panel. The major axis of the CO kinematics is shown as a black dashed line in the moment 1 maps. North is up and east to the left.

In the text

Fig. 6. PV diagrams along the CO kinematic major axis (PA = 290° = −70°) and minor axis (PA = 200° = 20°) of J1010 extracted with 3DBAROLO. We used a slit width of ∼0.8″. Blue contours correspond to the observed CO(2−1) emission above 2σ and red contours to the 3DBAROLO rotating disk model (inclination of i = 36°) shown in Fig. 5. Yellow dots are the average model velocities at different radii. Blue contours outside the boundaries of the model correspond to noncircular motions. Orange lines indicate the positive/negative pattern that, in the case of this QSO2, might be produced by either a coplanar inflow or a spatially extended vertical outflow.

In the text

	Fig. 7. Same as in Fig. 5 but for J1100. Continuum contours starting at 3σ are shown in pink in the moment 0 map (σ = 0.017 mJy beam−1), and beam size is 0.24″ × 0.20″.
In the text

	Fig. 8. Same as in Fig. 6 but for J1100. The CO kinematic major axis is PA = 69°, the minor axis PA ∼ 159° = −21° and the inclination i = 38°. We used a slit width of ∼0.2″.
In the text

	Fig. 9. Same as in Fig. 5 but for J1356 (N nucleus). Part of the W arm is also shown and labeled in orange in the moment 0 map. Continuum contours starting at 3σ are shown in pink in the moment 0 map (σ = 0.015 mJy beam−1), and beam size is 0.25″ × 0.24″.
In the text

	Fig. 10. Same as in Fig. 6 but for J1356. The CO kinematic major axis is PA = 110°, the minor axis PA = 200° = 20° and the inclination angle i = 52°. We used a slit width of ∼0.25″. The orange box indicates the region and velocities used to estimate the mass of the outflowing gas, which in the case of this QSO2 does not show the characteristic positive/negative pattern.
In the text

	Fig. 11. Same as in Fig. 5 but for J1430. Continuum contours starting at 3σ are shown in pink in the moment 0 map (σ = 0.015 mJy beam−1), and beam size is 0.21″ × 0.18″.
In the text

	Fig. 12. Same as in Fig. 6 but for J1430. The CO kinematic major axis is PA = 4°, the minor axis PA = 94° and the inclination angle i = 38°. We used a slit width of ∼0.2″. The orange boxes indicate the regions used to estimate the mass of the outflowing gas.
In the text

	Fig. 13. Same as in Fig. 5 but for J1509. Continuum contours starting at 3σ are shown in pink in the moment 0 map (σ = 0.013 mJy beam−1), and beam size is 0.26″ × 0.24″.
In the text

Fig. 14. Same as in Fig. 6 but for J1509. The CO kinematic major axis is PA = 262°= 82°, the minor axis PA = 352° = −8° and the inclination i = 43°. We used a slit width of ∼0.25″. The orange boxes indicate the regions and velocities used to estimate a lower limit to the mass of the outflowing gas. The gas shows the characteristic positive/negative pattern but we do not resolve the high-velocity gas within the outflow.

In the text

Fig. 15. Integrated spectra of the five QSO2s with CO(2−1) detection. J1010, J1100 and J1509 show double-peaked CO profiles characteristic of rotating disks. J1356 and J1430 show very broad CO profiles, indicating the presence of different emission-line components. For J1356 we also include the spectrum of the N nucleus only (J1356N). In all panels, v = 0 km s−1 corresponds to the CO frequencies listed in Table 5.

In the text

Fig. 16. IR (8–1000 μm) versus CO(2−1) luminosities. The QSO2s in our sample are shown with same colors and symbols as in Fig. 2. The solid line is the relation (log LIR = 1.05 log ) found for LIRGs and ULIRGs at z ≤ 0.1 and dusty star forming galaxies at z > 1 from Greve et al. (2014). The dashed lines indicate the correlation scatter, of 0.27 dex. The dot-dashed line is the fit of star-forming galaxies with z = 0–2 from Genzel et al. (2010). This relation (log LFIR = 1.15 log ) has a scatter of 0.31 dex and it corresponds to CO(1−0) instead of CO(2−1), so we assume r21 = 1 and use LIR = 1.3 × LFIR to estimate total IR luminosity as indicated in Genzel et al. (2010). Pink plus symbols are the PG quasars at z < 0.1 from Shangguan et al. (2020b), open diamonds the PUMA ULIRGs at z < 0.16 from Pereira-Santaella et al. (2021), and blue asterisks the QSO2s at z ∼ 0.1 from Jarvis et al. (2020). The APEX of J1010 is higher than the corresponding ALMA value because it includes two companion galaxies, it is lower in the case of J1100, equal for J1430, and for J1356 is an upper limit.

In the text

Fig. 17. SFR versus molecular gas masses for the QSO2s, shown as circles of the same colors as in Fig. 2. Average errors are indicated in the bottom right corner. Dotted, dashed and dot-dashed lines correspond to constant depletion times of 108, 109, and 1010 yr. Data from different samples are shown as follows: ATLAS3D are red diamonds, COLD GASS are blue/purple squares (disk/bulge-dominated galaxies), type-1 QSOs from Husemann et al. (2017) are triangles of the same colors as our QSO2s (orange=bulge-dominated, yellow=disk-dominated, light turquoise=major mergers), and QSO2s from Jarvis et al. (2020) are open circles.

In the text

Fig. 18. Molecular mass outflow rate as a function of AGN luminosity. Circles correspond to the QSO2s. Horizontal dotted lines indicate the position that they would occupy if we use the more conservative Lbol values derived from SED fitting from Jarvis et al. (2019). Blue squares are the CO outflows of Seyferts and ULIRGs compiled by Fiore et al. (2017), but dividing the outflow rates by a factor of 3 to match our outflow geometry (see Eq. (1)). Blue dashed line is the linear fit of the blue squares. Pink diamonds are the galaxies included in Table 7 that are not in Fiore et al. (2017) except NGC 1068 and IC 5063, for which they report different AGN luminosities and outflow rates.

In the text

	Fig. A.1. 1.48 mm (203 GHz) continuum map of J0232 with contours at 2, 3, 4, 5, and 6σ overlaid in black (σ = 0.013 mJy beam−1). The beam size (0.18″ × 0.15″) is shown in the bottom left corner. North is up and east to the left.
In the text

	Fig. B.1. 1.35 mm (222 GHz) continuum map of J1152 with contours at 2, 2.5, 3, and 3.5σ in black (σ = 0.023 mJy/beam). The beam size (0.23″ × 0.15″) is shown at the bottom left corner. North is up and east to the left.
In the text

	Fig. C.1. Large-scale CO(2−1) map of J1010 and two fainter companion sources at the same redshift. Contours correspond to 3, 5, 10, and 20σ, with σ = 0.09 Jy/beam. The beam size (0.80″ × 0.69″) is shown at the top left corner of the figure. North is up and east to the left.
In the text

Fig. C.2. CO(2−1) contours (0.2, 0.35, 0.5, 0.65, 0.8, 0.95 times the peak of CO emission) of J1010 overlaid on the R-I color map, shown in color scale, derived from continuum images obtained from VLT/MUSE data (see Table 3), and centered at 6200 and 8500 Å. The star corresponds to the peak of the ALMA continuum image and the color map corresponds to the central 6″ × 6″. North is up and east to the left.

In the text

Fig. D.1. CO(2−1) moment 0 maps of J1356 integrating over the emission detected in the “W arm,” westward of the N nucleus (left panel; black contours go from 2 to 7σ with σ = 0.0196 mJy beam−1), and in the S nucleus, ∼1.1″ south of the N nucleus (right panel; white contours go from 2 to 5σ with σ = 0.0156 mJy beam−1). The N and S nuclei are connected by a CO-emitting bridge or spiral arm. The beamsize is indicated in the bottom left corner of each panel.

In the text

	Fig. D.2. Same as in Figure C.2 but with the CO(2−1) contours of J1356 overlaid on the B-I color map derived from HST/WFC3 F438W and F814W images. The color map corresponds to the central 4″ × 4″.
In the text

	Fig. E.1. Same as in Figure C.2 but with the CO(2−1) contours of J1430 (0.15, 0.2, 0.35, 0.5, 0.65, 0.8, and 0.95 times the peak of CO emission), overlaid on the r-i color map derived from HST/WFC3 F621M and F763M images. The color map corresponds to the central 4″ × 4″.
In the text