Quasar feedback survey: molecular gas affected by central outflows and by ∼10-kpc radio lobes reveal dual feedback effects in ‘radio quiet’ quasars Open Access

ABSTRACT

We present a study of molecular gas, traced via CO (3–2) from Atacama Large Millimeter/submillimeter Array data, of four z < 0.2, ‘radio quiet’, type 2 quasars (L_bol ∼ 10^45.3–10^46.2 erg s⁻¹; L|$_{\mathrm{1.4\, GHz}}\sim 10^{23.7}\!-\!10^{24.3}$| W Hz⁻¹). Targets were selected to have extended radio lobes (≥ 10 kpc), and compact, moderate-power jets (1–10 kpc; P_jet ∼ 10^43.2–10^43.7 erg s⁻¹). All targets show evidence of central molecular outflows, or injected turbulence, within the gas discs (traced via high-velocity wing components in CO emission-line profiles). The inferred velocities (V_out = 250–440 km s⁻¹) and spatial scales (0.6–1.6 kpc), are consistent with those of other samples of luminous low-redshift active galactic nuclei. In two targets, we observe extended molecular gas structures beyond the central discs, containing 9–53  per cent of the total molecular gas mass. These structures tend to be elongated, extending from the core, and wrap-around (or along) the radio lobes. Their properties are similar to the molecular gas filaments observed around radio lobes of, mostly ‘radio loud’, brightest cluster galaxies. They have the following: projected distances of 5–13 kpc; bulk velocities of 100–340 km s⁻¹; velocity dispersion of 30–130 km s⁻¹; inferred mass outflow rates of 4–20 M_⊙ yr⁻¹; and estimated kinetic powers of 10^40.3–10^41.7 erg s⁻¹. Our observations are consistent with simulations that suggest moderate-power jets can have a direct (but modest) impact on molecular gas on small scales, through direct jet–cloud interactions. Then, on larger scales, jet-cocoons can push gas aside. Both processes could contribute to the long-term regulation of star formation.

1 INTRODUCTION

Active galactic nuclei (AGN) are observed as sites of growing black holes (Kormendy & Ho 2013) and are capable of converting the energy from accreted material into intense episodes of emitted energy in the form of radiation, accretion disc winds, and jets of relativistic particles. This energy can be extremely high, also exceeding the binding energy of the galaxy itself (Cattaneo & Best 2009; Bower, Benson & Crain 2012) and is theoretically capable of affecting the host galaxy through regulation of star-formation (Fabian 2012; McNamara & Nulsen 2012). Depending on how the available energy couples to the interstellar medium (ISM), the gas could be driven due to wide-angled accretion disc winds, radiation pressure on dust, and/or due to the acceleration by radio jets (e.g. Sijacki et al. 2007; Fabian 2012; King & Pounds 2015; Ishibashi & Fabian 2016; Mukherjee et al. 2016; Costa et al. 2018; Costa, Pakmor & Springel 2020; Tanner & Weaver 2022; Almeida, Nemmen & Riffel 2023). These processes can also influence the fuel available for feeding the black hole itself, thereby giving this process a self-limiting nature, and thus earning the name ‘AGN feedback’.

Direct evidence of the influence of AGN on the ISM comes from observations that have confirmed the presence of galactic-scale AGN outflows over different phases, including ionized, neutral, and molecular forms (e.g. Morganti, Tadhunter & Oosterloo 2005; Nesvadba et al. 2008; Alexander et al. 2010; Feruglio et al. 2010; Harrison et al. 2012; Liu et al. 2013; Rupke & Veilleux 2013; Cicone et al. 2014; Villar Martín et al. 2014; King & Pounds 2015; Fiore et al. 2017; Rupke, Gültekin & Veilleux 2017; Cicone et al. 2018; Harrison et al. 2018; Förster Schreiber et al. 2019; Davies et al. 2020; Roy et al. 2021; Venturi et al. 2021; Girdhar et al. 2022; Kakkad et al. 2022; Ramos Almeida et al. 2022; Kakkad et al. 2023). While each of these phases are crucial in forming a complete understanding of galaxy evolution, comprehending the impact of AGN on molecular gas is particularly popular because (i) molecular gas is the main reservoir for fuelling star-formation and the growth of supermassive black holes; (ii) most of the mass in galactic outflows is seen to reside in the molecular gas phase (e.g. Fiore et al. 2017 compiled literature measurements and found that for AGN with L|$_{\mathrm{bol}}\sim \, 10^{45}\!-\!10^{46}$| erg s⁻¹, the observed molecular outflows typically have 100 times more mass than the ionized outflows). It is hence important to understand the effect of powerful quasars on the molecular gas in their host galaxy (e.g. Feruglio et al. 2010; Alatalo et al. 2011; Mainieri et al. 2011; Cicone et al. 2014; Morganti et al. 2015; Fiore et al. 2017; Harrison 2017; Mainieri et al. 2021; Ward et al. 2022; Dall’Agnol de Oliveira et al.2023).

For the brightest AGN, with high accretion rates, the dominant feedback mechanism is typically expected to be due to accretion disc winds (which can propagate into the host galaxies) or directly due to radiation pressure. This can lead to the disturbance or removal of inter-stellar gas (e.g. Costa et al. 2018; Costa, Pakmor & Springel 2020). Many of the observational studies focusing on high accretion rate AGN have looked at starburst and highly luminous quasar targets. Specifically, there is a class of observational work searching for underlying wing components.¹ in CO emission-line profiles, as a tracer of molecular gas outflows, and then investigating these outflow properties as a function of star formation rates, stellar masses, and AGN luminosities (Cicone et al. 2014; Fiore et al. 2017; Fluetsch et al. 2019). Another class of observational studies have focused on massive, radio-luminous brightest cluster galaxies (BCGs), located in cool-core clusters, revealing molecular gas entrained in filamentary structures along with the large radio lobes and X-ray cavities in the systems (Salomé & Combes 2004; David et al. 2014; McNamara et al. 2014; Tremblay et al. 2016; Vantyghem et al. 2016; Russell et al. 2017, 2018; Tremblay et al. 2018; Olivares et al. 2019; Russell et al. 2019; Tamhane et al. 2022). Therefore, these classes of studies appear to investigate different types of feedback effects on the molecular gas, with the former assuming a dominant role of AGN winds/radiation (at least for driving the most powerful molecular outflows) and the latter finding a dominant role of radio jets.

One might conclude a simple overall picture of two feedback modes on the molecular ISM; one acting on larger scales, beyond the gas disc, and caused by powerful radio jets (e.g. in the BCGs) and one acting within the molecular gas discs, due to the radiative output of high accretion rate AGN. However, these different feedback mechanisms, acting on two scales are typically not investigated within the same objects. For example, potential radio-jet-related processes, are often assumed to be sub-dominant in radiatively luminous AGN with low to moderate radio powers (such as ‘radio quiet’ quasars). However, recent observational studies have come to highlight the importance of low- and moderate-power radio jets (P|$_{\mathrm{jet}}\le \, 10^{45}$| erg s⁻¹) in galaxies, which are traditionally classified as ‘radio quiet’ (because their radiative output dominates over that from jets). Low- and moderate-power jets in these systems have been observed to be driving turbulence, outflows, and excitation of the molecular gas (e.g. Morganti et al. 2015; Rosario et al. 2019; Girdhar et al. 2022; Audibert et al. 2023; Morganti et al. 2023), which have an impact which is, at least qualitatively, expected from the jets as seen in simulations (Mukherjee et al. 2016; Meenakshi et al. 2022; Tanner & Weaver 2022; Morganti et al. 2023). This all motivates an observational study to search for the impact on molecular gas, on multiple spatial scales, in systems that contain both radio jets and high luminosity AGN. With this goal in mind, we make use of the spatially resolved, multi-wavelength data from the Quasar Feedback Survey (QFeedS; Jarvis et al. 2021).

QFeedS² includes 42 quasars (L_bol ≳ 10⁴⁵ erg s⁻¹) that were selected from the parent population of AGN at z ≤ 0.2 by Jarvis et al. (2021). These luminosities are representative of the peak of the luminosity function, L_*, at the peak of the cosmic epoch of growth when quasar feedback is also expected to dominate (i.e. z ∼ 1). However, the low redshift provides the advantage to obtain spatially resolved, sensitive observations of such powerful quasars. While the galaxies studied in this sample are seen to be gas-rich and star-forming, it is a caveat that the conditions of the interstellar medium (ISM) may be different for the host-galaxies of quasars at z ∼ 1, and are not complete analogues of high-redshift AGN. The QFeedS data set is being used to extract information on the origin of radio emission in ‘radio quiet’ quasars; multi-phase outflows; and the impact of AGN on the host galaxies (Harrison et al. 2015; Lansbury et al. 2018; Jarvis et al. 2019, 2020, 2021; Girdhar et al. 2022; Silpa et al. 2022; Molyneux et al. 2023).

