Quasar Feedback Survey: molecular gas affected by central outflows and by ∼ 10 kpc radio lobes reveal dual feedback effects in ‘radio quiet’ quasars

We present a study of molecular gas, traced via CO (3–2) from ALMA data, of four z < 0.2, ‘radio quiet’, type 2 quasars (L bol ∼ 10 45 . 3 − 46 . 2 erg s − 1 ; L 1 . 4GHz ∼ 10 23 . 7 − 24 . 3 W Hz − 1 ). Targets were selected to have extended radio lobes ( ≥ 10 kpc), and compact, moderate-power jets (1–10 kpc; P jet ∼ 10 43 . 2 − 43 . 7 erg s − 1 ). All targets show evidence of central molecular outflows, or injected turbulence, within the gas disks (traced via high-velocity wing components in CO emission-line profiles). The inferred velocities (V out =250 – 440 km s − 1 ) and spatial scales (0.6 – 1.6 kpc), are consistent with those of other samples of luminous low-redshift AGN. In two targets, we observe extended molecular gas structures beyond the central disks, containing 9 – 53 % 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: projected distances of 5 – 13 kpc; bulk velocities of 100 –340 km s − 1 ; velocity dispersion of 30 – 130 km s − 1 ; inferred mass outflow rates of 4 – 20 M ⊙ yr − 1 ; and estimated kinetic powers of 10 40 . 3 − 41 . 7 erg s − 1 . 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.


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 disk 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 et al. 2012) and is theoretically capable of affecting the host galaxy through regulation of star-formation (McNamara & Nulsen 2012;Fabian 2012).Depending on how the available energy couples to the interstellar medium (ISM), the gas could be driven due to wideangled accretion disk winds, radiation pressure on dust, and/or due to the acceleration by radio jets (e.g., Sĳacki et al. 2007;Fabian 2012;King & Pounds 2015;Ishibashi & Fabian 2016;Mukherjee et al. 2016;Costa et al. 2018Costa et al. , 2020;;Tanner & Weaver 2022;Almeida et al. 2023).These processes can also influence the fuel available for feeding the black hole itself, thereby giving this process a selflimiting 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 et al. 2005;Nesvadba et al. 2008;Feruglio et al. 2010;Alexander et al. 2010;Harrison et al. 2012;Rupke & Veilleux 2013;Liu et al. 2013;Cicone et al. 2014;Villar Martín et al. 2014;King & Pounds 2015;Fiore et al. 2017;Rupke et al. 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;Ramos Almeida et al. 2022;Girdhar et al. 2022;Kakkad et al. 2022Kakkad 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 bol ∼ 10 45−46 erg s −1 , the observed molecular outflows typically have 100 × more mass than the ionised 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;Mainieri et al. 2011;Alatalo et al. 2011;Cicone et al. 2014;Morganti et al. 2015;Harrison 2017;Fiore et al. 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 disk 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. 2018Costa et al. , 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 1 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. 2017Russell et al. , 2018;;Tremblay et al. 2018;Russell et al. 2019;Olivares 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 disk, and caused by powerful radio jets (e.g., in the BCGs) and one acting within the molecular gas disks, 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 jet ≤ 10 45 erg s −1 ) 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).
QFeedS2 , includes 42 quasars (L bol ≳ 10 45 erg s −1 ) 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.,  ∼ 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 analogs of high redshift AGN.The QFeedS dataset 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. 2019Jarvis et al. , 2020Jarvis et al. , 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 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 molec-ular 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  0 = 70 km s −1 Mpc −1 , Ω  = 0.3 and Ω Λ = 0.7, throughout.In this cosmology, 1 arcsec corresponds to 2.47 kpc for the redshift of  = 0.14 (i.e., the average redshift of the galaxies studied here).When referred to, we define the radio spectral index, , using the relation,   ∝   , where   refers to the flux density at the corresponding frequency .

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]  5007 Å luminosities (L [O III] > 10 42.1 erg s −1 ).A moderate radio luminosity criteria of L 1.4GHz > 10 23.45 W Hz −1 is also applied to obtain the QFeedS sample; however, the sample still consists of 88 % 'radio quiet' sources, based on the criteria of Xu et al. (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 %; Zakamska et al. 2004).
Figure 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-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 deconvolved size of the core.The sample exhibits a wide range of radio sizes (∼0.1 kpc to ∼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).

