The Quasar Feedback Survey: characterising CO excitation in quasar host galaxies

We present a comprehensive study of the molecular gas properties of 17 Type 2 quasars at $z<$ 0.2 from the Quasar Feedback Survey (L$_{[OIII]}$>$10^{42.1}$ $\rm ergs^{-1}$), selected by their high [OIII] luminosities and displaying a large diversity of radio jet properties, but dominated by LIRG-like galaxies. With these data, we are able to investigate the impact of AGN and AGN feedback mechanisms on the global molecular interstellar medium. Using APEX and ALMA ACA observations, we measure the total molecular gas content using the CO(1-0) emission and homogeneously sample the CO spectral line energy distributions (SLEDs), observing CO transitions (J$_{up}$ = 1, 2, 3, 6, 7). We observe high $r_{21}$ ratios (r$_{21}$ = L'$_{CO(2-1)}$/L'$_{CO(1-0)}$) with a median $r_{21}$ = 1.06, similar to local (U)LIRGs (with $r_{21}$ $\sim$ 1) and higher than normal star-forming galaxies (with r$_{21}$ $\sim$ 0.65). Despite the high $r_{21}$ values, for the 7 targets with the required data we find low excitation in CO(6-5)&CO(7-6) ($r_{61}$ and $r_{62}$<0.6 in all but one target), unlike high redshift quasars in the literature, which are far more luminous and show higher line ratios. The ionised gas traced by [OIII] exhibit systematically higher velocities than the molecular gas traced by CO. We conclude that any effects of quasar feedback (e.g. via outflows and radio jets) do not have a significant instantaneous impact on the global molecular gas content and excitation and we suggest that it only occurs on more localised scales.


INTRODUCTION
A fundamental outstanding question of galaxy evolution is what impact active galactic nuclei (AGN) have on the interstellar medium (ISM) and star formation in their host galaxies.AGN can release energy into their host galaxies via processes known as AGN feedback, which are required by our current models of galaxy evolution to regulate star formation (Bower et al. 2006;Schaye et al. 2015), and believed to be the mechanism that regulates the co-evolution of accreting black holes (BH) and their host galaxies that is observed across cosmic time (e.g.Madau & Dickinson 2014;Kormendy & Ho 2013;Cresci & Maiolino 2018).Observational and theoretical studies have proposed both a suppression (Silk & Rees 1998;Hopkins et al. 2006;Booth & Schaye 2010;Feruglio et al. 2010; Cicone ★ Email: smolyneux.astro@gmail.comet al. 2014; King & Pounds 2015;Fiore et al. 2017;Costa et al. 2018;Ellison et al. 2021;Bertemes et al. 2023) and an enhancement of star formation in AGN host galaxies (Ishibashi & Fabian 2012;Silk 2013;Zubovas et al. 2013;Lacy et al. 2017;Fragile et al. 2017;Gallagher et al. 2019).However, these interactions are still not fully understood and the diversity of results stress the need for multiwavelength, multi-tracer studies to characterise the interplay between the central supermassive black hole and the host galaxy.
A natural assumption may be that the most powerful and luminous AGN and quasars will have the largest impact on their host galaxy.Studies suggest that they might be able to drive kpc-scale outflows across the entire galaxy, expelling the interstellar star-forming gas (e.g.Cicone et al. 2012;Harrison et al. 2014;Feruglio et al. 2015;Circosta et al. 2018;Longinotti et al. 2023).However, studies at low redshift () have also shown that AGN and quasars tend to reside in gas rich, star forming galaxies and find no instantaneous depletion of total gas content (Saintonge et al. 2017;Shangguan et al. 2020;Jarvis et al. 2020;Koss et al. 2021).These findings also agree with recent simulations (e.g.Piotrowska et al. 2022;Ward et al. 2022), therefore supporting the idea that large gas reservoirs are needed to fuel the accreting supermassive BHs in quasar and AGN hosts.It may therefore be the case that any impact from feedback is limited to a more localised scale and the global properties of the ISM are left largely unaffected.Indeed, there are works which show the possible impact of AGN feedback on the molecular gas content in central/localized part of the galaxy (e.g.Rosario et al. 2019;Feruglio et al. 2010;Ellison et al. 2021;Ramos Almeida et al. 2022;Audibert et al. 2023).
The molecular phase of the ISM, commonly traced by observing low transitions of carbon monoxide (CO), plays a critical role in galaxy evolution as it is this gas which is redistributed to both promote star-formation activity and fuel BH growth (e.g.McKee & Ostriker 2007;Carilli & Walter 2013;Vito et al. 2014;Tacconi et al. 2020).However, no consensus has yet been reached on the impact of AGN on the overall molecular gas content in the ISM (Kakkad et al. 2017;Perna et al. 2018;Kirkpatrick et al. 2019;Rosario et al. 2019;Circosta et al. 2021;Morganti et al. 2021).One possible reason for this might be due to the time scales over which any impact may take place (King et al. 2011;Mukherjee et al. 2018b;Zubovas 2018;Ward et al. 2022).There are also complexities due to the resolution of observations, biases in the sample selection, what tracers of the gas are used, and the uniformity in the observations.
While most studies of the AGN impact on molecular gas have focused on the total gas content, much is still unknown about other molecular gas properties such as molecular gas excitation.Knowledge of the ground state CO(1-0) line is a crucial reference that is often used to not only compare to higher transitions and measure the excitation, but to also convert to the total molecular gas content of the galaxy.However, there is discussion in the community about how reliable the ground state is in doing these calculations for different objects (e.g.star forming galaxies (SFGs) or (ultra) luminous infrared galaxies (U)LIRGs, see Leroy et al. 2022a;Montoya Arroyave et al. 2023).This therefore increases the importance in characterising the CO(1-0) across different samples.
Due to the expensive observations required to detect multiple emission lines for individual sources, most of our knowledge on molecular gas excitation is based on inhomogeneous coverage of few transitions, and limited to for the most luminous galaxies (e.g.Kakkad et al. 2017;Saintonge et al. 2017;Lamperti et al. 2020;Circosta et al. 2021;Boogaard et al. 2021;Valentino et al. 2021;Harrington et al. 2021;Leroy et al. 2022a).Further, studies have investigated the driving mechanism for the excitation of the molecular gas (e.g.Daddi et al. 2015;Pozzi et al. 2017;Mingozzi et al. 2018;Leroy et al. 2022b;Esposito et al. 2022), suggesting photodissociation regions (PDRs) and X-ray dominated regions (XDRs), of diverse temperature and gas densities, are the key physical components driving CO excitation.
Models of CO excitation suggest that AGN-related processes, such as X-ray emission (Meĳerink et al. 2007) and shock heating induced by AGN jets and outflows (Kamenetzky et al. 2016), would mainly affect the molecular gas excitation at the higher CO transitions.It is therefore crucial to study both low and high CO transitions as the impact of feedback may only be present at higher CO excitations (   > 5), whereas the bulk of the molecular gas content is still traced by the ground transition.However, this is challenging at low redshifts ( < 0.2), where even using the maximum frequency limit of the ALMA bands, only    ≤ 8 can be reached.Due to the observational difficulty of observing at higher frequencies, there are few examples of these critical higher CO transitions observed at low- (e.g.van der Werf et al. 2010;Greve et al. 2014;Rosenberg et al. 2015;Liu et al. 2015;Kamenetzky et al. 2016;Yang et al. 2017).Indeed, most    > 7 observations come from higher redshift ( > 1), highly luminous quasars, which are far more easily observed (Carilli & Walter 2013;Wang et al. 2019;Yang et al. 2019;Li et al. 2020;Pensabene et al. 2021;Decarli et al. 2022), which in turn lack observations of low-J transitions, consequently lacking a complete characterisation of the molecular gas content and excitation.
Previous studies have suggested a close relationship between AGN feedback diagnostics and the properties of the ionised phase of the ISM.For example, a study of optically selected AGN from SDSS found that those with higher radio luminosities were more likely to have larger FWHM [O III] (Mullaney et al. 2013), suggesting a relation between the radio emission and the kinematics of the ionised gas.Further work on the same sample showed that the most extreme ionised outflows (FWHM [O III] > 1000km s −1 ) were found to be more common when the radio emission was compact (Molyneux et al. 2019).With high resolution radio observations for a sample of 42 of these targets (presented in Jarvis et al. 2019) from the Karl G. Jansky Very Large Array (VLA), a prevalence of small scale radio jets (in the central few kpc) was found, leading to a suggestion that they could be the driver of these ionised outflows.Alternatively star formation driven outflows and quasar winds that shock the ISM may be responsible for producing the observed radio emission and correlation with outflow properties (e.g.Condon et al. 2013;Nims et al. 2015;Zakamska et al. 2016;Hwang et al. 2018;Panessa et al. 2019).
With the advent of deeper radio data, increasing evidence is being found of potentially ubiquitous low-level radio emission in radioquiet quasars (e.g.Mukherjee et al. 2018b;Jarvis et al. 2019Jarvis et al. , 2021;;Macfarlane et al. 2021).These observations suggest that radio jets are potentially an important feedback mechanism in radio quiet quasars.Indeed radio jets have been found to have an impact on the surrounding multi-phase ISM (e.g.Morganti et al. 2015;Oosterloo et al. 2017;Jarvis et al. 2019;Morganti et al. 2021;Girdhar et al. 2022).However, an outstanding question is how and when these jets can couple to the ISM, and have a positive and/or negative impact on the starforming molecular gas content (e.g.Silk 2013;Gabor & Bournaud 2014;Bieri et al. 2016;Costa et al. 2018).Studying the CO excitation of quasars with known outflows/jets is therefore key to solving these outstanding questions.
The Quasar Feedback Survey (QFeedS) is a multi-wavelength survey aiming to address these open questions in order to understand the co-evolution between quasars and their host galaxy, in particular in the context of ionised outflows and radio jets.These are luminous systems at  < 0.2 and so it is possible to study the impact that feedback (e.g. via radio jets) has on the multi-phase ISM on both resolved and global scales, whether it be driving outflows, disturbing the gas kinematics, affecting the molecular gas excitation, or impacting on star formation (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).
In this work we present a comprehensive study of the molecular gas properties of 17 quasars of the QFeedS sample which have multiwavelength data.We characterise the molecular excitation in these sources, presenting the CO(1-0), CO(2-1) and CO(3-2) for the entire sample and also  or CO(7-6) for 7 of the 17 targets.With the addition of ancillary multi-wavelength data, we will explore the impact of feedback, if any, on the total molecular gas content and molecular gas excitation within the quasar host galaxies.
In Section 2 we introduce the quasar sample presented in this work as part of the Quasar Feedback Survey.In Section 3 we describe the observations used and the data reduction.In Section 4 we describe the analysis techniques used to study our CO data, including spectral fitting, definition of detection, flux measurements and line profile characterisation.We also introduce the comparison samples from the literature that we utilise in our analysis.In Section 5 we present our results of the CO excitation, line profile properties and gas fractions, and at all times comparing to relevant samples from the literature.We then discuss our findings in the overall context of galaxy evolution and quasar feedback.Our final conclusions are presented in Section 6.

