CO in the ALMA Radio-source Catalogue (ARC): The molecular gas content of radio galaxies as a function of redshift

To evaluate the role of radio activity in galaxy evolution, we designed a large archival CO survey of radio galaxies (RGs) to determine their molecular gas masses at di ﬀ erent epochs. We used a sample of 120 RGs representative of the NVSS 1.4GHz survey, when ﬂux limited at 0.4Jy. Of those, 66 galaxies belonged to the ALMA Radio-source Catalogue (ARC) of calibrators and had spectral window tunings around CO (1–0), (2–1), (3–2), or (4–3). We reduced their ALMA data, determined their H 2 mass contents, and combined the results with similar results for the remaining 54 galaxies from the literature. We found that, while at all epochs the majority of RGs have undetectable reservoirs, there is a rapid increase in the H 2 mass content of the CO-detected RGs with z . At 1 < z < 2 . 5, one-fourth of the RGs have at least as much molecular gas as simulations would indicate for a typical halo mass of that epoch. These galaxies plausibly have “normal” or even starbursty hosts. Overall, reservoirs of 10 7 (cid:46) M H 2 (cid:46) 10 10 M (cid:12) are seen at z < 0 . 3, and 10 10 (cid:46) M H 2 (cid:46) 10 12 M (cid:12) at z > 1. Taking into account the completeness correction of the sample, we created the corresponding H 2 mass functions at 0 . 005 < z < 0 . 3 and 1 < z < 2 . 5. The local mass function reveals that the number density of low-z RGs with detectable molecular gas reservoirs is only a little lower (a factor of ∼ 4) than that of pure (or little star-forming) type 1 and 2 AGN in simulations. At 1 < z < 2 . 5, there is a signiﬁcant decrease in the number density of high-z RGs due to the rarity of bright radio galaxies. An estimate for the missing faint RGs would, nonetheless, bring populations close again. Finally, we ﬁnd that the volume density of molecular gas locked up in the brightest 1 / 5000–1 / 7000 RGs is similar in the examined z bins. This result likely indicates that the inﬂow rate on one hand and the star-formation depletion rate plus the jet-driven expulsion rate on the other hand counteract each other in the most luminous RGs of each epoch.


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), Reduced datacubes are only available in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr(130.79.128.5) or via http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/ 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 multiphase 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 M H 2 ∼ 10 10 − 10 11 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. 2011Emonts et al. , 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 Section 2 we describe the statistics used to reproduce a representative sample of NVSS radio galaxies.In Section 3 we describe the data reduction of the ALMA observations.Main results from the analysis of the literature and ALMA observations are presented in Section 4. The presentation and discussion of the mass and luminosity functions are described Section 4.3, followed by conclusions.In this paper we adopt a flat ΛCDM cosmological model with H 0 =70 km s −1 Mpc −1 , Ω M =0.3, and Ω Λ =0.7.

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 nonhomogeneous manner.
The ARC catalogue was downloaded from the European Southern Observatory (ESO) archival interface 1 .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 25827 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 q 24 = log(S 24µm /S 1.4GHz ), 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.Lim 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 q 24 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. .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.

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 151MHz (Hales et al. 1993), the Parkes Radio Sources Catalogue (PKS) at 2.7GHz (Wright & Otrupcek 1990) and at 4.8 GHz (Griffith et al. 1994), the MIT-Green Bank 5GHz Survey (Bennett et al. 1986), and the Australia Telescope 20GHz 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 correspond-  We also show the same analysis for the CO-ARC only using the ALMA calibrators (without including the literature sources, squares in Figure 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 Figure 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.5Gyr.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 fluxlimited 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 Figure 3.The top panel shows the z distribution of the two samples (with the fluxlimited 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 Figure 3 and is taken into account during the construction of the CO luminosity function in Section 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.
In Figure 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.
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 in the Appendix.

Galaxy
Line Beam size

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

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

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 continuumsubtracted, cleaned data cubes, keeping pixels down to 3σ noise levels.The resulting maps are shown in Figure 6   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. (2019bRuffa et al. ( , 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 noncircular 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;Papachristou et al. 2022 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 ; Figure 6 and  Figure 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 ringlike 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. (2019bRuffa et al. ( , 2022)).
The spatially-integrated CO spectra of the 17 CO-ARC detections are shown in the right panels of Figure 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 Fig. 6 and Fig. A.1).All of the eight above-mentioned sources showing clear signs of rotation in their moment 1 maps also exhibit the classic doublehorned 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, S CO ∆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 −1 channel 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.

Luminosity and Molecular mass
A common way to express the CO luminosity is via L CO in units of K km s −1 pc 2 , which is proportional to the surface brightness temperature, T B .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 L CO from the integrated fluxes S CO ∆V as: where ν obs is the observed frequency of the CO line, in GHz, D L is the luminosity distance of the galaxy, in Mpc, z is the redshift and S CO ∆V is the integrated flux in Jy km s −1 .The mid-J CO transitions can be converted to L CO(1−0) adopting a line ratio, called r J→1 = L CO (J→J−1) L CO (1→0) .For thermalized, optically thick gas emission, r J→1 =1.We adopted the r J→1 values reported by Carilli & Walter (2013), r 2→1 =0.99, r 3→1 =0.97 and r 4→1 =0.87, for typical gas excitation conditions of QSOs.The values derived for the CO-ARC sources are in the range 6.5<logL CO(1−0) <11.1, while the values reported for sources in the literature sample are 5.4<logL CO(1−0) <10.7.
Molecular masses can be calculated from L CO(1−0) using the CO-to-H 2 conversion factor M H 2 = α CO L CO(1−0) .We adopted the standard Galactic CO-to-H2, conversion factor of α CO = 4.36M (K km s −1 pc 2 ) −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×10 7 < M H 2 < 7.2 × 10 11 M for the 66 CO-ARC radio galaxies.For the 54 literature-based galaxies, including detections and upper limits, the reported M H 2 values vary from 1×10 6 to 4×10 11 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 H 2 gas content of galaxy halos from Popping et al. (2015), with halo masses of M vir = 10 13 and M vir = 10 14 .Interestingly, we find 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 M vir = 10 13 and M vir = 10 14 (Popping et al. 2015).
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×10 9 to ∼10 12 M , while at low z it ranges from 10 7 to 10 10 M .Still, at all z, the majority of galaxies have small (undetected) reservoirs than that of typical galaxy halos.

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/V max method of (Schmidt 1968), which provides the number density for each luminosity bin: where ∆ log L is the luminosity bin size (common in logarithmic scale), V max 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 Figure 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 H 2 luminosity (and consequently mass) function is shown in the left panel of Figure 8 (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)  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×10 9 K km s −1 pc 2 , and turn over at 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 arcmin 2 .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 AGNs.Surprisingly, even the number of RGs with large gas reservoirs at the bright-end of the COLF (>10 10 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 q 24 selection criteria described in Section 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 Figure 8, the bright sample properties are listed in Table C.2 in the Appendix).At low CO luminosities, the SPRITZ simulations do not have predictions for L < 10 7.25 K km s −1 pc 2 , 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 AGNs 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 Kcorrection, 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 Figure 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.
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 (Figure 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. (1998Condon et al. ( , 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.

Discussion: evolution of the M H 2 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 Mpc 3 of the Universe, ρ(M H 2 ), at different epochs (Figure 10).Indeed, this plot shows that the derived molecular gas content of bright enough galaxies to be detected at  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 ρ(M H 2 ) 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 AGNs 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 ρ(M H 2 ) 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).
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.4GHz luminosity function at 9×10 25 W/Hz, 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 Figure 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.

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 radioto-infrared ratio q 24µ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 × 10 7 < M H 2 < 7 × 10 10 M at low-z and 2 × 10 10 < M H 2 < 7 × 10 11 M 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 10 13 -10 14 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 H 2 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.Notes.The values for the molecular mass listed here are corrected by our adopted CO-to-H 2 conversion factor.References: Prandoni et al. (2010), Ruffa et al. (2019a), Olivares et al. (2019), Ocaña Flaquer et al. (2010), Leon et al. (2003), Taniguchi et al. (1990), Evans et al. (2005), Saripalli & Mack (2007), Nesvadba et al. (2010), Fotopoulou et al. (2019), Lo et al. (1999), Evans et al. (1996) and Emonts et al. (2014).

Fig. 1 .
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.

Fig. 2 .
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.

Fig. 3 .
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.

Fig. 4 .
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.

Fig. 5 .
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.
and in the Appendix in Figure A.1.Based on the integrated intensity maps (left panels of Fig. 6 and Fig. 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 (Fig. 6 and Fig. A.1, middle-

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

Fig. 7 .
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 M vir = 10 13 and M vir = 10 14(Popping et al. 2015).
for z <0.3 and in the right panel of the same figure for 1< z <2.5.

Fig. 8 .
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 ofKeres 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 byLagos 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.

Fig. 9 .
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.4GHz luminosity function from the radio AGN survey in the VLA-COSMOS field (Smolčić et al. 2017, open green squares) for each redshift bin.

Fig. 10 .
Fig. 10.Evolution of the cosmic molecular gas density, ρ(M H 2 ), 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 ρ(M H 2 ) 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 (L 1.4GHz > 3 × 10 25 W/Hz) 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.

Table C .
1. CO properties of 54 literature-based sources.

Table C .
2. CO properties of representative sample of the 6C catalogue.Notes.References for M H 2 are listed in the right column: Olivares et al. (2019), Ocaña Flaquer et al. (2010),