© ESO 2021

1. Introduction

The active phases supermassive black holes (SMBHs) go through while accreting material, when they are visible as active galactic nuclei (AGN), are thought to play a key role in galaxy evolution. During these phases, the central engine releases a huge amount of energy, which is injected into the surrounding interstellar medium (ISM). This energy, if coupled efficiently, could remove, heat, and/or dissociate the molecular gas, the fuel out of which stars form, as well as heat up the halo and suppress further gas accretion. Such a process, referred to as AGN feedback, may regulate the growth of the host galaxy (Fabian 2012; Kormendy & Ho 2013; King & Pounds 2015; Harrison 2017). Theoretically, AGN feedback is invoked to reproduce the observed properties of the galaxy population, for example the lack of very massive galaxies in the most massive galaxy haloes (Somerville et al. 2008), the bimodal color distribution of galaxies (Strateva et al. 2001), and the correlations between SMBH mass and host galaxy properties (Magorrian et al. 1998; Ferrarese & Merritt 2000). Although AGN feedback is a necessary ingredient in models of galaxy evolution, proving its role observationally remains a challenge. In particular, powerful outflows have been observed in many AGN host galaxies (e.g., Rupke & Veilleux 2013; Cicone et al. 2014; Genzel et al. 2014; Harrison et al. 2014; Carniani et al. 2015; Cresci et al. 2015; Kakkad et al. 2016; Nesvadba et al. 2017; Brusa et al. 2018; Davies et al. 2020; Veilleux et al. 2020), but their long-term impact on the global star-forming activity of the hosts is still an open issue (e.g., Cresci & Maiolino 2018; Gallagher et al. 2019; Scholtz et al. 2020).

To understand the link between AGN and star formation, several authors have investigated the star-formation rates (SFRs) of AGN host galaxies as a function of AGN luminosity (Lutz et al. 2010; Harrison et al. 2012; Page et al. 2012; Rosario et al. 2012; Rovilos et al. 2012; Delvecchio et al. 2015; Stanley et al. 2015; Harris et al. 2016). Many of the apparently discrepant results in these works can be explained by selection effects (e.g., controlling for stellar mass and redshift evolution; e.g., Stanley et al. 2017). However, while the observational results are broadly consistent with simple feedback models (e.g., Harrison 2017), they currently provide limited diagnostic power on the specifics of how AGN feedback works, at least in the absense of a wide range of theoretical predictions (e.g., Scholtz et al. 2018; Schulze et al. 2019).

Observations of cold molecular gas are more promising since any activity from the AGN (radiation, outflows, or jets) will have an impact on the molecular gas reservoir first, and then on the SFR, which previous studies have focused on. The molecular gas provides an instantaneous measure of the raw fuel from which stars form and can be used as a more direct tracer of potential feedback effects. Over the last decade, large observational efforts have been devoted to map the molecular gas reservoir of galaxies (e.g., Daddi et al. 2010a; García-Burillo et al. 2012; Bauermeister et al. 2013; Bothwell et al. 2013; Tacconi et al. 2013, 2018; Genzel et al. 2015; Silverman et al. 2015; Decarli et al. 2016; Cicone et al. 2017; Saintonge et al. 2017; Fluetsch et al. 2019; Freundlich et al. 2019). These studies are largely based on observations of carbon monoxide (CO) rotational emission lines, used as a tracer of cold molecular hydrogen H₂ (the ground-state rotational transition in particular traces the total molecular gas best), but there are also many studies based on dust emission (e.g., Tacconi et al. 2018 and references therein). The main targets of such CO campaigns have primarily been inactive¹ galaxies that mostly lie on the main sequence of star-forming galaxies (e.g., Noeske et al. 2007; Schreiber et al. 2015), where the majority of the cosmic star-formation activity occurs. The fundamental relation between SFR and the molecular gas content of galaxies, the Schmidt-Kennicutt relation, provides precious information about how efficiently galaxies turn their gas into stars (Schmidt 1959; Kennicutt 1989). This star-formation law is usually presented in terms of surface densities, therefore requiring spatially resolved measurements of galaxies. However, in high-redshift studies, an integrated form of this relation with global measurements of the SFR and molecular mass is normally used (e.g., Carilli & Walter 2013; Sargent et al. 2014).

A similar observational effort is needed to systematically characterize the cold gas phase of AGN host galaxies to understand whether this is different from inactive galaxies and quantify any potential effect of AGN feedback on the host galaxy ISM. Local studies agree in reporting no clear evidence for AGN to affect the ISM component of the host, by tracing the molecular gas phase (e.g., Saintonge et al. 2012, 2017; Husemann et al. 2017; Rosario et al. 2018; Jarvis et al. 2020; Shangguan et al. 2020b), the atomic one (e.g., Ellison et al. 2019), dust mass and dust absorption as a proxy of the gas mass (e.g., Shangguan & Ho 2019; Yesuf & Ho 2020). AGN appear to follow the same star-formation law as inactive galaxies. Interestingly, studies at redshift z > 1 present contrasting results. Since the redshift 1 < z < 3, the so-called cosmic noon, corresponds to the peak of accretion activity of SMBHs (Madau & Dickinson 2014; Aird et al. 2015), when the energy injected into the host galaxy may be maximized, this cosmic epoch is a crucial laboratory to look for AGN feedback effects. In particular, some studies find reduced molecular gas fractions (i.e., the molecular gas mass per unit stellar mass, f_gas = M_mol/M_*) and depletion timescales (i.e., the timescale needed to consume the available molecular gas content given the current SFR, t_dep = M_mol/SFR) of AGN compared with the parent population of inactive galaxies (Carilli & Walter 2013; Brusa et al. 2015; Carniani et al. 2017; Kakkad et al. 2017; Fiore et al. 2017; Perna et al. 2018; Talia et al. 2018; Loiacono et al. 2019, but see also Kirkpatrick et al. 2019; Herrera-Camus et al. 2019; Spingola et al. 2020). This has been interpreted as evidence for highly efficient gas consumption possibly related to AGN feedback affecting the gas reservoir of the host galaxies. For a few targets, fast ionized/molecular outflows were also detected (e.g., Brusa et al. 2015, 2018; Carniani et al. 2017). Complementary observations tracing the possible presence of outflows are therefore needed to generally confirm such conclusions (e.g., Vayner et al. 2017; Brusa et al. 2018; Herrera-Camus et al. 2019; Loiacono et al. 2019).

Nevertheless, these studies can be affected by assumptions and limitations. CO measurements at z > 1 are usually performed using various high-J transitions, and excitation corrections are needed to estimate the luminosity of the ground-state transition. The CO spectral line energy distribution (SLED) of a given target is rarely known (e.g., Weiß et al. 2007; Carilli & Walter 2013; Mashian et al. 2015) and therefore the excitation correction is typically highly uncertain. In addition, calculating gas masses from CO luminosities ($L_{\text{CO}}^{\prime }$) requires the assumption of a conversion factor α_CO, that depends on the conditions of the ISM and typically ranges between 0.8 and 4 M_⊙/(K km s⁻¹ pc²) for solar-metallicity galaxies (Carilli & Walter 2013; Accurso et al. 2017). When dealing with SFRs of AGN hosts, an additional complication is the difficulty to properly account for the AGN contribution to the far-infrared (FIR) luminosity. As recently shown by Kirkpatrick et al. (2019), different methods to estimate the AGN contribution can lead to completely different results and place the AGN population on the same star-formation law of inactive galaxies. Finally, AGN samples at high redshift are usually small and likely biased toward brighter objects (e.g., Brusa et al. 2018) or are heterogeneous (e.g., different selection criteria, CO lines observed) when assembled from literature data (e.g., Fiore et al. 2017; Perna et al. 2018; Kirkpatrick et al. 2019; Bischetti et al. 2021).

In this paper we present the first systematic and uniform analysis of the CO(3-2) emission of 27 X-ray selected AGN at z ≈ 2 to infer whether their activity affects the CO properties of the host galaxy. The targeted transition is the lowest-J transition accessible with ALMA at z ≈ 2. In this work, we controlled for many of the biases mentioned above by performing a uniform selection and multiwavelength characterization of a fairly representative sample of X-ray AGN. The sample includes less extreme sources, which cover a wider range in AGN bolometric luminosity and uniform distribution on the main sequence. We compared the CO emission properties of our AGN with those of inactive galaxies. The comparison sample was constructed as carefully as possible, taking into account stellar masses, SFRs, redshift and observed CO transition in order to minimize systematic differences. In particular, by considering targets observed in the same CO transition, we avoid assumptions on conversion factors, which could bias our conclusions.

The paper is structured as follows: in Sect. 2 we describe the sample selection criteria and the multiwavelength properties of our AGN; the sample of inactive galaxies used for comparison is presented in Sect. 3; in Sect. 4 we outline the ALMA observations and data analysis; we present and discuss our results in Sect. 5 and 6, respectively. Our conclusions are given in Sect. 7. In this work we adopted a WMAP9 cosmology (Hinshaw et al. 2013), H₀ = 69.3 km s⁻¹ Mpc⁻¹, Ω_M = 0.287 and Ω_Λ = 0.713.

2. The AGN sample

2.1. Target selection

As a tracer of AGN activity we used X-ray emission, which is very efficient thanks to the low contamination from the host galaxy. We combined catalogs from a deep and small-area as well as a shallow and wide-area X-ray survey, in order to include in our sample both high- and low-luminosity AGN and cover a wide range in AGN bolometric luminosity L_bol (Brandt & Alexander 2015; Circosta et al. 2018). The targets were drawn from the following surveys, by adopting as a threshold an absorption-corrected 2–10 keV X-ray luminosity L_X ≥ 10⁴² erg s⁻¹. (i) The COSMOS-Legacy survey (Civano et al. 2016; Marchesi et al. 2016), a 4.6 Ms Chandra observation of the COSMOS field, with a deep exposure over an area of about 2.2 deg² at a limiting depth of 8.9 × 10⁻¹⁶ erg cm⁻² s⁻¹ in the 0.5−10 keV band. (ii) The wide-area XMM-Newton XXL survey North (Pierre et al. 2016), a ∼25 deg² field surveyed for about 3 Ms by XMM-Newton, with a sensitivity in the full 0.5−10 keV band of 2 × 10⁻¹⁵ erg cm⁻² s⁻¹. The XMM-XXL subsample is culled from the X-ray source catalog presented by Liu et al. (2016) using the X-ray reduction pipeline described by Georgakakis & Nandra (2011), and with optical spectroscopy published by Menzel et al. (2016).