SAMPLE SELECTION
QFeedS, presented in Jarvis et al. (2021), is a multi-wavelength study of 42 quasars at  < 0.2.This main sample was selected from a parent sample of 24264 optically selected AGN from SDSS at  < 0.4 from Mullaney et al. (2013).These 42 quasars were selected to have L [O III] > 10 42.1 ergs s −1 and to cover the full range of FWHM Avg[O III] (a flux weighted average of the FWHM of the two Gaussian components present in the spectra) with velocities in the range = 339 -1289 km s −1 (see Figure 1).
Here we introduce a study as part of QFeedS to provide a detailed characterisation of molecular gas in 17 Type 2 quasars, studying properties such as molecular gas masses, gas fractions and CO excitation.The 17 targets were selected to be Type 2 quasars which are visible from the Very Large Telescope (VLT) and the Atacama Pathfinder EXperiment telescope (APEX) and that are representative of the parent population (see Figure 1).We selected Type 2 quasars in order to achieve a more robust characterisation of the host-galaxy stellar-emission properties (Jarvis et al. 2019).These 17 targets also have available optical (Multi Unit Spectroscopic Explorer, MUSE) and radio (VLA) data allowing us to perform a full multi-wavelength analysis of the quasar and host continuum emission, in addition to a multi-tracer characterisation of the ISM (these ancillary data are discussed further in Section 3.5).These 17 sources are representative of the survey sample as they cover the full range of QFeedS redshifts ( ∼ 0.1 -0.2) as well as [O III] and radio luminosities, L [O III] = 10 42.1 -10 43.2 erg s −1 and L 1.4GHz = 10 23.5 -10 24.4 W Hz −1 .
Based on the criteria of Xu et al. (1999) using the [O III] and radio luminosity division all 17 of our sample are defined as 'radio-quiet' (see also Jarvis et al. 2021).From previous work (Jarvis et al. 2019), we also know that at least 8 of these 17 sample galaxies are consistent with being luminous infrared galaxies (LIRGs, 10 11 L ⊙ ≲  IR,SF ≲ 10 12 L ⊙ , where  IR,SF is the far infrared luminosity associated with star formation).As only 9 of these targets have the required  IR,SF data then the number of sources consistent with being LIRGs is likely to be higher.This is an important consideration for when we make comparisons with samples in the literature.
The source selection for this study is shown in Figure 1.The colours for the 17 targets in this sample (shown in the legend of Figure 1) are used in all further figures in this work.Further, basic properties of these sources can be found in Table 1.