One benefit of QFeedS, for exploring different feedback mechanisms, is the availability of sensitive and high-resolution (HR) radio imaging provided by the Karl G. Jansky Very Large Array (VLA; Jarvis et al. 2021). In this work, we explore the feedback on the molecular gas of these quasar-host galaxies by comparing the radio emission with respect to the spatial distribution and kinematics of the molecular gas, traced with CO (3–2) transition, with data from the Atacama Large Millimeter/submillimeter Array (ALMA). The main focus is to compare with the prior feedback studies performed in BCGs, which look for molecular structures associated with radio lobes (Russell et al. 2019; Tamhane et al. 2022); and to also simultaneously search for the presence of CO emission-line wings (as a tracer of molecular outflows), as performed for a compilation of z < 0.2 AGN and starburst galaxies by Fluetsch et al. (2019).

This paper is structured as follows. In Section 2, we discuss the sample selection and the different observations and their reduction used for this analysis. In Section 3, we describe the approach for the emission-line fits to extract the molecular gas kinematics, followed by the stellar kinematics (using data obtained on the Very Large Telescope’s Multi Unit Spectroscopic Explorer; VLT/MUSE) and the methods used to extract morphological and kinematic properties of the molecular gas. In Section 4, we present the results and a discussion of these results, in the context of previous observations and simulation studies. Finally, in Section 5, we present our conclusions.

We have adopted the cosmological parameters to be H₀  = 70 km s⁻¹ Mpc⁻¹, Ω_M  = 0.3, and Ω_Λ  = 0.7, throughout. In this cosmology, 1 arcsec corresponds to 2.47 kpc for the redshift of z = 0.14 (i.e. the average redshift of the galaxies studied here). When referred to, we define the radio spectral index, α, using the relation, S_ν ∝ ν^α, where S_ν refers to the flux density at the corresponding frequency ν.

2 TARGETS, OBSERVATIONS, AND ANCILLARY DATA

We select our targets for this work from QFeedS (Jarvis et al. 2021), a survey of 42 sources that were originally selected from the parent population of emission-line AGN at z ≤ 0.2 (Mullaney et al. 2013), with quasar-like [O iii]|$\, \lambda$| 5007-Å luminosities (L_[O iii] > 10^42.1 erg s⁻¹). A moderate radio luminosity criteria of L_1.4 GHz > 10^23.45 W Hz⁻¹ is also applied to obtain the QFeedS sample; however, the sample still consists of 88  per cent ‘radio quiet’ sources, based on the criteria of Xu, Livio & Baum (1999) (see Jarvis et al. 2021 for full details). This is consistent with the ‘radio quiet’ fraction of the overall quasar population (i.e. ∼90  per cent; Zakamska et al. 2004).

Fig. 1 shows the [O iii] luminosities and the projected largest linear radio sizes, LLS_radio, of the 42 quasars from QFeedS. The LLS_radio measurements were calculated in Jarvis et al. (2021) from a set of 1.5–6 GHz VLA images, with a resolution ranging from 0.3 to 1 arcsec. LLS_radio is defined as the distance between the farthest radio emission peaks in the lowest resolution image where the source shows radio structures. If the source shows no morphological features beyond the core in any image, LLS_radio is defined as the beam de-convolved size of the core. The sample exhibits a wide range of radio sizes (∼0.1–∼60 kpc). Section 2.1 gives an overview of the sources selected from QFeedS for this study, followed by a description of the data used from ALMA, VLA, and VLT/MUSE (Section 2.2–2.4).

2.1 Sample Selection

Following our goal to analyse the molecular gas, we identified 9 out of 42 targets from the QFeedS sample that have available CO 12-m array ALMA data (highlighted by the empty black circles in Fig. 1). For all nine of these, CO(3–2) data are available in ALMA Band 7 observations; therefore, we decided to use this as our tracer of molecular gas for this work. We note two out of nine also have CO(2–1) data and one out of nine has CO(1–0) data that are presented in Ramos Almeida et al. (2022), Audibert et al. (2023), and Sun et al. (2014). To clearly separate CO emission related to galaxy discs from any extended emission outside the discs, which could be associated with extended radio lobes³, we further only selected the sources with LLS_radio ≥10 kpc (from VLA data; see Section 2.3). As shown in Fig. 1, four out of nine sources with ALMA data meet this criteria; namely, J0945+1439, J1000+1242, J1010+1413, and J1430+1339. Table 1 lists the basic properties of these four targets. The CO maps and radio images are shown for these four targets in Fig. 2 (centre panel; see Sections 2.2 and 2.3 for details).

All four targets are bright AGN with high, quasar-like bolometric luminosities of log  (L_bol/erg s⁻¹)  =  45.3–46.2 (from the fitting of the spectral energy distributions; Jarvis et al. 2019). These four targets are all classified as ‘radio quiet’ based on the L_1.4 GHz versus L_[O iii] criteria of Xu, Livio & Baum (1999). However, despite their modest radio luminosities of log (L_1.4 GHz/W Hz⁻¹)  =  23.7–24.3, all four of these targets have been confirmed to have an excess of radio emission over that expected from star-formation from their radio imaging (see Jarvis et al. 2019, 2021). HR VLA data at 6 GHz (see Jarvis et al. 2021) reveals collimated structures along with the presence of hotspots consistent with a jet morphology. Furthermore, the imaging and study of the polarization data of these four galaxies suggests a jet origin of the radio emission (see Silpa et al. 2022).

These targets also have known central AGN-driven outflows and/or high levels of turbulence identified in ionized gas, as traced via broad emission-line widths (≥ 600 km s⁻¹) of the [O iii] emission, extending over the central few kiloparsecs. In all cases, the interactions of radio jets with the ISM seem to be a significant driver with possible contributions from disc winds (Harrison et al. 2014; Jarvis et al. 2019; Venturi et al. 2023). Near-infrared spectroscopy of J1430+1339 and J0945+1737 further reveals multiple ionized outflow components through different gas tracers (Ramos Almeida et al. 2017; Speranza et al. 2022). Furthermore, for J1430+1339, there is evidence that the small-scale inner ∼1-kpc jet (Harrison et al. 2015) influences both the kinematics and excitation state of the cold molecular gas, as traced with CO (2–1) kinematics and CO (3–2)/CO (2–1) emission-line ratios (Ramos Almeida et al. 2022; Audibert et al. 2023).

In summary, these targets are well aligned with our goal to understand the impact of radio jets on the molecular gas on both small scales, within the molecular gas discs, (∼1 kpc) and on large scales (≳10 kpc), extended beyond the molecular gas discs, in ‘radio quiet’ quasars.

2.2 Observation and reduction of the ALMA data

We use 12-m array ALMA Band 7 observations to obtain the spatially resolved molecular gas emission, traced by the CO(3–2) transition. Three of the four targets (J0945+1737, J1000+1242, and J1010+1413) were observed in three 1-hr epochs under the ALMA project 2018.1.01767.S (PI: A.P. Thomson); using the C 43-4 array configuration. The chosen correlator setup comprises three spectral windows, with one spectral window covering the central frequency ν_obs  =  345.795 990 GHz, corresponding to the CO (3–2) line, and the other two windows partially overlapping the former spectral window for an enhanced signal. The fourth target, J1430+1339 was observed under the program code 2016.1.01535.S (PI: G. Lansbury) in the C40-3 configuration. The observation has a single pointing on-source integration time of 30.3 min. The spectral window has a bandwidth of 1.875 GHz and was centred at the CO (3–2) line, with the same frequency as mentioned above. The resulting angular resolutions of the observations are |$\mathrm{\theta _{res}}\sim 0.33\!-\!0.65$| arcsec (corresponding to linear scales of 0.75–1.09 kpc for the respective source redshifts) and a largest angular scale θ_LAS ∼ 4 arcsec (i.e. 4.6–6.4 kpc) for the former three and ∼ 19 arcsec (∼ 30 kpc), for J1430+1339.

