EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 668 (December 2022) A&A, 668 (2022) A67 Full HTML
Open Access
Issue
A&A
Volume 668, December 2022
Article Number A67
Number of page(s) 21
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202243666
Published online 05 December 2022

© A. Audibert et al. 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1. Introduction

Feedback from active galactic nuclei (AGN) is widely believed to be a crucial element of the mass assembly of galaxies: with feedback acting on gas and stars forming from gas, feedback regulates the stellar mass growth of galaxies (Bower et al. 2006; Croton et al. 2006; Di Matteo et al. 2008; Sijacki et al. 2007; Dubois et al. 2016; Cielo et al. 2018). By leading to the acceleration or heating of the gas in the interstellar medium (ISM), feedback is associated with winds and outflows that can even eject gas outside of galaxies.

It can also delay the deposition of new gas from the intracluster medium (ICM) or intergalactic medium (IGM) onto galaxies. Both effects help the quenching of star formation through the lack of fuel for accretion, simultaneously terminating the black hole activity (Alatalo et al. 2015; Lanz et al. 2016).

It is widely believed that AGN feedback operates in a bimodal manner. The quasar mode, also called radiative or wind mode, occurs through wide-angle winds driven by radiative pressure in luminous AGN with accretion rates close to Eddington. This mode can be met at any redshift z, and it is very common in high-z, young QSOs (King 2003; Harrison 2017; Bieri et al. 2017). The radio (or kinetic) mode involves the launch of collimated relativistic jet of particles (i.e. radio jet), and occurs mainly in low-luminosity AGN (< 0.01LEdd). Radio jets launched by supermassive black holes (SMBH) are seen to have an impact even on large scales, e.g., of tens of kpc, delaying cooling flows from the ICM/IGM and the formation of extremely massive galaxies (Fabian 2012; McNamara & Nulsen 2012). This is characteristically seen at the centers of fairly low z clusters and groups. Of course, they can also have a strong impact on (sub-)kpc scales, altering the properties of the ISM and driving molecular outflows (Morganti et al. 2015; Dasyra et al. 2016; Fotopoulou et al. 2019; Oosterloo et al. 2019; Aalto et al. 2020; Fernández-Ontiveros et al. 2020), and even in low-luminosity AGN (Combes et al. 2013; García-Burillo et al. 2014; Audibert et al. 2019; Ruffa et al. 2022). The observational findings are supported by 3D hydrodynamical simulations of relativistic jets coupling in inhomogeneous and clumpy ISM (Wagner & Bicknell 2011; Wagner et al. 2012; Mukherjee et al. 2018), showing that propagating radio jets can gradually disperse gas clouds and create a cocoon of shocked material associated with multi-phase outflows. As this jet-driven bubble expands, it also leads to the dispersal of molecular clouds. The final feedback impact depends on the jet kinetic energy and the distribution of the clumpy medium defining the jet’s path of minimum resistance.

Yet, a systematic and statistically meaningful observational measurement of the number of radio galaxies (RGs) with significant gas reservoirs, upon which the radio-mode feedback could be operating, is missing. So far, it has been shown that most of the local RGs have low gas content compared to spirals or infrared selected galaxies (e.g., Ruffa et al. 2022). High-z CO surveys of RGs started in the 1990s (Evans et al. 1996) and eventually led to the discovery of large reservoirs of molecular gas (with MH2 ∼ 1010 − 1011 M⊙) in some powerful RGs at z ≳ 2 or gas-depleted, dead hosts in other (Scoville et al. 1997; Papadopoulos et al. 2000; De Breuck et al. 2003, 2005; Greve et al. 2004; Nesvadba et al. 2008; Emonts et al. 2011, 2014; Castignani et al. 2019). The increased sensibility achieved with the advent of interferometric facilities such as the Atacama Large Millimeter Array (ALMA) and the Northern Extended Millimeter Array (NOEMA) nowadays offer a great opportunity to study the evolution of the molecular gas in low and high-z RGs.

In this paper, we present a new analysis of observations of CO molecular transitions from rotational numbers J = 1 − 0 to J = 4 − 3 performed with ALMA for RGs up to z < 2.5. To build a statistically significant sample, we analyzed the CO emission in a sub-sample of RGs from the ALMA Calibrator Source Catalog, which were selected to be representative of the NRAO/VLA Sky Survey (NVSS) catalogue (Condon et al. 1998), with respect to its z and 1.4 GHz flux distribution, down to a limit of 0.4 Jy. Complementing these with literature data, we carried out a large and complete archival CO survey of molecular gas in RGs. A main aim of our work was to analyze the evolution of the gas content of RGs with z, which is highly needed for an accurate benchmarking of the cosmological simulations. Such analysis also lead us to construct for the first time the CO luminosity function of RGs at low and high z. Another goal was to determine the occurrence of molecular outflows in the RGs in our sample.

The paper is structured as follows: in Sect. 2 we describe the statistics used to reproduce a representative sample of NVSS radio galaxies. In Sect. 3 we describe the data reduction of the ALMA observations. Main results from the analysis of the literature and ALMA observations are presented in Sect. 4. The presentation and discussion of the mass and luminosity functions are described Sect. 4.3, followed by conclusions. In this paper we adopt a flat ΛCDM cosmological model with H0 = 70 km s−1Mpc−1, ΩM = 0.3, and ΩΛ = 0.7.

2. Radio galaxies with ALMA observations

2.1. ALMA data mining

The selection of galaxies for our survey started from the ALMA Radio Source Catalogue (ARC), which comprises a list of mm and sub-mm bright objects that can serve as (flux, phase, bandpass, etc.) calibrators. This catalogue includes a compilation of catalogues from several facilities other than ALMA, including the Very Large Array Calibrator Manual, Submillimeter Array and Atacama Compact Array catalogues, the Parkes survey, and the Combined Radio All-Sky Targeted Eight-GHz Survey. As such, the ensemble of the ARC covers the full sky, but in a non-homogeneous manner.

The ARC catalogue was downloaded from the European Southern Observatory (ESO) archival interface1. Each catalogue entry corresponds to observations of a single object in a single band. More specifically, each entry shows the latest flux measurement of each object in a given band, with the exception of band 3. For band 3, two entries are provided, presenting the latest flux for each sideband. As of 01/08/2020, the downloaded catalogue contained 8679 entries. The number of unique objects (deduced by removing the information on the bands) was 3362. For some of the objects, the spectral windows of the calibrations covered the frequencies of molecular line transitions. Therefore, additional science on, e.g., the detection of CO lines, was feasible with the calibration data. Some of the sources in this catalogue were also observed as science targets. The combination of calibration and science related data led to a wealth of information that enabled our archival survey to be designed.

For each ARC source we used the “astroquery” tool (Ginsburg et al. 2019) to find the available ALMA observation searching within a radius of 3″ from the registered sky coordinates. Additionally, we iterated over all names given for each ARC catalogue entry in order to cross-check detection rates. In total, we found recorded 25 827 observations from 1562 unique sources.

To determine whether a CO transition was in the spectral window of an observation, we needed the redshift of the observed source. For this purpose, we first queried the NASA Extragalactic Database (NED), Simbad, and Vizier using “astroquery”. We performed an extra quality assessment by only keeping spectroscopically determined z values, and in fact, by manually examining whether each spectroscopic z was measured using two or more spectral lines. Once having the z information, we only kept observations of sources with CO(1–0), CO(2–1), CO(3–2), or CO(4–3) coverage in a spectral window. A total of 675 sources had at least one CO line with any integration time. We kept only sources with a minimum integration time of 5 min to ensure that the data analysis (e.g., the derivation of error bars) is meaningful.

As a next step, we had to ensure that only galaxies with radio emission associated with accretion onto black holes were kept. As the goal of our survey is to investigate the role of radio-mode feedback on galaxy evolution, all objects with radio emission associated with star formation activity had to be eradicated. For this reason, we imposed the 1.4GHz-to-24 μm flux criterion of Bonzini et al. (2013), identifying radio emission in excess of what supernovae can produce in star-forming galaxies via the quantity q24 = log(S24 μm/S1.4 GHz), which has to be < 0.52.

The 1.4 GHz fluxes of our sources were obtained directly from radio surveys, when available. Otherwise, we interpolated their values from the nearest available radio data, fitting the radio data with a power law in νFν. In total, this requirement left us with 204 sources. From here onward, we will refer to all radio galaxies with ALMA CO observations as the CO-ARC catalogue.

To serve as calibrators, the CO-ARC sources are point sources at millimeter and radio wavelengths (such as quasars or blazars). To also include RGs with spatially resolved, extended radio jets, we performed an extensive bibliographic search, which provided us with a pool of either single-object or dedicated samples with CO observations (e.g. Leon et al. 2000; Evans et al. 2005; Ocaña Flaquer et al. 2010; Ruffa et al. 2019a; Russell et al. 2019; Dabhade et al. 2020). The q24 criterion described above was then applied to the selected literature sources, leaving us with a total of 152 galaxies which were then considered suitable to complement our CO-ARC sample. Hereafter, we will refer to the joint CO-ARC and bibliographic super-sample of 356 galaxies as the extended CO-ARC.