OBSERVATIONS AND DATA REDUCTION
We use APEX and the Atacama Compact Array (ACA) to observe the carbon monoxide (CO) emission in the CO(1-0), CO(2-1), CO(3-2), CO(6-5) and CO(7-6) transitions for our sample of 17 Type 2 quasars (as detailed in Table A1 in the supplementary material and emission line properties are provided in Tables A2, A3, A4, A5 and A6 in the supplementary material.APEX is a single dish, 12 metre diameter telescope whereas the ACA is a subset of the Atacama Large Millimeter/submillimeter Array (ALMA) comprising of twelve 7 metre antennae.A description of all the observations used in this paper along with details of how the data was reduced is provided below.
The required sensitivity for the ACA observations were estimated based on two different approaches.In the case of sources with archival infrared data around the dust SED peak, a conversion was made from total IR (L IR ) from SED-fitting to CO luminosities L' CO .Otherwise, conversions were estimated based on the SED-inferred stellar mass and the average gas fraction value.
We image the CO(1-0) emission using the TCLEAN function in CASA and apply natural weighting with the Högbom deconvolver.Bin widths of 100 km s −1 were used for non-detections and 50 km s −1 bins were used if the S/N was high enough to see more structure in the line profile.In a few specific cases, slightly different bin sizes were used either to match to other available data or as a result of the data quality.
The beam size of the ACA observations ranged between 12 -14 arcsec.However, to obtain the CO(1-0) spectra we take an aperture equivalent to the APEX beam size when observing CO(2-1), which is ∼ 30 arcsec diameter at an observing frequency of ∼ 200 GHz (observation frequency of CO(2-1) at the samples median redshift of 0.14).Using this aperture consistently to extract the spectra allowed us to compare the fluxes obtained from the same regions, making calculations of line ratios and other properties more reliable.It further allowed us to investigate whether any extended diffuse gas was present, or at least detectable, when comparing to smaller apertures.The apertures used for the extraction and the contours of the CO(1-0) ACA data are shown in the Appendix (Figure C1 in the supplementary material), plotted over rgb images from the DESI Legacy Imaging Survey in the (z, r, g) bands.The 2 sigma contours in the highest S/N data show an extent of up to 27 arcsec (in J1108+0659).Further, those with low S/N (such as the cases of J0945+1737 and J1055+1102) show positional offsets extending out to the 30 arcsecond diameter aperture and slightly beyond.As positional uncertainty in the observations is proportional to beam size S/N , using the 30 arcsecond aperture also allows us to account for these potential offsets and ensure we accurately measure the fluxes.From the CO(1-0) data there are no signs of any companions that are spatially and spectrally aligned with our targets, such that they would impact upon the measured flux values, aside from the apparent mergers occurring in J1222-0007 and J1518+1403 (which are both treated as single systems in this work).Companions that are visible in the background rgb images do not appear in the ACA data and so either are not emitting at those frequencies, or our observations are not deep enough to observe the emission from them.Therefore we can be confident of the fluxes measured in our ACA observations and that these also represent the total CO fluxes in these galaxies.Since this is the only aperture we have control over for the flux/spectra extraction, we choose to match this to the APEX CO(2-1) aperture of 30 arcsec to be as consistent as we can be in the region we are calculating fluxes.
Recent work in the literature has shown evidence for extended, low surface brightness emission in quasars, with CO emission detected out to 100s kpc.(e.g.Cicone et al. 2021;Li et al. 2021;Scholtz et al. 2023).This provides further support to the approach taken here, where we extract spectra with an aperture diameter of ∼ 80 kpc at the median redshift.
To determine whether we obtain the total flux we plot the curves of growth of the ACA CO(1-0) (see Figure 2) where we indeed see that extracting the spectra at 30 arcsec is required to obtain a more accurate total flux value.Beyond 30 arcsec the flux density flattens off in almost all cases (note that in this figure we only plot those with an integrated signal-to-noise ratio (S/N) greater than 5).For a few sources we note that the curves of growth do continue to rise slightly after 30 arcsec, but by < 10% and within uncertainties.Given the larger uncertainties, we are still confident that we are consistent with obtaining the total flux.We also find that the curves of growth up to 30 arcsec follow the same trend and are consistent with each other so no conclusions can be drawn about any differences in morphology at these scales, with respect to galactic or feedback properties.
We note that if the CO(1-0) spectra were extracted using a 3 minimum level, we would measure on average 60% less flux when compared to the extraction at 30 arcsec, and in one case almost 90% less flux (values range from 25 -89%).Such differences would have a significant impact on the analysis of the excitation, stressing the importance of low-resolution data for a complete census of molecular gas content.
From available multi-wavelength data (Jarvis et al. 2021) we also note that J1518+1403 has a secondary source located 14.7 arcsec away to the north east and J1222-0007 has a secondary source 4.7 arcsec away which are both likely to be on-going mergers.Further, J1108+0659 and J1356+1026 show evidence of hosting two nuclei.In all these cases the flux from secondary sources is likely to be included in the flux calculations for each.However, since they will also be covered by the CO(2-1) and CO(3-2) APEX observations then this is unavoidable in the analysis.No other sources have known companions that would affect the total fluxes measured.

CO(2-1) Observations
The CO(2-1) emission of 8/17 sources of the sample were observed with SEPIA180 on APEX between 09 December 2020 and 21 June 2021, (proposal ID: E-0105.B-0713A-2020 [PI: Calistro-Rivera]).The remaining 9 targets had equivalent APEX archival CO(2-1) observations presented as a pilot sample by Jarvis et al. (2020) (Proposal ID: E-0100.B-0166, [PI: Jarvis]).For this work we have re-analysed the raw archival data from Jarvis et al. (2020) using the new analysis techniques presented here for consistency, however we note that we find the same results reported by Jarvis et al. (2020) within uncertainties.
The data was reduced using the standard procedures in the Continuum and Line Analysis Single-dish Software (CLASS) (Pety 2005).In all cases the reduction of the APEX data was done using a consistent strategy by modifying the template reduction script from APEX in CLASS.As with the CO(1-0) data we binned to 50 and 100 km s −1 where appropriate.The observing frequency range of 192.3 -212.5 GHz yields an APEX beam size in the range of 29 -32 arcsec, corresponding to a physical size of 68 -75 kpc at the median redshift of 0.14.
We fit Gaussians using standard procedures in python to obtain the integrated flux values and corresponding uncertainties for each target.We note that the results of fitting Gaussians in python match that of the Gaussian fits produced in CLASS.There is no spatial information for these data, however the beam size is large enough to cover the host galaxy, and so we are confident we measure the total flux, including any diffuse gas and do not over resolve.
Upper limits are calculated using line widths estimated from other CO transitions for the same targets where available.If none are available, the average CO  80 in that transition for all targets is used, for which these values are 397 km s −1 for CO(1-0), 477 km s −1 for CO(2-1) and 507 km s −1 for CO(3-2).We note that this strategy is different from the calculations used by Jarvis et al. (2020), where upper limits were calculated by taking the maximum CO line width from the sample in CO(2-1).The case by case approach based on information from the other CO lines for each target, and assuming we should see similar line widths between transitions, should therefore provide a more accurate and constraining upper limit.The two cases for which both methods have been applied are J0958+1439 and J1356+1026.For J0958+1439, the estimated CO(2-1) upper limit estimated with our method is 30 per cent of the value reported in Jarvis et al. (2020), whereas for J1356+1026 our estimate is 93 per cent of the value reported in Jarvis et al. (2020).
An important note to make is that during the period in which CO(2-1) and CO(3-2) observations were taken with APEX, the telescope was operating at significantly different efficiencies (at most by a factor of 40%).To account for this and to achieve accurate flux measurements, corrections have been made based on the following.Main beam characteristics have been determined from de-convolved continuum slews across Mars, Uranus and Jupiter.Using CO(3-2), this yielded a mean beam size   = 17.5 ± 0.2 ′′ , which we confirmed to be consistent with data on the CO line pointing sources (which are standard AGB stars used for pointing and focus calibration during observations).To determine the main beam efficiency, we used cross scans 1 obtained between December 2020 and December 2021 and cross-checked the result against the CO(3-2) flux of 8 line intensity monitoring sources2 .This analysis yielded main beam efficiencies (at 345 GHz) that depended on periods   = 0.63 ± 0.04 (Dec 2020), 0.55 ± 0.04 (May -Jun 2020), and 0.67 ± 0.04 (Aug -Dec 2020) and an antenna gain factors of Jy/K = 45 ± 4, 51 ± 4, 42 ± 4, respectively, which were converted to the science frequencies using the Ruze formula.In the observations we used the wobbler in symmetrical mode with an amplitude of 50 ′′ and frequency of 0.5 Hz.Pointing and focus were checked regularly against sources from the APEX line pointing catalog using the CO(3-2) emission line.We estimate the overall calibration uncertainty at 10% and that the pointing accuracy was typically within 2 ′′ .Baselines were stable and we only had to fit a first order baseline to each scan before averaging them.

CO(3-2) Observations
Of the 17 sources in our sample, 16 were observed in CO(3-2) with SEPIA345 between 09 December 2020 and 30 December 2021 (proposal ID: E-0105.B-0713B-2020 [PI: Calistro-Rivera]).The observing frequency range of 288.4 -318.7 GHz yields an APEX beam size in the range of 19 -22 arcsec, corresponding to a physical size of 44 -51 kpc at the median redshift of 0.14.These data were reduced, and the flux densities were extracted in the same way as CO(2-1), as described in Section 3.2.
For the one remaining target, J1430+1339, we utilise archival CO(3-2) ACA observations (proposal ID: 2016.1.01535.S [PI: Lansbury]) taken on 03 November 2016.These data were available on the ALMA archive with bin widths of 27 km s −1 .The coverage in the velocity space is not as wide as that of the equivalent APEX observations.Caution should be taken in this case as the difference in spatial resolution (here a beam size of 4.3 arcsec and a maximal recoverable scale of 23 arcsec) means that there is a possibility it is slightly over resolved and perhaps missing flux, but even with these caveats it provides a useful data point.

CO(6-5) and CO(7-6) Observations
We observed 7 targets in either CO(6-5) or CO(7-6) using SEPIA660 on APEX, which were selected for observation based on their brightness in the lower CO transitions.CO(6-5) and CO(7-6) were chosen to give an indication of the excitation in these higher transitions, with CO(6-5) preferred, but if it was not observable (due to the frequency range offered by SEPIA660) then we chose CO(7-6) instead.We also note that [C I](2-1) are also covered by our CO(7-6) observations, however, since all three have non-detections and there are no signs of detection across the obtained spectra, we do not perform any further analysis.
From the range of frequencies 613 -708 GHz, the corresponding beam size was 9 -10 arcsec, which relates to a physical size of 21 -23 kpc at the median redshift of 0.14.This beam size should still allow us to retrieve the full flux values for two reasons.Firstly, for targets in this sample that are observed at higher spatial resolution (0.2 arcsec) and presented by Ramos Almeida et al. (2022), the moment maps and position velocity diagrams show that the CO emission is confined within the APEX beam size.Furthermore, we would expect the CO(6-5) and CO(7-6) to be more compact than the emission of the lower transitions and we are therefore confident that we are measuring the total flux in these data.One caveat would be that if any extended, diffuse emission exists in these higher CO transitions we would potentially be resolving out some of the flux, but we consider this unlikely due to the reasons outlined above.

Ancillary Multiwavelength Data
To achieve a detailed characterisation of the AGN feedback processes in our sample, the 17 quasars studied in this work have ancillary radio and optical data from the VLA and MUSE on the VLT, respectively.Here we describe these data and all values used in this paper can be found in Table 1.
VLA radio data are available for all 17 quasars in the sample at 1-6 GHz and at a resolution of 0.3 -1 arcsec.For a full review and analysis of the radio data see Jarvis et al. (2019Jarvis et al. ( , 2021)).In this work we utilise knowledge of the 1.4 GHz radio data to aid in the interpretation of our findings.The quasars show a range of moderate radio luminosities of log(L 1.4 GHz /WHz −1 ) = 23.5 -24.4.Crucially, although the quasars in our sample are 'radio-quiet', according to widely-used radio-loudness definitions (Xu et al. 1999), the bulk of them exhibit extended radio structures on 0.2 -34 kpc scales, with evidence of jets and/or shocked winds being the dominant cause of the extended radio structures (Jarvis et al. 2019(Jarvis et al. , 2021)).
There is also strong evidence for these radio jets interacting with the ionised gas and driving outflows (Jarvis et al. 2019;Girdhar et al. 2022).We therefore know that within this sample we are observing a diversity of radio AGN emission and ionised outflows in quasars.Further, from these ancillary data we know that this sample is dominated by AGN which are driving outflows, host jets/winds and show interactions with the ISM.The question remains however, as to how these feedback mechanisms impact the molecular gas properties, which we aim to address in this work.
We have also obtained MUSE VLT observations for these 17 quasars (proposal ID: 0103.B-0071 [PI : Harrison]).In this work we use the MUSE data to extract spectra of the [O III]5007 emission line where possible, and otherwise [O III]4959 (1 case) or H (2 cases) if no [O III]5007 line were available.We use these lines as tracers of the ionised gas kinematics and to compare to the molecular CO gas presented here.The spectra were all extracted using the same aperture (diameter ∼ 30 arcsec) as our APEX data (details in Section 3) to make a comparison of ionised gas on the same scales.Specifically we use the [O III]  80 and the properties of the line profile to analyse the differences between the impact of feedback on the ionised and molecular gas properties.

ANALYSIS AND RESULTS
In this section we present the main analysis and results of this work.Firstly, in Section 4.1 we introduce the comparison samples that are used to put our results into context of the overall population of both AGN and non-AGN.We then present the analysis undertaken of the observed CO transitions in Section B. We present the calculations and results of the molecular gas masses and gas fractions in Section 4.3.Finally, we present our findings on the CO excitation via the use of CO Spectral Line Energy Distributions (CO SLEDs) and CO line ratios (Section 4.4).Further, an example of the spectra obtained can be found in Figure 3 and the remaining spectra, alongside tables of the line properties, are contained within the supplementary material.

Comparison samples
Throughout Section 4 we present comparison samples from the literature to put our work into context and aid in the interpretation of our analysis.These comparison samples are described below: We first utilise non-AGN and AGN from Tacconi et al. (2018) to put the gas fractions of our sources in context and show that they are consistent with both AGN and non-AGN.The comparison sample is a compilation of data from xCOLD GASS (Saintonge et al. 2017), EGNOG (Bauermeister et al. 2013) and GOALS (Armus et al. 2009) surveys as well as from the sample presented in Combes et al. (2011).We matched this sample to be within  ± 0.05 of the full range of redshifts spanned by our sample.AGN hosts for the sample were identified using BPT-based AGN classifications.The galaxies in this comparison sample also span the full range of stellar mass, sSFR, and ΔMS found for our sample (see Figure 6 in Jarvis et al. 2020) meaning that the dependency on the specific star-formation rate has been removed and we are focusing on any possible impact of having an active BH rather than the star-formation efficiency of the given galaxy.To ensure consistency in the comparison, the molecular gas masses presented for our sample are calculated using the same method as shown in Tacconi et al. (2018).This comparison is presented in Figure 4.For further information on this comparison sample also see Jarvis et al. (2020).
In our CO SLEDs (Section 4.4.1)we utilise the compilations by Valentino et al. (2021) and Carilli & Walter (2013) to compare to our CO SLEDs.From Valentino et al. (2020Valentino et al. ( , 2021) ) we present the CO SLEDs (both Figures 5a and 5b) of starburst and main sequence galaxies in the redshift range  ∼ 1 -2 with L   > 10 12 ⊙.These luminosities are similar to those in our sample (see Section 3.5) and therefore provide a useful comparison.From Carilli & Walter (2013) we utilise the compilation of high- quasars ( ∼ 1 -6), shown as "C&W13" in Figures 5a and 5b, to see how our low-z quasar sample compare to these more distant and more luminous objects (L bol > 10 47 erg −1 , compared to QFeedS with L bol < 10 46.5 erg −1 ).We also show these high- quasars in Figure 6 as individual points to show how they compare in the lower transitions to the range of line ratios in our sample and other comparison samples listed below.Finally, in Figure 7 we present the line ratios of our high J CO transitions compared to the high- quasar sample, as a function of bolometric luminosities.
As mentioned above, in Figure 6 we present the range of line ratios found in our sample compared to others in the literature in the form of violin plots.Montoya Arroyave et al. ( 2023) analysed a sample 40 local (U)LIRGs ( IR,SF ≳ 10 12 L ⊙ ) in the same redshift range as our sample ( < 0.2).The targets were selected based on OH absorption and not on the presence of radio jets, however this does not exclude radio jets being present.This sample also shows a range of AGN fractions, from 0 -0.92, with 50 % having an AGN fraction greater than 0.5.They find no correlation between AGN fraction or AGN luminosity within the sample of (U)LIRGs.Since 8 sources out of 9 in the QFeedS sample with the required measurements are known to be LIRGs, Montoya Arroyave et al. ( 2023) provides a useful comparison to determine whether the presence of radio jets or shocked winds and ionised outflows found in our sample makes a significant difference to the observed line ratios.
In Figure 6 we also perform a similar comparison to a sample of (U)LIRGs at  ≤ 0.1 from Greve et al. (2014).This sample of (U)LIRGs was selected against AGN, all with an AGN contribution of ≲ 0.3.With this, alongside the sample in Montoya Arroyave et al. ( 2023), we have comparisons samples with similar IR luminosities and a range of AGN contribution.Since the QFeedS sample is compiled of quasars with LIRG-like infrared luminosities, but with additional known radio jets and ionised outflows, any differences in the excitation of CO could potentially be attributed to the jet and outflow properties of our sample.
In order to also investigate how local AGN (with a median  ∼ 0.05) with lower luminosities (median L bol ∼ 10 44.8 ) compare with our sample, we utilise the sample by Lamperti et al. (2020).We can therefore test how our more luminous quasars are different in CO excitation.This comparison sample comprises of X-ray selected AGN, for which further information can also be found in Ricci et al. (2017); Koss et al. (2021).These data are also used in Figure 6.Finally, we also used a compilation of local star forming galaxies as a comparison Leroy et al. (2022a), which includes data from HERACLES (Leroy et al. 2009), the James Clerk Maxwell Telescope (JCMT) Nearby Galaxy Legacy Survey (NGLS Wilson et al. 2012), the CO Multiline Imaging of Nearby Galaxies (COMING) survey (Sorai et al. 2019), PHANGS ALMA (Leroy et al. 2021), IRAM 30m CO (2-1) observations, and Large APEX Sub-Millimetre Array (LASMA) CO (3-2) observations.

Spectral properties
We analyse the spectral properties across different CO transitions for all targets in our sample to investigate the integrated fluxes, line profiles, including line widths, velocity offsets and features within the lines (e.g.potential outflow components).These can then be compared to properties of the host galaxies and the line profiles of the ionised gas to search for any influence of AGN activity.
Spectra are mostly plotted for each target with the same bin widths across all transitions so that the line profiles of each transition can be easily compared (see Figure 3 and the remaining spectra in the supplementary material, Figures B1 -B17).In some exceptional cases, where we had enough S/N in some transitions to investigate the line profile in more detail, but not enough S/N in other transitions, we choose the bin widths accordingly.Central frequencies (where  = 0 km s −1 ) have been defined using the SDSS redshifts quoted in Table 1.We choose the SDSS redshifts as this has been used throughout the QFeedS survey work, and since we are only comparing CO and ionised gas lines within the same target, the specific reference velocity/redshift is not important.The aperture from which all the spectra are taken is consistent between CO(1-0) and CO(2-1) for each source, and for higher transitions will be slightly smaller due to the APEX beam size reducing as observing frequency increases (as mentioned in Section 3).However, our analysis is done in a way that is consistent between transitions and therefore any differences found are potential indications of the impact from feedback mechanisms being different on the different CO transitions.
As an example we present the spectra of all CO transitions (CO(1-0), CO(2-1), CO(3-2), and CO(6-5)) and the MUSE [O III] emission for one target in the main body of the paper (J1100+0846 shown in Figure 3) and all remaining spectra are then presented in the supplementary material.For all sources, [O III] line profiles are plotted in a separate panel below the CO spectra for comparison (in some cases the [O III] line was not available so the H emission line was used instead) extracted from MUSE data using the same aperture as that of the CO data.For some sources of our sample (J1010+0612, J1100+0846, J1356+1026 and J1430+1339), ALMA observations of the CO(2-1) at 0.2 arcsec resolution presented by Ramos Almeida et al. (2022) are available.We show these ALMA spectra plotted in orange over our APEX CO(2-1) spectra to compare.Making these comparisons required a velocity shift to the Ramos Almeida et al. (2022) data to match the zero velocity used here, which we determined using the SDSS redshift, as opposed to the approach taken in Ramos Almeida et al. (2022).Specifically, they used the SDSS redshift as the initial  = 0 km s −1 and then applied a small shift to make the peak (or centre of two peaks) at  = 0 km s −1 .Despite this small difference, we find consistent flux values and line profiles within errors when comparing to Ramos Almeida et al. (2022).We also note the availability of CO(1-0) and CO(3-2) data for J1356+1026 from Sun et al. (2014), however due to the differences in resolution (at 1.3 and 0.6 arcsec for CO(1-0) and CO(3-2) respectively) and therefore high chance of over-resolving the total CO emission compared to the QFeedS observations, these are excluded from any analysis.
From a first look at our CO spectra it is immediately clear that there is a large variety of line profiles, luminosities and detections for the different transitions within our sample.We also find a large diversity of broad line profiles, double peaked profiles as well as a blue wing and offsets from  = 0 km s −1 .A detailed comparison of the molecular and ionised gas line profiles is discussed in Section 5.3.
To analyse the line profiles of both the CO transitions and MUSE [O III] spectra, we fit either one, two, or three Gaussian components to the spectra where appropriate.In the CO data, since we have relatively low S/N, we find only two targets with more than one Gaussian component.However, in the MUSE data we see a wide range of line profiles, including many with multiple components.From the data we analyse the central velocity ( 50 ) and line width over which 80 percent of the flux is contained ( 80 ) which can either be done using the fits to the data or to the data itself.We choose to present the values calculated on the fits to the data for the following reasons.
There are 21 cases for which we are measuring  50 and  80 of spectra with S/N > 5 (used in future analysis) and of those 21, 17 of the  80 values from the data are within the uncertainties of the  80 from the fit, whilst 4 are outside the fit uncertainties.Three sources have higher S/N data published by Ramos Almeida et al. (2022); of these, two have line widths closer to the  80 values from the fit, while one is closer to the  80 of the data.Finally, of these 21 cases there are 16 where the  80 of the fit is larger, while there are five where the  80 of the data is larger.Given that both measures are mostly consistent with each other, that in two out of three cases the  80 of the fit is closer to higher S/N data in the literature and that in most cases  80 of the fit is larger, we decide that the  80 of the fits is more appropriate for this work.Since we are looking at relatively low S/N data across the sample, either measure used will give large uncertainties and so using the larger of the two means we are closer to the maximum value from the data, which is considered later in the analysis of the line profiles and further discussion in Section 5.3.
Therefore, in all cases we use the  50 and  80 of the fits to the data to analyse as the description of the line profiles.In the case of spectra which show a single Gaussian profile, the uncertainties on  50 and  80 are the errors on the fit.The uncertainties for the total flux is the uncertainty on the fit plus the uncertainty on the telescope efficiency.For those with multiple components, the uncertainty on  50 presented is the uncertainty on the narrowest components peak velocity.The uncertainty on  80 is the uncertainty on the  80 of the broadest component.The uncertainty on the flux for those with multiple components is the uncertainties on each component added in quadrature plus the uncertainty on the telescope efficiencies (∼ 5 -10% for all observations).The noise for each channel for APEX data was taken from the results of the CLASS reduction scripts.For ACA observations the rms was calculated by taking the median flux from the remaining spectra, whilst masking the line (where present).
In order to determine which lines and fits are robust and can be used for further analysis we measure the data quality based on the S/N of the lines within our sample.We then define different levels of detection: "detections", "low S/N detections" and "non-detections", which are classified in the following way: • "Detections" are defined as spectra which show lines with an integrated S/N ≥ 5.
• "Low S/N detections" are defined as those lines with 3 ≤ integrated S/N < 5.
• Anything with no clear line or a line with an integrated S/N < 3 is defined as a "non-detection".When referring to "detections" from this point onward we refer to both "detections" and "low S/N detections" as defined here unless otherwise stated.
For those with low S/N detections we present the fits to the spectra with a dashed black line to differentiate from detections shown by solid black lines.The integrated S/N values are also shown in the tables presenting CO data.For a list of the classifications of all detections, low S/N detections and non-detections see Table A1 in the supplementary material.
There are three sources which show non-detections in all CO transitions which are J0909+1052, J1222-0007 and J1356+1026.However, we note that J1356+1026 does have a CO(2-1) detection in the deeper ALMA data from Ramos Almeida et al. (2022) which we utilise in later analysis.Several other targets also show non-detections in at least one transition.To test whether the non-detections might show any low-brightness signal when combined, we stack these data by bringing the observations from different targets to the same velocity reference.However, from stacking non-detections we don't find any underlying flux and thus we can make no more conclusions based on these data.In all cases, solid black lines denote fits to the data.Here there are multiple components and as such, the dotted black lines denote the different components that make the total fit to the spectra.In the CO(2-1) spectra we also show the fit from higher resolution ALMA observations (solid orange line labelled RA+22, Ramos Almeida et al. 2022).Shaded grey regions represent the 1 level.Spectra for all targets are shown in the appendix, using the same presentation methods.MNRAS 000, 1-21 ( 2023) From these spectra we can calculate total fluxes, line widths ( 80 ), velocity offsets ( 50 ) and line luminosities (all data are presented in the appendix in Tables: A2, A3, A4, A5, A6 and A7.For those with non detections we choose to present the 3 upper limits for CO flux and luminosity.Further observations would be required to confirm any of these detections.Line luminosities (L ′  ) are calculated using the following equation from Solomon et al. (1997) (also used in the analysis of previous QFeedS work Jarvis et al. 2020): where  CO,rest is the rest frequency of the CO line, D L is the luminosity distance,  is the redshift and  is the velocity integrated line flux density measured in Jy km s −1 .
For those without CO(1-0) detections we use the CO(2-1) luminosity where possible and convert to CO(1-0).Conversions made from L ′ CO (2 − 1) use the median line ratios observed within this sample (presented in Table 3) as conversion factors (median line ratio of 1.06).Here the sources J1010+0612, J1010+1413, J1356+1753 and J1430+1339 are calculated using L ′ CO (2 − 1).For those with nondetections across all CO transitions we provide 3 upper limits of M CO based on the upper limits of the CO(1-0) flux.
The calculated values of  CO from our sample are within the range 4.0 -4.2 (shown in Table 2).These values are consistent with typical high redshift, high star-forming, quasar host galaxies (Bolatto et al. 2013;Carilli & Walter 2013).However,  CO may also be significantly lower in LIRGs, submillimeter galaxies, mergers, starbursts and AGN, with values as low as ∼ 0.6 -1 (e.g., Bolatto et al. 2013;Sargent et al. 2014;Calistro Rivera et al. 2018).Therefore there are uncertainties that arise in these values and the calculated gas masses.There can also be dependencies on the metallicity and SFR (see e.g.Bolatto et al. 2013;Sandstrom et al. 2013, and references therein), but for most galaxies a value of ∼ 4 is found, as is identified in our targets using the same method of calculation as the comparison sample in Tacconi et al. (2018) (which we follow for consistency).
The molecular gas masses in our sample range from 0.36 -5.50 ×10 10 M ⊙ and all molecular gas masses can be found in Table 2 along with stellar mass estimates obtained from SED fitting (see Jarvis et al. 2020, for a description of these calculations).There are cases with large uncertainties which result from poor constraints on the SED fitting.These uncertainties also follow through into Figure 4 where we present the stellar and CO gas masses.Gas fractions are also presented in Table 2 calculated as M CO /M ★ .We find values in the range 0.06 -1.4.These are used in Figure 4 to compare to AGN and non-AGN from the literature (Tacconi et al. 2018).See section 4.1 for further information about this comparison sample.We find that our sources are consistent with the comparison sample of non-AGN and AGN.There are a few cases where we see higher gas fractions, in particular J0945+1737, J1108+0659 and J0958+1439 but these are still consistent with both AGN and non-AGN from the literature (see Figure 4).In fact, AGN have been found to have similar, or higher gas fractions than non-AGN (e.g.Rosario et al. 2018;Kirkpatrick et al. 2019;Jarvis et al. 2020;Shangguan et al. 2020;Koss et al. 2021;Zhuang et al. 2021;Salvestrini et al. 2022).Further, in AGN CO gas can also be detected in a warm phase (e.g.Rosario et al. 2019) so the molecular gas levels here could also be considered a lower limit of the total molecular gas content.However, it is important to note that even with the presence of high gas fractions, it does not exclude the possibility of AGN feedback (including cases with the presence of outflows) as shown in comparisons of simulations (Ward et al. 2022).Jarvis et al. (2020) discuss gas fractions as well as analyse the star formation rate, specific star formation rate and distance from the main sequence for a subsample of these data.We note however, that there are differences in this work when comparing to the same targets in Jarvis et al. (2020), which used  21 of 0.8 to convert from CO(2-1) to CO(1-0).Given that we have the CO(1-0) measurements here we simply use those, and as discussed before we find higher values of  21 with a median of 1.06 within our sample.As such, the results of the total gas masses and gas fractions are also different here compared to Jarvis et al. (2020), with this work finding lower total gas masses, but within uncertainties and therefore still consistent.

CO excitation
To measure and analyse the excitation of the molecular gas we analyse both the observed shape of the CO SLEDs (Section 4.4.1) as well as the line ratios of the CO transitions (Section 4.4.2).Studying the excitation of the gas and comparing to literature samples of both AGN and non-AGN will again allow us to determine whether the excitation in the quasar host galaxies is systematically different to that of other relevant galaxy samples.

CO SLEDs
We present the CO SLEDs of our sources in two different sets.Firstly, in Figure 5a we show CO SLEDs from the ground state up to CO(3-2), making comparisons to literature samples.We plot only those with detections in CO(1-0) (including low S/N detections), so that we have a reliable normalisation to the ground state.For Name  2. We present   (calculated using equation 3), M CO (calculated using equations 2, 3 and 4), M ★ and gas fractions ( M CO /M ★ ) for each target in the sample.For J1116+2200 the stellar mass is unconstrained by the SED and so there is no value of M ★ (and therefore also  CO ). Calculated from CO(1-0) line luminosity. Converted to CO(1-0) from CO(2-1) using the median line ratio of those in the rest of the sample with detections. Median  CO value of 4.19, and M CO calculated from this median.
those without a detection in CO(1-0) deeper observations would be required to provide a reliable CO SLED.The reference SLEDs shown here are the Milky Way and thermalised SLEDs, shown by dotted lines (Carilli & Walter 2013).We further make comparison to starburst and main sequence galaxies at  ∼ 1 -1.7 taken from Valentino et al. ( 2021) (see Section 4.1 for further information on these comparison samples).
We report that the CO SLEDs of 2 out of 9 sources are consistent with being above the thermalised relation at the CO(2-1) level (excluding the upper limit of J1055+1102) and 2 out of 9 at the CO(3-2) level (see Figure 5a).
However, including the uncertainties on the CO(1-0) fluxes in addition to the CO(2-1) or CO(3-2) they are all still consistent with being thermalised, with the exceptions of  21 for J1100+0846 and  32 for J1010+1413 (see Table 3 for individual values).
For sources within our sample with available observations in higher CO transitions (J   = 6, 7) we plot the CO SLEDs extending to these higher transitions and normalise the SLED to CO(2-1) instead of CO(1-0) (see Figure 5b).Normalisation to the CO(2-1) transition was done for two reasons.Firstly, two of our targets with high CO transition data have not been observed in CO(1-0) and so to make comparisons within our own sample the next transition available with detections for all targets was CO(2-1).Secondly, based on our observations of several super-thermal SLEDs at CO(2-1) and the ongoing discussion in the community about the optical thickness of CO(1-0) (see Section 5.1.1),we argue that CO(2-1) might be more reliable transition for normalisation.Again, in Figure 5b we make comparisons to other samples from the literature.In addition to the literature samples mentioned previously we also include high- quasars ( ∼ 1 -6) compiled by Carilli & Walter (2013), to make comparisons to the low- counterparts in our sample.
Out of the 7 sources observed, 6 show non-detections and only 1 source, J1430+1339, has a detection in either CO(6-5) or CO(7-6), based on the same criteria stated in Section B. We find that this detection is relatively low in the CO SLED, suggesting that the peak of the SLED may be at J  < 6 (see Figure 5b).It appears to be more similar to main sequence and starburst galaxies as opposed to high- quasars (Carilli & Walter 2013).
For the remaining CO SLEDs in Figure 5b we only have upper limits, however two of these are very constraining namely J1100+0846 and J1108+0659, which are at the same excitation level as the already detected J1430+1339).5 targets are less excited that the high- quasar CO SLED, with 3 of these showing consistent excitation with starburst and main sequence galaxies.J1100+0846 shows signs of a detection, with a S/N of 2.2, but the 3 upper limit is at a similar excitation to our detected J1430+1339.J1108+0659 shows no sign of any detection at a very constraining level in the CO SLED, placing it already at a similar excitation to starburst and main sequence galaxies.Likewise, J1518+1403 shows no signs of a detection but is between the excitation of the high- and starburst CO SLED.J1010+0612 shows a very tentative signs of a line with a S/N of 1.1 at a higher excitation than starburst/MS galaxies, but still less than the high- QSO CO SLED.Further observations would be required to confirm any detection and place a proper constraint on the excitation.The remaining two targets (J0945+1737, J1116+2200) can't be analysed in much detail since the upper limits are not very constraining.However, the fact that these upper limits are already at the level of the thermalised relation and that of the high- quasars (see Figure 5b and Carilli & Walter 2013) shows we are likely observing systems with a significantly lower excitation.This along with the other targets clearly showing lower excitations shows a clear difference between our sample of quasars at  < 0.2 and those at  ∼ 1 -6.