The data for the four galaxies were reduced and calibrated using the Common Astronomy Software Applications (casa; McMullin et al. 2007). Using casa v6.4.3, the imaging of the cubes was made with the task tclean. The cleaning was performed in a mask centred in the peak luminosity pixel of each galaxy, with a radius varying between 5 and 8 arcsec to make sure all the resolved emission was included. A channel width of 25 km s⁻¹ was selected, and a pixel scale of 0.05 arcsec was used to sample all the synthesized beams. The Högbom CLEAN algorithm was run to a flux density threshold of 2 times the root mean square of each of the cubes. Different weightings of the Briggs mode were compared: a robustness = 2.0 (close to natural weighting), a robustness = 0.5 (between uniform and natural weighting), and applying a tapering of the visibilities in the u–v plane. We decided to use the robustness  = 2.0 to maximize the recovery of the extended emission without losing significant flux, or the central source small-scale structures. For this work, we use cubes with the continuum and allow a line component to fit the continuum (see Section 3.2). The final beam sizes of the observations were an average of 0.37 arcsec × 0.28 arcsec for the first three targets and 0.7 arcsec × 0.6 arcsec for J1430+1339.

2.3 Summary of the radio images

For our investigation of the relationship between the CO emission-line properties and radio morphology, we use the 6-GHz (C band) VLA radio images from Jarvis et al. (2019). We use both the ‘low-resolution’ (LR) and ‘high-resolution’ (HR), images described from Jarvis et al. (2019), which are constructed from a combination of A- and B-configuration VLA data. The 6-GHz (C band) VLA radio images were obtained by optimizing the imaging parameters and weighting schemes to enhance the extended morphological features for each source (i.e. using uniform, natural, or Briggs parameter = 0.5 weighting; see table 3 in Jarvis et al. 2019). In each of the figures in this work, we use these 6-GHz (C band) VLA radio images from Jarvis et al. (2019), to show the range of radio structures seen on the different spatial scales (e.g. Fig. 2). We note that this is why we prefer these over the radio images from Jarvis et al. (2021), where a simpler, but consistent, set of imaging parameters were applied to the whole QFeedS sample (but not optimized to show all morpohlogical structures). The LR images have major axis beam sizes of ∼1.0–1.2 arcsec (i.e. ∼ 2-kpc resolution for z = 0.1), whilst the HR images have beam sizes of ∼0.2–0.3 arcsec (i.e. ∼ 0.5-kpc resolution for z  =  0.1). For J0945+1737, we also make use of the 1.5-GHz e-MERLIN image from Jarvis et al. (2019), which has a beam size of 0.2–0.3 arcsec. This image highlights the ∼ 2.1-kpc bent jet-like structure in this source.

2.4 Stellar velocities and velocity dispersion from MUSE data

For this analysis, we are interested in measuring the stellar redshift (z_*) and stellar velocity dispersion (σ_*), integrated over the spatial extent of the galaxies. To do this, we use the available MUSE data for these targets and follow the procedure outlined in Girdhar et al. (2022) for another QFeedS target, making use of spectral fitting of stellar templates. Full details of the MUSE data and its reduction for these, and other QFeedS targets, are deferred to future works (e.g. Venturi et al. 2023). Therefore, we only provide brief details here.

The four targets have been observed with MUSE, in wide-field mode. This provides a field of view of 1 × 1 arcmin² and a pixel sampling of 0.2 arcsec. Observations of these targets were taken under proposal IDs 0103.B-0071 (PI: C. Harrison), 0102.B-107 (PI: Sartori), and 0104.B-0476 (PI: G. Venturi). We combine the data from these multiple programs to construct deep final stacked cubes, following the data reduction steps described in Girdhar et al. (2022).

We obtained stellar kinematics by employing the GIST pipeline (Bittner et al. 2019), following the detailed methodology and parameters described in Girdhar et al. (2022). GIST is a framework that inputs fully reduced MUSE cubes and prepares them for stellar continuum fitting to finally provide the stellar kinematics as per the following steps. Firstly, GIST performs a Voronoi tessellation routine (Cappellari & Copin 2003) to divide the galaxy into regions with a minimum signal-to-noise-ratio (SNR) in the continuum. We used a threshold of SNR  =  30. Following this, for each Voronoi region, GIST obtains the best fit to the stellar continuum, exploiting the pPXF routine (Cappellari & Emsellem 2004; Cappellari 2017), with a combination of stellar templates from XSL Library (Arentsen et al. 2019; Gonneau et al. 2020). This results in an accurate measure of the stellar velocity and stellar velocity dispersion in each Voronoi bin. We used flux-weighted averaging over all the Voronoi bins for the systemic redshift (z_*) and the stellar velocity dispersion (σ_*) of the galaxy. These systemic redshifts are used to shift all the molecular gas emission profiles to the rest frame (see Section 3.2). The errors in the values are determined as a median of the formal errors across all the Voronoi bins for each target, where formal errors are 1σ uncertainties as obtained by pPXF fitting. For the four targets, we obtain stellar velocity dispersion values between 170 and 270 km s⁻¹. All values for z_* and σ_* are listed in Table 2.

3 ANALYSIS OF THE CO EMISSION

In this section, we present our analysis steps to obtain the observed and derived properties of the molecular gas on different scales; for the whole galaxy, and for various spatially resolved scales. In Section 3.1, we formalize our approach to identify any molecular gas structures outside of the central galaxy discs (i.e. the extended molecular gas structures). In Section 3.2, we describe our emission-line fitting procedure to characterize the kinematics of the CO emission. In Section 3.3, we evaluate the properties of the molecular gas structures (velocity, velocity dispersion, projected extent, and estimated masses). Finally, in Section 3.4, we present a brief analysis of the broad CO emission-line wing components, as a tracer of central molecular outflows.

3.1 Identification of extended molecular gas structures

In Fig. 2, we show CO (3–2) emission-line images, collapsed over the full observed velocity width of the emission-line profiles (the velocity limits are indicated at the top right-hand side of each panel). Overlaid on these images are contours of the 6-GHz radio emission for both the LR and HR images (shown as white and red contours, respectively; see Section 2.3). In two of the galaxies, J0945+1737 (panel a1) and J1430+1339 (panel d1), we see that the CO emission is only observed in a central, contiguous region. However, for the other two galaxies, J1000+1242 (panel b1) and J1010+1413 (panel c1), in addition to the central molecular gas, we see molecular CO structures outside of the central regions.

To discern extended molecular gas structures in a systematic way, we formalized the following procedure, which is motivated by the qualitative methods used by Russell et al. (2019) and Tamhane et al. (2022) to search for extended molecular gas structures around BCGs. These two works will serve as our primary comparison sample (discussed in Section 4.2). Firstly, we define the central molecular gas discs as central, contiguous CO structures with smooth velocity gradients centred on the systemic galaxy velocities.

We used a visual inspection of the narrow-band images and individual velocity slices, as well as the kinematic maps, to identify extended molecular gas structures based on the following criteria:

• Emission with a clear morphological and/or kinematic separation from the central molecular gas disc.

• Emission detected at a ≥ 5σ significance level in individual velocity channels of the data cube, but also seen in more than one consecutive velocity channel.

• Structures extending to ≥ 1 kpc in projected size.

Using this systematic approach, for J0945+1737 and J1430+1339, we do not identify any extended molecular gas structures away from the central emission. This confirms the observation made from the total CO emission-line images shown in Fig. 2. These two quasars are hence not included in our analysis pertaining to the extended CO gas structures.

For J1000+1242 and J1010+1413, we identified five and three molecular gas structures, respectively, using the above method. We created narrow velocity slice CO images by collapsing over the consecutive velocity channels where any emission |$\ge \, 5\, \sigma$| was seen associated with these structures. These are shown in Fig. A1. A combined overview of the molecular gas structures is shown in Fig. 3. For this figure, to distinctly visualize each of these molecular gas structures, we performed a weighted combination of the narrow velocity slices of each filament shown in Fig. A1, with higher weights to the structures with lower surface brightness. Each filament is highlighted with a surrounding dashed grey box, and labelled following the notation as F1–5, respectively. These boxes cover the full observed extent (at |$\ge 3\, \sigma$|⁠) of the structures. In Section 3.3, we evaluate the properties of each of these structures.

For all eight identified gas structures, we estimated an axis ratio by fitting a 2D Gaussian over the surface brightness images of each. The uncertainty in the axis ratio was obtained by using the errors in the measurements of the major and minor axis of the 2D fitted Gaussian, and then propagating the errors to obtain the uncertainty on the axis ratio. These values are listed in Table 3. The axis ratios range from 1.2 to 6.5, with a median of 2.7, and all but two have an axis ratio >2. Therefore, we refer to each of the identified molecular gas structures as ‘filaments’. This follows the terminology adopted for the molecular gas structures seen in BCGs, which show a similar range in morphology (e.g. Russell et al. 2019; Tamhane et al.2022).