2.2. Selection of a radio-representative sample

Our next task was to find a sub-sample of the extended CO-ARC sample that is representative of a radio galaxy survey. As parent radio catalogue, we used the NVSS 1.4 GHz imaging survey, performed by the Very Large Array (VLA; Condon et al. 1998). The NVSS covers a large fraction of the sky for declinations above δ > −40°, and it is fairly deep in sensitivity, recovering 1.4 GHz flux densities as low as ∼2 mJy. This limit ensures completeness at the levels we are interested in examining, of a few hundred mJy, given the 1.4 GHz fluxes of the CO-ARC sources. In fact, we also examined a few other catalogues as potential parent surveys of ours: the Fourth Cambridge Survey (4C) at 178 MHz (Pilkington & Scott 1965), the Sixth Cambridge Survey (6C) at 151 MHz (Hales et al. 1993), the Parkes Radio Sources Catalogue (PKS) at 2.7 GHz (Wright & Otrupcek 1990) and at 4.8 GHz (Griffith et al. 1994), the MIT-Green Bank 5 GHz Survey (Bennett et al. 1986), and the Australia Telescope 20 GHz Survey Catalog (AT20G; Murphy et al. 2010). We found that NVSS is the optimal parent survey because of its similarity in the flux distribution and of its well characterized completeness3.

We then kept the largest possible sub-sample of the extended CO-ARC that is fully representative of the NVSS 1.4 GHz flux distribution. We used bootstrapping to iteratively and randomly remove more and more sources from our initial sample until we obtain a final sample consistent with the NVSS. Often, the convergence criterion of any two samples is checked with a Kolmogorov-Smirnov (KS) test: KS provides a significance level (p-value) based on the maximum distance between the two cumulative distributions (D-statistic), disproving the hypothesis that the two distributions are inconsistent. However, the p-value is rather arbitrary. Instead, we used our own precisely defined D-statistic criterion, so that bootstrapping can reliably provide the maximum possible CO-ARC sub-sample that is compatible with the NVSS. For this definition, we simulated new, mock sets of NVSS data, starting from the initial fluxes in the catalogue and adding to them Gaussian random noise with dispersion σf. Then, we measured how the D-statistic between the original NVSS fluxes and any N-sized subsample of simulated NVSS fluxes changes as a function of N. The curves are shown in Fig. 1, and they are calculated using the average D-statistic of 1000 implementations for each N. They indicate how close to the parent catalogue a child catalogue would be if its fluxes varied by 1σf and 3σf, respectively, using a typical flux uncertainty of σf = 10%. As the sample size is getting smaller, the corresponding D-statistic grows due to dominance of small number statistics: at small enough sample sizes any two samples can resemble each other and thus the statistics deteriorate.

thumbnail Fig. 1.

D-statistic, i.e., maximum distance between two cumulative distributions used for the sample selection versus the sample size. Solid lines show a theoretical comparison of the NVSS and a shuffled NVSS samples with 10 and 30% flux variations. Small black dots show how the sample D-statistic changes as a function of the number of bootstrapped sources (see text). The large circles show the final numbers for the CO-ARC (black) and for individual z ranges (color). The same analysis is shown for the CO-ARC only using the ALMA calibrators (without including the literature sources) is indicated as squares and equivalent color code for the redshift bin. We reach an optimal sample size of 120 (CO-ARC in black circle) for a D-statistics smaller than 10% and converging with the values computed to the theoretical shuffled NVSS sample down to 10% flux variation.

Having defined our theoretical convergence criterion, we proceeded with replacing the theorized child sample with the actual extended CO-ARC sample. Starting from the initial number of sources in the extended CO-ARC, we randomly drew consecutively smaller and smaller sub-samples and compared their D-statistic with the value established for the theorized NVSS sub-sample. Our bootstrapping again indicates how the D-statistic changes as a function of N: it decreases with decreasing N as (mainly bright) sources are excluded from the sample, until it reaches a minimum and starts rising again due to small number statistics. The two samples become fully compatible (down to ∼1σf or 10% flux variation) for 120 sources. These numbers are presented for a flux cut of 0.4 Jy. Looping over different flux thresholds and repeating the above procedure indicated that this choice led to the maximum possible sample size. The 1.4 GHz flux distribution of our final sample and the flux-limited NVSS parent sample are shown in Fig. 2.

thumbnail Fig. 2.

1.4 GHz flux distribution of NVSS and CO-ARC. The orange histogram shows the distribution for all NVSS parent sample, while the blue histogram is for the NVSS with redshift information. The thin black lines indicate the re-shuffled samples of NVSS with flux variations up to 10%. The CO-ARC in green is well represented by both parent samples within the 10% errors.

We also show the same analysis for the CO-ARC only using the ALMA calibrators (without including the literature sources, squares in Fig. 1) and the CO-ARC with ALMA calibrators and literature sources (circles) for the total sample and for the respective redshift bins, indicated by different color codes in Fig. 1. The inclusion of literature sources helps to improve the statistics and to converge the CO-ARC sample to the theoretical D-statistic distribution of the shuffled NVSS.

It is noteworthy that our final sample is representative of the NVSS flux distribution also in each of the 3 individually examined z ranges. We have chosen to use the ranges 0.005 < z < 0.3, 0.3 < z < 1,1 < z < 2.5, to study the evolution of the galaxies in comparable time steps of about 3.5–4.5 Gyr. The CO-ARC sources of the z slices are from 1σf to 3σf away from the NVSS parent catalogue at the same z. For the sake of completion, we note that this also holds separately for the sources from the ARC and the sources from the literature. The only exception comes from the few CO-ARC sources at intermediate z, where the construction of the luminosity function is unfeasible due to small number statistics.

Having built a final sample that is representative of the flux-limited NVSS, we then calculated its completeness with respect to the parent survey number of sources in the sky as a function of z. For this purpose, we split both samples in even smaller z steps (of 0.77 Gyr) and present the comparison in Fig. 3. The top panel shows the z distribution of the two samples (with the flux-limited NVSS being divided by a factor of 25 for comparison purposes). It shows that the final CO-ARC is over-represented at z < 0.2 and under-represented at z ∼ 0.8, but otherwise follows the parent catalogue’s z distribution. The completeness correction, i.e., the z-dependent number with which we need to multiply the CO-ARC galaxies to make them as frequent (in z and in the sky) as the flux-limited NVSS, is shown in the bottom panel of Fig. 3 and is taken into account during the construction of the CO luminosity function in Sect. 4.3. There is one more correction to be taken into account in the luminosity function construction: the fact that not all NVSS sources have known z. On average, 1 on every 4 NVSS sources have z information. This correction is, in principle, also z-dependent, but its behaviour with z is unknown. We, thus, also multiply the number density of sources times 4, except for the local Universe (at z < 0.3), in which we assume the redshift acquisition to be complete.

thumbnail Fig. 3.

Redshift completeness of CO-ARC. Top: number of sources as a function of z in our final CO-ARC sample and in the NVSS parent catalogue (divided by 25 for comparison purposes). The NVSS distribution was scaled to a factor of 25 for displaying purposes. Bottom: completeness correction so that our final sample reproduces the number of sources in the flux-limited NVSS with z.

In Fig. 4, we show the distribution of the 120 sources of our survey in the sky plane projection. The symbol sizes are indicative of the full width at half power (FWHP) of the primary beam of the observations. The total survey area covered by our sample was estimated by summing the FWHP of all individuals observations, taking into account the dependency of the FWHP on the frequency of each observed CO line. This is feasible as the sources are distributed randomly in the sky with almost no overlap between them.

thumbnail Fig. 4.

Distribution of the sample in the sky in the aitoff projection. The color code corresponds to the redshift and the sizes are according to the FWHP of the primary beam of the observations. The dotted curve and the star indicate the Galactic plane and center.

Our final sample, consisting of 120 sources and being representative of the NVSS down to 0.4 Jy, contains 66 sources from the CO-ARC plus 54 literature sources. Details on the CO-ARC sources are given in Tables B.1 and 1, including the spectral line of interest per galaxy, the ID of the ALMA observations that we found and reduced per line, as well as the redshift and its source of origin per galaxy. For the bibliographic sources, all pertinent information is summarized in Table C.1.

Table 1.

CO properties of the 66 CO-ARC sample sources.

3. ALMA observations

3.1. Data reduction

The ALMA observations that were analyzed as part of the CO-ARC survey were carried out between Cycle 0 and Cycle 5. The sources were observed as flux, phase, bandpass, or polarization calibrators. The relevant properties of each observing run (e.g. scan intent, on-source integration time, CO transition, etc.) are reported in Table B.1.