These fields are covered by a rich multiwavelength set of ancillary data spanning from the X-rays to the radio regime, which are essential to obtain robust measurements of the target properties. Our AGN were then selected to have spectroscopic redshift in the range z = 2.0−2.5, whose quality was flagged as “secure” in the respective catalogs as well as a coverage of AGN and galaxy properties as wide and uniform as possible. Overall, our sample consists of 27 AGN, namely 7 from XMM-XXL and 20 from COSMOS. IDs, coordinates and redshift of the sources are reported in Table 1.

The ALMA program presented in this paper is part of SUPER (SINFONI Survey for Unveiling the Physics and Effect of Radiative feedback; Circosta et al. 2018, PI: Mainieri), which aims at studying AGN feedback at cosmic noon. SUPER started as an ESO’s VLT/SINFONI Large Programme to perform the first systematic investigation of ionized outflows in a sizeable and blindly-selected sample of 39 X-ray AGN at z ≈ 2, which reaches high spatial resolutions (∼2 kpc) thanks to the adaptive optics-assisted integral field spectroscopy observations. The goal is to have a reliable overview of the ionized gas component as well as AGN-driven outflows traced by [O III]λ5007Å. In parallel, we ran an ALMA campaign to trace the molecular gas properties of the host galaxies. The target sample was initially similar for the SINFONI and ALMA programs. However, due to scheduling constraints on the SINFONI observations the two programs evolved differently, but 43% of the ALMA targets feature good-quality SINFONI observations. Moreover, SINFONI and ALMA targets share the same properties and selection criteria (see Circosta et al. 2018). In Table 1, the targets with good-quality SINFONI data from our Large Program are flagged (Kakkad et al. 2020; Vietri et al. 2020; Perna et al. in prep.). This distinction is made for future reference since we defer the comparison of outflow and molecular gas properties obtained from the SINFONI and ALMA datasets, respectively, to a future paper. In this work, we focus on the analysis of the CO(3-2) properties of our AGN.

2.2. Multiwavelength properties of the sample

We characterized the physical properties of our sources by exploiting the multiwavelength coverage available from the X-ray to the radio regime. In particular, we followed the procedure described in Circosta et al. (2018) to collect the multiwavelength data and derive the properties of the targets that were not analyzed in our previous work, through X-ray spectral analysis and broad-band spectral energy distribution (SED)-fitting. We refer the reader to Circosta et al. (2018) for an extensive description, but we summarize in the following some key pieces of information.

The counterparts to the X-ray sources in COSMOS are provided along with the optical-to-MIR multiwavelength photometry in the COSMOS2015 catalog (Laigle et al. 2016). We complemented this dataset with FIR data from Herschel/PACS and SPIRE, when available, using a positional matching radius of 2″ (we note that we used 24 μm-priored catalogs, which in turn are IRAC-3.6 μm priored). Spitzer/MIPS photometry at 24 μm and Herschel/PACS at 100 and 160 μm was taken from the PEP DR1(Lutz et al. 2011). Herschel/SPIRE photometry at 250, 350 and 500 μm was retrieved from the data products presented in Hurley et al. (2017)². We also added ALMA continuum data in Band 7 and Band 3 available from the COSMOS A3 catalog (Liu et al. 2019), Scholtz et al. (2018) and our ALMA dataset (this work and Lamperti et al., in prep.).

As for XMM-XXL, the counterparts to the targets are known in the SDSS optical images and the corresponding associations to the X-ray sources, as well as the photometry from UV-to-mid-IR, are given in Fotopoulou et al. (2016), Menzel et al. (2016) and Georgakakis et al. (2017). Herschel/PACS and SPIRE data are those released by the HerMES collaboration in the Data Release 4 and 3 respectively (Oliver et al. 2012)³.

Mid- and far-IR photometric datapoints from 24 to 500 μm were considered as detections when the signal-to-noise ratio S/N > 3, where the total noise was given by the sum in quadrature of both the instrumental and the confusion noise (Lutz et al. 2011; Oliver et al. 2012; Magnelli et al. 2013). The detections below this threshold were converted to 3σ upper limits. All the data used in this work were corrected for Galactic extinction (Schlegel et al. 1998).

The datasets described above were modeled by using the Code Investigating GALaxy Emission (CIGALE⁴; Boquien et al. 2019; Noll et al. 2009), a publicly available state-of-the-art galaxy SED-fitting technique (we used version 0.11.0). This code disentangles the AGN contribution from the emission of the host by adopting a multicomponent fitting approach and includes attenuated stellar emission, dust emission heated by star formation, AGN emission (both primary accretion disk emission and dust-heated emission), and nebular emission. Moreover, it takes into account the energy balance between the UV-optical absorption by dust and the corresponding reemission in the FIR. The parameters of interest (namely stellar masses, FIR luminosities and AGN bolometric luminosities) and their uncertainties are determined by the code through a Bayesian statistical analysis, by building up a probability distribution function (PDF) that takes into account the different models used during the fit. The fiducial results correspond to the mean value of the PDF and the associated error is the standard deviation (Noll et al. 2009). To reproduce the stellar component we used: Bruzual & Charlot (2003) stellar population models with solar metallicity, a delayed star-formation history, the modified Calzetti et al. (2000) attenuation law, and a Chabrier (2003) initial mass function (IMF). The contributions from star-forming dust and AGN were modeled with the libraries presented by Dale et al. (2014) and Fritz et al. (2006), respectively. Finally, the nebular emission was reproduced with the models of Inoue (2011). Tables 1 and 2 in Circosta et al. (2018) provide a list of the photometric filters and input parameters of the models used for the SED-fitting procedure.

The results of the SED-fitting analysis are reported in Table A.1, together with the optical spectroscopic classification in broad-line (BL) and narrow-line (NL) AGN, depending on the presence of broad (FWHM > 1000 km s⁻¹) or narrow (FWHM < 1000 km s⁻¹) permitted emission lines in their optical spectra, respectively. From now on, we will refer to type 1 and type 2 AGN according to the optical spectroscopic classification. We refer the reader to Circosta et al. (2018) for the SEDs of most of our sample. The SEDs of the targets that were not previously presented (flagged in Table 1), are shown in Appendix C. When possible, SFRs were derived by assuming the Kennicutt (1998) calibration, converted to a Chabrier (2003) IMF by subtracting 0.23 dex, from the IR luminosity integrated in the rest-frame wavelength range 8–1000 μm, after removing the AGN contribution (11 targets). For non-detections at observed λ > 24 μm, we provide a 3σ upper limit on the SFR, derived as the 99.7th percentile from the PDF of the FIR luminosity (11 targets). When data were not available at observed λ > 24 μm, the dust templates were not included in the fitting routine. We therefore provide the SFR as derived from the modeling of the stellar component in the UV-to-near-IR(NIR) regime with SED fitting, averaged over the last 100 Myr of the galaxy star-formation history (two targets, flagged in Table A.1). Three targets without data coverage at λ > 24 μm are also dominated by the AGN emission at UV-to-NIR wavelengths (namely X_N_128_48, X_N_102_35, and X_N_104_25). For these targets we could not derive an estimate of SFR and stellar mass and therefore only the AGN bolometric luminosity is reported in Table A.1⁵. These three targets are excluded from the analysis where estimates of stellar masses and SFRs are required (see Sect. 5). Moreover, measurements of stellar mass for X_N_53_3, X_N_6_27 and cid_1605 could not be recovered because of the strong AGN contamination in the optical/NIR regime, but we report their SFRs. We measured stellar masses in the range log(M_*/M_⊙) = 9.6−11.2, FIR luminosities log (L_FIR / erg s⁻¹) < 45.0 − 46.4, SFRs < 25−686 M_⊙ yr⁻¹, and AGN bolometric luminosities mainly in the range log (L_bol / erg s⁻¹) = 44.7 − 46.9 covering two orders of magnitude. In Table A.1 we also report black hole masses, collected from the literature (Menzel et al. 2016) and measured from our SINFONI dataset (Vietri et al. 2020) as well as spectra available from SDSS and the COSMOS survey. The results range between 10^7.64 and 10^10.04 M_⊙. Overall, our sample is composed of 17 type 1 and 10 type 2 AGN.

From the analysis of the X-ray spectra we derived obscuring column densities N_H and X-ray luminosities L_X. Again, we followed the same procedure described in Circosta et al. (2018) to extract and fit the spectra. The energy ranges considered are 0.5−7 keV band for Chandra and 0.5−10 keV band for XMM-Newton. All the fits were performed by using XSPEC v. 12.9.1⁶ (Arnaud 1996) and adopting the Cash statistic (Cash 1979) and the direct background option (Wachter et al. 1979). The spectra were binned to 1 count per bin to exclude empty channels. For sources with more than 30 (50) net counts (reported in Table A.1) for Chandra (XMM-Newton), we used a simple spectral model, consisting of a power law modified by both intrinsic and Galactic absorption, as well as a secondary power law to reproduce any excess in the soft band (0.5–2 keV), due to scattering or partial covering in obscured sources. The photon index was left free to vary (typical values within Γ = 1.5−2.5) for spectra with more than ∼100 net counts, otherwise we fixed it to the canonical value of 1.8 (e.g., Piconcelli et al. 2005). For targets with less than 30 (50) counts N_H values were derived at the source redshift from hardness ratios ($HR=\frac{H-S}{H+S}$, where H and S are the number of counts in the hard 2−7 keV and soft 0.5−2 keV bands, respectively; Lanzuisi et al. 2009). In both cases, we propagated the uncertainty on N_H when deriving the errors on the intrinsic luminosity. The results derived for column densities and 2−10 keV absorption-corrected luminosities are listed in Table A.1. Our sample spans a range of obscuration properties including unobscured targets (N_H < 10²² cm⁻²) up to Compton-thick ones (N_H > 10²⁴ cm⁻²) and X-ray luminosities in the range log (L_{[2−10 keV]} / erg s⁻¹) = 43.0 − 45.4.