Line ratios
We can also investigate the excitation of the gas by calculating the ratios of line luminosities of different CO transitions.We calculate  (5) where L ′ CO is the CO luminosity of a given CO transition line.
All measured line ratios for each individual target are shown in Table 3.We present the observed line ratios found in our sample of quasars and compare to literature values in Figure 6.In the lower panel of Figure 6 we present histograms of the line ratios from our sample, as well as violin plots for different reference samples in the upper panels.We choose violin plots as these show more information of the distribution of the data, with a wider section showing a larger number of data, as well as the maximum and minimum values in the range of the data plus a defined median and 16/84th quartiles.For targets within our sample which show multiple components in their CO spectra, we calculate the line ratios for these individual components and these are also presented in Table 3.
The overall line ratios observed in this sample (only including those targets with detections, and ignoring non-detections) are as follows: For those in our sample with detections in both CO(1-0) and CO(2-1) (6 sources), we calculate the line ratios using equation 5 and find a median  21 of 1.06 +0.53  −0.18 , where the negative and positive uncertainties represent the 16th and 84th quartiles respectively.For  31 we find a median of 0.77 +0.31  −0.20 from 6 sources with detections in both CO(3-2) and CO(1-0).We find a median  32 of 0.61 +0.43  −0.21 from 8 sources with detections in both CO(3-2) and CO(2-1).We find that in those targets with multiple components, the excitation of the different, individual components and of the total emission across the entire spectral line from the galaxy are consistent (within uncertain-ties).Therefore, with the data we have available, we cannot measure any difference between the excitation levels in these different components.The median line ratios in our sample along with literature comparisons can be found in Table 4.
Making comparisons to lower redshift, less luminous samples of AGN and star forming galaxies however, we find higher line ratios (see Figure 6).For example, Lamperti et al. (2020) present a study of 36 Hard X-ray selected AGN at z = 0.002 -0.04, conducted as part of the BASS AGN sample.They present twelve targets with the requisite data to calculate the  21 values and find a median of 0.72 +0.17 −0.17 .These are lower than those measured in our sample.Further, the  32 of the same sample also shows lower excitation with a median of 0.50 +0.20  −0.19 .
As expected the line ratios of our sample are higher than for normal, star forming galaxies (e.g.Saintonge et al. 2017;Yajima et al. 2021;den Brok et al. 2021;Leroy et al. 2022a).For example, a compilation of low-z samples found a median  21 of 0.65 Leroy et al. (2022a).Further, we analyse the higher CO line ratios, which are mostly upper limits for our sample, and utilise data of high- quasars which have CO line fluxes in at least one higher transition (taken from the compilation by Carilli & Walter 2013) as well as corresponding bolometric luminosities (taken from Trentham 1995;Lewis et al. 1998;Lutz et al. 2007;Aravena et al. 2008;Bradford et al. 2009;Wang et al. 2010).We also note that in comparing bolometric luminosities we convert our [O III] luminosities to bolometric luminosities via the following equation: L bol /L [O III] = 3500 (from Heckman et al. 2004).We present these data in Figure 7 by also comparing the bolometric luminosities of the targets.