We note that our approach of selecting these structures may not uniquely select physically distinct ‘filaments’. For example, the identified structures may be part of a larger connected ‘flow’, and/or each ‘filament’ can contain sub-structures. Indeed, the contours in Fig. 3 do show sub-structure. However, we have used the systematic approach described above to identify these structures, without laying much emphasis on the sub-structures. We note that this does not affect our scientific goals, which are primarily to compare the properties of these overall structures to those seen in BCGs. Further, the filaments analysed in the BCGs were also identified using a very similar approach and definitions. We discuss the origin and properties of these structures in Section 4.

3.2 Emission-line fitting procedure

We evaluate the molecular gas velocity and velocity dispersion by performing fitting to the CO emission-line profiles. As described below, we studied the properties on various spatial scales from the datacubes: (1) in individual spatial pixels to produce maps; (2) from a region covering the entire CO emitting region for the full galaxies (shown through dashed grey boxes in Fig. 2); (3) regions covering each of the individual filaments (shown through dashed white boxes in Fig. 3); and (4) central ‘outflow’ regions where we identify broad CO emission-line wing components (shown through dashed purple boxes in Fig. 2).

We used the scipy curve fit routine (Virtanen et al. 2020b) to obtain the best fit to the data, within the velocity range of ± 700 km s⁻¹. We modelled the emission-line profiles using one and two Gaussian components; in addition to a linear component for characterizing any underlying continuum emission. To statistically select the best fit to the data, we used the difference in the BIC values (Bayesian information criterion; Schwarz 1978), i.e. the model with the lowest BIC value was selected.

While we use multiple Gaussian components, to characterize the emission-line profiles, we adopt a non-parametric approach for most of our analysis (following e.g. Harrison et al. 2014; Girdhar et al. 2022). The bulk velocity is measured in terms of the median velocity of the line profile, V₅₀; and the velocity dispersion is measured through velocity width in terms of W₈₀, i.e. the width containing 80  per cent of the emission-line flux. For a single Gaussian, W₈₀ is approximately related to the full-width at half-maxima (FWHM) as W|$_{80}\, =\, 1.088\, \times$| FWHM; where the FWHM itself can be related to σ as FWHM  =  2.35σ.

All line profiles and velocity maps presented in this work are shifted from the observed to the rest-frame using the stellar systemic redshift as listed in Table 2 (following Section 2.4). The application of this emission-line fitting process to different spatial scales is explained in detail as follows:

• For the entire galaxy scale:

  For the total CO emission, we extract the spectrum over the region shown through a grey-dashed box in Fig. 2 for all four targets. The obtained emission-line spectra along with the best fit following the fitting procedure above are shown for each galaxy in the panels (a2), (b2), (c2), and (d2), respectively.

• For individual spatial pixel fits:

  To apply the emission-line fitting routine on individual spaxels, we first re-grided the ALMA cubes from an initial spatial resolution of 0.05 × 0.05 to 0.15 × 0.15 arcsec², to increase the SNR of each spatial-unit. We first checked that the SNR ≥ 3 for the emission line to be considered as detected. For a further conservative check, we then compared the single Gaussian fit, with a simple straight line fit using ΔBIC. If the line fit had a lower BIC value, the spaxel was discarded from further kinematic analysis. After these checks to confirm a detected emission line, the line fitting routine continued as above. For the two targets that show the filamentary molecular gas structures, the resulting kinematic maps are illustrated in Figs 4 and 5 for J1000+1242 and J1010+1413, respectively.

• For individual molecular gas structures:

  To quantify the kinematic properties of the molecular gas structures (the ‘filaments’), we fit the integrated spectra from the region shown through dashed, gray boxes in Fig. 3. Kinematic maps for the individual structures, along with the extracted CO emission-line profiles and fits, are shown in the Appendix Figs B1 and B2.

3.3 Molecular gas properties

In this section, we measure the properties of the CO filaments seen in our sources. To compare the properties to similar structures seen in the BCGs at comparable redshifts (z ≲ 0.2); we follow methods motivated by the techniques used in Russell et al. (2019) and Tamhane et al. (2022). These works have performed an extensive study of the morphology and kinematics of the filaments for 14 unique BCGs (across both samples). The measured CO properties for the whole galaxy measurements are listed in Table 2. For the individual filaments, the properties are provided in Tables 3 and 4.

3.3.1 Molecular gas velocity

To measure the bulk velocities in the molecular gas, we utilize the median velocity V₅₀ (see Section 3.2) derived from the emission-line profiles depending on the spatial level in consideration, i.e. at the galaxy level, the velocity will be V_{50, gal} (or V_gal for simplicity). This median velocity is calculated from the fit to the galaxy-level spectra as shown in the top right-hand panels of Fig. 2. The velocity maps for J1000+1242 and J1010+1413 are shown in the left-hand panels of Figs 4 and 5, respectively. For the filamentary molecular gas structures, the velocity of the filaments (or V_fil) is obtained from the spectra extracted from the filamentary regions (see appendix Figs B1 and B2). To estimate the uncertainties for our velocity values, we performed Monte Carlo (MC) simulations. For this purpose, multiple simulated representations of the CO emission-line data were obtained by adding random Gaussian noise to our best-fitting model (on the scale of the residual noise in the continuum). A fit was obtained for each of these simulated spectra with V₅₀ measured each time. The standard deviation of the distribution of the measured V₅₀ values was then used as the uncertainty. We also compared our method with that used by Tamhane et al. (2022), where filament velocities are defined as ‘flow velocities’ (V_flow). They obtain their value from a flux-weighted average, over the filament regions, from a velocity map, and always use a single Gaussian component for their fits. Following the approach used by Tamhane et al. 2022, the resulting velocities are similar (i.e. within ≤ 5  per cent), as compared to the the previously described V₅₀ method. We also note that in the case of more than one filament for a galaxy, Tamhane et al. (2022) only provides an average value over all the filaments.

3.3.2 Molecular gas velocity dispersion

To obtain the velocity dispersion, we refer to the analysis of Russell et al. (2019). They define the molecular velocity dispersion σ_gal (σ_mol in Russell et al. 2019) as the σ width of a single Gaussian component, fitted to the CO line emission over the entire galaxy. They also separately estimate the CO line velocity dispersion of the filaments, σ_fil, by fitting a single Gaussian for the emission-line profiles obtained only over the individual filamentary regions. We follow the same process to obtain the CO velocity dispersion for the total galaxy spectrum and for each of the extended filamentary structures in our sample. While we allow for multiple Gaussian fits, for consistency, we obtain an equivalent σ value from the non-parametric values of velocity width, following |$\sigma \, =\, \mathrm{\it W}_{\mathrm{80}}/2.6$|⁠. To obtain the uncertainty in the obtained values, we perform MC simulations as described for estimating the molecular gas velocity (see Section 3.3.1).

3.3.3 Radial extent of filaments

We quantify the maximum projected radial extent, R_fil, of the filamentary structures as the maximum projected distance from the nucleus to the most distant part of each filament. For this, we used the images collapsed over the narrow velocity ranges in which the individual filaments were detected (see Fig. A1). We measured the spatial separation between the centre of the galaxy which was identified using the position of the radio core; and the farthest point of the filaments. This was measured for a few spaxels in the filaments and the maximum projected distance was used as a measure of the R_fil for each filament. This was done following the same method as Tamhane et al. (2022) (defined as R_flow in their work). The uncertainty on the R_fil values are estimated to be the equivalent size of the beam’s major axis.

3.3.4 Molecular gas mass estimates

To obtain the total molecular gas mass in the galaxy (M_tot), we first measured the integrated line flux for the CO (3–2) emission from the galaxy spectra. Likewise, for the filament mass (M_fil), we used the integrated flux from the spectra extracted from filamentary regions. The uncertainty in flux values is obtained from the emission-line fitting routine as the square root of the diagonal of the covariance matrix of each free parameter used for the emission-line fit.

To estimate the molecular mass, we used the same conversion factors as Tamhane et al. (2022), for a consistent comparison. We first converted CO (3–2) fluxes to CO (1–0), by using integrated line flux ratios of 7.2 (Vantyghem et al. 2016). We then converted the integrated flux density of CO(1–0) line (S|$_{\mathrm{CO}}\, \Delta \,$|v) into molecular gas mass (M_mol) using the following relation (Solomon & Vanden Bout 2005; Bolatto, Wolfire & Leroy 2013):