The data of the sources observed as flux, phase, and polarization calibrators were reduced using the Common Astronomy Software Application (CASA; McMullin et al. 2007) pipeline, which automatically processes the data by performing standard calibration and flagging. Various pipeline versions were used, according to those adopted for the corresponding archival data reduction. When one of our targets was observed in more than one project, once calibrated separately according to the appropriate CASA version, the visibilities were then combined using the CASA task concat. In a minority of cases the projects were observed in the first cycles (0, 1 and early 2) and were calibrated using casa version < = 4.2.2, which computes the visibility weights in a different way and therefore they had to be corrected to the same definition using the CASA task statwt4. Bright quasars were typically used as flux calibrators, resulting in a standard 10% flux calibration uncertainty.

A different procedure was required for sample sources observed as bandpass calibrators, as any line emission in the source would be flattened by the default pipeline corrections. This is because the bandpass calibration corrects the observed visibilities for frequency-dependent amplitude and phase corruptions. Such corrections are calculated using as model the bandpass calibrator, which is typically a bright quasar (i.e. a point-like source in the sky). Due to the Fourier Transform properties, a point-like sky brightness distribution ideally has flat and constant amplitudes in the visibility plane. The bandpass calibration calculates the corrections needed to make the data as similar as possible to such ideal model. Therefore, any line emission in the spectrum of the bandpass calibrator would be treated as instrumental fluctuation and thus flattened. In such cases, the data were then manually calibrated using customised data reduction scripts, which were modified so that the bandpass calibrator was substituted by the phase (or flux) calibrator in the relevant calibration steps. A standard data reduction process was then carried out also in these cases.

3.2. Imaging

Where detected, the CO line emission was identified and isolated in the uv-plane using the CASA task uvcontsub. This latter forms a continuum model from linear fits in frequency to line-free channels, and then subtracts this model from the visibilities.

Deconvolution and imaging were performed using the CASA task tclean with Briggs weighting and two robust parameters: 0.5 and 2. The former was chosen as it provides a good trade-off between angular resolution and signal-to-noise ratio (S/N); the latter the best sensitivity. For each type of weighting, we produced three data cubes with different channel widths: 20, 60 and 100 km s−1. For low spectral resolution observations in which the raw channel width was larger than 20 or 60 km s−1, we adopted the lowest possible channel binning. The continuum-subtracted dirty cubes were iteratively cleaned in regions with line emission (first identified through auto-masking) down to selected thresholds, and then primary-beam corrected.

The achieved synthesised beams and root mean square (rms) noise levels of both detections and non-detections are listed in Table 1. The median rms noise obtained for the emission line data cubes of ALMA calibrators is 0.48 mJy and some scatter is seen along this value, however, overall the rms noises per ALMA band are fairly homogenous and below 1 mJy for most of the targets, as we can see in Fig. 5 the distribution of σrms per frequency (or per ALMA band). We note that even though some of the CO-ARC sources have been presented in the literature before, we opted to re-image them, so that our archival survey has homogeneous imaging procedure.

thumbnail Fig. 5.

Distribution of achieved rms sensitivities per ALMA band. The black circles represent the rms noise levels (in log scale and mJy units) achieve for the emission line datacubes of the ALMA calibrators sample in function of the respective observed frequencies (in GHz) or ALMA bands highlighted in colors. The dashed line indicates the median value of 0.482 mJy for the 66 sources drawn from the CO-ARC.

4. Results

4.1. Moment maps and spectra of CO-ARC detections

Overall, we report the detection of CO emission in 17 out of the 66 CO-ARC sources listed in Table 1. A marginal detection for J2341+0018 (that most likely reflects a potential flux loss/bad calibration of its data). The remaining 49 analyzed sources are undetected in CO. Together with the 25 CO detections in the 54 literature sources listed in Table C.1, we find that 35% of radio galaxies in the extended CO-ARC sample contain detectable gas reservoirs. The fraction differs for local and high-z sources, varying from 31/60 sources in the local Universe to 2/23, and 9/37 sources at intermediate and high-z, respectively.

Integrated intensity (moment 0), mean line-of-sight velocity (moment 1), and velocity dispersion (moment 2) maps of the CO-ARC detections were created from the continuum-subtracted, cleaned data cubes, keeping pixels down to 3σ noise levels. The resulting maps are shown in Fig. 6 and Fig. A.1.

thumbnail Fig. 6.

Example of integrated intensity (moment 0; left panels), mean line-of-sight velocity (moment 1; middle-left panels), velocity dispersion (moment 2; middle-right) maps and spectra (right panels) of 2 detections in the CO-ARC sample. The moments maps for all the 17 detections are shown in the Appendix in Fig. A.1. The synthesized beam is shown as a blue ellipse in the bottom-left corner of each moment 0 map. The bar to the top of each moment map shows the color scales (in Jy beam−1 km s−1 and km s−1 for moment 0 and moment 1/2 maps, respectively). East is to the left and north to the top. The observed CO transition is indicated in the top-right corner of each spectral profile. The aperture within which the illustrated spectra have been extracted are overlain in green in each moment 0 maps.

Based on the integrated intensity maps (left panels of Figs. 6 and A.1), spatially resolved extended gas distributions are observed in six targets (NGC 315, IC 4374, NGC 3100, J0623–6436, NGC 6328, and J2009–4349), most of which presenting disc- or ring-like structures. Marginally resolved or unresolved gas structures are instead observed in all the other cases. The mean line-of-sight velocity maps (Figs. 6 and A.1, middle-left panels) show clear evidence of gas rotation in eight cases, that are NGC 315, IC 4296, NGC 3557, NGC 3100, J0623–6436, NGC 6328, J2009–4849, and NGC 7075. Among these, clear kinematic distortions (i.e. s-shaped iso-velocity contours) are observed in IC 4296, NGC 3100, NGC 6328. Such disturbances usually trace the presence of unrelaxed sub-structures in the gas distribution, possibly caused by either warps or non-circular motions, or a combination of both.

Accurate 3D kinematic modelling can help differentiating between the two. This has been carried out in details for IC 4296 and NGC 3100 by Ruffa et al. (2019b, 2022). In the former case, such an analysis led to the conclusion that the presence of a position angle warp best reproduces the observed kinematic asymmetries. The additional presence of non-circular motions, however, cannot be fully excluded nor fully established due to the fact that the gas distribution is only marginally resolved in IC 4296. The features observed in NGC 3100 are instead found to be the result of a more complex kinematics, including the presence of both a position angle and inclination warp and non-circular motions. These latter have been recently partly identified into a mild jet-induced outflow with Ṁ = 0.12 M⊙ yr−1 (Ruffa et al. 2022).

Another case of full kinematic modeling was that of NGC 6328, which has been carried out in the context of the CO-ARC project and can be found in separate papers (Papachristou et al. 2021, in prep.). In short, our detailed analysis suggests the presence of a jet-induced of outflow as likely explanation for the kinematics observed in NGC 6328, in addition to a highly warped gas distribution. A massive molecular outflow with Ṁ = 2300 M⊙ yr−1 has been also recently reported by Vayner et al. (2017) in the galaxy 3C 298.

The line-of-sight gas velocity dispersion (σgas; Figs. 6 and A.1, middle-right panels) varies significantly among different detections and across the gas sky distribution observed in each source. In most of the detections with disc- or ring-like shapes (e.g. NGC 315, J1505+0326, IC 4296, NGC 3557, NGC 3100, etc.), the gas is observed to be dynamically cold towards outskirts of CO distribution, with σgas ∼ 10 − 20 km s−1 in agreement with the typical ranges expected in unperturbed cold gas clouds (e.g. van de Voort et al. 2018). σgas then progressively increases towards the central gas regions, with the largest values (≥50 km s−1) typically observed around the location of the line peak(s). In first instance, such large values would imply the presence of turbulence within the gas clouds. Observed line-of-sight velocity dispersions, however, can be significantly overestimated due to observational effects. Among these, the limited spatial resolution and/or low detection significance can introduce severe beam smearing (i.e. contamination from partially resolved velocity gradients within the host galaxy), which usually dominates in areas of large velocity gradients (such as the central regions of galaxies). Full kinematic modelling is again the only way to derive reliable estimates of the intrinsic gas velocity dispersion: this has been carried out in details for NGC 315 by Boizelle et al. (2021) and for IC 4296, NGC 3557, NGC 3100, and NGC 7075 by Ruffa et al. (2019b, 2022).

The spatially integrated CO spectra of the 17 CO-ARC detections are shown in the right panels of Fig. A.1. These have been extracted from the cleaned data cubes within circular or elliptical apertures enclosing all the observed CO emission (as illustrated in the left panels of Figs. 6 and A.1). All of the eight above-mentioned sources showing clear signs of rotation in their moment 1 maps also exhibit the classic double-horned spectral shape expected from a rotating disc. In addition, in IC 4296 a strong absorption feature is observed. This has been analyzed in detail by Ruffa et al. (2019a) and mainly interpreted as absorption against a bright core continuum. In all the other detections, the spectral profiles show Gaussian-like, centrally peaked shapes. The integrated CO flux of each target, SCOΔv, has been estimated from the obtained spectral profiles and, together with the line full width at half-maximum (FWHM), is listed in Table 1.