Gas temperature and density
We can utilise the measured line ratios to help understand the physical conditions of the molecular gas in individual targets.In particular, comparing the measured line ratios can give us an indication of the gas temperature and density within these quasar host galaxies (Peñaloza et al. 2017;Leroy et al. 2022a).Specifically, the ratio of  32 with  21 gives hints as to these properties (see Figure 8).We present these data plotted over simulations of the expected parameter space covered for these line ratios based on variable temperature, densities and optical depths (Leroy et al. 2022a).Figure 8 shows this effect of temperature, density, and optical depth on the observed line ratios (grey points from Leroy et al. 2022a).The expected line ratios for starbursts and AGN would be close to 1 (e.g.Mao et al. 2010;Lamperti et al. 2020;Yajima et al. 2021).These values are consistent with having both higher densities and hotter gas (Leroy et al. 2022a).
To further investigate the temperature and density of gas in these quasars we analysed the 5 targets with detections in all of the first three transitions (J1000+1242, J1100+0846, J1108+0659, J1116+2200 and J1518+1403), which provide best opportunity to test these properties.Using the Dense Gas Toolbox (van der Tak et al. 2007;Leroy et al. 2017;Puschnig 2021), the density of the gas and the dense gas fraction (fraction of gas with density > 10 5 cm −3 ) were calculated by fixing the temperatures in 5 degree increments in the range of 10 -50 K.In 4 out of 5 sources (except J1000+1242) a temperature of less than 35K lead to the highest density provided in the model (10 5 cm −3 ) and dense gas fractions greater than 90%.From 40K -50K these four targets give densities in the range 10 5 -5000 cm −3 .J1000+1242 is the only target showing different properties of the temperature and density.Only at 20K does the model give the highest density with a dense gas fraction greater than 90%.From 50 -25K a gas density ranging between 600 -4000 cm −3 was calculated.This is significant as J1000+1242 has the lowest line ratios amongst the 5 and is the only one with  21 < 1.This therefore shows that within our sample, those with high line ratios would require higher temperatures (> 35K) and densities.For those in this sample with lower line ratios, temperatures as low as 25K can still provide realistic scenarios.
From analysing the line ratios in Figure 8 we find that the molecular gas in J1116+2200 seems consistent with being optically thicker and J1100+0846 seems consistent with being optically thinner.J1000+1242 and J1518+1403 are most likely somewhere in the middle and J1108+0659 is more difficult to determine due to the larger uncertainties.Further studies with higher S/N observations across the entire sample would be required to make any further conclusions.However, from the large variety identified in these limited data we can say that there is not a particular tendency in optical depth valid for all our sources.