where z is the redshift of the galaxy, D_L is the luminosity distance, and X_CO is the CO-to-H₂ conversion factor, with X_{CO, Gal}  =  2 × 10²⁰ cm⁻² (K km s⁻¹)⁻¹ (Solomon et al. 1987; Solomon & Vanden Bout 2005). This relation is the corollary of the relation |$\mathrm{\it M}_{\mathrm{mol}}\, =\, \alpha \, \mathrm{\it L}^{\prime }_{\mathrm{CO(1-0)}}$|⁠, where X_CO and α_CO are both referred to as ‘CO-to-H₂’ conversion factor. For X_CO =  2 × 10²⁰ cm⁻² (K km s⁻¹)⁻¹, the corresponding α_CO  = 4.3 M_⊙ (K km s⁻¹ pc²)⁻¹. It is a caveat that this exact X_{CO, Gal} factor may not apply to our galaxies (and also neither to BCGs), and there may be significant variation, pertaining to environmental variations (reviewed by Bolatto, Wolfire & Leroy 2013). However, we use these factors for a consistent comparison between these different studies, and we assume a systematic uncertainty of ∼ 0.5 dex on any derived quantities related to molecular gas masses, throughout (following Tamhane et al. 2022).

3.4 Spatially mapping central outflows in molecular gas phase

We also aim to characterise any central molecular outflows, which are often traced with underlying wings in the CO emission-line components. Although a detailed kinematic analysis of molecular outflows is beyond the scope of this work (following e.g. Ramos Almeida et al. 2022), we compare the CO properties to previous works that investigate such CO components, under the assumption that they are tracing outflows. We focus our comparison to Fluetsch et al. (2019), who study the CO outflow kinematics for 45 active galaxies (starburst and AGN) at z<0.2, with L|$_{\mathrm{AGN}}\, \sim \, 10^{40}\!-\!10^{46}$| erg s⁻¹. Therefore, we are motivated by the methods of Fluetsch et al. (2019) and consequently, we measure the spatial extent of the region over which CO emission-line wings are identified.

Following Section 3.2, we mapped the CO emission in the central regions, identifying pixels where two emission-line components were required. The velocity width maps within the central regions reveal broad velocity widths across all four galaxies (i.e. ≳ 400 km s⁻¹). This motivated us to undertake a more comprehensive analysis to identify any disturbed gas in the central regions of the four targets. The BIC-based selection was effective in selecting the required number of Gaussian components for obtaining the fits. However, acknowledging the complexities of emission-line kinematics, we also visually inspected the fits to identify the regions that show clear signs of a secondary, underlying high-velocity wing component (as opposed to two narrow components). The regions over which broad CO wings are clearly identified are shown as dashed purple boxes in the panels a1, b1, c1, and d1 of Fig. 2. Using these regions, we extracted the spectral profile for studying the properties of the central outflows (shown in the respective a3, b3, c3, and d3 panels of Fig. 2). We measured R_out as the projected distance between the farthest spaxel from the central spaxel over these regions. The uncertainty in the projected distance was taken as the major axis of the respective beams.

We use the CO emission-line profiles from these central regions to obtain the velocities of the outflowing gas as follows: V_out = FWHM_broad/2+V_broad, i.e. the same definition as Fluetsch et al. (2019). For the uncertainty in the velocity values, we combined the uncertainties in FWHM and V_broad (see Section 3.3). The central outflow properties are listed for all four targets in Table 5 and plotted in Fig. 6. We discuss these properties later in Section 4.3. We note that when we employ the same methods used by Fluetsch et al. (2019) (which involves a simplified approach of producing CO images over the high-velocity wings of the profiles), we obtain very similar values, and our derived values are also close to the previous studies of the CO emission for the case of J1430+1339 (see Audibert et al. 2023). In summary, our values are sufficient for the simple parameter-space comparison of CO emission-line profile properties presented in Fig. 6.

4 RESULTS AND DISCUSSION

In this section, we present our results, and discuss their interpretation, from our analysis of the molecular gas (traced via CO (3–2) emission) of four quasars from the QFeedS (Fig. 1). Specifically, in Section 4.1, we summarize the properties of the identified extended molecular gas structures in terms of their morphology, radial extent, and kinematics. In Section 4.2, we make a comparison to similar structures found in BCGs. Hence, in the first two sections, we discuss the interaction of the radio lobes with the molecular gas at larger scales. This is followed by Section 4.3, where we present the observations of the central molecular gas outflows. Finally, in Section 4.4, we discuss the evidence for two feedback mechanisms acting on the molecular gas, in the same targets, and discuss possible implications for an evolutionary sequence of feedback via low- and moderate-power radio jets in ‘radio quiet’ quasars.

4.1 Properties of the extended molecular gas structures

Fig. 2 reveals molecular gas in the form of extended filamentary structures for two of the four galaxies. As presented in Section 3.1, these structures have morphologies that are typically elongated (with a median axis ratio of 2.7; see Table 3). Following the terminology used for similar morphological structures seen in BCGs, we refer to these gas structures as filaments. As shown in Fig. 3, we identify five filaments for J1000+1242 and three for J1010+1413. Figs 4 and 5 show the velocity and velocity-width maps (in terms of W₈₀) over the entire CO-emitting regions. A zoomed-in version of these maps, and corresponding CO (3–2) emission-line profiles extracted from the regions of the filaments, are provided for the individual filaments in Appendix  B. The observed filament properties are listed in Table 3.

We present values of filament velocity and radial extent in Fig. 6, as teal-coloured triangles for J1000+1242 and stars for J1010+1413. For J1000+1242, across the five filaments, there is a radial extent range of R_fil  =  5–11 kpc and velocities of V_fil  =  | 100–190 | km s⁻¹. In case of the three filaments in J1010+1413, we see a radial extent range of 12–13 kpc with comparatively higher velocities of 220–340 km s⁻¹. This gives us an average radial extent of 8 and 12 kpc; and an average velocity of ∼ 150 and ∼ 280 km s⁻¹ for J1000+1242 and J1010+1413, respectively.

In Fig. 7, we compare the molecular velocity dispersion of the filaments (σ_fil) with the stellar velocity dispersion (σ_*) of the host galaxies. The filaments show velocity dispersion values in the range of 30–130 km s⁻¹ for J1000+1242 and 47–90 km s⁻¹ for J1010+1413. In general, the velocity dispersion of the filaments is much lower than the stellar velocity dispersion values, with a median ratio of σ_fil/σ_*  =  0.32 across all eight filaments.

In the left-hand panel of Fig. 8, we present the fraction of total molecular gas mass located in the filaments, with respect to the radio luminosity (L_1.4 GHz), where the colour-scaling corresponds to the galaxy’s stellar velocity dispersion. These ratios are simply the ratio of CO (3–2) flux across all filaments, divided by the total CO (3–2) flux for each galaxy. That is, we are assuming the same CO flux to mass conversion factor for both the filaments and total gas mass. J1000+1242 is observed to have the highest value, with 53  per cent of the gas located within these structures. J1010+1413 has only 9  per cent of the CO(3–2) emitting gas located in these structures. For J0945+1737 and J1430+1339, where no filaments were detected, we estimated upper limits for the molecular gas mass of the filaments of 4  per cent, and 9  per cent respectively. For estimating the upper limits, we used the flux ratio of the faintest detected filament (i.e. filament 3 of J1010+1413) to total galaxy flux and scaled it by the noise in the respective cubes of J0945+1737 and J1430+1339. We note that these mass ratio measurements can be affected by the sensitivity to structures on different scales, depending on the distribution of the CO (3–2) emitting gas. For example, we may be missing low-surface-brightness CO (3–2) emission (either contained in filaments or the central galaxy). Towards this, we compared our total CO flux measurements from the 12 m ALMA observations with single-dish observations of the three targets detected in CO (3–2) in APEX data (i.e. all but J0945+1737; Molyneux et al. 2023). We found that ALMA/APEX flux ratios range from 0.53–1.3. Although this adds some additional uncertainty (at the ∼0.3-dex level on these mass ratios), our measurements are sufficient for a broad comparison to the values for similar structures, obtained using similar data sets, seen in BCGs (Section 4.2).

The CO (3–2) emitting molecular gas seen as elongated structures in the two galaxies, appear to be entrained along or around the radio bubbles seen in these targets (see Figs 3, 4, and 5). The surface brightness and velocity maps of the filaments reveal some clumpy sub-structures, but the velocity gradients across the filaments are relatively smooth, and are typically small (30–100 km s⁻¹), following the major axes of the structures. Further, the more clumpy molecular gas is seen to be coincident with bends in the radio bubbles. We note that for J1000+1242, the velocity structures of the northern filaments seen in Fig. 4, could be consistent with seeing both the blueshifted (filaments 4 and 5) and redshifted parts (filament 1) of an expanding bubble. Qualitatively, all of these morphological and kinematic structures, and their spatial connection to expanding bubbles, are similar to those we see associated with BCGs, hosted in cool core clusters that are rich in molecular gas (e.g. see Balmaverde et al. 2018; Tremblay et al. 2018; Russell et al. 2019; Capetti et al. 2022; Tamhane et al. 2022). Therefore, it is warranted to make a more quantitative comparison between the structures observed in our new observations of quasars, with those seen in BCGs, which are not classed as ‘radio quiet’ quasars (Fig. 8).