For non-detections, upper limits were calculated in a customary Poissonian way multiplying 3σrms (i.e. three times the root mean square noise level determined from line-free channels of each data cube) by the square root of the number of channels contained in an assumed line of FWHM = 300 km s−1. We used the σrms value of the data cubes made with robustness parameter equal to 2 and 100 km s−1channel width. This is valid if we consider that the whole CO emission is concentrated within the synthesized beam, which is likely the case for the most distant sources. For nearby targets, to take into account possibly extended emission, we assumed a typical galaxy size of 10 kpc and corrected the previous noise by multiplying it by the square root of the number of synthesized beams inside the assumed galaxy size.

4.2. Luminosity and molecular mass

A common way to express the CO luminosity is via in units of K km s−1 pc2, which is proportional to the surface brightness temperature, TB. Its advantage comes from its definition: since it originates from the Rayleigh-Jeans flux multiplied by the beam area, the frequency squared in the Rayleigh-Jeans law gets cancelled out with 1/ν2 of the beam area. More specifically, the equation of Solomon & Vanden Bout (2005) provides from the integrated fluxes SCOΔV as:

(1)

where νobs is the observed frequency of the CO line, in GHz, DL is the luminosity distance of the galaxy, in Mpc, z is the redshift and SCOΔV is the integrated flux in Jy km s−1. The mid-J CO transitions can be converted to adopting a line ratio, called . For thermalized, optically thick gas emission, rJ → 1 = 1. We adopted the rJ → 1 values reported by Carilli & Walter (2013), r2 → 1 = 0.99, r3 → 1 = 0.97 and r4 → 1 = 0.87, for typical gas excitation conditions of QSOs. The values derived for the CO-ARC sources are in the range 6.5 < log  < 11.1, while the values reported for sources in the literature sample are 5.4 < log  < 10.7.

Molecular masses can be calculated from using the CO-to-H2 conversion factor . We adopted the standard Galactic CO-to-H2, conversion factor of αCO = 4.36 M⊙ (K km s−1 pc2)−1 (Tacconi et al. 2013; Bolatto et al. 2013) for consistency throughout the sample. Based on these assumptions, we derived molecular masses in the range of 1.5 × 107 < MH2 < 7.2 × 1011 M⊙ for the 66 CO-ARC radio galaxies. For the 54 literature-based galaxies, including detections and upper limits, the reported MH2 values vary from 1 × 106 to 4 × 1011 M⊙, converting to the same αCO = 4.36 used in this work for consistency.

Figure 7 illustrates the molecular masses as a function of z for all the sources analyzed in this work. We compare our mass measurements from semi-empirical predictions of the typical H2 gas content of galaxy halos from Popping et al. (2015), with halo masses of Mvir = 1013 and Mvir = 1014. Interestingly, we find that even though at the very local Universe all radio galaxies have molecular gas contents significantly (1-2 orders of magnitude) below those predicted for typical galaxy halo values, the picture changes at high-z. At z > 1, several radio galaxies have molecular gas reservoirs that are very similar to (and even higher than) those of the typical galaxy halo contents. The detected molecular gas reservoirs of radio galaxies at high-z ranges from 5 × 109 to ∼1012 M⊙, while at low z it ranges from 107 to 1010 M⊙. Still, at all z, the majority of galaxies have small (undetected) reservoirs than that of typical galaxy halos.

thumbnail Fig. 7.

Molecular masses derived for the CO-ARC sample. Detection are shown with filled circles and upper limits are indicated with open circles and arrows. For comparison, we show the predicted evolution of the molecular gas content with redshift from semi-empirical models of normal star-forming galaxies with halo masses of Mvir = 1013 and Mvir = 1014 (Popping et al. 2015).

4.3. CO luminosity function

Previous studies reported luminosity functions for normal and starburst galaxies in the local Universe (Keres et al. 2003; Obreschkow & Rawlings 2009; Saintonge et al. 2017), however this is the first time that the CO luminosity function of radio galaxies is presented. To construct the luminosity function, we adopted the 1/Vmax method of (Schmidt 1968), which provides the number density for each luminosity bin:

(2)

where Δlog L is the luminosity bin size (common in logarithmic scale), Vmax is the maximum comoving volume that any galaxy i could be observable given its measured 1.4 GHz flux and the survey limit and ωi is the weight shown in Fig. 3 corresponding to the completeness correction with respect to the NVSS.

We present the luminosity function for two different redshift ranges. For the range 0.005 < z < 0.3 (average ⟨z⟩ = 0.07), comprising 60 radio galaxies, and for the range 1 < z < 2.5 (with average ⟨z⟩ = 1.53), containing 37 radio galaxies.

A luminosity function at intermediate redshifts is not presented due to the small number of detections. For either redshift range, we present both the luminosity function measured using detections only, and an alternative version using the combination of detections and upper-limits for non-detections. The H2 luminosity (and consequently mass) function is shown in the left panel of Fig. 8 for z < 0.3 and in the right panel of the same figure for 1 < z < 2.5.

thumbnail Fig. 8.

Luminosity functions for 0.005 < z < 0.3 in the left panel and for 1 < z < 2.5 in the right panel. For the local mass and luminosity functions we compare our findings with those of different surveys: the XCOLD GASS (Saintonge et al. 2017; Fletcher et al. 2021, open diamonds) as the gray dots the FIR-selected CO survey of Keres et al. (2003) from the IRAS bright galaxy sample, the empirical CO luminosity function derived from the Herschel IR luminosity (Vallini et al. 2016, black line) and the light orange boxes are results from the ASPECS LP (Decarli et al. 2019). The SHARK semi-analytic models by Lagos et al. (2020) for the CO luminosity function of sub-mm selected galaxies are shown in dotted lines. Predictions from the SPRITZ simulation (Bisigello et al. 2021) are show in orange for all galaxies, and for the different AGN populations in yellow for SF-AGN, in light blue for SB-AGN, and for AGN1 and AGN2 are indicated in red and dark blue, respectively.

For the local luminosity function, we have enough number statistics and bins to fit the datapoints in logarithmic scale with an analytical Schechter (1976) function, as defined in Riechers et al. (2019):

(3)

In this equation, Φ(L′) is the number of galaxies per comoving volume per luminosity bin, Φ* is the scale density of galaxies per unit volume; L* is the scale luminosity defining the knee of the LF and α is the slope of the faint end of the luminosity function. The best Schechter function fit to our data is shown in Figure 8 as a dashed line, with 95% confidence intervals shown as shaded areas. We found a negative slope of α = −0.6, scale luminosity L* of 5.83 × 109 K km s−1 pc2, and turn over at Φ* = 2.39  ×  10−7 Mpc−3dex−1.

For the local Universe, we compared our results with those from the eXtended CO Legacy Database for GASS (XCOLD GASS) survey (Saintonge et al. 2017; Fletcher et al. 2021), the CO survey of FIR-selected galaxies from the IRAS bright galaxy sample (Keres et al. 2003), and the empirical CO luminosity function derived from the Herschel IR luminosity (Vallini et al. 2016). Figure 8 shows that the gas mass content locked up in local radio galaxies is ∼2–4 orders magnitude less than that in normal star-forming galaxies, with the gap closing at high masses. The same applies when comparing the CO-ARC results to those of the CO observations in the blind ALMA Spectroscopic Survey of the Hubble Ultra Deep Field (ASPECS LP; Decarli et al. 2019; orange shaded boxes in Fig. 8) for the average z of each bin. ASPECS presents CO luminosity functions for galaxies of any kind, from the local Universe to z ≲ 4 and covers an area of 4.6 arcmin2. Deviations in the (bright-end) shape are seen for the sub-millimeter selected galaxies from the SHARK semi-analytical model predictions by Lagos et al. (2020).

The comparison with the spectro-photometric Realisations of Infrared-selected Targets at all-z (SPRITZ) simulation predictions (Bisigello et al. 2021) turns out to be very instructive, as SPRITZ galaxies are separated in different populations including: unobscured and obscured AGN-dominated systems (hereafter AGN1 and AGN2, respectively) with little star formation, and composite systems of AGN and star-forming galaxies. One population of composite systems is similar to spirals hosting a low luminostity AGN (SF-AGN) and the other population resembles starburst galaxies hosting an obscured AGN (SB-AGN).

The CO luminosity function in the local Universe interestingly shows that the number of radio galaxies with H2 gas detection is only little lower than that of pure type 1 or 2 AGN. The comparison between our findings and the SPRITZ simulation at low-z shows that one on every ∼four RGs (i.e. a comparable number of RGs) have similar gas reservoirs to pure type 1 or type 2 AGN. Surprisingly, even the number of RGs with large gas reservoirs at the bright-end of the COLF (> 1010 M⊙) is not considerably lower, so RGs are not entirely red and dead as it is often believed to be the case. The agreement between the molecular gas content in RGs and the AGN-dominated populations with little star formation (AGN1 and AGN2) from SPRITZ point out that we are comparing AGN populations where the star formation does not dominated the SED, since our q24 selection criteria described in Sect. 2.1 assures that we selected galaxies whose radio emission is associated to black hole accretion instead of star formation. A deficit of gas at high masses is only seen for the very bright comparison sample, made to be representative of the 6C radio survey, flux limited to 3 Jy (green points of Fig. 8, the bright sample properties are listed in Table C.2). At low CO luminosities, the SPRITZ simulations do not have predictions for L′< 107.25 K km s−1 pc2, although the shape of our Schechter function fit appears to agree with an extrapolated curve of the SPRITZ AGN1 and AGN2 predictions at lower luminosities.