Sample Selection
Following our goal to analyze the molecular gas, we identified 9/42 targets from the QFeedS sample which have available CO 12-m array ALMA data (highlighted by the empty black circles in Figure 1).For all 9 of these, CO(3-2) data is available in ALMA Band 7 observations; therefore, we decided to use this as our tracer of molecular gas for this work.We note 2/9 also have CO(2-1) data and 1/9 has CO(1-0) data which are presented in Ramos Almeida et al. ( 2022 the disks, which could be associated with extended radio lobes3 , we further only selected the sources with LLS radio ≥10 kpc (from VLA data, see Section 2.3).As shown in Figure 1, 4/9 sources with ALMA data meet this criteria; namely, J0945+1439, J1000+1242, J1010+1413, and J1430+1339.All four targets are bright AGN with high, quasar-like bolometric luminosities of log (L bol /erg s −1 ) = 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.4GHz versus L [OIII] criteria of Xu et al. (1999).However, despite their modest radio luminosities of log(L 1.4GHz /W Hz −1 ) = 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. 2019Jarvis et al. , 2021)).High-resolution VLA data at 1.4 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 −1 ) 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 disk 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  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 disks, (∼1 kpc) and on large scales (≳10 kpc), extended beyond the molecular gas disks, in 'radio quiet' quasars.

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, one-hour 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.795990GHz, 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 minutes.The spectral window has a bandwidth of 1.875 GHz and was centered at the CO (3-2) line, with the same frequency as mentioned above.The resulting angular resolutions of the observations are  res ∼ 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 centered in the peak luminosity pixel of each galaxy, with a radius varying between 5 arcsec and 8 arcsec to make sure all the resolved emission was included.A channel width of 25 km s −1 was selected, and a pixel scale of 0.05" 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.