3. Comparison sample of inactive galaxies

We built a comparison sample of inactive galaxies from the Plateau de Bure high-z Blue Sequence Survey (PHIBSS; Tacconi et al. 2018). PHIBSS aims at investigating the gas properties of galaxies across cosmic time by using a sample of 1444 targets spanning a redshift range z = 0−4.4 with estimates of molecular gas mass given by CO measurements, FIR SED dust measurements or 1mm dust measurements. In particular, in the redshift range of our interest (z = 2−2.5) PHIBSS galaxies were selected to have stellar mass ≥2.5 × 10¹⁰ M_⊙ and SFR ≥ 30 M_⊙ yr⁻¹, in order to ensure a relatively uniform coverage of the SFR-stellar mass relation. At z ≈ 2, the PHIBSS sample contains CO(3-2) measurements obtained with IRAM PbBI and the updated IRAM NOEMA (Tacconi et al. 2013; Freundlich et al. 2019), as well as CO measurements of main sequence star-forming galaxies and submm galaxies (SMGs) collected from the literature (Tacconi et al. 2006; Bothwell et al. 2013; Saintonge et al. 2013; Silverman et al. 2015; Decarli et al. 2016), performed with both ALMA and NOEMA. From the PHIBSS parent sample, we selected for our purposes galaxies with z = 2−2.5 (the same redshift range of our targets), whose molecular gas mass was derived through CO(3-2) observations. By requiring that also the control sample is observed in CO(3-2), we can directly compare the observed CO luminosity and we are free from assumptions on conversion factors (e.g., the excitation correction and α_CO), which could bias the results. All stellar masses and SFRs assume a Chabrier IMF. The typical systematic uncertainties on both the stellar mass and the SFR were assumed to be 0.2 dex (Freundlich et al. 2019). As for the SFRs, it is worth noting that they were mainly derived from IR measurements, or a combination of IR and UV. However, a small subsample of galaxies feature SFR estimates from Hα fluxes (Tacconi et al. 2013). Since we want to compare CO and FIR luminosities of these galaxies, we relied on the agreement between SFRs derived from Hα and IR in main-sequence galaxies at cosmic noon (Rodighiero et al. 2014; Puglisi et al. 2016; Shivaei et al. 2016) and converted all SFRs into FIR luminosities, using the Kennicutt (1998) relation corrected for a Chabrier (2003) IMF (i.e., by subtracting 0.23 dex). We used CO fluxes retrieved from the corresponding papers, when data were taken from the literature, or provided by the PHIBSS collaboration (L. Tacconi, private communication), since the survey paper by Tacconi et al. (2018) provides only the final estimates of gas masses. We looked for the presence of AGN in the sample, by checking the corresponding papers the targets are from (see Tacconi et al. 2018). They were identified by using power-law NIR component fraction (i.e., rest-frame NIR excess in their SED), color selection, X-ray emission and line ratios. Six targets were therefore excluded. We did not include these targets in our AGN sample because (i) the methods used to identify AGN were not the same as ours and (ii) stellar masses and SFRs were not derived by accounting for the AGN contribution hence these parameters may be overestimated. The inactive galaxies satisfying the selection criteria mentioned above cover the same range of stellar mass of our AGN sample within the uncertainties (log M_* ≈ 9.5−11.5). From this set of objects we further excluded those with extreme SFRs, that is the four targets with SFRs larger than the highest value in our sample per bin of stellar mass (width of 0.5 dex; Fig. 1), in order to have a similar distribution on the main sequence. We note that even without excluding such targets from the comparison sample, the results presented in Sect. 5 would not change significantly. The final sample selected from PHIBSS is made of 42 objects. The distribution on the main sequence of our AGN and the comparison sample is shown in Fig. 1. Although we are unavoidable limited by 50% SFR upper limits in our AGN sample considered in Fig. 1, the comparison sample has a broadly similar distribution on the main sequence.

Fig. 1. Distribution of host galaxy properties in the SFR-M* plane for the 21 AGN (type 1s marked by triangles and type 2s marked by circles) in our sample with an estimate of both stellar mass and SFR. The comparison sample of inactive galaxies is depicted by blue squares. The two data points with orange edges represent the targets with SFR derived through modeling of the stellar emission using SED fitting. Black arrows represent 3σ upper limits. The color coding indicates the AGN bolometric luminosity for our sample, which covers two orders of magnitude in Lbol. The black solid line reproduces the main sequence of star-forming galaxies from Schreiber et al. (2015) at the average redshift of our target sample (i.e., ∼2.3). The dashed lines mark the scatter of the main sequence (equal to 0.3 dex). The histograms show the projected distribution of the two quantities along each axis, in red for our sample and blue for the comparison one.

4. ALMA data: Observations and analysis

The ALMA program was carried out in Cycle 4 and 5 (Project codes: 2016.1.00798.S and 2017.1.00893.S; PI: V. Mainieri) using 42–47 antennas and maximum baselines between 704 m and 1.1 km. Observations were taken in Band 3 between November 2016 and May 2017 (Cycle 4) and in March 2018 (Cycle 5). In each observation, one spectral window, with a bandwidth of 1.875 GHz, was centered at the expected redshifted frequency of the CO(3-2) emission line (the rest-frame frequency is 345.8 GHz) based on the spectroscopic redshift of the target (Table 1), while the other three spectral windows were used to sample the continuum emission. Each spectral window was divided into 240 channels and had a native spectral resolution of 7.8 MHz, corresponding to 21–23 km s⁻¹ at the line frequency. The on-source exposure time varied between ∼9 and ∼80 min.

ALMA visibilities were calibrated using the CASA software⁷ version 4.7.0 for Cycle 4 data and 5.1.1 for Cycle 5, as originally used for the reduction with the pipeline. The presence of continuum emission was estimated from the continuum maps generated with TCLEAN by averaging all spectral windows, and subtracted in the uv plane by using the UVCONTSUB task within CASA, by excluding the spectral range where we expected the line. Continuum emission at ∼100 GHz was detected for five targets, as shown in Appendix E. For these targets, a fitting procedure using a two-dimensional Gaussian was performed on the continuum maps in the image plane to estimate the flux densities, reported in Table B.1. The final continuum-subtracted datacubes were generated with the CASA task TCLEAN in velocity mode, using “natural” weighting, which optimizes the point-source sensitivity in the image plane, the “Hogbom” cleaning algorithm, a cellsize of 0.2″, and a velocity bin width of 24 km s⁻¹. As a threshold, we used 3 times the rms derived from the dirty cubes (reported in Table B.1). The emission-line cubes have angular resolutions in the range 0.7–1.7″, and 1σ rms in the range 0.18–0.72 mJy beam⁻¹ per 24 km s⁻¹ velocity bin, as listed in Table B.1.

To determine the flux and size of the CO(3-2) line we proceeded as follows. For each target, we considered the continuum-subtracted cube and created moment 0 maps (i.e., velocity-integrated flux maps) by collapsing, with the task IMMOMENTS, spectral channels around the expected CO(3-2) frequency over increasingly larger velocity widths, from 100 km s⁻¹ to 1000 km s⁻¹ in steps of 100 km s⁻¹ (in case of successful line detection, the procedure was repeated by collapsing spectral channels around the peak of the emission line). The rms of the collapsed CO maps was estimated over an area approximately half that of the primary beam. We then assessed the significance of the emission line in the moment 0 maps. If the line was detected (i.e., S/N ≳ 3), we extracted the one-dimensional spectrum from the emission-line cube using as extraction area the region above 2σ significance in the moment 0 map (∼1–2 arcsec, centered on the target; see contours in Appendix D). Among the different integrated maps produced, we considered the one where the line had the highest significance. The one-dimensional line spectrum was fitted using a single Gaussian profile (python function CURVE_FIT in scipy.optimize) to estimate parameters such as line centroid (from which z_CO, Table B.1), FWHM and integrated flux density. The uncertainties on fluxes, FWHM, z_CO, and peak fluxes were estimated by adopting a Monte Carlo approach. We created 100 mock spectra by adding to the model spectrum random noise proportional to the noise measured per channel. To estimate the noise in each channel of the spectrum we extracted 100 spectra from random regions within the cube. Such regions had an area similar to the source extraction region and were not located too close to the source. We considered the standard deviation of the flux densities, in each spectral channel, of the 100 random spectra as an indication of the rms per channel, and used these values in the Monte Carlo routine. We note that the rms resulted to be quite uniform throughout the spectrum. We performed the fit for each of our mock spectra and considered the standard deviation of the distribution of the resulting parameters as uncertainties. As for the integrated fluxes, the final uncertainties are given by adding in quadrature a typical ALMA flux calibration error equal to 5 per cent of the flux, as presented by Bonato et al. (2018), to the uncertainties obtained from the Monte Carlo routine. We detect the CO(3-2) line in 11 targets out of 27 observed. The line is considered detected if the emission in the line-integrated map is significant at a level ≳3σ (the integrated rms is given in Table B.1). For the sources without CO detection, we provide a 3σ upper limit, calculated from the rms of the velocity-integrated maps assuming a line width of 348 km s⁻¹ (i.e., the mean of the FWHM values measured for our detections). The results of our analysis, namely the line flux, FWHM and redshift of the line, are reported in Table B.1. We measured FWHM of the CO line in the range 97–732 km s⁻¹ and flux densities in the range 76–1942 mJy km s⁻¹. The redshift of the CO line is in good agreement with that reported in Table 1. Moment 0 maps and spectra of the targets with a detected line are reported in Appendix D.