At high-z, the gas mass function of detected RGs with CO emission is lower (by up to 2.5 orders of magnitude) than that of local RGs with a corresponding mass content. It is even lower compared to other AGN from the SPRITZ survey, and all galaxies in the ASPECS LP (Decarli et al. 2019). This discrepancy largely originates from the fact that only a small number of RGs from the 1.4 GHz luminosity function are detectable in radio wavelengths at high-z. One of possible reason is due to the fact that the continuum spectrum of RGs at ∼GHz frequencies decreases with frequency and therefore it suffers from a strong K-correction, contrary to the CO lines and continuum at mm wavelength, whose fluxes increase with frequency, leading to a negative K-correction and the possibility to be detected more easily at high-z. Another explanation comes from our survey flux-limit of 0.4 Jy, implying that only the most extremely luminous galaxies at the bright-end of the radio luminosity function can fit this criteria (see Fig. 9 and text below). In other words, bright RGs are rarer to find, so our result needs to be roughly corrected with respect to the number density of undetectable RGs.

thumbnail Fig. 9.

Radio luminosity function at 1.4 GHz for the low (left), intermediate (middle) and high (right) redshift bins. We also report the luminosity function of Condon et al. (2002) in the local Universe for star-forming and AGN population as open squares and stars, respectively. The SPRITZ predictions for radio-loud AGN are shown in gray diamonds (Bisigello et al. 2021), as well as the 1.4 GHz luminosity function from the radio AGN survey in the VLA-COSMOS field (Smolčić et al. 2017, open green squares) for each redshift bin.

A rough attempt to correct this bias is made from the fraction of RGs in the inaccessible part of the 1.4 GHz luminosity function, assuming a similar gas mass distribution throughout. For this reason, we created the 1.4 GHz luminosity function of each z bin (Fig. 9) and identified what part of it our observations cover. At low-z, our survey covers a large part of the luminosity range of the SPRITZ simulations or of the NVSS and UGC galaxy samples Condon et al. (1998, 2002). However, at high-z, the integral of the covered 1.4GHz luminosity function is much lower than that of the SPRITZ simulations, indicating that we only see the brightest of every ∼7000 galaxies. Therefore ∼7000 galaxies are missed, and based on our previous results, 35% of them would, on average, contain detectable gas fractions. When accounting for the ∼2000 missing faint RGs in this way, we find that the gas mass content would be significantly higher than that in the local Universe, and comparable to that of other populations at high-z.

5. Discussion: Evolution of the MH2 content of gas in radio galaxies with z

By integrating L′Φ(L′) over all luminosity function bins and converting it into mass, we obtain the molecular gas content per comoving Mpc3 of the Universe, ρ(MH2), at different epochs (Fig. 10). Indeed, this plot shows that the derived molecular gas content of bright enough galaxies to be detected at 1.4 GHz (circles) shows a decrease with z due to the Malmquist bias. When correcting for the density of missing RGs as indicated above, i.e., by multiplying the derived high-zρ(MH2) by ∼2000, then the corrected value is less than an order of magnitude below that of AGN1 plus AGN2, as in the local Universe. So, under the rough assumption that detectable and undetectable RGs have a similar gas mass distribution, the amount of gas locked-up in RGs and AGN is similar at z > 1 – as the local Universe. Likewise, in the intermediate z, we only see the brightest of every 5000 galaxies, from which we deduce the appropriate ρ(MH2) correction. However, we draw no further conclusions from this bin, as it suffers from small number statistics (only 2 detections out of 23 sources there).

thumbnail Fig. 10.

Evolution of the cosmic molecular gas density, ρ(MH2), with redshift. Our results derived directly from the detection of the COLF are shown as the purple points and corresponding uncertainties boxes. The extrapolation of the number density of missing sources in the intermediate and high z bins due to our flux limit is shown in the filled triangles. The ρ(MH2) for the representative sample of the 6C catalog is indicated in the green box, and for the very bright end of our radio luminosity function (L1.4 GHz > 3 × 1025 W Hz−1) is shown in the purple star. The light orange boxes indicate the results from the ASPECS LP survey (Decarli et al. 2019) and the dotted box from the COLDz (Riechers et al. 2019). The dark and light blue lines represent the predictions from SPRITZ from all AGN hosts, including SF-AGN and SB-AGN, and the AGN-dominated (only AGN1+AGN2) populations.

We further investigated our findings by only examining the most luminous objects at different epochs. Due to the z evolution of luminosities (Pracy et al. 2016), we avoided directly comparing equal-L galaxies at different z. Instead we compared the gas mass content of the brightest of every few thousand RGs at any given z. For comparison, we cut the local 1.4 GHz luminosity function at 9 × 1025 W Hz−1, thus using only the brightest of every 7000 local galaxies. We then reconstructed the CO luminosity function and local Universe gas mass density using these galaxies. The outcome of this computation is the open star in Fig. 10. We find that the same amount of molecular gas is locked up in the brightest 1/5000–1/7000 radio galaxies at any given redshift. This result indicates that the net effect of intergalactic medium inflows on one hand and black hole feedback and star formation on the other hand, integrated over time, leads to the same amount of residual molecular gas in the brightest radio galaxies. It is, thus, plausible that radio-mode feedback indeed acts as a “maintenance” mode. To further check this, an independent measurement of the star formation rate in these galaxies is needed. This will be presented in a future paper from the fitting of the IR data.

6. Summary

To quantify the gas mass content of radio galaxies at different epochs, we examined a sample of 120 radio galaxies from the local Universe up to z ∼ 2.5. Our sample was meticulously constructed to be complete, representative of the 1.4 GHz NVSS survey when flux limited to 0.4 Jy. All galaxies in it have radio emission associated with black hole activity, as indicated by a radio-to-infrared ratio q24μm < 0.5. Of the sources, 66 were ALMA calibrators, originating from the ALMA Radio-Source Catalogue, and were often observed as such. Another 54 sources, often comprising more extended radio emission, were compiled from the literature. For all galaxies, CO (1–0) up to (4–3) data were used for the gas mass content derivation. Our main findings can be summarized as follows.

  • We detected CO emission in 35% of the sample galaxies, with molecular masses ranging from 1 × 107 < MH2 < 7 × 1010 M⊙ at low-z and at high-z. Two sources with outflows were seen among the 17 sources with newly reduced ALMA data.

  • One quarter of the 1 < z < 2.5 radio galaxies contain at least as much molecular gas as a typical galaxy halo of 1013–1014 M⊙. Therefore, a large number of radio galaxies at high-z would reside at a “normal” or even starbursty for that epoch host. Still, at any z, most radio galaxies are more depleted and thus evolved than the typical galaxy of a simulated halo.

  • The CO luminosity function in the local Universe interestingly shows that the number of radio galaxies with H2 gas detection is only little lower than that of pure type 1 or 2 AGN.

  • At 1 < z < 2.5, the CO luminosity function is significantly lower than that predicted by simulations for other galaxy types, because only hard-to-find bright radio galaxies are detectable. A rough correction of the fraction of missing sources of the 1.4 GHz luminosity function implies that the CO luminosity function would be 2000 times higher. This would close the gap between the various high-z populations and it would imply that the gas mass locked up in high-z galaxies is considerably higher than in the local Universe.

  • The gas mass content of the brightest 1/5000-1/7000 radio galaxies at any z is similar. It is plausible that these molecular gas reservoirs are maintained at similar levels thanks to both the consumption of gas due to star formation and the delay of inflows from the intergalactic medium due to radio activity.

This work demonstrate a potential case of exploring the wealth of the ALMA archival observations, including calibrators, for scientific purposes. The CO luminosity function of radio galaxies derived here might be useful for benchmarking of cosmological simulations and constraining the gaseous content of radio galaxies at the local Universe up to the star formation/activity history peak at z ≲ 2.5, where feedback could be more crucial for galaxy evolution.

2

The “classical” definition of radio loudness via the R indicator, the ratio between the rest-frame radio-to-optical flux density with typical values of R ∼ 10 characterizing RGs (Kellermann et al. 1989), is often insufficient to identify radio-quiet (RQ) objects, because both star-forming and RQ galaxies can have low R values.

3

Despite using the NVSS as radio parent sample, a derivation of a luminosity function with another sample representative of the 6th Cambridge Survey of radio sources (6C survey) is also discussed in Sect. 4.3, for comparison purposes.

Acknowledgments