Summary of the radio images
For our investigation of the relationship between the CO emissionline properties and radio morphology, we use the 6 GHz (C-band) VLA radio images from Jarvis et al. (2019).We use both the 'lowresolution' (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 (for example, using uniform, natural or Briggs pa-rameter=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., Figure 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 optimised 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.

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 and properties listed in Table 5.The purple curves shows the combined fits while the blue, and orange dashed curves represent the individual Gaussian components.In the case of central spectra, we use the orange curves as the "broad" Gaussian component.3. The identified filament regions are highlighted with the dashed-grey boxes and are labeled respectively as F1 -5.We overlay as green contours, within each box, the CO(3-2) emission at the 3 and 5 levels from the associated images presented in Figure A1.The radio contours are the same as in Figure 2 but shown here in blue and white colours for high-and low-resolution, respectively.A legend is shown in the right panel.(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 -270 km s −1 .All values for z * and  * are listed in Table 2.

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 disks (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 emissionline wing components, as a tracer of central molecular outflows.

Identification of extended molecular gas structures
In Figure 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 of each panel).Overlaid on these images are contours of the 6 GHz radio emission for both the low-resolution and high-resolution 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 disks as central, contiguous CO structures with smooth velocity gradients centered 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.
(i) Emission with a clear morphological and/or kinematic separation from the central molecular gas disk.(ii) 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.(iii) Structures extending to ≥ 1 kpc in projected size.
Table 2. Measured global galaxy properties using the spectra shown in Figure 2 (panels a2 -d2) and extracted from the dashed-grey regions (shown in the panels a1 -d1): (1) quasar name; (2-3) stellar redshift and stellar velocity dispersion, respectively, measured from the stellar kinematics following Section 2.4; (4-7) measurements following Section 3.2, namely, (4) mean velocity (V 50 ); (5) velocity width (W 80 ); (6) velocity dispersion estimated as  =  80 /2.56; (7) flux; (8) molecular gas mass; (9) jet kinetic power obtained from the median of the total low-resolution flux density and core flux density at 5.2 GHz (Jarvis et al. 2019).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 Figure 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 5 and 3 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 ≥ 5  was seen associated with these structures.These are shown in Figure A1.A combined overview of the molecular gas structures is shown in Figure 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 Figure A1, with higher weights to the structures with lower surface brightness.Each filament is highlighted with a surrounding dashed grey box, and labeled following the notation as F1 -5, respectively.These boxes cover the full observed extent (at ≥ 3 ) of the structures.In Section 3.3, we evaluate the properties of each of these structures.
For all 8 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 2-D 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 -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 Figure 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 analyzed 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.

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 Figure 2); (3) regions covering each of the individual filaments (shown through dashed-white boxes in Figure 3); and (4) central 'outflow' regions where we identify broad CO emission-line wing components (shown through dashed-purple boxes in Figure 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 −1 .We modeled 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 50 ; and the velocity dispersion is measured through velocity width in terms of W 80 , i.e., the width containing 80% of the emission-line flux.For a single Gaussian, W 80 is approximately related to the full-width-at-half-maxima (FWHM) as W 80 = 1.088 × 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 below.
• For the entire galaxy-scale: For the total CO emission, we extract the spectrum over the region shown through a grey-dashed box in Figure 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 arcsec 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 Figure 4 and Figure 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 Figure 3. Kinematic maps for the individual structures, along with the extracted CO emission-line profiles and fits, are shown in the Appendix Figures B1 and B2.

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.

Molecular gas velocity
To measure the bulk velocities in the molecular gas, we utilize the median velocity V 50 (see Section 3.2) derived from the emissionline 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 panels of Figure 2. The velocity maps for J1000+1242 and J1010+1413 are shown in the left panels of Figure 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 Figures 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-fit model (on the scale of the residual noise in the continuum).A fit was obtained for each of these simulated spectra with V 50 measured each time.The standard deviation of the distribution of the measured V 50 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 fluxweighted 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%), as compared to the the previously described V 50 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.

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 fila-ments,  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  = W 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).

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 Figure 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.

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 CO Δ v) to molecular gas mass (M mol ) using the following relation (Solomon & Vanden Bout 2005;Bolatto et al. 2013): where z is the redshift of the galaxy, D L is the luminosity distance, and X CO is the CO-to-H 2 conversion factor, with X CO,Gal = 2×10 20 cm −2 (K km s −1 ) −1 (Solomon et al. 1987;Solomon & Vanden Bout 2005).This relation is the corollary of the relation M mol =  L ′ CO(1−0) , where X CO and  CO are both referred to as "CO-to-H 2 " conversion factor.For X CO = 2×10 20 cm −2 (K km s −1 ) −1 , the corresponding  CO = 4.3 M ⊙ (K km s −1 pc 2 ) −1 .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 et al. 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).

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 AGN ∼ 10 40−46 erg s −1 .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 −1 ).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 Figure 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 Figure 2).We measured R out as the projected distance between the farthest 0 100 200 300 400 500 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: 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 Figure 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 parameterspace comparison of CO emission-line profile properties presented in Figure 6.

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 (Figure 1).Specifically, in Section 4.1, we summarise 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.