Finally, we find that all but one target are spatially unresolved and/or the S/N is not enough to firmly constrain their CO(3-2) sizes in the current observations. cid_346 is spatially resolved with deconvolved CO(3-2) FWHM size (0.61 ± 0.15)″ × (0.31 ± 0.28)″ as derived from a two-dimensional Gaussian fit in the image plane. As a consistency check, we also analyzed the CO maps in the uv plane, and compared the results with those obtained from the analysis of the CO maps in the image plane. The analysis in the uv plane was performed with GILDAS (UVFIT). We fitted the velocity-integrated visibilities of the CO lines with either a point-source (for unresolved sources) or Gaussian functions (for resolved sources), and found that RA, Dec positions, fluxes and sizes are consistent with those obtained from the analysis in the image plane. In the maps of cid_1253 and cid_971 we found a companion at 2.19 and 5.61 arcsec, respectively. Instead, a tentative emission-line detection (S/N ≈ 3) was found in the field of cid_1215 at 3.6 arcsec. As for cid_1253 and cid_1215, we note that the 2σ contours of the targets include their companion and therefore the whole system was considered when extracting the CO(3-2) spectra. Presumably the large beam characterizing Herschel observations, used to derive the FIR luminosities, cannot separate the central target from its companion. Since cid_971 lacks FIR data and the SFR was estimated from the UV-to-NIR SED (see Sect. 2.2), it was sensible to exclude the companion from the extraction region.

The number of detections found (11 out of 27 targets) corresponds to a detection rate of ∼40%. This means that CO luminosities for non detections are lower than what we were expecting from estimates based on samples of star-forming galaxies (Genzel et al. 2010). We note that most (70%) of the non-detections are also characterized by upper limits on the FIR luminosity, which may have contributed with additional uncertainty to the predicted CO luminosity.

We finally derived CO(3-2) luminosities (in Table B.1) as given by Solomon & Vanden Bout (2005):

$$\begin{array}{c}\hfill L_{\mathrm{CO}}^{\prime }=3.25\times {10}^7{I}_{\mathrm{CO}}{\nu }_{\mathrm{obs}}^{-2}{D}_{\mathrm{L}}^2{(1+z)}^{-3},\end{array}$$(1)

where I_CO is the velocity-integrated flux, D_L is the luminosisty distance, ν_obs is the observed frequency of the line and z is the redshift. We measured CO(3-2) luminosities in the range $log(L_{\mathrm{CO}}^{\prime }/{\mathrm{K}\mathrm{km}\mathrm{s}}^{-1}{\mathrm{pc}}^2)=9.33-10.80$. For completeness, we also derived gas masses for both our AGN sample and the comparison one (in the range log M_mol/M_⊙ = 10.19−11.66) by adopting uniform assumptions, that is the excitation factor $r_{31}=L_{\text{CO}(3\text{-}2)}^{\prime }/L_{\text{CO}(1\text{-}0)}^{\prime }=0.5$ and α_CO = 3.6 M_⊙/(K km s⁻¹ pc²) commonly used in the literature for star-forming galaxies at similar redshifts as our AGN sample (Daddi et al. 2010b; Tacconi et al. 2013; Kakkad et al. 2017). However, literature results show that AGN and inactive galaxies may feature different excitation factors (e.g., Kirkpatrick et al. 2019, but see also Riechers et al. 2020). Moreover, α_CO prescriptions dependent on metallicity and star-formation conditions can be adopted (Bolatto et al. 2013). The decision to use the same assumptions for both samples follows our goal of minimizing systematic uncertainties, which could produce artificial offsets between our samples. Molecular gas masses were derived to provide a reference to a physically motivated parameter in the plots presented in Sect. 5, besides the CO(3-2) luminosity, which is the main focus of our study and is adopted as a consistent proxy for the molecular gas mass that allows us to avoid systematics involved when assuming conversion factors.

5. Results

In this section we perform comparisons between CO and FIR luminosities as well as stellar masses of AGN and inactive galaxies. In our analysis we focused on the CO(3-2) luminosity for both AGN and inactive galaxies to limit systematic uncertainties and avoid assumptions on conversion factors (e.g., excitation correction and α_CO) needed to estimate the CO(1-0) luminosity and the gas mass. Therefore, the CO(3-2) luminosity can be considered as a proxy for the molecular gas mass. Similarly, we used the FIR luminosity, from which we subtracted the AGN contribution (Sect. 2.2), as a proxy for the SFR. Stellar masses were derived from broad-band SED fitting for both samples. Bayesian methods were used in order to properly take into account the upper limits on both CO and FIR luminosities, which especially concern the AGN sample. Throughout our analysis, we relied on the lack of assumptions on conversion factors, the uniform characterization of the multiwavelength properties of our sample as well as the Bayesian methods including upper limits adopted.

5.1. CO versus FIR luminosities

In Fig. 2 we study the correlation between CO and FIR luminosities for our sample and the comparison one. This is the integrated form of the Schmidt-Kennicutt relation, although we use the CO(3-2) transition instead of CO(1-0). In order to quantify the distribution of the two samples in this plane we fitted a linear model to the data by using the ordinary least-square (OLS) bisector fit method (Isobe et al. 1990), that is, by taking into account uncertainties on L_FIR and $L_{\text{CO}}^{\prime }$ separately and then considering the bisector of the two lines. To derive the best-fit parameters of the individual fits we adopted a Bayesian framework. When constructing the likelihood function, we assumed that uncertainties are Gaussian-distributed for detections, while upper limits were taken into account by integrating the Gaussian likelihood from minus infinity to the value of the upper limit (i.e., the 3σ flux upper limits derived in Sect. 4), which resulted in the error function (e.g., see Sawicki 2012). The likelihood function was sampled using EMCEE (Foreman-Mackey et al. 2013), a python implementation of the affine invariant MCMC (Monte Carlo Markov Chain) ensemble sampler of Goodman & Weare (2010). We sampled the posterior distribution in the parameter space to derive the marginalized posterior probability distribution. The initial guesses were given by the maximum likelihood estimates obtained with the python module scipy.optimize. We also included an intrinsic scatter to the relation as third free parameter in the fit that accounts for underestimated uncertainties, given the presence of luminosities with uncertainties spanning a large range of values, which could bias the fitting results. Best-fit parameters were estimated by taking the median of the sampled marginalized posterior distribution of the OLS bisector fit parameters and the 16th and 84th percentiles as uncertainties. We then used these best-fit parameters to derive slope and intercept of the bisector as given in Isobe et al. (1990). The results of the analysis are shown in Fig. 2, where we plot the AGN sample in red and inactive galaxies in blue, with the corresponding fits and their dispersions. The dispersion shown is obtained by taking 500 realizations of the bisector fit by considering the values of slope and intercept within one sigma of the sampled marginalized posterior distributions of the two fits (along the x and y axes). The fit ($log{L}_{\mathrm{CO}}^{\prime }/{\mathrm{K}\mathrm{km}\mathrm{s}}^{-1}{\mathrm{pc}}^2=mlog{L}_{\mathrm{FIR}}/L_{\odot }+b$) to the AGN sample has slope m = 1.35 ± 0.32 and intercept $b=-6.{90}_{-3.96}^{+3.90}$, while for inactive galaxies m = 1.14 ± 0.19 and $b=-3.{91}_{-2.28}^{+2.16}$. Because of the scatter characterizing the two samples and the small luminosity range probed, the uncertainties are quite large. Within the uncertainties, our fits indicate an almost linear logarithmic slope for both samples. Previous analyses in the literature show different slopes; for example, Sargent et al. (2014) find a value of 0.81 ± 0.03 over a large range of redshifts and galaxy types, while Sharon et al. (2016) report a linear relation for z ≈ 2 AGN and SMGs that becomes superlinear (m ≈ 1.15−1.2) when low-redshift IR-bright galaxies are included. However, a detailed analysis of the slope of the relation and/or trends across a wide range of luminosities and redshift is beyond the scope of this paper (e.g., see Yao et al. 2003; Sargent et al. 2014; Kamenetzky et al. 2016; Sharon et al. 2016). Instead, we aim at quantifying any potential shift between the distributions of AGN and inactive galaxies. As can be seen in Fig. 2, the dispersions of the two fits are quite large and the two distributions do not seem to be statistically different within the uncertainties. We provide an estimate of the shift between the two samples at log(L_FIR/L_⊙) = 12.2, where the dispersion of the fits in the y direction is better constrained (∼0.15 dex). AGN have CO luminosities a factor of 2.7 lower (0.43 dex), which are different at the ∼2σ level. To verify how different the two distributions are, we also performed a Kolmogorov-Smirnov (KS) test for two-dimensional datasets⁸. Following Press & Teukolsky (1988) we implemented a Monte Carlo approach and 100 sets of synthetic data were generated from our observed CO and FIR luminosities. As done before, a Gaussian distribution was assumed for detections. As for upper limits, we included them by using the same value (i.e., the value of the upper limit) in each synthetic dataset. We then considered the fraction of synthetic datasets showing a two-dimensional KS statistic (often indicated as D) larger than the value measured for the observed data. This fraction represents how significant the difference between the two samples is. The test was first performed by excluding upper limits, and the two distributions result to be different at the ∼28% level, although the number of suitable AGN is just eight. When including upper limits, the significance increases to 99%. Although $L_{\text{CO}}^{\prime }$ upper limits could contribute in making this result a lower limit, L_FIR upper limits may go in the opposite direction. Overall, this is an indication that taking upper limits into account is crucial to understand the difference between the two samples.