We thank the referee for the constructive review and comments. KMD, AA, MP, JAFO, and IR acknowledge financial support by the Hellenic Foundation for Research and Innovation (HFRI), under the first call for the creation of research groups by postdoctoral researchers that was launched by the General Secretariat For Research and Technology under project number 1882 (PI Dasyra, Title: “Do massive winds induced by black-hole jets alter galaxy evolution? Evidence from galaxies in the ALMA Radio-source Catalog (ARC)”). AA also acknowledges the projects “Quantifying the impact of quasar feedback on galaxy evolution”, with reference EUR2020-112266, funded by MICINN-AEI/10.13039/501100011033 and the European Union NextGenerationEU/PRTR, and from the Consejería de Economía, Conocimiento y Empleo del Gobierno de Canarias and the European Regional Development Fund (ERDF) under grant “Quasar feedback and molecular gas reservoirs”, with reference ProID2020010105, ACCISI/FEDER, UE. JAFO acknowledges the financial support from the Spanish Ministry of Science and Innovation and the European Union – NextGenerationEU through the Recovery and Resilience Facility project ICTS-MRR-2021-03-CEFCA. IR also acknowledges support from the UK Science and Technology Facilities Council through grant ST/S00033X/1. This paper makes use of the ALMA data with project IDs listed in Table B.1. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. We made use of the NASA/IPAC Extragalactic Database (NED), and of the VizieR catalogue access tool, CDS, Strasbourg, France. This research made use of Astropy, a community developed core Python package for Astronomy.