Comparing CO and ionised gas line profiles
To investigate any differences between the molecular and ionised phase of the ISM we analyse the CO and [O III] line profiles.We do this by comparing both the velocity offsets ( 50 ) and the line width ( 80 ) for the CO lines compared to MUSE observations of the ionised gas (in most cases traced by [O III] and for two cases H, where the [O III] was not available).We present these results in Figure 9.We only perform this analysis on those with an integrated velocity integrated S/N greater than 5 so that the line profiles we compare are more reliable.This threshold will ensure that there is a lower uncertainty on the CO line profile, allowing us to make more reliable comparisons between the ionised and molecular phase.From these results, in all cases we observe consistent or broader [O III] line profiles than in the CO transitions (as measured by  80 , shown in Figure 9).
As well as the line widths, we can also analyse the velocity offsets of CO compared to [O III] (see Figure 9).From this we find that a number of targets are consistent in terms of the  50 measurements, whilst others do show significant differences, both positive and negative.While sources which show more redshifted [O III] lines than CO also exhibit broader [O III] line widths, the few cases with more blueshifted [O III] lines present equivalent [O III] line widths to those observed in CO, which are among the most similar line profiles that we observe.
For our sample overall we also find more complex line profiles in the ionised gas than in the molecular phase.Specifically we identify 2 out of 17 sources with CO spectra that show multiple components in CO whereas in the MUSE ionised gas data we find 13 out of 17 sources with multiple components.This is limited by S/N in our CO data but interestingly, in the 2 cases where we find multiple components in the CO (J1100+0846 & J1108+0659) we identify the following: J1100+0846 shows a clear double peaked line profile in the first three CO transitions.This double peaked profile was also identified in Ramos Almeida et al. (2022) and we confirm the similarity by plotting them together onto our CO(2-1) spectra.This double peak however is not present in the [O III] where instead a broad line is identified.
Another target within our sample (J1010+0612) also shows a double peak CO(2-1) profile (identified in higher resolution ALMA data by Ramos Almeida et al. 2022).In our APEX spectra we also see tentative signs for a double peak profile in both our CO(2-1) and CO(3-2) data (see spectra in Figure B5 in the supplementary material).However, like J1100+0846, J1010+0612 shows a broad line in the MUSE [O III] data, without signs of two narrower components as identified in CO.
J1108+0659 shows a prominent blue wing component, identified in the first three CO transitions.The spatial resolution in our CO(1- 0) data is not enough to spatially locate where this potential outflow component is located.Interestingly, the H line profile also shows the same blue wing component, but with a more obviously present blue wing in H than the CO.The peaks of the two components are also at very similar velocities across the three CO transitions and the H data.
Looking at the ionised gas in isolation, we extracted [O III] spectra from the MUSE cubes at both 30 arcsec (to match our CO observations) as well as at 3 arcsec (to match the SDSS observations).We find that the observed [O III] line profiles in SDSS and extracted from 3 arcsec apertures from MUSE are consistent within uncertainties (Figure 10a).However, we find a larger scatter when comparing SDSS line profiles to those extracted from 30 arcsec apertures in MUSE (Figure 10b).There are 8 cases for which they are not consistent, 7 where the SDSS lines are broader, and 1 case where MUSE at 30 arcsec is broader.Cases where SDSS are broader suggest a larger impact on velocities of ionised gas close to the core but further work would be needed to confirm this.Interestingly, the only case which is broader in MUSE at 30 arcsec than SDSS is J1016+0028, with the largest radio size, known radio lobes extending to distances of ∼ 15 arcsec, suggesting that impact on velocity of ionised gas is present out to these larger distances.There is no clear overall picture or trends from these results and so any differences likely depend on a case by case basis, but these observations do hint at potential differences in velocities of ionised gas in the halos, and that the integrated velocity is dominated by the kinematics in the core, since 7 out of 8 are broader at 3 arcsec.However, further and more focused studies, are needed to investigate these effects in more detail.

DISCUSSION
Here we discuss the findings and an interpretation of the results in relation to the literature.In Section 5.1 we discuss the observed CO excitation within our sample using CO SLEDs and line ratios, including comparisons to the literature.In Section 5.3 we discuss the differences in the line profiles seen in our CO data compared to ionised gas observations from MUSE.Finally, in section 4.5 we discuss the implications of our observed line ratios on the molecular gas temperature and density in our sources, based on theoretical models.

No observed impact on total CO excitation
Here we will discuss the findings of the CO excitation, in the context of the selected comparison samples from the literature (presented in Section 4.1).Overall, these sources which have evidence for AGNdriven ionised outflows and jets/winds do not have exceptional global CO content or excitation, beyond that seen in relevant comparison samples of (U)LIRGS and less luminous AGN.However, we discuss further that a local impact on CO excitation is still a possible scenario, and if present is likely to be driven by AGN feedback processes such as via radio jets.

Consistent CO excitation to local (U)LIRGs
Previous studies into the molecular excitation have found no consensus on whether AGN have systematically different r 21 line ratios (see Ocaña Flaquer et al. 2010;Papadopoulos et al. 2012;Xia et al. 2012;Husemann et al. 2017;Shangguan et al. 2020), or the r 31 line ratios (Sharon et al. 2016).With this work we looked for any correlations with these line ratios and AGN or galactic properties.
Our  21 values (median  21 of 1.06 +0.53 −0.18 ) are consistent with those of recent studies of (U)LIRGs (Greve et al. 2014;Montoya Arroyave et al. 2023) with median values of 0.95 and 1.05 respectively.The  31 and  32 values are also found to be consistent with these (U)LIRG studies (as shown in Figure 6 and Table 4).As mentioned previously, 9 quasars in our sample of 17 have L IR,SF data, of which 8 are consistent with being LIRGs (see Section 2 and Jarvis et al. 2019).The ninth source, J1430+1339, is just below the threshold with log[L IR,SF ] = 44.32+0.06 −0.07 erg −1 .The fact that we see similar line ratios to these studies of local (U)LIRGs is an indication that the presence of a quasar and the observed radio jets and ionised outflows has no significant impact on the excitation of the global molecular gas content at the mentioned transitions.
As mentioned above, we find slightly higher excitation in our sample than found in lower redshift AGN (Lamperti et al. 2020), with the distribution of QFeedS line ratios extending to higher values.These targets have a lower luminosity compared to our sample with a median L bol ∼ 10 44.8 which may explain the small differences observed in the CO excitation, despite the presence of AGN.However, it is important to note that the differences are not large (particularly when considering the smaller sample size in QFeedS) and so there is no clear impact from the presence of AGN-driven ionised outflows and jets/winds.
With the diversity of lines ratios observed for different populations, depending on the type of galaxy and properties, we note that extreme care should be taken when converting from higher J  lines to the ground state to make sure that correct line ratio is used based on these factors.

Optical depth effects
As presented in Section 4.4.1, we find cases where the CO SLEDs are above the thermalised relation in CO(2-1) and CO(3-2).The most straight-forward explanation might be that some of the CO(1-0) emission may be resolved out given the ACA spatial resolution, particularly if there is extended diffuse emission beyond the beam size of the observations, therefore obtaining a lower flux value than the true value.The larger uncertainties on the CO(1-0) flux values may also be a factor in this.Interestingly two out of the four non-detections in the ACA CO(1-0) data are known to be ongoing mergers/dual AGN (J1222-0007 and J1316+1753, see Jarvis et al. 2021).However, apart from this there are no special properties of these sources that would differentiate them from the rest of the sample, and as mentioned in Section 3.1, the sensitivity calculations for the observations were done using the same method across the sample.
For those that are above the thermalised level in the CO SLEDs, and assuming the thermalised level should be the maximum flux, then this would imply that between 7 and 50 per cent of the flux would be in extended emission at scales greater than 188kpc (median recoverable scale of the observations).Since this would be an unrealistic scenario, we favour a physical interpretation rather than an observational one.In addition, we can compare to Montoya Arroyave et al. ( 2023) which contains more nearby objects (at a median  = 0.09 compared to a median  = 0.14 for the QFeedS sample) and would therefore more subject to over-resolution effects with ACA CO(1-0).Despite this, they find very similar line ratios to ours (Figure 6), again lending credence to the interpretation that this is a physical and not an observational effect.
Although we cannot exclude the influence of the above mentioned observational effects we note that other investigations on ULIRG samples have reported similar trends, indicating that these superthermal  21 ratios may have a physical explanation.For example, the high  21 values in our sample could be due to optical depth effects, with highly excited gas in combination with a low opacity in the CO(1-0) transition (Zschaechner et al. 2018).It has been argued that low opacities can be driven by large velocity gradients and would require the presence of turbulent or outflowing gas, perhaps also in a diffuse, warm phase (Cicone et al. 2018;Montoya Arroyave et al. 2023).In the case of our sample this could mean kinematic disturbances as a result of quasar feedback (e.g.Jarvis et al. 2019;Ramos Almeida et al. 2022;Girdhar et al. 2022).However, in the data presented here we only see tentative indications of such outflows, but this may be the result of limited S/N and further studies with deeper observations at a higher spatial resolution are required to investigate this further.
It therefore may well be the case that the CO(2-1) is a more reliable line, especially when investigating the excitation through CO SLEDs.For example Montoya Arroyave et al. ( 2023) claimed that one cannot determine the CO excitation from the low-J CO lines, most likely because the optical depth effects add too much noise to the calculation.They also identified that weak correlations are found only from ratios involving the CO(1-0) line and L IR and SFR.This therefore supports the idea that CO(1-0) is the only line affected by optical depth effects and not the CO(2-1) or CO(3-2).

Low excitation at CO J 𝑢 𝑝 = 6, 7
As shown in Figure 5b, we observe relatively low excitation in the CO(6-5) and CO(7-6) transitions, for at least 4 out of 7 targets with the appropriate data.We note that the differing beam sizes between the observations at CO(2-1) compared to CO(6-5) and CO(7-6) may play a role in the low excitation observed, especially if flux is resolved out in the smaller beam for observations in higher CO transitions.However, similarly low excitation is found when compared to the CO(1-0) ground state which has similar beam sizes to the higher CO transition APEX observation.From previous work into the CO(2-1) kinematics for some of these targets (Ramos Almeida et al. 2022) and the fact that higher CO transitions are expected to be more centrally located we argue that these line ratios are reliable and still very constraining.
One explanation (as discussed previously in relation to the  21 values) could be due to a difference in bolometric luminosity.The high- quasars that show much higher line ratios do indeed have significantly higher bolometric luminosities (see Figure 7), with luminosities of 10 47 -10 48 ergs −1 compared to our sample which are all L bol < 10 46.5 ergs −1 .The high excitation levels observed in quasar host galaxies with extreme luminosities, like those presented in Carilli & Walter (2013), and in particular unobscured sources (Banerji et al. 2018;Wang et al. 2019;Li et al. 2020;Bischetti et al. 2021), can often be associated with depleted gas reservoirs (e.g.Brusa et al. 2015;Perna et al. 2018;Förster Schreiber et al. 2019;Circosta et al. 2021).These depleted gas reservoirs, potentially as a result of quasar activ-ity, could therefore lead to lower observed CO fluxes at J  < 3, and consequently high excitation at J  > 5. Since we observe gas-rich systems in our quasar host galaxies, the lower excitations that we measure might therefore be expected.Further studies of both low and high- quasars covering much of the CO SLED are required to test this further (e.g.Novak et al. 2019;Pensabene et al. 2021).
It is also possible that the effect of AGN may only be detected at J CO > 10.CO is most commonly excited by photodissociation regions (PDRs) from the UV photons emitted from young stars.However, excitation at these higher CO transitions requires shocks and/or X-ray emission (through X-ray-dominated region models, XDR), both of which can be powered by AGN or jets (Pereira-Santaella et al. 2013;Carniani et al. 2019).The effects of both XDR and PDR regions play an important role but observations with parsec scale resolution are required to disentangle and analyse these regions (Wolfire et al. 2022).Therefore, observations at even higher CO transitions and high resolution may be needed to detect the XDR-dominated CO lines and provide a more complete constraint on the influence of AGN (see e.g.van der Werf et al. 2010;Mashian et al. 2015;Carniani et al. 2019).
Differences in the overall star formation and ISM conditions at low and high redshifts could also contribute to the differences we observe in CO excitation between our sample and high- quasars.Estimates of the cosmic molecular gas density indeed suggest that the molecular gas fractions peak at redshifts of  = 1 -3 (see Péroux & Howk 2020, for a review), roughly mirroring the cosmic star formation rate density (Madau & Dickinson 2014) and black hole accretion density (e.g., Aird et al. 2015).As a result, we might therefore expect higher gas temperatures and densities at nuclear scales in quasar host galaxies at redshifts 1 -3 (higher than our sample), excited by stronger radiation fields from star formation and the AGN.However, there are cases where high- quasars have also been found to have lower than expected CO excitation.For example, a recent study of nine  ∼ 3 quasars (Muñoz-Elgueta et al. 2022) finds that their CO SLEDs peak in the range J up = 5 -7 compared to the expected J up > 6 -8.Further studies with coverage across a large range of the CO SLED is required for both low and high- quasars to investigate these findings further.
Another factor could be the potential impact of obscuration and line-of-sight effects, when comparing Type 1, Type 2, and red quasars.For example, the CO emission in unobscured quasars might be dominated by gas within the ionising cone through which we observe Type 1 quasars (e.g.Vayner et al. 2021;Stacey et al. 2022).Despite this, no significant differences in the molecular gas content and star formation efficiencies have been reported between obscured and unobscured AGN in the context of AGN unification and line-of-sight effects (e.g.Perna et al. 2018).However, the increased incidence of high-velocity [O III] outflows and radio emission in red quasars has been associated with higher nuclear dust reddening at high redshifts (e.g., Perrotta et al. 2019;Klindt et al. 2019;Fawcett et al. 2020;Calistro Rivera et al. 2021;Andonie et al. 2022).This connection might suggest an increased amount of obscuring material (dust and gas) at nuclear scales for luminous quasars with ongoing outflows and jets or winds.