Fig. 2. $L_{\text{CO}}^{\prime }$ and molecular gas masses as a function of LFIR for our AGN sample (red circles) and inactive galaxies (blue squares; Tacconi et al. 2018). To display the distribution of the two samples in the plot, we performed a bisector fit by adopting a Bayesian framework (see text). The dispersion of the fits is given by plotting 500 realizations of the bisector fit. Where the fits are better constrained (log LFIR/L⊙ = 12.2) AGN show CO luminosities 0.43 dex lower than inactive galaxies at the ∼2σ level.

We performed the linear fit by removing from the comparison sample the two outliers with low CO luminosities ($log{L}_{\mathrm{CO}}^{\prime }/{\mathrm{K}\mathrm{km}\mathrm{s}}^{-1}{\mathrm{pc}}^2<9$), which are lensed galaxies from Saintonge et al. (2013) and, given the known uncertainties on the magnification correction (e.g., Sharon et al. 2016), could be a source of bias in our analysis. Without these two targets, the slope of the fit is $0.{98}_{-0.16}^{+0.18}$ and the intercept $-1.{86}_{-2.17}^{+1.96}$, which is consistent with the previous result within the uncertainties. As for the two-dimensional KS test, we obtained a significance of ∼35%(99%) without (with) upper limits.

5.2. CO luminosities versus stellar masses

We compare $L_{\text{CO}}^{\prime }$ with the stellar mass of AGN and inactive galaxies in Fig. 3. To quantify the differences between the two samples, we performed the bisector fit as described in Sect. 5.1. The results of the fits ($log{L}_{\mathrm{CO}}^{\prime }/{\mathrm{K}\mathrm{km}\mathrm{s}}^{-1}{\mathrm{pc}}^2=mlog{M}_{*}/M_{\odot }+b$) to the AGN sample is $m=0.{98}_{-0.44}^{+0.47}$ and $b=-1.{00}_{-5.10}^{+4.73}$, while for inactive galaxies $m=1.{11}_{-0.17}^{+0.19}$ and $b=-2.{08}_{-1.99}^{+1.82}$, as displayed in the top panel of Fig. 3. We also note that the statistics of the AGN sample is slightly reduced (21 objects) with respect to the previous fitting analysis, since three targets do not have an estimate of the stellar mass (see Sect. 2.2 and Table A.1) and therefore are excluded from the fit. Our results show large uncertainties and we do not find any clear difference between the two samples. Similarly, the KS test for two-dimensional datasets implementing a Monte Carlo approach provides a significance of ∼32% without upper limits. Instead, the significance including upper limits is ∼65%. This result can be taken as a lower limit, since CO luminosities for the targets without a detection could be lower than the value assigned in the Monte Carlo simulation (i.e., the value of the upper limit). We performed the analysis by removing from the fit the two outliers in the comparison sample with low CO luminosities. The slope is $0.{96}_{-0.15}^{+0.17}$ and the intercept $-0.{48}_{-1.79}^{+1.58}$, which is consistent with the previous result given the large uncertainties. The two-dimensional KS test provides a significance of ∼30% (∼54%) without (with) upper limits.

Fig. 3. $L_{\text{CO}}^{\prime }$ and molecular gas masses as a function of M* for our AGN sample (red circles) and inactive galaxies (blue squares; Tacconi et al. 2018). Top: bisector fit performed by adopting a Bayesian framework to display the distribution of the two samples (see text). The dispersion of the fits is given by plotting 500 realizations of the bisector fit. We do not identify a clear difference between the two samples, since the fits are consistent within the uncertainties. Bottom: mean $L_{\text{CO}}^{\prime }$ per bin of stellar mass obtained by using a hierarchical Bayesian approach plotted in darker colors (see text). In the background, with lighter colors, we plot the distribution of the observed quantities, as in the top panel. Due to the large number of upper limits in our sample and the small statistics, error bars are quite large therefore it is difficult to assess in a robust way the true distribution of AGN. The difference is better constrained in the high-mass bin, where AGN show a mean CO luminosity lower than inactive galaxies by 0.72 dex at the ∼3σ level.

We performed an additional analysis in this plane, by dividing the targets in bins of stellar mass and for each we computed the mean $L_{\text{CO}}^{\prime }$. We joined the targets of the low-mass bins (log M_*/M_⊙ = 9.5−10 and 10−10.5) given the very low number of targets for both samples. The other mass bins have a width of 0.5 dex. We also note that this type of analysis does not take into account errors on the stellar mass, since each target is assigned to a bin based on its M_* measured value. To derive the mean $L_{\text{CO}}^{\prime }$, we adopted a hierarchical Bayesian approach, by including the prior assumption that CO luminosities follow a common distribution, that is the prior distribution from now on. In our framework, we assumed that this prior distribution is Gaussian, described by two parameters (called hyper-parameters), mean μ and standard deviation σ. In particular, μ is the mean $L_{\text{CO}}^{\prime }$ per stellar mass bin that we want to estimate. The key point of this approach is that μ and σ are not given as inputs, but we infer their distributions directly from the data, by fitting the CO luminosities of our targets simultaneously in each bin. Given the consistent number of CO upper limits in our sample, we used the hierarchical framework to provide better constraints on the distribution of upper limits, by exploiting the prior distribution followed by CO luminosities.

The validity of the assumption that CO luminosities follow a Gaussian distribution was tested on the xCOLD-GASS reference survey

(Saintonge et al. 2017), which provides a complete sampling of the molecular gas content in galaxies across the main sequence in the

nearby Universe, purely selected by mass. We considered galaxies with M_* > 10^9.5 M_⊙ on the main sequence and the

distribution of their CO luminosities results to be Gaussian. In order to check for redshift effects, we performed the same test on

the PHIBSS survey, from which main-sequence galaxies at 2 < z < 2.5 with M_* > 10^9.5 M_⊙ featuring CO observations were selected. Although the statistics is much worse than xCOLD-GASS, the CO luminosities of the sample roughly follow a Gaussian

distribution. We used uniform priors (only positive values) on the hyper-parameters μ and σ of the prior distribution. Similarly to the linear fit described in Sect. 5.1, the likelihood function was constructed by assuming that the uncertainties are Gaussian-distributed for detections. For non detections, we used the error function. The likelihood function was sampled with EMCEE and the mean CO luminosity per bin of stellar mass is given by the 50th percentile of the sampled marginalized posterior distribution of μ for both the AGN and inactive galaxy sample. The 16th and 84th percentiles of the distribution are taken as uncertainties. Our results are plotted in darker colors in Fig. 3 (bottom panel), while the background points display the observed quantities. Overall, it seems that AGN show mean CO luminosities 0.3–1.0 dex lower than inactive galaxies. However, due to the large number of upper limits in our sample and the small statistics, error bars are quite large and therefore it is difficult to assess in a robust way the true distribution of AGN CO luminosities. Indeed, upper limits dominate the distribution of μ in the low stellar mass bins, which are characterized by a tail at low CO luminosities as proved by the low mean value (see Fig. 3, bottom panel). However, the difference between the two samples is better constrained at stellar masses higher than 10¹¹ M_⊙, given the larger number of datapoints and especially detections. In this mass bin the mean is lower by 0.72 dex at the ∼3σ level, with values of $9.{43}_{-0.60}^{+0.21}$ for AGN and 10.15 ± 0.11 for inactive galaxies. Since the two outliers of the inactive-galaxy sample reside in the low-mass bin, their removal from the analysis does not produce significant differences, given the uncertainties.

5.3. Gas fractions

To gain further insight into the comparison of CO properties and stellar masses of our samples, we finally analyzed the ratio between $L_{\text{CO}}^{\prime }$ and stellar mass, which is a proxy for the gas fraction. A similar analysis on the ratio between $L_{\text{CO}}^{\prime }$ and L_FIR (which is a proxy for the depletion time) was not carried out because 40% of the targets suitable for this investigation are characterized by an upper limit on both quantities. We applied the same hierarchical Bayesian method described in Sect. 5.2 with the aim of quantifying the mean gas fraction of the two samples by taking upper limits into account. Again, we assumed that gas fractions in the two samples follow a common Gaussian distribution and we aim at estimating the mean of the distribution (μ). Through the analysis of the whole set of AGN or inactive galaxies without any binning, we want to increase the statistics and obtain better constraints on μ. The other advantage of the hierarchical approach is that, by fitting the sample of gas fractions simultaneously, we derive the posterior distribution of the hyper-parameters of the prior (μ and σ) as well as the posterior distribution of the gas fraction for each target, used below to obtain the total marginalized posterior distribution.

As done before, the assumption that gas fractions follow a Gaussian distribution was verified on the xCOLD-GASS survey. We selected galaxies with M_* > 10^9.5 M_⊙ on the main sequence and the distribution of their gas fractions results to be Gaussian and slightly asymmetric toward low values. This could be due to a combination of decreasing gas fractions with increasing stellar mass and the fact that above M_* = 10^10.5 M_⊙ xCOLD-GASS starts to have a more important contribution from galaxies at the bottom-end of the main sequence and below, featuring lower gas fractions. The same check was performed on main-sequence galaxies with M_* > 10^9.5 M_⊙ featuring CO observations from the PHIBSS survey. The distribution of gas fractions follows a Gaussian trend. We used uniform priors (only positive values) on the standard deviation of the prior distribution.