Properties of the extended molecular gas structures
Figure 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 Figure 3, we identify five filaments for J1000+1242 and three for J1010+1413.Figures 4  and 5 show the velocity and velocity-width maps (in terms of W 80 ) 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 Figure 6, as teal-coloured triangles for J1000+1242 and stars for J1010+1413.For J1000+1242, across the 5 filaments, there is a radial extent range of R fil = 5 -11 kpc and velocities of V fil = | 100 -190 | km s −1 .In case of the 3 filaments in J1010+1413, we see a radial extent range of 12 -13 kpc with comparatively higher velocities of 220 -340 km s −1 .This gives us an average radial extent of 8 kpc and 12 kpc; and an average velocity of ∼ 150 km s −1 and ∼ 280 km s −1 for J1000+1242, and J1010+1413, respectively.
In Figure 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 −1 for J1000+1242 and 47 -90 km s −1 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 8 filaments.
In the left panel of Figure 8, we present the fraction of total molecular gas mass located in the filaments, with respect to the radio luminosity (L 1.4GHz ), 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% of the gas located within these structures.J1010+1413 has only 9% 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%, and 9% 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 suf-ficient for a broad comparison to the values for similar structures, obtained using similar datasets, 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 Figure 3; Figure 4; and Figure 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 −1 ), 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 Figure 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 Tremblay et al. 2018;Balmaverde et al. 2018;Russell et al. 2019;Tamhane et al. 2022;Capetti 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 (Figure 8).

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. (2015Hogan et al. ( , 2017) ) and Pulido et al. (2018)  4 .
The BCGs are typically massive galaxies (with  ★ ≈200 -500 km s −1 ; see Figure 7; and M ★ = 10 10.6 − 12.5 M ⊙ ) and 'radio loud' galaxies (see Figure 8; right panel).In contrast, our sample uniquely consists of 'radio quiet' quasars (see Section 2.1; Figure 8) and has comparatively lower stellar masses (i.e., M ★ = 10 9.9 − 11 M ⊙ ; Jarvis et al. 2020).The BCGs, with known strong 'radio-mode' feedback and large reservoirs of molecular gas (10 9 − 10 11 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 (Figure 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 (Figure 6).Further, we see similar velocity dispersion values of ∼ 10 -160 km s −1 , shown in the right panel of Figure 7, across both samples 5 .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 panel of Figure 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% in Hydra-A (disk-dominated) to ∼90% in AS1101 (filamentdominated).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%.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 Figure 8, we look for trends in the observed filamentto-total molecular gas fractions, with respect to the galaxy properties, such as, radio luminosity (L 1.4GHz ), [O iii] luminosity (L [OIII] ), and stellar velocity dispersion ( * ) for the combined sample of BCGs plus our sample.In the left panel of Figure 8, we see no clear trends with 4 When radio flux and [O iii] luminosities were not available in these work, we obtained these values using NASA/IPAC Extragalactic Database (Helou et al. 1991).However, for 5/15 BCGs we didn't obtain either the radio flux or the [O iii] luminosities, and hence these are excluded from the Figure 8 (right panel).For 4/15 BCGs, we could not recover the stellar dispersion values and hence these are represented as empty symbols in Figure 8 (left panel), and discussed later. 5For consistency, only the 10/12 galaxies common between Russell et al. radio luminosity, nor with stellar velocity dispersion (represented by the colour-scaling).In the right panel of Figure 8, we compare our targets to the BCG sample in the L 1.4GHz vs. L [OIII] 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 [OIII] ≥ 10 42.1 erg s −1 ; in the quasar regime) have a typically lower mass-fraction of gas in the filaments (∼ 22%), compared to the higher average mass-fraction (∼ 52%) 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%).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.

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. 2014Harrison et al. , 2015;;Ramos Almeida et al. 2017;Speranza et al. 2022;Venturi et al. 2023), with velocities reaching 1000 km s −1 , 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 highvelocity 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 −1 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, 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 40−46 erg s −1 ; in comparison to our sample that lies in the higher quasar luminosity regime, i.e., 10 45−46 erg s −1 .Further, both the samples lie in the same redshift range, z ∼ 0.2.In Figure 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.The rate of kinetic energy vs. the jet kinetic power (P jet ).The teal symbols show the individual filaments for our targets.The yellow symbols show the total rates over all the filaments for each galaxy and the red circles show the values for the BCG sample (which are averaged over multiple filaments).The black error bar shows the variation in the jet power depending on whether the core or total radio luminosity was used.The dashed olive lines on the right correspond to coupling efficiencies of 5%, 0.5%, and 0.05%.In the bottom-right the average systemic error assumed for each derived value is shown.

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  ≤ 0.2 BCGs with the properties of highvelocity wings observed in CO emission lines (as a tracer of central molecular outflows) in the sample of  ≤ 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 Section 4.1 and Section 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 4 of the sample.In this section, we discuss the implications for the observed feedback effects on multiple scales.

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 Figure 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. 2014McNamara et al. , 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 (Figures 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 shockheated (McNamara et al. 2000;Fabian 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 1-3 orders of magnitude larger than the central outflows found in Fluetsch et al. (2019).In two of our sample, 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 disks.
Again, following Tamhane et al. (2022), we compute the mass flow rate for the filaments ( M 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: E kin = 1/2 M fil V 2 fil .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 Figure 9, we compare the filaments in our targets to BCGs in terms of their mass outflow rates (left panel) and kinetic powers (in the right 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 dex 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 Figure 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 −1 ) are comparable to the BCGs.At the same time, the kinetic power in the filaments (10 40−42 erg s −1 ) 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 % percent of the jet kinetic power (consistent with all but 2 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 % 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.,  kin / bol ∼ 10 −6 i.e., 10 −4 % (see Table 4).
This brings us to the conclusion that Figure 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.

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: a) ultra-fast outflows/small-scale winds (Costa et al. 2014(Costa et al. , 2020)), (b) jets (Talbot et al. 2022a,b), (c) radiation pressure (Costa et al. 2018) and even the more phenomenological AGN feedback models used in state-of-theart 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 et al. 2012;Mukherjee et al. 2016Mukherjee et al. , 2018;;Mandal et al. 2021;Talbot et al. 2022aTalbot et al. , 2023;;Tanner & Weaver 2022;Talbot et al. 2022b) 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 et al. 2023).On small spatial scales (and shorter timescales), 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 jetcocoon which 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 timescales, 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.

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−46.2erg s −1 ; L 1.4 GHz ∼ 10 23.7−24.3W Hz −1 ), 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 Figure 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 Figure 2 for a data overview).Our main findings are summarised below: 1. We identify filamentary molecular gas structures in and around ∼10 kpc radio lobes in two out of the four 'radio quiet' quasars (Figure 3, 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 −1 and velocity dispersion values of  mol,filament = 30 -130 km s −1 .We observe that ∼ 53% and ∼ 9% 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% of the total CO (3-2) emitting gas could be contained in such structures (see Figure 8; left 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 Figure 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 (Figures 6, 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 Figures 6 and 7).
The inferred mass outflow rates (4 -20 M ⊙ yr −1 ) and kinetic powers (10 40.3−41.7 erg s −1 ) are also comparable to those seen in the BCGs (Figure 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 [OIII] ≤ 10 42 erg s −1 ) tend to show higher fractions of the molecular mass in filaments, with an average filament mass-fraction of 55 %, while the higher-[O iii] luminous sources, including our targets, have an average filament mass-fraction of ≤16 % (see Figure 8; right 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 (Figure 4, 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 −1 ) of the CO(3-2) emission-line profiles (see Figure 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 moderatepower 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 et al. 2022a;Morganti et al. 2023) that suggest that two feedback effects can take place due to low-and moderate-power radio jets (P jet ≲ 10 44 erg s −1 ).On the smaller scales, jet-ISM interactions can drive turbulence and central outflows.On larger scales, radio lobes have penetrated beyond galaxy disks 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(Mukherjee et al. , 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.

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 ≥ 5  was seen associated with these structures.These are shown in Figure A1.For J1000+1242, we identified 5 filaments and for J1010+1413, we identified 3 filaments, as shown in each of the panels below.The velocity windows used for collapsing each image are labeled on the top-right of each panel, and each filament is highlighted with a surrounding dashed white box.These boxes cover the full observed extent (at ≥ 3 ) of the structures, where the 3  contours are also shown.A combined overview figure for each of the targets can be seen in Figure 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 Figure B1 and Figure B2 for each of the filaments of J1000+1242 and J1010+1413, respectively.From left to right, each of the panels in the figure represents

)Figure 1 .
Figure 1.Largest linear size of radio structures (LLS radio ) versus [O iii] luminosity for the parent sample of 42 QFeedS targets, colour-coded by their 1.4 GHz radio luminosity (Jarvis et al. 2021).The 9/42 QFeedS targets with the required CO(3-2) ALMA data are highlighted with a black circle.A selection criteria of LLS radio ≥ 10 kpc was then applied (dashed line) to select the four targets for this work (labeled with their names, see Section 2.1).

Figure 2 .
Figure 2. The Panels a1, b1, c1, d1: show CO (3-2) flux (moment 0) maps for the four targets selected following Section 2.1.The CO(3-2) images are created by integrating over the velocity ranges annotated in the top-right of each panel.The green contours show CO emission corresponding to levels indicated in the individual legends, and the dashed green contours show the −3  level of CO emission.The red and the white contours show the radio emission from the high-and low-resolution 6 GHz images, respectively.The contour levels are mentioned in the individual panel legends and are chosen to highlight the important structures in each target.Zoom-ins of the central emission are shown at the top for J1430+1339 and J0945+1737 (where the e-MERLIN image at 1.5 GHz, is used and shown through blue contours).A 5 kpc scale bar is shown in each panel at the bottom-right and the ALMA beams are shown at the bottom-left of the panel.The dashed-grey boxes show the region over which the galaxy-integrated spectra are extracted and are shown in the Panels a2, b2, c2, d2: with the properties listed in Table 2.The dashed purple boxes show the central outflow regions with the spectra shown in Panels a3, b3, c3, d3:

Figure 3 .
Figure 3. CO(3-2) emission-line maps, produced using a weighted combination of the individual narrow velocity-range images, highlighting each filament for J1000+1242 and J1010+1413, as shown in Figure A1 (see Section 3.1).The individual velocity ranges over which the filaments were detected are listed in Table3.The identified filament regions are highlighted with the dashed-grey boxes and are labeled respectively as F1 -5.We overlay as green contours, within each box, the CO(3-2) emission at the 3 and 5 levels from the associated images presented in FigureA1.The radio contours are the same as in Figure2but shown here in blue and white colours for high-and low-resolution, respectively.A legend is shown in the right panel.

Figure 4 .
Figure 4. Kinematic maps of the CO(3-2) emission line for J1000+1242 (see Section 3.2).The left panel is a velocity map (V 50 ), the right panel is a map of velocity width (W 80 ).Overlaid on each map, the black and olive-green contours correspond to the high-and low-resolution 6 GHz radio images, same as in Figure2.The black, dashed boxes highlight the filamentary regions identified in Figure3, and labeled respectively as F1 -5.A legend is shown in the left panel and a 5 kpc scale bar in the right panel.

Figure 5 .
Figure 5. Kinematic maps of the CO(3-2) emission line for J1010+1413 (see Section 3.2).The left panel is a velocity map (V 50 ), the right panel is a map of velocity width (W 80 ).The black, dashed boxes highlight the filamentary regions identified in Figure 3, and labeled respectively as F1 -3.The rest is the same as in Figure 4.

Figure 8 .
Figure8.Left: Ratio of molecular gas in filaments to the total molecular gas as a function of 1.4 GHz radio luminosity of the respective target for both the samples: symbols as in the legend, and the 14 BCGs fromTamhane et al. 2022, represented as circles.Hydra A is also shown which has no identified filaments, through a plus symbol(Rose et al. 2019).The data points are colour-coded by stellar velocity dispersion of the galaxies, except where these data were not found (shown as empty circles).Right: 1.4 GHz radio luminosity versus [O iii] luminosity for the targets as in the left panel, but for the entire galaxy (excluding 4 BCGs for which L[OIII]  was not available), and colour-coded by the ratio of molecular mass found in filaments to the total.The dot-dashed line separates 'radio loud' and 'radio quiet' sources (followingXu et al. 1999) and the vertical dashed line is the 'quasar' luminosity threshold used to select our sample (i.e., L [OIII] ≥ 10 42.1 erg s −1 ; see Section 2.1).
(2019) andTamhane et al. (2022) have been shown in this plot.The two remaining BCGs fromRussell et al. (2019) that are not shown are A262 and A2052, which do not affect the overall scientific interpretation in this context.

Figure 9 .
Figure 9. Left panel: Rate of Molecular Mass Outflow vs. the Jet Kinetic Power (P jet ) for our filaments and the 14 comparison BCG.Right panel:The rate of kinetic energy vs. the jet kinetic power (P jet ).The teal symbols show the individual filaments for our targets.The yellow symbols show the total rates over all the filaments for each galaxy and the red circles show the values for the BCG sample (which are averaged over multiple filaments).The black error bar shows the variation in the jet power depending on whether the core or total radio luminosity was used.The dashed olive lines on the right correspond to coupling efficiencies of 5%, 0.5%, and 0.05%.In the bottom-right the average systemic error assumed for each derived value is shown.
(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 80 ); 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 Figure 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.

Figure A1 .
Figure A1.CO(3-2) emission-line maps, integrated over specific velocity ranges (quoted in the top-right of each panel), to highlight the identified filamentary structures (see Section 2.1).The filament regions are highlighted with dashed-white boxes, and labeled respectively as F1 -5 for J1000+1242 and as F1 -3, for J1010+1413.The green contours are the CO(3-2) emission inside these boxes, at the 3 and 5 level.The radio contours are the same as in Figure 2 just shown here in blue and white colours for high-and low-resolution at 6 GHz respectively.A legend is shown on the bottom-right panel.