Line ratios as tracers of temperature and density
In Section 4.5 we presented the results of using the line ratios  21 and  32 to give an indication of the temperature and density of gas within the QFeedS sample.We are limited by the lack of detections across the CO SLED, with only 5 sources having detections in all of the first three J  transitions, and by the large uncertainties in the line ratios (see Figure 8).Further, the dynamic range of the line ratios observed in AGN, quasars, (U)LIRGs and star forming galaxies is not large and therefore, more accurate, high S/N observations are required to analyse these to a high degree (i.e.placing accurately on Figure 8).
Despite this, we can place constraints and upper/lower limits on the temperature and density as shown in Section 4.5.Indeed, we have shown a difference between those with line ratios ≥ 1 and those < 1.The 4 out of 5 sources with higher line ratios ≥ 1 in our sample required temperatures > 35K to show dense gas fractions less than 90%.On the other hand, for J1000+1242 with both  21 and  32 < 1, this was reduced to 25K.
Observations to constrain the CO properties of other targets in QFeedS would be an important next step, as well as observing other lines such as CI, HCN or HCO+ which would also help further constrain the the gas temperature and densities.

A localised impact on CO?
Although our sample consist of luminous quasars with known ionised outflows, radio jets and/or large radio structures, we report no signs of enhanced CO excitation on the global scales in which we are measuring when compared to local (U)LIRGs of similar luminosities.
There has previously been found a positive relation with L IR and the line ratios  21 and  31 (Montoya Arroyave et al. 2023), as well as a positive relation between  31 and SFR.However, as mentioned in Section 5.1.2they claim that one cannot determine the CO excitation from the low-  lines because of optical depth effects.
We performed some basic analysis on the line ratios compared with galactic and AGN feedback diagnostics such as SFR, radio luminosities, FWHM [O III] , finding no correlations within our sample.However, when working with such small sample (less than 10 sources for which we have the requisite data) and small parameter space, it is hard to take strong conclusions from these investigations.
Despite this, a lack of any correlation may not be surprising.For example, a study on local Seyferts found no clear evidence for a systematic reduction in the molecular gas reservoir at galactic scales with respect to SFGs (Salvestrini et al. 2022).Previous studies have also found weak or no correlation with properties such as stellar mass, AGN fraction and SFR offset to the main sequence (Liu et al. 2021).Recent work has also found no correlation between the cold molecular gas properties and AGN properties (Molina et al. 2023).Additionally, studies at both low and high redshifts have found no differences in low-  excitation between samples of star forming galaxies (taken form the xCOLD GASS Survey Saintonge et al. 2011) and AGN host galaxies (Sharon et al. 2016;Lamperti et al. 2020).This finding again indicates a lack of influence from the AGN on the total molecular gas content.Further support for a lack of impact on global scales is the evidence of extended CO(1-0) emission as shown in Figure C1 in the supplementary material These data may be an indication that for those with detections we still observe an extended molecular gas reservoir.On the other hand, those with non-detections may be an indication of disruption on these global scales, but with the data that is available we cannot make any conclusions about why a few of the targets have gone undetected in CO(1-0).
The typical theoretical prediction is that AGN outflows do not efficiently disrupt disc systems, because the outflow is deflected into the halo (Costa et al. 2020), therefore supporting the hypothesis of a lack of impact on the total molecular gas content.Some simulations predict that an outflow will carve out a small cavity (on scales of ∼ 1 kpc) in the galactic nucleus.As we are far from resolving on these scales here then it is not surprising that we do not see a significant impact on galactic scales.This would also mean that any immediate impact on star formation is also likely to be modest (see also Gabor & Bournaud 2014;Ward et al. 2022;Piotrowska et al. 2022).However, this does not rule out any long-term impact through the effect of outflows on halo gas.
An explanation for our findings could therefore be that the impact of feedback and/or galactic properties on the excitation of the molecular gas may occur on a more localised scale, and once looking at the total molecular gas content, this effect is no longer observed.Indeed, higher spatial resolution studies have found differences in the excitation on scales < 1 kpc (e.g.Dasyra et al. 2016;Oosterloo et al. 2017;Rosario et al. 2019;Zhang et al. 2019;Ellison et al. 2021;Audibert et al. 2023), in some cases also localised next to radio jets.
Further supporting this idea, two of this sample have been studied in a resolved way in which the CO velocity dispersion was observed to be effected perpendicular to radio jets (J1316+1737 Girdhar et al. 2022) and CO temperature ratios were enhanced perpendicular to the radio jet (J1430+1339 Audibert et al. 2023).This provides further motivation for resolved studies of the CO excitation and kinematics around small scale radio-jets.One interesting scenario would also be to determine any dependence on the inclination of any radio jet with respect to galaxy plane (Mukherjee et al. 2018b;Venturi et al. 2021;Girdhar et al. 2022;Meenakshi et al. 2022;Audibert et al. 2023), something which needs to be studied further within this survey.
It therefore may be the case that any impact of AGN feedback on the excitation of the molecular ISM seems to only occur on localised scales (Morganti et al. 2021), but the impact does not take effect over the whole galaxy.Therefore, since we are observing the total molecular gas content in the host galaxies, these smaller scale effects are likely to be lost in the full picture.
One reason for this limited impact may be due to the power of the radio jets, being too weak to penetrate throughout the entire galaxy and they are deflected by interactions with the ISM and are contained within the central region of the galaxy.Another potential for these small scale jets is the potential time scales involved, and that what we observe are young jets which have not yet made their way to have an influence over the whole galaxy (e.g.O'Dea et al. 1991;Morganti 2017;Bicknell et al. 2018).An alternative scenario may be quasar driven winds that drive ionised outflows and simultaneously shock the ISM to produce radio emission in the same region of the galaxy (e.g.Wagner et al. 2013;Zakamska & Greene 2014;Nims et al. 2015;Zakamska et al. 2016;Hwang et al. 2018).

Comparing CO and ionised gas line profiles
Across our sample we see differences between the line profiles of the ionised gas and the CO (see individual spectra in the appendix), suggesting potential differences in the impact of AGN feedback on the different gas phases.
As shown in Section 4.6 we identify broader [O III] line profiles than in the CO transitions indicating a larger impact from feedback (e.g. via known radio jets) on the ionised gas kinematics than in the molecular gas.As discussed in Section B, the CO  80 measure based on the fit that can be considered as an upper limit as it is mostly higher than that measured on the data.
Despite this, we still see that the CO widths are less than, or consistent with, the [O III] line widths.This difference in molecular and ionised gas velocities could be attributed to the different densities, with the denser molecular gas being naturally more difficult to drive to higher velocities (Nayakshin & Zubovas 2012;Mukherjee et al. 2016Mukherjee et al. , 2018a;;Girdhar et al. 2022).
The presence of double-peaked CO lines (e.g. in the case of J1100+0846 and J1010+0612) can be indicative of jet-gas interactions whereby jets are pushing the gas in opposite directions (Kharb et al. 2021).Alternatively, these profiles can indicate that the gas is in a disk, or that binary black holes with individual broad and narrow line regions are present.Binary black holes could have resulted from galaxy mergers.
However, as mentioned the lack of S/N in our CO data means that interpreting these multiple components is limited.On the other hand, we should nonetheless be sensitive to the overall gas kinematics, and so our finding of broader total [O III] line profiles compared to CO is reliable.It is worth mentioning that another explanation for this is that what we now observe as ionised gas was originally molecular gas that became ionised and heated in an outflow.In this case, the fact that [O III] is broader might reflect a shorter survival time of cold, dense gas in the outflow (see e.g.Costa et al. 2015Costa et al. , 2018)), rather than the impact of AGN feedback.
Differences in the kinematics of CO compared to the ionised gas could also indicate that the molecular gas is not mixed in the outflowing ionised medium.This may be as a result of cold gas clouds being unable to survive in hot winds (e.g.Farber & Gronke 2022).