The sampled marginalized posterior distribution of the mean log($L_{\text{CO}}^{\prime }$/M_*) for both AGN and inactive galaxies is shown by the filled histograms in Fig. 4. As data points, we plot individual detections and upper limits in the bottom part of the Figure. The unfilled histograms represent the total sampled marginalized posterior distributions and were obtained by joining the sampled posterior distributions of each target in the AGN and inactive galaxy sample. The total posterior distributions cover the same range of values spanned by the individual datapoints, but they also incorporate the information on upper limits. The filled histograms represent the distribution of μ, the mean of the prior distribution that describes the trend followed by gas fractions. We took the 50th, 16th and 84th percentiles from the distributions of μ as best value and uncertainties. AGN show a mean log($L_{\text{CO}}^{\prime }$/M_*) of $-1.{39}_{-0.30}^{+0.23}$, while for inactive galaxies we find $-0.{82}_{-0.11}^{+0.14}$. The AGN log($L_{\text{CO}}^{\prime }$/M_*) is lower than the inactive galaxy one by a factor ∼3.7 (0.57 dex), and they are different at the 2.2σ level. By performing a two-sample KS test, the Student t-test and the Anderson-Darling test on the posterior distributions (both the mean and the total) of the two samples we can reject the null hypothesis that they are drawn from the same distribution (P > 99%). The logrank and Gehan’s tests (Schmitt 1985) performed on the observed data provide a modestly significant result (p-value ≈ 0.02–0.05).

Fig. 4. Distribution of the ratio between $L_{\text{CO}}^{\prime }$ and stellar mass (proxy of gas fraction) for AGN (red) and inactive galaxies (blue). The filled histograms show the sampled posterior distribution of the mean (μ) of the hierarchical Gaussian prior adopted in our Bayesian analysis. In the bottom part of the plot, individual detections and upper limits are displayed. The unfilled histograms represent the total distributions and were obtained by joining the sampled posterior distributions of each target in the AGN or inactive galaxy sample. AGN show a mean log($L_{\text{CO}}^{\prime }$/M*) of $-1.{39}_{-0.30}^{+0.23}$, while for inactive galaxies we find $-0.{82}_{-0.11}^{+0.14}$. The log($L_{\text{CO}}^{\prime }$/M*) ratio of AGN is lower than inactive galaxies by a factor ∼3.7 (0.57 dex), at the 2.2σ level. The two distributions are significantly different, as confirmed by statistical tests (see text).

By removing the two low-CO luminosity outliers from the inactive-galaxy sample, the mean log($L_{\text{CO}}^{\prime }$/M_*) of inactive galaxies does not change substantially. Similarly, the results of the statistical tests mentioned above remain almost unchanged, although the result of the logrank test is now significant (p-value ≈ 0.004).

6. Discussion

In this work we search for possible signatures of AGN feedback. We compared CO and FIR luminosities as well as stellar masses of X-ray selected AGN and inactive galaxies, and quantified their differences by: 1. performing a linear fit in the $L_{\text{CO}}^{\prime }$ − L_FIR plane; 2. performing a linear fit in the $L_{\text{CO}}^{\prime }$ − M_* plane; 3. dividing our samples in bins of stellar mass and computing mean CO luminosities for each bin; 4. deriving the mean distribution of $L_{\text{CO}}^{\prime }$/M_* (a proxy of gas fraction) for both samples; 5. assessing how significant the differences are by performing statistical tests, such as the KS test (for both one- and two-dimensional datasets), the Student t-test, the Anderson-Darling test, and the Gehan’s and logrank tests. Differently from previous work, we used the same CO transition for a fairly representative sample of AGN and inactive galaxies, controlled for assumptions and sources of bias and treated upper limits statistically. The two samples show statistically consistent trends in the $L_{\text{CO}(3\text{-}2)}^{\prime }-L_{\text{FIR}}$ and $L_{\text{CO}(\text{3-2})}^{\prime }-M_{\ast }$ planes. However, when we focus on the subset of parameters where the results are better constrained (i.e., L_FIR ≈ 10^12.2 L_⊙ and M_* > 10¹¹ M_⊙) and on the distribution of the mean $log(L_{\text{CO}(\text{3-2})}^{\prime }/M_{\ast })$, there are indications that AGN are underluminous in CO (0.4–0.7 dex) with respect to inactive galaxies at the 2–3σ level.

Since we presented a homogeneous characterization of the multiwavelength properties of our targets (Sect. 2.2), we exploited this set of physical parameters by investigating potential trends of $L_{\text{CO}}^{\prime }$/L_FIR and $L_{\text{CO}}^{\prime }$/M_* with respect to AGN bolometric luminosity, obscuring column density, X-ray luminosity and Eddington ratio⁹. We do not find any significant correlation with these properties (e.g., Shangguan et al. 2020a). As for the bolometric luminosity, there is a known limitation due to different timescales of the physical processes probed by L_FIR and L_bol. Star-forming activity as traced by FIR observations has a timescale of ∼100 Myr while AGN activity can vary on much shorter time intervals, < 10 Myr (e.g., Hickox et al. 2014; Stanley et al. 2015). As for the column density, evolutionary scenarios of powerful AGN predict that the unobscured phase follows the obscured one, after the AGN has removed some gas and dust from the galaxy because of feedback mechanisms such as outflows (e.g., Di Matteo et al. 2005; Hopkins et al. 2006). However, previous studies in the literature looked at possible differences between obscured and unobscured AGN in the $L_{\text{CO}}^{\prime }$ − L_FIR plane but no trend with obscuration is found (Perna et al. 2018; Shangguan & Ho 2019). Identifying the evolutionary phase of an AGN based on the amount of obscuration is challenging because of the nonuniform distribution of the absorbing material and rapid variations of the AGN duty cycle. Therefore, potential trends may be washed out in the scatter. Lastly, given the smaller subset (17/28) of targets with a black hole mass measurement available, the reduced statistics combined with upper limits could hamper a detailed analysis.

AGN are normally assumed to have higher CO excitation (e.g., Weiß et al. 2007; Carilli & Walter 2013), hence higher excitation correction r₃₁ (common values are 0.92 and 0.5 for AGN and inactive galaxies, respectively; Daddi et al. 2015; Tacconi et al. 2018; Kirkpatrick et al. 2019, but see also Riechers et al. 2020). We find on average lower CO(3-2) luminosities in AGN than in inactive galaxies (at given stellar masses and FIR luminosities), therefore even by conservatively assuming the same r₃₁ we would find lower molecular gas masses in AGN, and assuming higher r₃₁ would exacerbate this difference. Lastly, different assumptions on α_CO, based for example on metallicity and star-formation properties (Bolatto et al. 2013), can further increase, but also decrease, the difference. Measurements of the CO(1-0) line, for example with the JVLA, would be necessary to investigate excitation properties and obtain a reliable estimate of the total molecular gas reservoir in our AGN.

A deficit of CO(3-2) emission in AGN has also been reported by Kirkpatrick et al. (2019), who performed a study of CO emission properties of galaxies as a function of AGN contribution in the MIR, for a heterogeneous sample collected from the literature. Although the results are not statistically robust because of substantial uncertainties and small sample sizes, they find systematically lower CO fluxes in AGN. The ratio between CO(3-2) fluxes of AGN and star formation-dominated galaxies is 0.58 ± 0.2, for which they have a larger statistics with ∼30 detections in total across the two samples. Such a result is in agreement with what we find for our samples and would translate in lower gas masses in high-redshift AGN.

6.1. Comparison with previous work

Previous work on 1 < z < 3 AGN often found systematically reduced molecular gas reservoirs compared to inactive galaxies. For example, Carilli & Walter (2013) and Perna et al. (2018) present lower CO luminosities for an heterogeneous collection of AGN from the literature (but see also Kirkpatrick et al. 2019); Kakkad et al. (2017) analyzed a sample of 10 AGN with homogeneous ALMA CO(2-1) observations; Carniani et al. (2017), Brusa et al. (2018) and Loiacono et al. (2019) focus on a few targets with detected ionized outflows. Sharon et al. (2016) looked for differences between AGN and SMGs and find the two samples to be consistent in their CO and FIR properties for both CO(3-2) and CO(1-0) transitions. Conversely, Spingola et al. (2020) observed two lensed AGN at z ≈ 2−3 and find larger CO(1-0) luminosities with respect to their FIR emission when compared with star-forming galaxies, meaning that they may be less efficient at forming stars.

Potential sources of bias in previous work and the way we addressed them are as follows.

CO luminosities. When samples are assembled from the literature by joining observations of different transitions, excitation factors have to be assumed in order to estimate the CO ground transition for all the targets. These correction factors depend on the SLED and, given the steeper SLED often found for AGN (Carilli & Walter 2013; Kirkpatrick et al. 2019), the differences between AGN and inactive galaxies are systematically amplified when deriving the CO(1-0) luminosity. However, the CO flux at the various transitions can vary in individual targets based on the physical conditions of the ISM (Carilli & Walter 2013). Different assumptions on α_CO, for example a starburst-like value adopted for AGN in some cases (e.g., Perna et al. 2018), further contribute to the differences. In this work this issue was addressed by comparing directly the observed CO(3-2) luminosity for both samples.

FIR luminosities and stellar masses. A complication in deriving reliable FIR luminosities of AGN host galaxies is that the AGN can provide a significant contribution to the FIR luminosity. This issue is magnified when different methods to estimate FIR luminosities are used in the same analysis (e.g., in literature samples). For example, Kirkpatrick et al. (2019) and Perna et al. (2018) show how different methods to derive this quantity bring to different conclusions. A homogeneous SED-fitting decomposition for our AGN sample was performed in our analysis in order to derive FIR luminosities associated with star-formation activity. As for the comparison sample of inactive galaxies, we note that SFRs are taken from the literature and derived with different tracers in a few cases, although we rely on the observational evidence that main-sequence galaxies show agreement among Hα and IR star-formation tracers (Shivaei et al. 2016, see also Sect. 3). Uniformly measuring FIR luminosities for the comparison sample would allow us to quantify the entity of this bias, if any. Moreover, stellar masses are not available for many bright AGN (often the main targets of previous observations), which is a limitation toward assessing where they lie on the main sequence and performing a meaningful comparison with matched samples of inactive galaxies. Estimates of the stellar mass are available for the majority (∼78%) of our targets, which include more common AGN.