References

  1. Aalto, S., Falstad, N., Muller, S., et al. 2020, A&A, 640, A104 [EDP Sciences] [Google Scholar]
  2. Adelman-McCarthy, J. K., Agüeros, M. A., Allam, S. S., et al. 2008, ApJS, 175, 297 [Google Scholar]
  3. Alatalo, K., Lacy, M., Lanz, L., et al. 2015, ApJ, 798, 31 [Google Scholar]
  4. Albareti, F. D., Allende Prieto, C., Almeida, A., et al. 2017, ApJS, 233, 25 [Google Scholar]
  5. Allington-Smith, J. R., Peacock, J. A., & Dunlop, J. S. 1991, MNRAS, 253, 287 [Google Scholar]
  6. Audibert, A., Combes, F., García-Burillo, S., et al. 2019, A&A, 632, A33 [Google Scholar]
  7. Baker, J. C., Hunstead, R. W., Kapahi, V. K., & Subrahmanya, C. R. 1999, ApJS, 122, 29 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bennett, C. L., Lawrence, C. R., Burke, B. F., Hewitt, J. N., & Mahoney, J. 1986, ApJS, 61, 1 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bieri, R., Dubois, Y., Rosdahl, J., et al. 2017, MNRAS, 464, 1854 [Google Scholar]
  10. Bisigello, L., Gruppioni, C., Feltre, A., et al. 2021, A&A, 651, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Boizelle, B. D., Walsh, J. L., Barth, A. J., et al. 2021, ApJ, 908, 19 [Google Scholar]
  12. Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207 [Google Scholar]
  13. Bonzini, M., Padovani, P., Mainieri, V., et al. 2013, MNRAS, 436, 3759 [Google Scholar]
  14. Bower, R. G., Benson, A. J., Malbon, R., et al. 2006, MNRAS, 370, 645 [Google Scholar]
  15. Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105 [Google Scholar]
  16. Castignani, G., Combes, F., Salomé, P., et al. 2019, A&A, 623, A48 [Google Scholar]
  17. Cielo, S., Bieri, R., Volonteri, M., Wagner, A. Y., & Dubois, Y. 2018, MNRAS, 477, 1336 [Google Scholar]
  18. Combes, F., García-Burillo, S., Casasola, V., et al. 2013, A&A, 558, A124 [Google Scholar]
  19. Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693 [Google Scholar]
  20. Condon, J. J., Cotton, W. D., & Broderick, J. J. 2002, AJ, 124, 675 [Google Scholar]
  21. Croom, S. M., Smith, R. J., Boyle, B. J., et al. 2004, MNRAS, 349, 1397 [Google Scholar]
  22. Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11 [Google Scholar]
  23. Dabhade, P., Combes, F., Salomé, P., Bagchi, J., & Mahato, M. 2020, A&A, 643, A111 [Google Scholar]
  24. Dasyra, K. M., Combes, F., Oosterloo, T., et al. 2016, A&A, 595, L7 [Google Scholar]
  25. De Breuck, C., Neri, R., & Omont, A. 2003, New Astron. Rev., 47, 285 [Google Scholar]
  26. De Breuck, C., Downes, D., Neri, R., et al. 2005, A&A, 430, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Decarli, R., Walter, F., Gónzalez-López, J., et al. 2019, ApJ, 882, 138 [Google Scholar]
  28. Di Matteo, T., Colberg, J., Springel, V., Hernquist, L., & Sijacki, D. 2008, ApJ, 676, 33 [NASA ADS] [CrossRef] [Google Scholar]
  29. Drinkwater, M. J., Webster, R. L., Francis, P. J., et al. 1997, MNRAS, 284, 85 [NASA ADS] [Google Scholar]
  30. Dubois, Y., Peirani, S., Pichon, C., et al. 2016, MNRAS, 463, 3948 [Google Scholar]
  31. Emonts, B. H. C., Feain, I., Mao, M. Y., et al. 2011, ApJ, 734, L25 [Google Scholar]
  32. Emonts, B. H. C., Norris, R. P., Feain, I., et al. 2014, MNRAS, 438, 2898 [Google Scholar]
  33. Evans, A. S., Sanders, D. B., Mazzarella, J. M., et al. 1996, ApJ, 457, 658 [NASA ADS] [CrossRef] [Google Scholar]
  34. Evans, A. S., Mazzarella, J. M., Surace, J. A., et al. 2005, ApJS, 159, 197 [NASA ADS] [CrossRef] [Google Scholar]
  35. Fabian, A. C. 2012, ARA&A, 50, 455 [Google Scholar]
  36. Fernández-Ontiveros, J. A., Dasyra, K. M., Hatziminaoglou, E., et al. 2020, A&A, 633, A127 [Google Scholar]
  37. Fletcher, T. J., Saintonge, A., Soares, P. S., & Pontzen, A. 2021, MNRAS, 501, 411 [Google Scholar]
  38. Fosbury, R. A. E., Mebold, U., Goss, W. M., & van Woerden, H. 1977, MNRAS, 179, 89 [Google Scholar]
  39. Fotopoulou, C. M., Dasyra, K. M., Combes, F., Salomé, P., & Papachristou, M. 2019, A&A, 629, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Fricke, K. J., Kollatschny, W., & Witzel, A. 1983, A&A, 117, 60 [NASA ADS] [Google Scholar]
  41. García-Burillo, S., Combes, F., Usero, A., et al. 2014, A&A, 567, A125 [Google Scholar]
  42. Gattano, C., Lambert, S. B., & Le Bail, K. 2018, A&A, 618, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Ginsburg, A., Sipőcz, B. M., Brasseur, C. E., et al. 2019, AJ, 157, 98 [Google Scholar]
  44. Glikman, E., Helfand, D. J., White, R. L., et al. 2007, ApJ, 667, 673 [NASA ADS] [CrossRef] [Google Scholar]
  45. Greve, T. R., Ivison, R. J., & Papadopoulos, P. P. 2004, A&A, 419, 99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Griffith, M. R., Wright, A. E., Burke, B. F., & Ekers, R. D. 1994, ApJS, 90, 179 [NASA ADS] [CrossRef] [Google Scholar]
  47. Hales, S. E. G., Baldwin, J. E., & Warner, P. J. 1993, MNRAS, 263, 25 [NASA ADS] [Google Scholar]
  48. Harrison, C. M. 2017, Nat. Astron., 1, 0165 [Google Scholar]
  49. Healey, S. E., Romani, R. W., Cotter, G., et al. 2008, ApJS, 175, 97 [Google Scholar]
  50. Hewett, P. C., & Wild, V. 2010, MNRAS, 405, 2302 [NASA ADS] [Google Scholar]
  51. Hewitt, A., & Burbidge, G. 1989, ApJS, 69, 1 [NASA ADS] [CrossRef] [Google Scholar]
  52. Hook, I. M., Shaver, P. A., Jackson, C. A., Wall, J. V., & Kellermann, K. I. 2003, A&A, 399, 469 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Jauncey, D. L., Batty, M. J., Wright, A. E., Peterson, B. A., & Savage, A. 1984, ApJ, 286, 498 [NASA ADS] [CrossRef] [Google Scholar]
  54. Jones, D. H., Read, M. A., Saunders, W., et al. 2009, MNRAS, 399, 683 [Google Scholar]
  55. Katgert, P., Mazure, A., den Hartog, R., et al. 1998, A&AS, 129, 399 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Kellermann, K. I., Sramek, R., Schmidt, M., Shaffer, D. B., & Green, R. 1989, AJ, 98, 1195 [Google Scholar]
  57. Keres, D., Yun, M. S., & Young, J. S. 2003, ApJ, 582, 659 [NASA ADS] [CrossRef] [Google Scholar]
  58. King, A. 2003, ApJ, 596, L27 [Google Scholar]
  59. Lagos, C. d. P., da Cunha, E., Robotham, A. S. G., et al. 2020, MNRAS, 499, 1948 [Google Scholar]
  60. Landoni, M., Falomo, R., Paiano, S., & Treves, A. 2020, ApJS, 250, 37 [NASA ADS] [CrossRef] [Google Scholar]
  61. Landt, H., Padovani, P., Perlman, E. S., et al. 2001, MNRAS, 323, 757 [Google Scholar]
  62. Lanz, L., Ogle, P. M., Alatalo, K., & Appleton, P. N. 2016, ApJ, 826, 29 [NASA ADS] [CrossRef] [Google Scholar]
  63. Leon, S., Lim, J., Combes, F., & Trung, D. V. 2000, ApJ, 545, L93 [NASA ADS] [CrossRef] [Google Scholar]
  64. Leon, S., Lim, J., Combes, F., & Trung, D. V. 2003, in ASP Conf. Ser., 290, 525 [NASA ADS] [Google Scholar]
  65. Lim, J., Leon, S., Combes, F., & Trung, D. V. 2003, in ASP Conf. Ser., 290, 529 [NASA ADS] [Google Scholar]
  66. Lo, K. Y., Chen, H. W., & Ho, P. T. P. 1999, A&A, 341, 348 [NASA ADS] [Google Scholar]
  67. McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, ASP Conf. Ser., 376, 127 [Google Scholar]
  68. McNamara, B. R., & Nulsen, P. E. J. 2012, New J. Phys., 14, 055023 [Google Scholar]
  69. Morganti, R., Oosterloo, T., Oonk, J. B. R., Frieswijk, W., & Tadhunter, C. 2015, A&A, 580, A1 [Google Scholar]
  70. Mukherjee, D., Bicknell, G. V., Wagner, A. Y., Sutherland, R. S., & Silk, J. 2018, MNRAS, 479, 5544 [Google Scholar]
  71. Murphy, T., Sadler, E. M., Ekers, R. D., et al. 2010, MNRAS, 402, 2403 [Google Scholar]
  72. Nesvadba, N. P. H., Lehnert, M. D., De Breuck, C., Gilbert, A. M., & van Breugel, W. 2008, A&A, 491, 407 [Google Scholar]
  73. Nesvadba, N. P. H., Boulanger, F., Salomé, P., et al. 2010, A&A, 521, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Obreschkow, D., & Rawlings, S. 2009, MNRAS, 394, 1857 [Google Scholar]
  75. Ocaña Flaquer, B., Leon, S., Combes, F., & Lim, J. 2010, A&A, 518, A9 [Google Scholar]
  76. Olivares, V., Salome, P., Combes, F., et al. 2019, A&A, 631, A22 [Google Scholar]
  77. Oosterloo, T., Morganti, R., Tadhunter, C., et al. 2019, A&A, 632, A66 [Google Scholar]
  78. Owen, F. N., Ledlow, M. J., & Keel, W. C. 1995, AJ, 109, 14 [NASA ADS] [CrossRef] [Google Scholar]
  79. Papachristou, M., Dasyra, K. M., Fernández-Ontiveros, J. A., et al. 2021, Astron. Nachr., 342, 1160 [NASA ADS] [CrossRef] [Google Scholar]
  80. Papadopoulos, P. P., Röttgering, H. J. A., van der Werf, P. P., et al. 2000, ApJ, 528, 626 [NASA ADS] [CrossRef] [Google Scholar]
  81. Perlman, E. S., Padovani, P., Giommi, P., et al. 1998, AJ, 115, 1253 [Google Scholar]
  82. Peterson, B. A., Wright, A. E., Jauncey, D. L., & Condon, J. J. 1979, ApJ, 232, 400 [NASA ADS] [CrossRef] [Google Scholar]
  83. Peterson, S. D. 1979, ApJS, 40, 527 [CrossRef] [Google Scholar]
  84. Pietsch, W., Bischoff, K., Boller, T., et al. 1998, A&A, 333, 48 [NASA ADS] [Google Scholar]
  85. Pilkington, J. D. H., & Scott, J. F. 1965, MmRAS, 69, 183 [Google Scholar]
  86. Popping, G., Behroozi, P. S., & Peeples, M. S. 2015, MNRAS, 449, 477 [Google Scholar]
  87. Pracy, M. B., Ching, J. H. Y., Sadler, E. M., et al. 2016, MNRAS, 460, 2 [Google Scholar]
  88. Prandoni, I., Laing, R. A., de Ruiter, H. R., & Parma, P. 2010, A&A, 523, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Richter, P., Wakker, B. P., Fechner, C., et al. 2016, A&A, 590, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Riechers, D. A., Pavesi, R., Sharon, C. E., et al. 2019, ApJ, 872, 7 [Google Scholar]
  91. Ruffa, I., Davis, T. A., Prandoni, I., et al. 2019a, MNRAS, 489, 3739 [Google Scholar]
  92. Ruffa, I., Prandoni, I., Laing, R. A., et al. 2019b, MNRAS, 484, 4239 [Google Scholar]
  93. Ruffa, I., Prandoni, I., Davis, T. A., et al. 2022, MNRAS, 510, 4485 [Google Scholar]
  94. Russell, H. R., McNamara, B. R., Fabian, A. C., et al. 2019, MNRAS, 490, 3025 [Google Scholar]
  95. Saintonge, A., Catinella, B., Tacconi, L. J., et al. 2017, ApJS, 233, 22 [Google Scholar]
  96. Saripalli, L., & Mack, K. H. 2007, MNRAS, 376, 1385 [Google Scholar]
  97. Sbarufatti, B., Ciprini, S., Kotilainen, J., et al. 2009, AJ, 137, 337 [NASA ADS] [CrossRef] [Google Scholar]
  98. Schechter, P. 1976, ApJ, 203, 297 [Google Scholar]
  99. Schmidt, M. 1968, ApJ, 151, 393 [Google Scholar]
  100. Schmidt, M. 1977, ApJ, 217, 358 [NASA ADS] [CrossRef] [Google Scholar]
  101. Scoville, N. Z., Yun, M. S., Windhorst, R. A., Keel, W. C., & Armus, L. 1997, ApJ, 485, L21 [NASA ADS] [CrossRef] [Google Scholar]
  102. Shaw, M. S., Romani, R. W., Cotter, G., et al. 2013, ApJ, 764, 135 [NASA ADS] [CrossRef] [Google Scholar]
  103. Sijacki, D., Springel, V., Di Matteo, T., & Hernquist, L. 2007, MNRAS, 380, 877 [Google Scholar]
  104. Smith, R. J., Lucey, J. R., Hudson, M. J., Schlegel, D. J., & Davies, R. L. 2000, MNRAS, 313, 469 [Google Scholar]
  105. Smolčić, V., Novak, M., Delvecchio, I., et al. 2017, A&A, 602, A6 [Google Scholar]
  106. Solomon, P. M., & Vanden Bout, P. A. 2005, ARA&A, 43, 677 [Google Scholar]
  107. Son, D., Woo, J.-H., Kim, S. C., et al. 2012, ApJ, 757, 140 [NASA ADS] [CrossRef] [Google Scholar]
  108. Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74 [Google Scholar]
  109. Taniguchi, Y., Kameya, O., Nakai, N., & Kawara, K. 1990, ApJ, 358, 132 [NASA ADS] [CrossRef] [Google Scholar]
  110. Titov, O., Jauncey, D. L., Johnston, H. M., Hunstead, R. W., & Christensen, L. 2011, AJ, 142, 165 [NASA ADS] [CrossRef] [Google Scholar]
  111. Vallini, L., Gruppioni, C., Pozzi, F., Vignali, C., & Zamorani, G. 2016, MNRAS, 456, L40 [Google Scholar]
  112. van de Voort, F., Davis, T. A., Matsushita, S., et al. 2018, MNRAS, 476, 122 [Google Scholar]
  113. Vayner, A., Wright, S. A., Murray, N., et al. 2017, ApJ, 851, 126 [Google Scholar]
  114. Véron-Cetty, M. P., Woltjer, L., Staveley-Smith, L., & Ekers, R. D. 2000, A&A, 362, 426 [Google Scholar]
  115. Wagner, A. Y., & Bicknell, G. V. 2011, ApJ, 728, 29 [Google Scholar]
  116. Wagner, A. Y., Bicknell, G. V., & Umemura, M. 2012, ApJ, 757, 136 [CrossRef] [Google Scholar]
  117. White, G. L., Jauncey, D. L., Savage, A., et al. 1988, ApJ, 327, 561 [NASA ADS] [CrossRef] [Google Scholar]
  118. Wilkes, B. J. 1986, MNRAS, 218, 331 [Google Scholar]
  119. Wills, D., & Lynds, R. 1978, ApJS, 36, 317 [NASA ADS] [CrossRef] [Google Scholar]
  120. Wisotzki, L., Christlieb, N., Bade, N., et al. 2000, A&A, 358, 77 [NASA ADS] [Google Scholar]
  121. Wright, A., & Otrupcek, R. 1990, PKS Catalog [Google Scholar]
  122. Wright, A. E., Ables, J. G., & Allen, D. A. 1983, MNRAS, 205, 793 [NASA ADS] [Google Scholar]
  123. Zafar, T., Popping, A., & Péroux, C. 2013, A&A, 556, A140 [Google Scholar]