CONCLUSIONS
We present a molecular gas excitation survey, observing a range of CO transitions (J = 1, 2, 3, 6, 7) for a sample of 17 quasars at  < 0.2.Our goal is to measure the molecular gas properties such as molecular gas masses, fractions, and CO excitation, as well as gas kinematics in order to identify any impact due to the presence of radio jets and ionised outflows on a global scale.
-From all the evidence presented here we suggest that the presence of ionised outflows and radio jets in these LIRG type systems does not significantly impact the CO excitation on a global scale, but that given evidence from the literature, localised effects are likely, and do not extend to the scales of the entire galaxy.
-We find no differences between the molecular gas fractions of our sample of quasars as compared to non-AGN in the literature (see Figure 4), in agreement with previous works.
-We suggest that optical depth effects may contribute to the high line ratios involving CO(1-0) that are observed, in agreement with previous studies.
-From analysing the CO SLEDs in 7 targets of our sample (see Figures 5a and 5b), we observe lower excitation in CO(6-5) and CO(7-6) as compared to a sample of quasars at higher redshift ( = 1 -6).We suggest this difference is due to higher bolometric luminosities in the higher redshift quasars (see Figure 7).We conclude that we detect no evidence of impact of AGN feedback on the CO SLEDs up to J ≤ 7 for our quasar sample, despite the strong feedback signatures that characterise them (i.e., a sample with prevalent radio jets and/or shocked winds and ionised outflows).
-We observe differences between the CO and [O III] line profiles, both in the line widths and velocity offsets, finding systematically broader [O III] line profiles than CO.The median difference in  80 between [O III] and CO is ∼ 200 km s −1 , with a maximum difference of ∼ 650 km s −1 .This suggests a larger impact of feedback on the ionised gas than on the molecular gas (see Figure 9).Alternatively this can indicate cold gas clouds are unable to survive in hot winds.
-We identify consistent [O III] line profiles in SDSS data compared to MUSE data extracted at a 3 arcsec aperture.However, differences in the line profiles are identified when extracted from MUSE data at a larger aperture of 30 arcsec.This suggests that the effects of feedback processes (such as outflows, radio jets or winds) are likely more dominant at smaller scales, closer to the central AGN/quasar (see Figures 10a and 10b. Overall we conclude that in these sample of quasars at  < 0.2 the impact of these quasars on the total molecular gas content, both in excitation and velocities, is likely to be minimal.On a global scale we see no real divergences from ULIRGs.This work therefore adds to the growing body of evidence that on global scales there is a minimal impact on CO excitation and total gas content, even in the extreme cases of luminous quasars with ionised outflows and extended radio structures.However, we note that on smaller scales an increased velocity dispersion (Girdhar et al. 2022) and increased line ratios (Audibert et al. 2023) for two targets in our sample plus displaced molecular gas in another two targets (Girdhar et al, submitted), have been previously identified with a spatial relation to the observed radio jets.The question remains as to whether this impact is seen across the entire sample and further resolved studies will shed light on the the impact on the multi-phase ISM, in particular further investigation into targets with well characterised radio emission.

Figure 1 .
Figure 1.Sample Selection: [O III] luminosity versus emission line widths, demonstrating how the 42 quasars in the QFeedS sample (star symbols) are selected from  < 0.2 AGN sample (Mullaney et al. 2013, blue points and contours).The 17/42 selected for this work are highlighted by the larger coloured stars and the remaining 25/42 of the QFeedS sample are shown by small grey stars.The colours for the 17 targets in this sample (shown in the legend) are carried through in figures throughout the rest of this paper.The dashed line represents the selection criteria in QFeedS of L [O III] > 10 42.1 ergs s −1 .

Figure 2 .
Figure 2.For the 7 targets with S/N > 5 in CO(1-0) we plot the flux extracted from various diameter apertures of the ACA data up to 40 arcsec.The flux presented here is normalised to 30 arcsec which is where we extract our flux to match with APEX data (dashed lines).Uncertainties on each flux estimate are presented with grey error bars.

Figure 3 .
Figure 3. Example of multiple CO spectra and MUSE data obtained, here showing J1100+0846.Top panel: CO(1-0) from ACA. Second panel CO(2-1) from APEX (in this case the data was from Jarvis et al. 2020).Third panel: CO(3-2) data from APEX.Fourth panel: CO(6-5) data from APEX.Fifth panel: MUSE spectra for the [O III] line extracted from a 30 arcsec diameter aperture.In all cases, solid black lines denote fits to the data.Here there are multiple components and as such, the dotted black lines denote the different components that make the total fit to the spectra.In the CO(2-1) spectra we also show the fit from higher resolution ALMA observations (solid orange line labelled RA+22, RamosAlmeida et al. 2022).Shaded grey regions represent the 1 level.Spectra for all targets are shown in the appendix, using the same presentation methods.MNRAS 000,1-21 (2023)

Figure 4 .
Figure 4.For all targets with data for both stellar mass (M ★ ) and CO gas masses (M CO ) we present the stellar mass versus the gas fraction (M CO /M ★ ).Large coloured stars show targets from this work, with colours as in Figure1.Background small green circles and pink squares show a compiled literature sample of AGN and non-AGN fromTacconi et al. 2018.Here 16/17 targets are presented, with J1116+2200 missing due to a lack of stellar mass information.
(a) For the 9 targets which have detections in CO(1-0) we present the CO SLEDs, normalised to CO(1-0) with Milky Way and thermalised relations for comparison (Carilli & Walter 2013).Also shown are  = 1 -2 Starburt and main sequence (MS) CO SLEDs from Valentino et al. (2021) (V21).The CO SLEDs are colour coded by target (as in Figure 1).3 sigma upper limits are indicated by empty circles with downward arrows.No uncertainties are shown for CO(1-0) as these are incorporated into the uncertainties on CO(2-1) and CO(3-2) to best represent the uncertainty on the excitation level.(b) For the 7 targets which have either CO(6-5) or CO(7-6) data we present the CO SLEDs up to these higher transitions.Here we are normalising to the CO(2-1).Literature comparison data showing high- quasars Carilli & Walter (2013) (C&W13) and  = 1 -2 Starburst and main sequence (MS) Valentino et al. (2021) (V21).3 sigma upper limits are indicated by empty circles with downward arrows.No uncertainties are shown for CO(2-1) as these are incorporated into the uncertainties on CO(1-0), CO(3-2) and CO(6-5) to best represent the uncertainty on the excitation level.

Figure
Figure 5. CO SLEDs Table ofmedian line ratios along with 16th and 84th quartiles (indicated by plus and minus values) for all targets and comparisons to the literature used in Figure 6 along with the  range for each sample.Montoya Arroyave et al. 2023, Leroy et al. 2022a, Greve et al. 2014, Lamperti et al. 2020, Carilli & Walter 2013.

Figure 7 .
Figure 7.For the 7 targets which have either CO(6-5) or CO(7-6) data we compare the bolometric luminosity and line ratios of a sample of high- quasars (from Carilli & Walter 2013, shown by red diamond markers) to our sample of low-z quasars (shown by star markers and arrows indicating those that are upper limits).The numbers in each marker indicates the line ratio being plotted, e.g."61" indicates r 61 .Markers of the QFeedS targets are also coloured in the same way as in all previous plots (see Figure 1).

Figure 8 .
Figure 8.For the 5 targets which have detections in all three transitions r 21 plotted vs r 32 to give an indication of temperature, density and optical thickness.Background grey points taken from Leroy et al. 2022a.

Figure 9 .
Figure 9.For the 11 targets which have lines with S/N > 5 we present an analysis of the velocity offsets ( 50 ) and line widths ( 80 ) of the CO and [O III] spectra.Colours of markers for the individual targets are the same as throughout this work (colours shown in the legend of Figure 1).Some targets have multiple transitions shown, which are differentiated by: Upward triangles indicate CO(1-0), squares are CO(2-1) and downward triangles are CO(3-2) data (see legend).Each region of the plot is shaded to highlight the corresponding properties and the differences between the [O III] and the CO.Note: J1108+0659 is H not [O III].Also note that the CO data for J1356+1026 (yellow square) is taken from Ramos Almeida et al. 2022.
(a) Difference between SDSS [O III]  80 and [O III]  80 extracted from MUSE at 3 arcsec, as a function of [O III]  80 from SDSS.(b) Difference between SDSS [O III]  80 and [O III]  80 extracted from MUSE at 30 arcsec, as a function of [O III]  80 from SDSS.

Figure 10 .
Figure 10.Comparing line widths of the [O III] spectra extracted from SDSS and MUSE data cubes at two different extraction radii.

Figure B1 .
Figure B1.J0909+1052, non-detections across all three CO transitions.All shown with 100km s −1 bins.MUSE [O III] observations show a single with Gaussian.Figure B2.J0945+1737 spectra across four CO transitions.All with 100 km s −1 bins.Detection in CO(2-1) and low S/N detection in CO(1-0), nondetection in CO(3-2) (although tentative signs of a line at S/N = 2.7) or CO(6-5).MUSE [O III] observations show a two component fit with a narrow and broad component.

Figure B5 .
Figure B5.J1010+0612 spectra across three CO transitions, with no CO(1-0) observations available.All with 100 km s −1 bins.Detection in CO(3-2) and a low S/N detection in CO(2-1) which matches the two component line profile identified in Ramos Almeida et al. 2022 (spectra shown in orange, labelled RA+22) for which we see tentative signs here.MUSE [O III] observations show a two component fit with a narrow and broad component.

Figure B9 .
Figure B9.J1100+0846 spectra across four transitions.CO(1-0) with 50 km s −1 bins and CO(2-1) and CO(3-2) with 25 km s −1 bins.Detection in all four transitions, all with same double peak profile.The same double peaked profile was found by Ramos Almeida et al. 2022 CO(2-1) observations (spectra shown in orange, labelled RA+22).MUSE [O III] observations show a two component fit with a narrow and broad component.

Figure B10 .
Figure B10.J1108+0659 spectra across four CO transitions.First three CO transitions plotted with 50 km s −1 bins and CO(7-6) plotted with 100 km s −1 bins due to a lack of detection.Blue wing in the first three CO transitions indicating a potential molecular outflow component.There is also potentially a small velocity shift in this outflow component across the first three CO transitions.We find a non-detection in CO(7-6).MUSE [O III] observations showing a similar blue wing to the CO data.

Figure B12 .
Figure B12.J1116+2200 spectra across four CO transitions.The first three transitions with 50 km s −1 bins and the CO(7-6) with 100 km s −1 bins.Detection of a broad line in the first three transitions.Non-detection in CO(7-6).MUSE [O III] observations show a single with Gaussian line profile.

Figure B13 .
Figure B13.J1222-0007 non-detections in all three CO transitions.MUSE [O III] observations show a two component fit to the data.

Figure B14 .
Figure B14.J1316+1753, low S/N detection in CO(2-1) and non-detection in CO(1-0) and CO(3-2).MUSE [O III] observations show a three component fit to the data.The low S/N detection in our CO(2-1) data seems to match the middle component present in the MUSE [O III] spectra.

Figure B15 .
Figure B15.J1356 spectra across two CO transitions with no CO(1-0) data available.All set to 100 km s −1 bin widths.Non-detection in both transitions.Also shown is the detection from Ramos Almeida et al. 2022 detection in CO(2-1) (labelled RA+22 N in orange and RA+22 in green, where orange represents the nucleus only and green is the full line).We use the RA+22 data in our analysis as well lack a detection in CO(2-1) .MUSE [O III] observations show a blue wing component.

Figure B17 .
Figure B17.J1518+1403 spectra across four CO transitions.All with 50 km s −1 bins (except for CO(7-6) which is 100 km s −1 due to lack of a detection).Detections in both CO(1-0) and CO(3-2) and a low S/N detection in CO(2-1).MUSE [O III] observations show a two component fit with a narrow and broad component.

Figure C1 .
Figure C1.Contour plots for the 9 targets with an integrated S/N > 3 detection in the ACA CO(1-0) data, where red contours show the 2  levels and upwards.The yellow circle denotes the beam size of the CO(2-1) observations, which is also the aperture used to extract the CO(1-0) spectra.The background images are rgb images from the DESI Legacy Imaging Survey in the (z, r, g) bands.