Sample size. Previous conclusions on the CO properties of AGN were sometimes biased toward bright objects and therefore the more extreme ones of the AGN population (e.g., Carniani et al. 2017; Brusa et al. 2018). In analyses of larger samples, albeit more representative, upper limits were excluded in some cases (e.g., Kirkpatrick et al. 2019). The study carried out in this work, focused on X-ray selected AGN, presents more common sources spanning a wide range of properties such as AGN bolometric luminosities, stellar masses and SFRs, as well as a Bayesian analysis to take upper limits into account.

We performed a careful selection of the comparison sample, with the aim of building a set of targets with similar properties to our AGN on the main sequence and looking for differences in CO properties. By taking into account the biases mentioned above, we find that the two samples do not show statistically significant differences although there are hints of AGN being deficient in their CO(3-2) luminosities at the 2–3σ level. Therefore, caution must be exercised when mixing targets with observations tracing different CO transitions, since a priori assumptions on excitation corrections can artificially produce larger offsets between AGN and inactive galaxies. Finally, we note that our sample selection is limited by 50% SFR upper limits of AGN on the main sequence (Fig. 1). Indeed, if the SFRs of these targets, and in turn their FIR luminosities, are actually much lower than those of the comparison sample we cannot rule out that the lower CO luminosities of AGN are consistent with their potentially lower FIR luminosities, according to the L_CO − L_FIR relation (Sargent et al. 2014). Future improvements of our analysis rely on ALMA as the key facility to provide stronger constraints on the results presented in Sect. 5 and reduce the dispersion of our fits. More sensitive CO observations are needed to detect the CO(3-2) emission of our undetected AGN as well as dust-continuum observations to allow a precise characterization of the FIR emission for the targets with just upper limits from Herschel. These observations, combined with larger and more uniformly selected samples, will allow us to assess whether there is a true difference between AGN and inactive galaxies.

6.2. The effects of AGN activity

Our results point toward a difference in CO luminosity between AGN and inactive galaxies, which could be ascribed to the presence of the central engine that may have a role in heating, exciting, dissociating, and/or ejecting the gas. However, understanding the physical mechanism producing such a difference requires tracing different gas phases and several CO transitions. With only the CO(3-2) data in hand, we provide possible scenarios to interpret our findings.

AGN radiation field. Interesting examples of the impact of AGN activity on spatially resolved CO properties are presented by Feruglio et al. (2020) and Rosario et al. (2019). Thanks to high-resolution observations of different gas phases (both molecular and ionized) in two nearby Seyfert galaxies, they find CO(2-1) cavities in the inner 200 pc around the AGN, which are filled with ionized and warm molecular gas (traced by H₂ NIR emission lines). They conclude that CO may be excited higher up in the rotational ladder in the CO cavities, giving rise to the apparent lack of CO emission at low J. It is also possible that CO is dissociated by the AGN radiation field and/or shocks. In particular, Kawamuro et al. (2020) propose that the AGN X-ray field can be responsible for weakening the CO(2-1) emission. In those Seyfert galaxies the action of the central engine is limited to the inner regions within a few hundred parsec from the AGN. It is unlikely that this cavity-effect would be detected in integrated measurements of CO transitions as it would be diluted in the beam (e.g., see Lamperti et al. 2020 for a study of integrated CO(3-2)/CO(1-0) properties of z < 0.05 AGN and star-forming galaxies). At high redshift, however, when the average AGN population is more active and more luminous, the AGN could have an impact out to larger spatial scales in the host galaxy, possibly out to kpc scales, and in this case it would be possible to detect an overall deficit of CO emission. Some indications of this scenario at higher redshift are provided by lensed AGN. Spingola et al. (2020) and Paraficz et al. (2018) observed a decrease in the surface brightness of the CO(1-0) and (2-1) emission (at z ≈ 2 and 0.6), respectively, at the location of the AGN emission and the center of host galaxies. This is interpreted as a potential decrement in the molecular gas at low excitation close to the AGN possibly due to AGN radiative feedback. As argued by Fogasy et al. (2020), the lack of emission at a given CO transition does not necessarily mean that the source is gas poor. Indeed, they analyzed integrated ALMA observations of a z ≈ 2.8 AGN and find that the target is undetected at low-J transition lines but appears bright at high-J transitions, indicating the presence of warm and highly excited molecular gas. For our AGN, observations of more CO transitions to probe the SLED and CO observations at higher resolution would be necessary to verify the impact of AGN activity on CO excitation and smaller physical scales.

AGN-driven outflows. Another possible effect of the AGN activity on the molecular gas is through outflows. This possibility is supported by observations of individual objects: For example, Carniani et al. (2017), Brusa et al. (2018) and Loiacono et al. (2019) find low gas fractions in powerful AGN at cosmic noon hosting high-velocity molecular and ionized outflows (but see also Herrera-Camus et al. 2019). AGN feedback in action in these targets could be depleting the molecular gas reservoir (Brusa et al. 2015). Förster Schreiber et al. (2019), studying outflows in a large sample of 0.6 < z < 2.7 galaxies through integral field spectroscopy of the Hα emission line, find that incidence, strength, and velocity of AGN-driven winds are strongly correlated with the stellar mass. In particular, they find that high-velocity (∼1000–2000 km s⁻¹) AGN-driven outflows are commonly detected at masses above log(M_*/M_⊙) = 10.7, and present in up to 75% of the population for log(M_*/M_⊙) > 11.2. Interestingly, above this stellar-mass threshold we find a significant CO luminosity deficit in our AGN sample with respect to inactive galaxies (Fig. 3, bottom). Moreover, our AGN show on average gas fractions 0.57 dex (by using uniform assumptions, Sect. 4) lower than inactive galaxies at the 2.2σ level. Quantitatively, this translates into M_gas, mol/M_* ≈ 0.3 for AGN (0.16 if we use r₃₁ = 0.92; Kirkpatrick et al. 2019) and ≈1 in inactive galaxies. This representative value for our AGN is in line but not as low as previous work targeting extremely powerful sources (e.g., M_gas, mol/M_* < 0.05 in Brusa et al. 2018). Our team is performing a systematic investigation of ionized gas outflows with SINFONI as part of the SUPER survey, and 11 targets of our ALMA sample have complementary good quality SINFONI data (Kakkad et al. 2020; Perna et al., in prep.). For some of them we measured [O III] line widths larger than 600 km s⁻¹, interpreted as a clear signature of the presence of an AGN-driven outflow in these objects (Kakkad et al. 2020). A detailed comparison between outflow and CO properties for these targets will be presented in a future work. Distinguishing among the scenarios described above is challenging with the current dataset. AGN feedback could proceed in different ways and different mechanisms likely overlap in shaping the properties of the molecular gas reservoir. For example, AGN radiation could both heat and/or dissociate CO molecules. In this case, AGN would produce a feedback mechanism that does not require outflows but would potentially work toward inhibiting further star formation. As for AGN-driven outflows, they could impact the gas content by ejecting material out of the galaxy (e.g., Travascio et al. 2020), or they could produce CO heating or dissociation due to shocks. Additionally, numerical simulations predict that AGN-driven outflows may heat via shocks a significant quantity of the gas in the ISM, reaching the high temperatures required for the excitation of high-J CO transitions (Costa et al. 2018). To reach a deeper understanding of the impact of AGN on the molecular gas reservoir, also on longer timescales, predictions from simulations providing the spatial scales and effects of AGN activity on CO properties as a function of cosmic time are needed.

7. Conclusions

In this work we presented the first systematic investigation of the CO(3-2) emission of AGN at z ≈ 2 by using ALMA observations for a sample of 27 X-ray selected AGN (Sect. 2.1). We characterized their AGN and host galaxies properties through SED-fitting and X-ray spectral analysis (Sect. 2.2) and measured stellar masses (log M_*/M_⊙ = 9.6−11.2), FIR luminosities (log L_FIR / erg s⁻¹ < 45.0 − 46.4), AGN bolometric luminosities (log L_bol / erg s⁻¹ = 44.7 − 46.9), hydrogen column densities (N_H < 2 × 10²⁴ cm⁻²) and X-ray luminosities (log L_{[2−10 keV]} / erg s⁻¹ = 43.0 − 45.4). CO emission was detected in 11 out of 27 targets and we find CO luminosities in the range $log(L_{\mathrm{CO}}^{\prime }/{\mathrm{K}\mathrm{km}\mathrm{s}}^{-1}{\mathrm{pc}}^2)=9.33-10.80$, line FWHM 97–732 km s⁻¹ and molecular gas masses log(M_mol/M_⊙) = 10.19−11.66 (Sect. 4). To infer whether AGN activity affects the CO emission of the host galaxy, we compared the CO properties of our sample with those of inactive galaxies with similar redshift, stellar masses and SFRs selected from the PHIBSS survey (Sect. 3, Fig. 1). We further controlled for systematic differences possibly introduced by conversion factors (e.g., excitation corrections, α_CO) by comparing directly the CO(3-2) luminosities for both samples. In order to properly account for upper limits we adopted a Bayesian approach throughout our analysis. Our findings can be summarized as follows.

• We compared CO and FIR luminosities of AGN and inactive galaxies (Sect. 5.1) and quantified the distribution of the two samples in the $L_{\text{CO}}^{\prime }$ − L_FIR plane by fitting a linear model through the ordinary least-square bisector fit method (Fig. 2). The resulting fits are characterized by similar (almost linear) slopes and large dispersions, and do not show a significant shift within the uncertainties. However, where the results are best constrained (i.e., log(L_FIR/L_⊙) ≈ 12.2), we find that AGN have CO luminosities 0.43 dex lower than inactive galaxies, different at the 2σ level. By applying a KS test for two-dimensional datasets implementing a Monte Carlo approach, the two samples result to be different at the 99% (∼28%) level including (excluding) upper limits.