4.2 Comparison with BCGs

We compare our observations of molecular gas structures with a sample of 14 BCGs (at z ≤ 0.2) compiled from Tamhane et al. (2022) in terms of their radial extent, velocity, and mass. For comparing in terms of velocity dispersion, we use the velocity dispersion data for only the 10 BCGs from Russell et al. (2019), that are in common over both the samples. We also note that the filament properties from Russell et al. (2019) are provided as individual filaments, whilst for Tamhane et al. (2022), they are an average over all the filaments. When referred to, we make this clear for each comparison. These works make use of ALMA to measure the properties of gas structures observed via low CO transitions (i.e. CO (3–2), CO (2–1), and CO (1–0)), similar to our observations and approaches. The stellar velocity dispersion, the radio fluxes, and the [O iii] luminosities for the BCGs are taken from Hogan et al. (2015, 2017) and Pulido et al. (2018)⁴

The BCGs are typically massive galaxies (with σ_⋆ ≈ 200–500 km s⁻¹; see Fig. 7; and M_⋆  =  10^10.6–10^12.5M_⊙) and ‘radio loud’ galaxies (see Fig. 8; right-hand panel). In contrast, our sample uniquely consists of ‘radio quiet’ quasars (see Section 2.1; Fig. 8) and has comparatively lower stellar masses (i.e. M_⋆  =  10^9.9–10¹¹ M_⊙; Jarvis et al. 2020). The BCGs, with known strong ‘radio-mode’ feedback and large reservoirs of molecular gas (10⁹–10¹¹M_⊙) serve as an interesting comparison for our targets. This helps to explore the feedback processes across different populations, and over an extended parameter space in terms of radiative and radio luminosities (Fig. 8).

In comparison to the BCGs, the CO filaments observed in our targets have comparable properties in the V_fil vs. R_fil parameter space (Fig. 6). Further, we see similar velocity dispersion values of ∼ 10–160 km s⁻¹, shown in the right-hand panel of Fig. 7, across both samples.⁵ Furthermore, the vast majority of BCG filaments fall significantly lower than half of the stellar velocity dispersion, as is also seen for those in our sample.

In the left-hand panel of Fig. 8, a comparison sample of 15 BCGs is used to observe the spread in filament mass fractions. Along with the 14 BCGs comparison sample, we also add Hydra-A to this comparison list (studied in Rose et al. 2019), only for this plot, for an overall representation of the parameter space covered by the BCGs. They show a significant spread in the filament mass fraction from 0  per cent in Hydra-A (disc-dominated) to ∼90  per cent in AS1101 (filament-dominated). Three of the four galaxies from our sample lie towards the lower end of this filament mass-fraction, with J1000+1242 lying towards the middle at 53  per cent. Due to the archival and inhomogeneous nature of the BCG sample, and the small sample of ‘radio quiet’ quasars investigated here, it is not yet possible to rigorously assess if the distribution of filament-to-total molecular mass fractions of the two populations is consistent. A more complete, systematic survey of the two populations is required.

Overall, in Fig. 8, we look for trends in the observed filament-to-total molecular gas fractions, with respect to the galaxy properties, such as, radio luminosity (L_1.4 GHz), [O iii] luminosity (L_[O iii]), and stellar velocity dispersion (σ_*) for the combined sample of BCGs plus our sample. In the left-hand panel of Fig. 8, we see no clear trends with radio luminosity, nor with stellar velocity dispersion (represented by the colour-scaling). In the right-hand panel of Fig. 8, we compare our targets to the BCG sample in the L_1.4 GHz vs. L_[O iii] plane; with the colour-scaling corresponding to the filament mass fraction. Our sources uniquely lie in the radio-quiet but radiative quasar regime. However, in terms of radio-luminosity, we see no obvious trend combining this sample and the BCG sample. In contrast, in terms of the radiative luminosity, the BCG sources with a higher [O iii] luminosity (L_[O iii] ≥ 10^42.1 erg s⁻¹; in the quasar regime) have a typically lower mass-fraction of gas in the filaments (∼ 22  per cent), compared to the higher average mass-fraction (∼ 52  per cent) seen for those with lower [O iii] luminosities. The average filament mass-fraction for luminous [O iii] sources reduces further when our targets are also included (≤ 16  per cent). This appears to indicate that higher radiative power does not result in higher fractions of mass involved in these filaments. Nonetheless, we reiterate that a more homogeneous and complete investigation across both samples is required to confirm any such trend (or lack thereof) between radiative power and the fraction of molecular gas located in filamentary structures.

4.3 Central multi-phase outflows

All four quasars are already known to contain central ionized gas outflows (traced by broad emission-line components) from previous work (see Harrison et al. 2014, 2015; Ramos Almeida et al. 2017; Speranza et al. 2022; Venturi et al. 2023), with velocities reaching 1000 km s⁻¹, and extending to spatial extents of ≳1–10 kpc. In all cases, the jet–ISM interactions have been proposed as an important driving mechanism of outflows and turbulence, with possible additional contributions from quasar-driven winds.

As reported in Section 3.4, we have found evidence of high-velocity wings in the CO (3–2) emission-line profiles in the central regions of the galaxies. Such signatures are attributed to outflowing molecular gas in Fluetsch et al. (2019). Following their definition (also see Section 3.4), we observe the velocities of the central disturbed gas to be V_out  =  249, 376, 441, and 331 km s⁻¹ for J0945+1737, J1000+1242, J1010+1413, and J1430+1339, respectively. We measure the projected radial extent of this outflowing phase to be R_out  =  0.60, 0.65, 1.64, and 0.80 kpc, for these same targets. In Table 5, we summarise all the properties, along with the measured uncertainties of the central outflows. We note that for J1430+1339, similar properties of outflowing molecular phase using the CO (3–2) and CO (2–1) emitting gas were presented in Ramos Almeida et al. (2017) and Audibert et al. (2023), which is attributed to the inner radio jet seen in this source (Harrison et al. 2015; Jarvis et al. 2019). These molecular outflow components are less extreme in both velocity and spatial extent than seen in the corresponding ionized gas, in agreement with the multi-phase study of Girdhar et al. (2022) for a different QFeedS target.

The sample studied by Fluetsch et al. (2019) covers an AGN luminosity range of 10⁴⁰–10⁴⁶ erg s⁻¹; in comparison to our sample that lies in the higher quasar luminosity regime, i.e. 10⁴⁵–10⁴⁶ erg s⁻¹. Further, both the samples lie in the same redshift range, z|$\, \leq \,$| 0.2. In Fig. 6 , we show that the velocities and radial extents of molecular outflow properties estimated by Fluetsch et al. (2019) (shown as orange squares) lie closely in the parameter space to our observed central outflow properties (yellow symbols), despite the differences in AGN luminosities.

4.4 Feedback on two spatial scales

With a goal to understand the relative importance of feedback processes in different populations, a recent study by Tamhane et al. (2022) compared the properties of molecular filaments (located around radio bubbles) in z ≤ 0.2 BCGs with the properties of high-velocity wings observed in CO emission lines (as a tracer of central molecular outflows) in the sample of z ≤ 0.2 AGN and starburst galaxies in Fluetsch et al. (2019). They conclude that radio feedback is generally more effective at lifting the gas in galaxies compared to the AGN and starburst winds. However, as acknowledged by Tamhane et al. (2022), there is not a systematic investigation of possible ‘radio feedback’ in the Fluetsch et al. (2019) sample. Nor is the same CO broad wing analyses, as performed by Fluetsch et al. (2019), applied to the BCG sample. In our study, we have searched for both types of molecular gas features (extended filaments and central outflows) in our sample of four ‘radio quiet’ quasars. In Sections 4.1 and 4.2, we showed that two of the four targets show molecular gas filaments located around radio lobes, with similar properties to those seen in BCGs. Further, in Section 4.3, we show the presence of central outflows in all four of the samples. In this section, we discuss the implications for the observed feedback effects on multiple scales.

4.4.1 Radio lobes impact on ∼10-kpc scale molecular gas