Appendix A: Moment maps and spectra

We present in Figure A.1 the integrated intensity, mean line-of-sight velocity and velocity dispersion maps of the 17 detections out of the 66 CO-ARC sources listed in Table 1 and discussed in Section 4.1.

thumbnail Fig. A.1.

Integrated intensity (moment 0; left panels), mean line-of-sight velocity (moment 1; middle-left panels), velocity dispersion (moment 2; middle-right) maps and spectra (right panels) of the 17 detections in the CO-ARC sample. The synthesized beam is shown as a blue ellipse in the bottom-left corner of each moment 0 map. The bar to the top of each moment map shows the color scales (in Jy beam−1 km s−1 and km s−1 for moment 0 and moment 1/2 maps, respectively). East is to the left and north to the top. The observed CO transition is indicated in the top-right corner of each spectral profile. The aperture within which the illustrated spectra have been extracted are overlain in green in each moment 0 maps.

thumbnail Fig. A.1.

continued.

thumbnail Fig. A.1.

continued.

thumbnail Fig. A.1.

continued.

Appendix B: Summary of the ALMA data

In Table B.1 we report main properties (name, coordinates, redshift and 1,4 GHz fluxes) of the the 66 ALMA CO-ARC sources and the respective ALMA projects, integration times and scan intents used in this work.

Table B.1.

CO-ARC sample ancillary information and information on the reduced ALMA observations.

Appendix C: Properties of the molecular gas in the literature data

We present in Table C.1 the properties of the molecular gas measurements for the 54 literature sources that are part of our CO-ARC sample.

Table C.1.

CO properties of 54 literature-based sources.

The properties of the bright sample representative of the 6C catalogue when flux limited to 3 Jy at 151 MHz are listed in Table C.2.

Table C.2.

CO properties of representative sample of the 6C catalogue.

All Tables

Table 1.

CO properties of the 66 CO-ARC sample sources.

In the text
Table B.1.

CO-ARC sample ancillary information and information on the reduced ALMA observations.

In the text
Table C.1.

CO properties of 54 literature-based sources.

In the text
Table C.2.

CO properties of representative sample of the 6C catalogue.

In the text

All Figures

thumbnail Fig. 1.

D-statistic, i.e., maximum distance between two cumulative distributions used for the sample selection versus the sample size. Solid lines show a theoretical comparison of the NVSS and a shuffled NVSS samples with 10 and 30% flux variations. Small black dots show how the sample D-statistic changes as a function of the number of bootstrapped sources (see text). The large circles show the final numbers for the CO-ARC (black) and for individual z ranges (color). The same analysis is shown for the CO-ARC only using the ALMA calibrators (without including the literature sources) is indicated as squares and equivalent color code for the redshift bin. We reach an optimal sample size of 120 (CO-ARC in black circle) for a D-statistics smaller than 10% and converging with the values computed to the theoretical shuffled NVSS sample down to 10% flux variation.
In the text
thumbnail Fig. 2.

1.4 GHz flux distribution of NVSS and CO-ARC. The orange histogram shows the distribution for all NVSS parent sample, while the blue histogram is for the NVSS with redshift information. The thin black lines indicate the re-shuffled samples of NVSS with flux variations up to 10%. The CO-ARC in green is well represented by both parent samples within the 10% errors.
In the text
thumbnail Fig. 3.

Redshift completeness of CO-ARC. Top: number of sources as a function of z in our final CO-ARC sample and in the NVSS parent catalogue (divided by 25 for comparison purposes). The NVSS distribution was scaled to a factor of 25 for displaying purposes. Bottom: completeness correction so that our final sample reproduces the number of sources in the flux-limited NVSS with z.
In the text
thumbnail Fig. 4.

Distribution of the sample in the sky in the aitoff projection. The color code corresponds to the redshift and the sizes are according to the FWHP of the primary beam of the observations. The dotted curve and the star indicate the Galactic plane and center.
In the text
thumbnail Fig. 5.

Distribution of achieved rms sensitivities per ALMA band. The black circles represent the rms noise levels (in log scale and mJy units) achieve for the emission line datacubes of the ALMA calibrators sample in function of the respective observed frequencies (in GHz) or ALMA bands highlighted in colors. The dashed line indicates the median value of 0.482 mJy for the 66 sources drawn from the CO-ARC.
In the text
thumbnail Fig. 6.

Example of integrated intensity (moment 0; left panels), mean line-of-sight velocity (moment 1; middle-left panels), velocity dispersion (moment 2; middle-right) maps and spectra (right panels) of 2 detections in the CO-ARC sample. The moments maps for all the 17 detections are shown in the Appendix in Fig. A.1. The synthesized beam is shown as a blue ellipse in the bottom-left corner of each moment 0 map. The bar to the top of each moment map shows the color scales (in Jy beam−1 km s−1 and km s−1 for moment 0 and moment 1/2 maps, respectively). East is to the left and north to the top. The observed CO transition is indicated in the top-right corner of each spectral profile. The aperture within which the illustrated spectra have been extracted are overlain in green in each moment 0 maps.
In the text
thumbnail Fig. 7.

Molecular masses derived for the CO-ARC sample. Detection are shown with filled circles and upper limits are indicated with open circles and arrows. For comparison, we show the predicted evolution of the molecular gas content with redshift from semi-empirical models of normal star-forming galaxies with halo masses of Mvir = 1013 and Mvir = 1014 (Popping et al. 2015).
In the text
thumbnail Fig. 8.

Luminosity functions for 0.005 < z < 0.3 in the left panel and for 1 < z < 2.5 in the right panel. For the local mass and luminosity functions we compare our findings with those of different surveys: the XCOLD GASS (Saintonge et al. 2017; Fletcher et al. 2021, open diamonds) as the gray dots the FIR-selected CO survey of Keres et al. (2003) from the IRAS bright galaxy sample, the empirical CO luminosity function derived from the Herschel IR luminosity (Vallini et al. 2016, black line) and the light orange boxes are results from the ASPECS LP (Decarli et al. 2019). The SHARK semi-analytic models by Lagos et al. (2020) for the CO luminosity function of sub-mm selected galaxies are shown in dotted lines. Predictions from the SPRITZ simulation (Bisigello et al. 2021) are show in orange for all galaxies, and for the different AGN populations in yellow for SF-AGN, in light blue for SB-AGN, and for AGN1 and AGN2 are indicated in red and dark blue, respectively.
In the text
thumbnail Fig. 9.

Radio luminosity function at 1.4 GHz for the low (left), intermediate (middle) and high (right) redshift bins. We also report the luminosity function of Condon et al. (2002) in the local Universe for star-forming and AGN population as open squares and stars, respectively. The SPRITZ predictions for radio-loud AGN are shown in gray diamonds (Bisigello et al. 2021), as well as the 1.4 GHz luminosity function from the radio AGN survey in the VLA-COSMOS field (Smolčić et al. 2017, open green squares) for each redshift bin.
In the text
thumbnail Fig. 10.

Evolution of the cosmic molecular gas density, ρ(MH2), with redshift. Our results derived directly from the detection of the COLF are shown as the purple points and corresponding uncertainties boxes. The extrapolation of the number density of missing sources in the intermediate and high z bins due to our flux limit is shown in the filled triangles. The ρ(MH2) for the representative sample of the 6C catalog is indicated in the green box, and for the very bright end of our radio luminosity function (L1.4 GHz > 3 × 1025 W Hz−1) is shown in the purple star. The light orange boxes indicate the results from the ASPECS LP survey (Decarli et al. 2019) and the dotted box from the COLDz (Riechers et al. 2019). The dark and light blue lines represent the predictions from SPRITZ from all AGN hosts, including SF-AGN and SB-AGN, and the AGN-dominated (only AGN1+AGN2) populations.
In the text
thumbnail Fig. A.1.

Integrated intensity (moment 0; left panels), mean line-of-sight velocity (moment 1; middle-left panels), velocity dispersion (moment 2; middle-right) maps and spectra (right panels) of the 17 detections in the CO-ARC sample. The synthesized beam is shown as a blue ellipse in the bottom-left corner of each moment 0 map. The bar to the top of each moment map shows the color scales (in Jy beam−1 km s−1 and km s−1 for moment 0 and moment 1/2 maps, respectively). East is to the left and north to the top. The observed CO transition is indicated in the top-right corner of each spectral profile. The aperture within which the illustrated spectra have been extracted are overlain in green in each moment 0 maps.
In the text
thumbnail Fig. A.1.

continued.
In the text
thumbnail Fig. A.1.

continued.
In the text
thumbnail Fig. A.1.

continued.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.