• A similar investigation was performed in the $L_{\text{CO}}^{\prime }$ − M_* plane (Sect. 5.2), where the results show large uncertainties and different slopes, and we do not find any clear difference between the two samples (Fig. 3, top). We further divided our samples in bins of stellar mass and computed mean CO luminosities for each bin. In order to provide constraints on the distribution of upper limits, we implemented a hierarchical Bayesian framework by adopting the prior assumption that CO luminosities follow a Gaussian distribution. Despite the error bars, the difference results to be best constrained in the high-mass bin log(M_*/M_⊙) > 11 (Fig. 3, bottom), where AGN are ∼0.72 dex less luminous than inactive galaxies at the 3σ level.

• We finally performed our hierarchical Bayesian analysis on the total distribution of log($L_{\text{CO}}^{\prime }$/M_*), a proxy of gas fraction (Sect. 5.3). AGN show a mean log($L_{\text{CO}}^{\prime }$/M_*) 0.57 dex lower than inactive galaxies at the 2.2σ level (Fig. 4). By applying statistical tests such as the two-sample KS, the Anderson-Darling and the Student t-test to the sampled marginalized posterior distributions, and the logrank and Gehan’s tests to the observed data we assessed that the two samples are different at the 99% level.

Summarizing, in our analysis we controlled for assumptions and sources of bias, focused on a sample of targets covering a wide range of AGN bolometric luminosity and treated upper limits statistically. When averaged over a fairly representative sample of targets including also less extreme AGN, the differences between AGN and inactive galaxies are not very large, but when data are best constrained AGN result to be underluminous in CO as a function of FIR luminosities and stellar masses at the 2–3σ level. Overall, our AGN feature CO(3-2) luminosities lower than inactive galaxies (at given stellar masses and FIR luminosities), therefore even by conservatively assuming the same r₃₁ we would find lower molecular gas masses in AGN, and assuming higher r₃₁ would exacerbate this difference. We interpreted our result as a hint toward the potential effect of AGN activity, which may be able to heat, excite, dissociate, and/or deplete the gas reservoir of the host galaxies. To help establish whether the driving mechanisms are AGN radiation and/or outflows, observations tracing the outflowing multiphase gas are needed, as well as observations of more CO transitions to probe the SLEDs and CO observations at higher resolution to verify if the impact of AGN activity on molecular gas takes place on smaller scales.

Acknowledgments

We thank the anonymous referee for carefully reading the manuscript and providing helpful comments. C. Circosta thanks Linda Tacconi for useful advice on the PHIBSS survey and providing the table with measured CO fluxes; Vinod Arumugam and Emanuele Daddi for helpful discussions about ALMA data analysis; Sotiria Fotopoulou and Antonis Georgakakis for providing the photometry of the XMM-XXL targets. C. Circosta and AS acknowledge support from the Royal Society. A. P. gratefully acknowledges financial support from STFC through grants ST/T000244/1 and ST/P000541/1. M. P. is supported by the Programa Atracción de Talento de la Comunidad de Madrid via Grant 2018-T2/TIC-11715. We acknowledge support from PRIN MIUR project “Black Hole winds and the Baryon Life Cycle of Galaxies: the stone-guest at the galaxy evolution supper”, contract #2017PH3WAT. GV acknowledges financial support from Premiale 2015 MITic (PI B. Garilli). This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.00798.S and ADS/JAO.ALMA#2017.1.00893.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This research has made use of the following data: data based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under ESO programme ID 179.A-2005 and on data products produced by TERAPIX and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. Data from HerMES project. HerMES is a Herschel Key Programme utilising Guaranteed Time from the SPIRE instrument team, ESAC scientists and a mission scientist. The HerMES data was accessed through the Herschel Database in Marseille (HeDaM – http://hedam.lam.fr) operated by CeSAM and hosted by the Laboratoire d’Astrophysique de Marseille. HerMES DR3 was made possible through support of the Herschel Extragalactic Legacy Project, HELP (http://herschel.sussex.ac.uk), HELP is a European Commission Research Executive Agency funded project under the SP1-Cooperation, Collaborative project, Small or medium-scale focused research project, FP7-SPACE-2013-1 scheme.

References

Appendix A: Multiwavelength properties of the sample

Appendix B: Results of the ALMA CO data analysis

Appendix C: Spectral energy distributions of the sample

Fig. C.1. Rest-frame SEDs of the targets that were not presented in Circosta et al. (2018). The black dots represent the observed multiwavelength photometry, while the empty dots indicate 3σ upper limits. The black solid line is the total best-fit model, the orange curve represents the stellar emission attenuated by dust, the green template reproduces the AGN emission, the red curve accounts for dust emission heated by star formation. Emission lines in the black curves are part of the nebular emission component, included in the overall SED.

Appendix D: CO(3-2) emission

Fig. D.1. CO line emission maps (left) and spectra (right) extracted from the region above 2σ significance. Left: black contours in steps of 1σ, starting from 2σ. Dashed contours correspond to −1σ. The beam of each observation is shown as a gray ellipse on the bottom-left of the maps. Right: observed spectrum plotted in black, Gaussian best-fit model depicted in red.

Fig. D.1. continued.

Fig. D.1. continued.

Appendix E: Continuum emission

Fig. E.1. Continuum emission maps at ∼100 GHz. Black contours are in steps of 1σ, starting from 2σ, for all targets but cid_451 (σ × [2, 4, 8, 16, 24]) and X_N_102_35 (σ × [2, 4, 8, 16, 32, 64]). Dashed contours correspond to −1σ. The beam of each observation is shown as a gray ellipse on the bottom-left of the maps.

All Tables

All Figures

Fig. 1. Distribution of host galaxy properties in the SFR-M* plane for the 21 AGN (type 1s marked by triangles and type 2s marked by circles) in our sample with an estimate of both stellar mass and SFR. The comparison sample of inactive galaxies is depicted by blue squares. The two data points with orange edges represent the targets with SFR derived through modeling of the stellar emission using SED fitting. Black arrows represent 3σ upper limits. The color coding indicates the AGN bolometric luminosity for our sample, which covers two orders of magnitude in Lbol. The black solid line reproduces the main sequence of star-forming galaxies from Schreiber et al. (2015) at the average redshift of our target sample (i.e., ∼2.3). The dashed lines mark the scatter of the main sequence (equal to 0.3 dex). The histograms show the projected distribution of the two quantities along each axis, in red for our sample and blue for the comparison one.

In the text

Fig. 2. $L_{\text{CO}}^{\prime }$ and molecular gas masses as a function of LFIR for our AGN sample (red circles) and inactive galaxies (blue squares; Tacconi et al. 2018). To display the distribution of the two samples in the plot, we performed a bisector fit by adopting a Bayesian framework (see text). The dispersion of the fits is given by plotting 500 realizations of the bisector fit. Where the fits are better constrained (log LFIR/L⊙ = 12.2) AGN show CO luminosities 0.43 dex lower than inactive galaxies at the ∼2σ level.

In the text

Fig. 3. $L_{\text{CO}}^{\prime }$ and molecular gas masses as a function of M* for our AGN sample (red circles) and inactive galaxies (blue squares; Tacconi et al. 2018). Top: bisector fit performed by adopting a Bayesian framework to display the distribution of the two samples (see text). The dispersion of the fits is given by plotting 500 realizations of the bisector fit. We do not identify a clear difference between the two samples, since the fits are consistent within the uncertainties. Bottom: mean $L_{\text{CO}}^{\prime }$ per bin of stellar mass obtained by using a hierarchical Bayesian approach plotted in darker colors (see text). In the background, with lighter colors, we plot the distribution of the observed quantities, as in the top panel. Due to the large number of upper limits in our sample and the small statistics, error bars are quite large therefore it is difficult to assess in a robust way the true distribution of AGN. The difference is better constrained in the high-mass bin, where AGN show a mean CO luminosity lower than inactive galaxies by 0.72 dex at the ∼3σ level.

In the text

Fig. 4. Distribution of the ratio between $L_{\text{CO}}^{\prime }$ and stellar mass (proxy of gas fraction) for AGN (red) and inactive galaxies (blue). The filled histograms show the sampled posterior distribution of the mean (μ) of the hierarchical Gaussian prior adopted in our Bayesian analysis. In the bottom part of the plot, individual detections and upper limits are displayed. The unfilled histograms represent the total distributions and were obtained by joining the sampled posterior distributions of each target in the AGN or inactive galaxy sample. AGN show a mean log($L_{\text{CO}}^{\prime }$/M*) of $-1.{39}_{-0.30}^{+0.23}$, while for inactive galaxies we find $-0.{82}_{-0.11}^{+0.14}$. The log($L_{\text{CO}}^{\prime }$/M*) ratio of AGN is lower than inactive galaxies by a factor ∼3.7 (0.57 dex), at the 2.2σ level. The two distributions are significantly different, as confirmed by statistical tests (see text).

In the text

Fig. C.1. Rest-frame SEDs of the targets that were not presented in Circosta et al. (2018). The black dots represent the observed multiwavelength photometry, while the empty dots indicate 3σ upper limits. The black solid line is the total best-fit model, the orange curve represents the stellar emission attenuated by dust, the green template reproduces the AGN emission, the red curve accounts for dust emission heated by star formation. Emission lines in the black curves are part of the nebular emission component, included in the overall SED.

In the text

	Fig. D.1. CO line emission maps (left) and spectra (right) extracted from the region above 2σ significance. Left: black contours in steps of 1σ, starting from 2σ. Dashed contours correspond to −1σ. The beam of each observation is shown as a gray ellipse on the bottom-left of the maps. Right: observed spectrum plotted in black, Gaussian best-fit model depicted in red.
In the text

	Fig. D.1. continued.
In the text

	Fig. D.1. continued.
In the text

	Fig. E.1. Continuum emission maps at ∼100 GHz. Black contours are in steps of 1σ, starting from 2σ, for all targets but cid_451 (σ × [2, 4, 8, 16, 24]) and X_N_102_35 (σ × [2, 4, 8, 16, 32, 64]). Dashed contours correspond to −1σ. The beam of each observation is shown as a gray ellipse on the bottom-left of the maps.
In the text