Possible explanations of the molecular gas structures seen around the radio bubbles (which contain radio jets; see Fig. 3), is either a thin cover of clumpy molecular gas, expanding along with the expanding radio bubbles, or molecular gas that is in situ condensed in the updrafts (e.g. McNamara et al. 2014, 2016; Russell et al. 2019; Zanchettin et al. 2023). The gas is then expected to appear the brightest around the edges of the bubbles, aligned with the line of sight, thus giving a filamentary appearance. In general, the filaments are observed to have slow velocities and narrow velocity widths (Figs 6, 7; and also in Russell et al. 2019; Tamhane et al. 2022). This could be because they retain the velocity structure of the rising bubbles which themselves may be relatively cool, and not shock-heated (Fabian et al. 2000; McNamara et al. 2000) as opposed to the typically energetic jet–ISM interactions.

Whilst we can not be conclusive about the origin of the molecular gas structures observed, following Tamhane et al. (2022), we assume the off-nuclear molecular gas structures as a flow, noting that some of this gas may be flowing towards the central galaxies as opposed to a pure outflow (see Russell et al. 2016; Balmaverde et al. 2022). Tamhane et al. (2022) noted that these molecular flows in BCGs are one to three orders of magnitude larger than the central outflows found in Fluetsch et al. (2019). In two of our samples, we find evidence of both types of molecular gas structures in the same sources. We find that those molecular gas structures identified around the radio lobes are roughly an order of magnitude larger in size. However, it is important to note the approaches taken to search for these types of flows are very different and are somewhat biased by the requirement to have large radio bubbles outwith the central molecular gas discs.

Again, following Tamhane et al. (2022), we compute the mass flow rate for the filaments (⁠|$\dot{\mathrm{\it M}}_{\mathrm{fil}}$|⁠) as the molecular mass in the filament (M_fil) divided by the time (t_fil) it would have taken the filament to reach the projected radial extent (R_fil) at the velocity of the filament (V_fil), which is computed as t_fil  =  R_fil/V_fil. We compute the kinetic power of the filament as follows: |$\dot{\mathrm{\it E}}_{\mathrm{kin}}$|  =  1/2 |$\dot{\mathrm{\it M}}_{\mathrm{fil}}$| V|$_{\mathrm{fil}}^2$|⁠. For all these derived quantities, we assume a systemic error of 0.5 dex on our computed values (following Tamhane et al. 2022). These derived values for the filaments are listed in Table 4.

In Fig. 9, we compare the filaments in our targets to BCGs in terms of their mass outflow rates (left-hand panel) and kinetic powers (in the right-hand panel) in relation to estimated jet powers. We obtained the jet power P_jet using the Merloni & Heinz (2007) relation. For each target, a range of jet powers were calculated using the 5.2-GHz radio luminosity corresponding to the radio core (component HR: A, in Jarvis et al. 2019) and across all the radio structures combined (LR: Total, in Jarvis et al. 2019). The variation in jet powers depending on whether the core luminosity was used or the total radio luminosity was used, varies by 0.08 and 0.15 dex for J1000+1242 and J1010+1413, respectively. We use the median of these two values for our primary calculations and data points and the range for an error bar in the figures (also quoted in Table 2). In Fig. 9, individual filaments are shown as triangles for J1000+1242 and as stars for J1010+1413, and the total for all filaments is represented using yellow symbols.

It can be seen that the mass outflow rates in the filaments for our targets (4–20 M_⊙ yr⁻¹) are comparable to the BCGs. At the same time, the kinetic power in the filaments (10⁴⁰–10⁴² erg s⁻¹) also seems to be comparable to the BCGs. However, similar to the case of the BCGs, the kinetic energy transferred in the filaments is typically a small fraction of the energy available in the radio jet. The rate of energy transferred to the filaments, under these assumptions, is observed to be lower than 5  per cent per cent of the jet kinetic power (consistent with all but two of the BCGs). Summing over all filaments, for our targets, this suggests a jet coupling efficiency of 0.0005–0.04, i.e. 0.05–4  per cent for our sample, which is significantly higher than if we were to assume the AGN bolometric luminosity was responsible for driving these flows, i.e. |$\dot{E}_{\rm kin}/L_{\rm bol}\sim$| 10⁻⁶, i.e.  10⁻⁴ per cent (see Table 4).

This brings us to the conclusion that Fig. 6 does not necessarily always correspond to two different galaxy populations for different feedback mechanisms on the molecular gas, i.e. where radiative energy (‘quasar’) drives central CO outflows in luminous, high accretion rate AGN, and radio jets drive molecular flows in typically ‘radio loud’, low accretion rate sources. In our sample, we see both central outflows and turbulence (in all four galaxies) and large-scale filamentary structures (in two). Furthermore, despite the high radiative output, it appears that radio jets and lobes also have a significant role in the impact on the molecular gas on multiple scales.

4.4.2 Dual feedback effects and evolution of moderate power radio jets

We have found two different types of impact on molecular gas in ‘radio quiet’ quasars, acting on different scales. Indeed, theoretical studies do show that specific AGN physical mechanisms are expected to result in distinct concomitant forms of AGN feedback, operating on different spatial and temporal scales, for example: (i) ultra-fast outflows/small-scale winds (Costa, Sijacki & Haehnelt 2014; Costa, Pakmor & Springel 2020), (ii) jets (Talbot, Sijacki & Bourne 2022), (iii) radiation pressure (Costa et al. 2018) and even the more phenomenological AGN feedback models used in state-of-the-art cosmological boxes (e.g. Zinger et al. 2020).

In the case of radio jets, even low- and moderate-power jets, have gained recognition as potentially causing significant disturbance through direct jet–ISM interactions from observations (Alatalo et al. 2011; Tadhunter et al. 2014; Morganti et al. 2015; Venturi et al. 2021; Girdhar et al. 2022; Morganti et al. 2023; Nandi et al. 2023). Hydrodynamic simulations of jets (see Sutherland & Bicknell 2007; Wagner, Bicknell & Umemura 2012; Mukherjee et al. 2016, 2018; Mandal et al. 2021; Tanner & Weaver 2022; Talbot, Sijacki & Bourne 2022, 2023) have studied the progression of a jet through a clumpy interstellar medium to understand its impact on the ISM.

Through these works comes a possible evolutionary sequence of jet–ISM interactions as motivated by both simulations and observations (see Morganti 2023). On small spatial scales (and shorter time-scales), the jet directly interacts with the ISM causing turbulence, and as the jet propagates, it grows along with the cocoon of shocked and heated gas and ISM plasma, that takes over at larger scales. This may cause the molecular gas to couple and rise in the wake of this growing radio bubble to greater radial extents. On larger scales (∼ 5–10 kpc), the feedback may hence be moderated by the jet-cocoon that heats or disperses the molecular gas causing cavities and pushing it aside, which is then observed as filaments. However, we cannot be certain about the origin of the filaments and it could also be possible they were cooled in situ around the radio bubbles (see Russell et al. 2019). Furthermore, a single radio event can continue to drive gas outwards for a long time after a quasar of similar power has shut down (see discussion in Tamhane et al. 2022). We suggest the presence of two ongoing effects in our targets, i.e. (i) effects on central scales due to ongoing accretion activity; and (ii) effects due to the impact of the larger scale radio bubbles. We suspect the latter may take over as the dominant mechanism on longer time-scales, as is also observed in the case of more evolved BCGs as the dominant feedback mechanism.

Our observations stem from a small number of objects and need to be confirmed by a larger sample. One significant step towards this would be to quantify the feedback from jets depending on their properties of inclination, power, and evolutionary stage. Additionally, it will be important to establish the relative importance of jets, quasar-driven winds, and radiation pressure for driving multi-phase outflows and turbulence across a homogenous sample. Further, obtaining observations of other resolved CO transitions could be crucial for deriving their excitation and physical properties and will also help unveil the full ‘population’ of filaments.

5 CONCLUSIONS

We present the study of the molecular gas properties in four z < 0.2, ‘radio quiet’ type 2 quasars from QFeedS (L_bol ∼ 10^45.3–10^46.2 erg s⁻¹; L|$_{\mathrm{1.4\, GHz}}\sim 10^{23.7}\!-\!10^{24.3}$| W Hz⁻¹), namely, J0945+1737, J1000+1242, J1010+1413, J1430+1339. These targets were selected based on the availability of high-spatial-resolution ALMA data, to trace the CO(3–2) emission, and their projected radio linear sizes of LLS_radio ≥ 10 kpc (see Fig. 1, and Table 1 for an overview of properties). We explored the kinematics and spatial distribution of the CO (3–2) emission with ∼ 0.33–1.09 kpc spatial resolution. This was compared to the morphology seen in 6-GHz radio images, previously obtained from the VLA (see Fig. 2 for a data overview). Our main findings are summarized as follows:

1. We identify filamentary molecular gas structures in and around ∼10-kpc radio lobes in two out of the four ‘radio quiet’ quasars (Figs3, 4, and 5). Both J1000+1242 and J1010+1413 show filamentary molecular gas structures that appear to wrap around the radio lobes. They have maximal radial extents of 5–13 kpc, velocities of V_fil  =  | 100–340 | km s⁻¹ and velocity dispersion values of σ_{mol, filament}  =  30–130 km s⁻¹. We observe that ∼ 53 and ∼ 9  per cent of the total molecular gas mass is contained within these structures for J1000+1242 and J1010+1413, respectively. For J0945+1737 and J1430+1339, we do not see any such structures, but estimate a maximum of 4 and 9  per cent of the total CO (3–2) emitting gas could be contained in such structures (see Fig. 8; left-hand panel). Our observations of radial extents and mass fractions are in close agreement to simulations of low- and moderate-power radio jets predicting the relation between the radial extent and the mass of molecular gas being driven (see fig. 20 in Mukherjee et al.2016).

2. The molecular gas filaments in these ‘radio quiet’ quasars have properties comparable to those seen driven by radio jets in, predominantly radio-loud, BCGs (Figs6, 8, 7, and 9). The velocities, velocity dispersion and maximal radial extent of the molecular filaments from our ‘radio quiet’ quasars show very similar values to those seen in BCGs, for which similar analyses have been performed (see Figs 6 and 7).

   The inferred mass outflow rates (4–20 M_⊙ yr⁻¹) and kinetic powers (10^40.3–10^41.7 erg s⁻¹) are also comparable to those seen in the BCGs (Fig. 9). Combining our sample with the BCGs, we observe no obvious trends with stellar velocity dispersion or radio luminosity and the fraction of mass found in filaments. Although limited by source statistics, there is tentative evidence that less [O iii] luminous sources (i.e. L_[O iii] ≤ 10⁴² erg s⁻¹) tend to show higher fractions of the molecular mass in filaments, with an average filament mass-fraction of 55  per cent, while the higher [O iii] luminous sources, including our targets, have an average filament mass-fraction of ≤ 16  per cent (see Fig. 8; right-hand panel). This may tentatively suggest kinetic power to be a more compelling driver of these molecular filaments than radiative power at larger scales (∼ 10 kpc).

3. Evidence for both central molecular outflows and large-scale radio feedback on the molecular gas in ‘radio quiet’ quasars (Figs4, 5, and 6). In all four quasars, we see evidence for central (0.6–1.6 kpc), outflows traced by high-velocity wings (V_out  =  250 – 440 km s⁻¹) of the CO(3–2) emission-line profiles (see Figs 4 and 5). These have properties comparable to those seen for the archival AGN and starburst galaxies (of typically lower L_bol; Fluetsch et al. 2019). This adds to the evidence for central multi-phase outflows in these systems, likely caused by an interaction between moderate-power radio jets and the ISM.

   We have shown that both central molecular outflows, typically associated with luminous AGN, and molecular gas filaments around radio lobes on ∼10 kpc-scales (analogous to those found in BCGs), can also be found in ‘radio quiet’ quasars. This implies that both feedback mechanisms can act within the same systems. Our observations are consistent with recent simulations and observations (e.g. Talbot, Sijacki & Bourne 2022; Morganti et al. 2023) that suggest that two feedback effects can take place due to low- and moderate-power radio jets (P_jet ≲ 10⁴⁴ erg s⁻¹). On the smaller scales, jet–ISM interactions can drive turbulence and central outflows. On larger scales, radio lobes have penetrated beyond galaxy discs and cause a more gentle pushing aside of molecular gas. Unlike high-power jets that escape swiftly, low- and moderate-power jets are seen to be trapped for longer in simulations. Their effect is amplified due to the development of an energy bubble which is basically a cocoon of a shocked ISM and plasma. This interaction leads to the constant stirring of the ISM with the energy bubble, thus inhibiting the star formation (see Mukherjee et al. 2016).

   Our results underscore that the availability of higher radiative energy as we see in quasars, does not necessarily imply that it would also be the dominant feedback mechanism on all spatial scales. We should therefore take caution in assuming different dominant mechanisms, based on if a quasar is ‘radio quiet’ or ‘radio loud’. A radio jet, even if low in power, has considerable potential to couple with the galaxy ISM and lead to a significant impact on the host galaxy gas (Mukherjee et al. 2016, 2018; Meenakshi et al. 2022). We acknowledge that all four studied targets are different to some extent and the two targets showing the presence of filaments are themselves quite unalike. To draw firm conclusions about the population as a whole, calls for a larger study, across a wider sample, to understand: how common each of these different mechanisms are; the relative importance of radio jets compared to radiative processes; and what AGN or galaxy properties determine the amount of molecular mass associated with both central outflows and larger scale radio lobe interactions.

ACKNOWLEDGEMENTS

We thank the referee for their valuable comments. We thank Helen R. Russell for helping us obtain the values for the CO velocity widths.

CMH and AN acknowledge funding from the United Kingdom Research and Innovation grant (code: MR/V022830/1). ACE acknowledges support from Science and Technology Facilities Council (STFC) grant ST/P00541/1. EPF is supported by the international Gemini Observatory, a program of National Science Foundation's (NSF) NOIRLab (formally named the National Optical-Infrared Astronomy Research Laboratory), which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the NSF, on behalf of the Gemini partnership of Argentina, Brazil, Canada, Chile, the Republic of Korea, and the USA. MB acknowledges funding support from program JWST-GO-01717, which was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA (National Aeronautics and Space Administration) contract NAS 5–03127. PK and SS acknowledges the support of the Department of Atomic Energy, Government of India, under the project 12-R&D-TFR-5.02-0700. SS acknowledges financial support from Millenium Nucleus NCN19_058 (TITANs).

This paper makes use of the following ALMA (Atacama Large Millimeter Array) data: ADS/JAO.ALMA#2016.1.01535.S and ADS/JAO.ALMA#2018.1.01767.S. ALMA is a partnership of ESO (European Southern Observatory; 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. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology. This research mainly uses the python packages: astropy,⁶ a community-developed core python package for astronomy (Astropy Collaboration et al. 2013, 2018); scipy (Virtanen et al. 2020a); numpy (Harris et al. 2020); matplotlib (Hunter 2007).

DATA AVAILABILITY

The MUSE and ALMA data presented in this analysis were accessed from the ESO and ALMA archives under the proposal IDs, namely 0103.B-0071, 0102.B-107, and 0104.B-0476, for the MUSE data and ADS/JAO.ALMA#2016.1.01535.S and ADS/JAO.ALMA#2018.1.01767.S for the ALMA data. The VLA images used in this work are available at Newcastle University’s data repository (https://data.ncl.ac.uk) and can also be accessed through our Quasar Feedback Survey website.

Footnotes

References

APPENDIX A: INDIVIDUAL FILAMENT NARROW-BAND IMAGES

To identify the filaments following our approach in Section 3.1, we created narrow velocity slice CO images by collapsing over the consecutive velocity channels where any emission |$\ge \, 5\, \sigma$| was seen associated with these structures. These are shown in Fig. A1. For J1000+1242, we identified five filaments, and for J1010+1413, we identified three filaments, as shown in each of the panels below. The velocity windows used for collapsing each image are labelled on the top right-hand side of each panel, and each filament is highlighted with a surrounding dashed white box. These boxes cover the full observed extent (at |$\ge 3\, \sigma$|⁠) of the structures, where the |$3\, \sigma$| contours are also shown. A combined overview figure for each of the targets can be seen in Fig. 3.

APPENDIX B: FIGURES OF THE FILAMENTARY MOLECULAR GAS STRUCTURES

The velocity and velocity width maps obtained following the methods outlined in Section 3.3, are shown in Figs B1 and B2 for each of the filaments of J1000+1242 and J1010+1413, respectively. From the left- to right-hand side, each of the panels in the figure represents (a) CO(3–2) emission narrow band image; (b) a map of the median velocity values in a group of (3 × 3) spaxels; (c) a map of the velocity width (W₈₀); and (d) emission-line profile extracted over the entire filamentary regions. The CO (3–2) emission map in panel (a) shows a zoom-in of each of the filaments identified following the definition in Section 3.1 and shown in Fig. 3. The green and cyan contours show the emission at 3σ and 5σ, respectively. The two subsequent panels (b) and (c) represent the spatially resolved kinematics for each of the filaments with the 5σ contour highlighting the edges and the colour bars at the bottom. The final panel (d) shows the integrated emission-line profile over the entire filamentary region enclosed within the 5σ contour.