The VLA-COSMOS 3~GHz Large Project: AGN and host-galaxy properties out to z$\lesssim$6

TL;DR

This study analyzes the properties of radio-selected AGN and their host galaxies up to redshift 6, revealing distinct evolutionary behaviors and host characteristics for different AGN classes over cosmic time.

Contribution

It provides the largest high-redshift radio source sample with detailed multiwavelength analysis, distinguishing between radiative and mechanical AGN modes and their host galaxy properties.

Findings

HLAGN are mainly radiative, MLAGN show significant mechanical luminosity

Host galaxies of MLAGN are more massive and less star-forming at z<1.5

Reversal in stellar mass distribution of AGN hosts at z>1.5

Abstract

We explore the multiwavelength properties of AGN host galaxies for different classes of radio-selected AGN out to z$\lesssim$6 via a multiwavelength analysis of about 7700 radio sources in the COSMOS field. The sources were selected with the Very Large Array (VLA) at 3 GHz (10 cm) within the VLA-COSMOS 3 GHz Large Project, and cross-matched with multiwavelength ancillary data. This is the largest sample of high-redshift (z$\lesssim$6) radio sources with exquisite photometric coverage and redshift measurements available. We constructed a sample of moderate-to-high radiative luminosity AGN (HLAGN) via spectral energy distribution (SED) decomposition combined with standard X-ray and mid-infrared diagnostics. Within the remainder of the sample we further identified low-to-moderate radiative luminosity AGN (MLAGN) via excess in radio emission relative to the star formation rates in their host galaxies. We show that AGN power in HLAGN occurs predominantly in radiative form, while MLAGN display a substantial mechanical AGN luminosity component. We found significant differences in the host properties of the two AGN classes, as a function of redshift. At z$<$1.5, MLAGN appear to reside in significantly more massive and less star-forming galaxies compared to HLAGN. At z$>$1.5, we observed a reversal in the behaviour of the stellar mass distributions with the HLAGN populating the higher stellar mass tail. We interpret this finding as a possible hint of the downsizing of galaxies hosting HLAGN, with the most massive galaxies triggering AGN activity earlier than less massive galaxies, and then fading to MLAGN at lower redshifts. Our conclusion is that HLAGN and MLAGN samples trace two distinct galaxy and AGN populations in a wide range of redshifts, possibly resembling the radio AGN types often referred to as radiative- and jet-mode (or high- and low-excitation), respectively.

Tables

Table 1. Table 1: Comparison between our classification method and that presented in Bonzini et al. ( 2013 ) . For this check, we cut our sample at total 1.4 GHz flux S 1.4 > subscript 𝑆 1.4 absent S_{\rm 1.4}> 37 μ 𝜇 \leavevmode\nobreak\ \mu Jy to match the radio selection adopted by the authors.

Classifications			(B13)
(this work,		RQ AGN	RL AGN	SFGs
$S_{1.4}>$ 37$\mu$Jy)	total	(609)	(865)	(3099)
HLAGN	(1044)	609	169	266
MLAGN	(1032)	0	569	463
Rest of the sample	(2497)	0	127	2370

Table 2. Table 2: Number of 3 GHz radio sources studied in this work as a function of redshift and AGN class: MLAGN, HLAGN, and HLAGN with radio excess (in brackets). For each redshift bin we report the mean redshift ⟨ z ⟩ delimited-⟨⟩ 𝑧 \langle z\rangle of the corresponding population.

redshift bin	$⟨z⟩$	MLAGN	HLAGN
			(radio-excess)
0.01 0.30	0.21	22	36 (9)
0.30 0.70	0.51	221	232 (66)
0.70 1.20	0.94	375	416 (135)
1.20 1.80	1.48	350	350 (98)
1.80 2.50	2.08	225	307 (84)
2.50 3.50	2.89	111	217 (74)
3.50 5.70	4.21	29	46 (15)

Table 3. Table 3: Results from the two-sample Kolmogorov-Smirnov (K-S) test. The table shows, in different redshift bins, the probability P(K-S) that the distributions of a given parameter ( M ⋆ subscript 𝑀 ⋆ M_{\star} , SFR or [NUV- r 𝑟 r ] ) for HLAGN and MLAGN (shown in Figs. 8 , 9 and 10 ) are drawn from the same parent distribution. A lower probability indicates a more significant difference. Values in brackets report the significance of the K-S test in units of σ 𝜎 \sigma .

		P(K-S)
Redshift bin	$M_{\star }$	SFR	[NUV-$r$]
	% ($\sigma$)	% ($\sigma$)	% ($\sigma$)
0.01 0.30	0.76 (2.67)	2.97 (2.17)	0.11 (3.25)
0.30 0.70	2.07$\times {10}^{-3}$ (4.26)	($>$10)	($>$10)
0.70 1.20	3.84$\times {10}^{-4}$ (4.62)	($>$10)	($>$10)
1.20 1.80	34.3 (0.95)	($>$10)	3.50$\times {10}^{-14}$ (8.16)
1.80 2.50	4.51$\times {10}^{-7}$ (5.86)	($>$10)	1.06$\times {10}^{-4}$ (4.88)
2.50 3.50	1.17 (2.52)	1.52$\times {10}^{-13}$ (7.97)	1.99$\times {10}^{-8}$ (6.36)
3.50 5.70	13.5 (1.49)	4.44$\times {10}^{-3}$ (4.08)	10.7 (1.61)

Table 4. Table 4: Median stacked fluxes derived in each redshift bin through a bootstrapping procedure. Lower and upper errors (in mJy) correspond to the 16 th and 84 th percentiles of the cumulative flux distribution. In case the error is larger than the median flux, we report the 1 σ 𝜎 \sigma upper flux. The number of stacked sources N stack subscript 𝑁 stack N_{\rm stack} is reported for each redshift bin.

redshift bin	$N_{\mathrm{stack}}$	PACS 100	PACS 160	SPIRE 250	SPIRE 350	SPIRE 500
		[mJy]	[mJy]	[mJy]	[mJy]	[mJy]
0.01 0.30	37	0.42	1.90${}_{-0.69}{}^{}{}_{}{}^{+0.44}$	2.20${}_{-0.66}{}^{}{}_{}{}^{+1.92}$	2.24${}_{-0.47}{}^{}{}_{}{}^{+1.09}$	3.02
0.30 0.70	366	1.09${}_{-0.11}{}^{}{}_{}{}^{+0.15}$	2.22${}_{-0.27}{}^{}{}_{}{}^{+0.30}$	1.64${}_{-0.37}{}^{}{}_{}{}^{+0.20}$	0.89${}_{-0.62}{}^{}{}_{}{}^{+0.68}$	1.09
0.70 1.20	642	1.02${}_{-0.09}{}^{}{}_{}{}^{+0.08}$	1.91${}_{-0.14}{}^{}{}_{}{}^{+0.20}$	1.98${}_{-0.27}{}^{}{}_{}{}^{+0.28}$	1.47${}_{-0.37}{}^{}{}_{}{}^{+0.44}$	0.70${}_{-0.43}{}^{}{}_{}{}^{+0.30}$
1.20 1.80	620	1.44${}_{-0.08}{}^{}{}_{}{}^{+0.10}$	3.03${}_{-0.24}{}^{}{}_{}{}^{+0.20}$	4.21${}_{-0.49}{}^{}{}_{}{}^{+0.47}$	4.25${}_{-0.50}{}^{}{}_{}{}^{+0.52}$	2.75${}_{-0.37}{}^{}{}_{}{}^{+0.20}$
1.80 2.50	306	0.95${}_{-0.11}{}^{}{}_{}{}^{+0.13}$	2.49${}_{-0.30}{}^{}{}_{}{}^{+0.24}$	3.33${}_{-0.38}{}^{}{}_{}{}^{+0.25}$	3.77${}_{-0.77}{}^{}{}_{}{}^{+0.45}$	3.10${}_{-0.58}{}^{}{}_{}{}^{+0.50}$
2.50 3.50	175	1.00${}_{-0.13}{}^{}{}_{}{}^{+0.10}$	2.25${}_{-0.30}{}^{}{}_{}{}^{+0.34}$	5.44${}_{-1.36}{}^{}{}_{}{}^{+0.99}$	5.73${}_{-0.54}{}^{}{}_{}{}^{+0.88}$	4.55${}_{-1.13}{}^{}{}_{}{}^{+1.13}$
3.50 5.70	57	0.79${}_{-0.37}{}^{}{}_{}{}^{+0.22}$	2.15${}_{-0.40}{}^{}{}_{}{}^{+0.88}$	5.19${}_{-0.98}{}^{}{}_{}{}^{+1.19}$	7.68${}_{-1.35}{}^{}{}_{}{}^{+0.94}$	6.80${}_{-1.35}{}^{}{}_{}{}^{+0.60}$

Table 5. Table 5: Table listing the properties and classification for the first 20 radio sources (sorted by ID 3 GHz) used in this work. The full Table will be made available through the IPAC/IRSA database.

(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14-17)	(18)
ID	RA 3 GHz	Dec 3 GHz	Redshift	Redshift	$S_{3\mathrm{GHz}}$	$L_{3\mathrm{GHz}}$	$L_{1.4\mathrm{GHz}}$	$L_{\mathrm{IR},\mathrm{SF}}$	$\ge$3$\sigma$ Herschel	$M_{\star }$	SFRIR	[NUV-$r$]	AGN	Class
3 GHz	[deg]	[deg]		type	[$\mu$Jy]	log[W Hz-1]	log[W Hz-1]	log[$L_{\odot }$]	detection	log[$M_{\odot }$]	$M_{\odot }$ yr-1	AB mag	criteria
1	149.64771	2.09546	1.546	phot	16147.0	25.84	25.72	11.66	true	11.22	46.0	2.42	0 0 01	MLAGN
3	150.33360	2.57880	1.555	phot	11369.2	25.80	25.77	12.12	true	10.21	131.4	-0.29	1 1 11	HLAGN
5	150.72035	1.93047	2.446	phot	7746.0	26.57	26.90	12.01	false	11.45	102.5	1.84	0 0 01	MLAGN
6	150.47405	2.83167	1.259	spec	9212.9	25.90	26.20	11.26	false	11.23	18.3	3.14	1 0 01	HLAGN
8	150.00256	2.25863	2.450	spec	6749.5	25.79	25.68	11.48	true	11.08	30.3	2.07	1 1 11	HLAGN
16	149.51338	2.23267	2.238	phot	13523.3	26.62	26.88	11.42	false	11.03	26.2	1.34	1 0 01	HLAGN
17	150.72133	1.58238	1.051	phot	6278.1	25.04	24.80	11.41	false	11.22	25.8	2.49	1 0 11	HLAGN
18	149.95944	1.80146	0.684	spec	3235.1	24.35	23.99	10.49	false	10.92	3.1	2.85	1 0 01	HLAGN
19	149.42624	2.07387	1.081	phot	8626.8	25.75	26.10	11.91	true	10.76	81.3	0.63	0 0 11	HLAGN
23	150.56024	2.58613	2.099	phot	3151.0	25.85	26.07	11.40	false	11.24	25.3	2.26	0 0 01	MLAGN
25	150.06906	2.44399	2.436	phot	3178.2	26.05	26.30	12.70	true	10.61	501.8	-0.18	0 1 11	HLAGN
29	150.44730	2.05394	0.323	spec	4133.0	24.12	24.36	9.98	false	11.45	0.9	4.90	1 0 0 1	HLAGN
30	149.96638	2.09516	1.356	phot	2277.0	24.84	24.67	12.27	true	11.52	187.5	1.17	0 0 01	MLAGN
31	150.61989	2.28940	2.625	spec	2372.8	25.88	26.07	11.82	false	11.23	65.5	1.17	1 0 01	HLAGN
33	150.37977	2.49021	0.349	spec	2969.6	24.06	24.31	10.27	false	11.30	1.9	4.42	1 0 01	HLAGN
37	149.78699	1.60178	1.680	phot	2405.4	25.75	26.15	11.64	false	10.74	43.8	1.22	0 0 01	MLAGN
39	149.61919	1.91632	0.913	spec	1868.1	24.87	25.17	10.87	false	10.67	7.4	2.42	0 0 01	MLAGN
40	149.57238	2.26267	0.706	spec	1825.5	24.52	24.71	11.49	true	11.31	30.9	1.97	1 0 11	HLAGN
41	149.55937	1.63104	0.901	spec	2102.8	24.71	24.77	11.10	false	11.32	12.6	2.32	1 0 01	HLAGN
43	150.67227	1.71565	0.609	phot	1510.3	23.92	23.54	11.35	true	11.53	22.4	2.33	0 0 01	MLAGN

Equations

$$\log {\left(\frac{L_{1.4\mathrm{GHz}}}{\mathrm{SF}{\mathrm{R}}_{\mathrm{IR}}}\right)}_{\mathrm{excess}}=21.984\times (1+z{)}^{0.013}.$$

Full text

1 Introduction

Multiwavelength observations of galaxies hosting active galactic nuclei (AGN) in the past 20 years have extensively shown that supermassive black holes (SMBHs) are one of the key ingredients in shaping the evolution of galaxies through cosmic time. In particular, it is now well-established that AGN activity and star formation in their hosts are related processes, which are likely driven by a common fuelling mechanism such as accretion of cold gas supplies (e.g. Vito et al. 2014). In the local Universe, hints of such a connection have been suggested by the empirical correlations found between black hole mass and galaxy properties (e.g. Magorrian et al. 1998; Gebhardt et al. 2000; Ferrarese 2002; Gültekin et al. 2009). In the distant Universe, this co-evolution scenario is supported by the similarity between volume-averaged cosmic star formation history and black hole accretion history, which both peak at z$\sim$2 and decline towards the local Universe (e.g. Madau & Dickinson 2014 for a review). Semi- analytic models (e.g. Bower et al. 2006) and numerical simulations interpret these correlations as originated from a long-lasting ($\gtrsim$ Gyr) self-regulation process between the SMBH and its host (e.g. Bower et al. 2006; Croton et al. 2006; Hopkins et al. 2008; Lagos et al. 2008; Menci et al. 2008), which occurs in two flavours: quasar mode (QSO mode) and radio mode.

On the one hand, the quasar mode is usually associated with a radiatively efficient phase of SMBH accretion through isotropically distributed ionising winds and molecular outflows as means to prevent the runaway growth of SMBHs (e.g. Booth & Schaye 2009; Dubois et al. 2014). On the other hand, a subsequent radio-mode phase is usually invoked to prevent further episodes of galaxy star formation through mechanical feedback, such as collimated and relativistic jets (Monaco et al. 2000; Dubois et al. 2014), in order to regulate the galaxy stellar mass (Croton et al. 2006; Marulli et al. 2008; Hopkins & Quataert 2010) and to reproduce the galaxy colours observed in the local Universe (Strateva et al. 2001).

Though it is today widely accepted that the evolution of active SMBHs is connected to the evolution of their hosts, the underlying mechanisms explaining the transition between these two key stages of the AGN and galaxy life cycles are still poorly constrained. Testing this paradigm is challenging as this is supposed to be a quick transition (on timescales of $\sim$100 Myr) from highly accreting to fading AGN. Promising studies of individual smoking-guns, mostly X-ray and optically detected AGN (e.g. Farrah et al. 2012; Cicone et al. 2014; Perna et al. 2015; Brusa et al. 2015; Brusa et al. 2016), have shown compelling evidence of ongoing AGN feedback, but are currently limited to a small number of candidates. A combined study with large AGN samples is necessary to constrain the role of AGN feedback in a more statistical sense, and to shed light on the connection between AGN and their hosts at different cosmic epochs.

A fundamental, complementary perspective in the framework of the AGN-host evolution comes from radio observations. Indeed, radio observations are essential to capture possible signatures of relativistic jets powered by a central SMBH (e.g. Hogan et al. 2015), which are detected via the synchrotron emission of the jet (e.g. Miller et al. 1993). In addition, radio continuum emission may arise from the diffusion of cosmic ray electrons produced in supernovae and their remnants in high-mass star-forming regions, and this emission has been calibrated on star-forming galaxies to provide an almost dust-unbiased star formation rate (SFR) indicator (Condon 1992; Yun et al. 2001; Bell 2003). This underlines the great potential of radio observations in unveiling a mixture of AGN and star-forming galaxies. Nevertheless, radio observations need to be supplemented by multiwavelength data to fully characterise the nature of the radio sources.

Outstanding progress has been made through the analysis of large samples of radio-selected AGN in the local Universe (Smolčić 2009; Best & Heckman 2012). For instance, Smolčić (2009) identified a twofold population of radio-emitting AGN, namely high-excitation and low-excitation radio galaxies (HERGs and LERGs, respectively), on the basis of the presence of high- or low-excitation lines in their optical spectra taken from the Sloan Digital Sky Survey (SDSS; York et al. 2000). Interestingly, the author found that HERGs preferentially live in galaxies within the green valley (in terms of optical [NUV-$r$] colours and stellar mass, $M_{\star}$), while LERGs usually populate the red sequence of massive and passive systems. Such a dichotomy observed in the host-galaxy properties between HERGs and LERGs may reflect physically different modes of SMBH accretion and presumably different stages of AGN-galaxy evolution (Hardcastle et al. 2006; see Heckman & Best 2014 for a review). While SMBHs in HERGs are thought to accrete via cold gas inflows from galaxy mergers or secular processes (Sijacki et al. 2007; Di Matteo et al. 2008; Bournaud et al. 2012), AGN activity in LERGs is probably induced by a continuous gas inflow coming from the atmosphere of the hot halo (Bower et al. 2006; Ellison et al. 2015). In the former case, accretion is radiatively efficient and covers a wide range of the electromagnetic spectrum (up to X-ray frequencies), while for the latter scenario the feedback is predominantly mechanical and does not outshine the host galaxy in most bands except radio.

In their comprehensive study, Hickox et al. (2009) have thoroughly investigated the AGN, host galaxy, and environmental properties of X-ray, mid-IR (MIR), and radio-selected AGN at 0.25$<$z$<$0.8 in the Böotes field (Jannuzi & Dey 1999). In particular, Hickox et al. (2009) defined as “radio AGN” those sources with (rest-frame) 1.4 GHz luminosity $L_{\rm 1.4\leavevmode\nobreak\ GHz}>$ 1024.8 W Hz*-1* to minimise the contamination from star-forming galaxies. These authors found that most radio-selected AGN have very low accretion rates (Eddington ratio $\lambda_{\rm Edd}\lesssim$ 10*-3*) and populate overdense regions similarly to the most massive galaxies (e.g. Georgakakis et al. 2007; Silverman et al. 2009; Coil et al. 2009). In contrast, X-ray and MIR selected AGN are characterised by active star formation and less dense environments. These results, which have been corroborated through a similar analysis up to intermediate redshifts ($z\lesssim 1.4$; see Goulding et al. 2014), strongly suggest that various AGN selection criteria might be sensitive to physically distinct classes of AGN and galaxies. In particular, the peculiarity of black hole and galaxy properties observed in radio-selected AGN stands out more than in any other AGN sample. For these reasons, it is now widely recognised that a multiwavelength investigation of radio-selected sources is essential to constrain the AGN-galaxy properties in a key stage of their cosmic evolution.

In this paper, we exploit the largest compilation of high-redshift (z$\lesssim$6) radio-selected galaxies in the Cosmic Evolution Survey (COSMOS; Scoville et al. 2007) field. The Karl G. Jansky Very Large Array (VLA) observations were conducted at 3 GHz (10 cm) over the entire COSMOS field, in the framework of the VLA-COSMOS 3 GHz Large Project (PI: V. Smolčić, Smolcic et al. 2017b), reaching a 1$\sigma$ sensitivity of 2.3$\leavevmode\nobreak\ \mu$Jy beam*-1*. The rich multiwavelength (X-ray to radio) data set of photometry and redshifts available in the COSMOS field allows us to investigate the physical properties of these sources from a panchromatic perspective. The main goals of the present work are twofold: first, to provide a value-added catalogue that includes classification and physical properties for each 3 GHz VLA-selected source in the COSMOS field, and, second, to explore the multiwavelength properties of AGN hosts for different classes of radio-selected AGN out to z$\lesssim$6.

The paper is structured as follows. In Sect. 2 we describe our sample selection and the cross-match with ancillary photometry. In Sect. 3 we decompose the multiwavelength spectral energy distribution (SEDs), while the classification of our sample is discussed in Sect. 4. A brief description of the value-added catalogue is given in Sect. 5. Sect. 6 illustrates the average radio-selected AGN host-galaxy properties out to z$\lesssim$6, while the interpretation of our results are presented and discussed in detail in Sect. 7. We list our concluding remarks in Sect. 8. In Appendix A we show the results of infrared stacking, while Appendix B shows a portion of the value-added 3 GHz radio catalogue including some physical parameters used in this work. Throughout this paper, magnitudes are given in the AB system (Oke 1974). We assume a Chabrier (2003) initial mass function (IMF) and a flat cosmology with $\Omega_{\rm m}$ = 0.30, $\Omega_{\rm\Lambda}$ = 0.70, and H0 = 70 km s*-1* Mpc*-1* (Spergel et al. 2003).

2 Sample selection

2.1 3 GHz radio sources

Radio data at 3 GHz were collected from 384 hours of observations with VLA over 2.6 deg2, reaching an average rms sensitivity of 2.3 $\mu$Jy beam*-1* and an angular resolution of about $0\aas@@fstack{\prime\prime}75$. A detailed description of the survey strategy, data reduction, and radio source catalogue is given in Smolcic et al. (2017b). The catalogue includes 10 830 radio sources, identified at peak surface brightness $\geq$5$\sigma$, out of which 67 are multi-components. The present catalogue represents the deepest compilations of radio sources available to date across an area of 2.6 deg2. Our sample covers a wide redshift range (0$<$z$\lesssim$6, see Sect. 2.3) and is around three times larger than the 1.4 GHz sample taken from the Westerbork Synthesis Radio Telescope (WSRT, 3 172 sources) in the NOAO Deep Wide-Field Survey (NDWFS, de Vries et al. 2002). Moreover, our sample outnumbers the previous 1.4 GHz VLA-COSMOS survey by a factor of about four (2 865 sources; see Schinnerer et al. 2007, 2010) and by more than one order of magnitude the 1.4 GHz VLA survey in the Extended Chandra-Deep Field South (E-CDFS; 883 sources, Miller et al. 2013) survey.

We derived (rest-frame) 3 GHz radio luminosity ($L_{\rm 3\,GHz}$) for radio sources with multiwavelength counterparts and redshifts (see Sects. 2.2 and 2.3). Under the assumption of purely synchrotron emission, the radio spectrum behaves like a power law $S_{\nu}\propto\nu^{\alpha}$, where the spectral index $\alpha$ is set to the observed 1.4–3 GHz slope for sources detected also at 1.4 GHz (about 30%) in the 1.4 GHz VLA-COSMOS survey (Schinnerer et al. 2010); the spectral index is set to –0.7, which is consistent with a non-thermal synchrotron index, (e.g. Condon 1992; see also Smolcic et al. 2017b) if the sources are detected at 3 GHz alone. In Fig. 1 (bottom panel) we show $L_{\rm 3\,GHz}$ as a function of redshift with respect to the 5$\sigma$ luminosity limit. For comparison, we show the corresponding 1.4 GHz luminosity $L_{\rm 1.4\,GHz}$ on the right $y$ axis. Our sample clearly spans a wide luminosity range (up to 4–5 dex), which allows us to investigate the multiwavelength properties of our sample in different radio luminosity regimes.

In particular, we are able to detect 2.5 times intrinsically fainter sources (under the assumption $\alpha$=–0.7), at a given redshift, compared to the previous 1.4 GHz VLA-COSMOS survey (Schinnerer et al. 2007, 2010).

2.2 Optical to (sub)millimetre photometry

The COSMOS field benefits from an exquisite photometric data set, covered from the X-rays to the submillimetre domain111An exhaustive overview of the COSMOS field and multiwavelength data products is available at: http://cosmos.astro.caltech.edu/ . Cross-matching our 3 GHz selected sample to existing ancillary data is essential to derive physical properties of galaxies. The multiwavelength photometry is taken from the COSMOS2015 catalogue (Laigle et al. 2016), which combines optical photometry222Optical photometry is taken from Subaru Hyper-Suprime Cam observations over the full 2 deg2 (Capak et al. 2007), and also from the Canada-France-Hawaii Telescope Legacy Survey (CFHT-LS; McCracken et al. 2001) in the central 1 deg2., the most recent UltraVISTA (DR2333DR2 replaces the previous DR1 by McCracken et al. (2012). A detailed description of the survey and data products can be retrieved at: http://ultravista.org/release2 ) data over the central 1.5 deg2 in the near-infrared (NIR) bands $Y$, $J$, $H$, and $K_{s}$444Outside the UltraVISTA coverage, NIR photometry includes CFHT $H$ and $K_{s}$ observations obtained with the WIRCam (McCracken et al. 2001)., and MIR photometry obtained from the Infrared Array Camera (IRAC), which is recently complemented by deeper IRAC 3.6 and 4.5$\leavevmode\nobreak\ \mu$m observations with the Spitzer Large Area Survey with Hyper-Suprime-Cam (SPLASH; Steinhardt et al. 2014; P. Capak et al. in prep.). In addition, this data set has been cross-matched with 24$\leavevmode\nobreak\ \mu$m photometry (Le Floc’h et al. 2009) from the Multi-Band Imaging Photometer for Spitzer (MIPS). Laigle et al. (2016) provides further details.

The cross-match to associate a possible optical-NIR counterpart with each radio source is fully described in Smolcic et al. (2017a) (see their Sect. 3). First, they excluded stars and masked regions in the COSMOS2015 catalogue because of the less accurate optical photometry, which reduces the effective area of the COSMOS field to 1.77 deg2 and our 3 GHz selected sample to 8 696 radio sources. Secondly, they performed a nearest-neighbour matching, by selecting for each radio source only candidate counterparts within $0\aas@@fstack{\prime\prime}8$ searching radius and, at the same time, requiring a false match probability (i.e. probability of being a spurious association) lower than 20%. This approach yields an average expected fraction of spurious associations of about 1% (see Smolcic et al. 2017a). The percentage of radio sources with multiple optical-NIR counterparts within $0\aas@@fstack{\prime\prime}8$ is around 1%, for which the cut in false-match probability ensures the selection of the most probable counterpart. After this cut, our final sample consists of 7 729 radio sources with optical-NIR counterparts, corresponding to about 89% of our radio-selected sample within the common 1.77 deg2.

To enrich the spectral coverage of our analysis and derive robust star formation rates (SFRs) for as many sources as possible, we used also Herschel photometry at far-infrared and submillimetre wavelengths provided in the COSMOS2015 catalogue. Herschel imaging covers the entire COSMOS field with the Photoconductor Array Camera and Spectrometer (PACS; 100 and 160$\leavevmode\nobreak\ \mu$m, Poglitsch et al. 2010) and Spectral and Photometric Imaging Receiver (SPIRE; 250, 350, and 500$\leavevmode\nobreak\ \mu$m, Griffin et al. 2010) data, as part of the PACS Evolutionary Probe (PEP; Lutz et al. 2011) and the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012). Herschel fluxes were extracted and de-blended by using 24$\leavevmode\nobreak\ \mu$m positional priors and unambiguously associated with the corresponding optical-NIR counterpart via 24$\leavevmode\nobreak\ \mu$m sources listed in both catalogues. In total, the number of radio sources with ($\geq$3$\sigma$) Herschel detection in at least one PACS or SPIRE band are 4 836/7 729 (63%). This percentage decreases with redshift, being 87% at z$<$0.3 and 45% at z$>$3.5.

To obtain reliable dust-based SFRs also in potential high-redshift candidates (z$>$3), where Herschel observations are incomplete even towards ultraluminous infrared galaxies (ULIRGs, i.e. having rest-frame 8-1000$\leavevmode\nobreak\ \mu$m infrared luminosity $L_{\rm IR}\geq$ 1012 L⊙, e.g. Sanders & Mirabel 1996), photometry at longer wavelengths is essential. For around 115 radio sources, we retrieved additional (sub)millimetre photometry from at least one of the following data sets: JCMT/SCUBA-2 at 450 and 850$\leavevmode\nobreak\ \mu$m (Casey et al. 2013), LABOCA at 870$\leavevmode\nobreak\ \mu$m (F. Navarrete et al. priv. comm.), Bolocam (PI: J. Aguirre), JCMT/AzTEC (Scott et al. 2008) and ASTE/AzTEC (Aretxaga et al. 2011) at 1.1 mm, MAMBO at 1.2 mm (Bertoldi et al. 2007), and interferometric observations at 1.3 mm with ALMA (PI: M. Aravena, M. Aravena et al. in prep.) and PdBI (Smolčić et al. 2012; Miettinen et al. 2015). The (sub)mm positions were cross-matched to the COSMOS2015 positions via a nearest neighbour matching, using $1\aas@@fstack{\prime\prime}$ searching radius (the smallest beam width of the (sub)mm data we collected). A thorough visual inspection of the counterpart associations has been performed for the 1.3 mm detected ALMA sources (68% of the (sub)mm photometry we collected), which is detailed in Brisbin et al. (2017) and Miettinen et al. (2017).

We also collected X-ray data from the Chandra-COSMOS (Elvis et al. 2009; Civano et al. 2012) and COSMOS-Legacy catalogues (Civano et al. 2016). The optical-NIR counterparts of X-ray sources were matched via a maximum likelihood algorithm and are presented in Marchesi et al. (2016). We matched their catalogue to our 3 GHz selected sample of 7 729 optical-NIR counterparts via COSMOS2015 IDs. This match yields 903 X-ray sources, corresponding to 12% (903/7 729) of our radio sample, and to 32% (903/2 804) of the X-ray sources with optical-NIR association in unmasked areas.

2.3 Spectroscopic and photometric redshifts

We collected photometric redshifts for the 7 729 radio sources with a counterpart in the COSMOS2015 catalogue. Photometric redshift estimates are included in the catalogue and were derived using the Le Phare SED-fitting code (Arnouts et al. 1999; Ilbert et al. 2006) following the procedure detailed in Ilbert et al. (2009, 2013). Based on the comparison with the spectroscopic redshifts available in the COSMOS field, Laigle et al. (2016) report an average photometric redshift accuracy of $\left\langle|\Delta z/(1+z)|\right\rangle=$ 0.021 for $K_{s}>$22, which becomes less than 0.010 for brighter sources.

For X-ray sources, we used a different set of photometric redshifts from M. Salvato et al. (in prep.), which are more suitable for AGN-dominated sources as they account for AGN variability and adopt additional AGN templates (Salvato et al. 2009, 2011).

An exhaustive list of spectroscopic redshifts was compiled (April 2015, Salvato et al. in prep.) and made internally accessible to the COSMOS team. Most of the spectroscopic redshifts used in this paper were taken from the $z$COSMOS survey (Lilly et al. 2007, 2009), either the public $z$COSMOS-bright or the proprietary $z$COSMOS-deep database, the DEep Imaging Multi-Object Spectrograph (DEIMOS, Capak et al. in prep.), and the FOcal Reducer and low dispersion Spectrograph (FORS2, Comparat et al. 2015). M. Salvato et al. (in prep.) provide for a full reference list.

For each radio source with a multiwavelength counterpart in this spectroscopic compilation, we replaced the photometric redshifts with new spectroscopic values only in case of secure or very secure measurements555The reliability of each spectroscopic redshift is determined by its quality flag. In case of spectroscopic redshift from the zCOSMOS survey (Lilly et al. 2007, 2009), we followed the prescription recommended on the zCOSMOS IRSA webpage: https://irsa.ipac.caltech.edu/data/COSMOS/spectra/z-cosmos/Z-COSMOS_INFO.html For the other surveys we selected quality flag $Qf\geq$ 3 and discarded less reliable spectroscopic redshifts from our analysis.. In addition, we included the latest spectroscopic redshifts from the VIMOS Ultra Deep Survey (VUDS; Le Fèvre et al. 2015; Tasca et al. 2016), from which we found 25 associations to our radio sources.

After these checks, the number of radio sources with spectroscopic redshift is 2 734/7 729 (around 35%). Every radio source with multiwavelength counterpart has its own redshift estimate. Fig. 1 (top panel) shows the final redshift distribution for our 7 729 radio sources. The number of spectroscopic and photometric redshifts are comparable out to z$\sim$1, while photometric redshifts become more numerous at higher redshift. We tested the accuracy of the photometric redshifts in our sample based on the spectroscopic measurements available for 2 734 sources. We found a median $\left\langle|\Delta z/(1+z)|\right\rangle=$ 0.010, which becomes as high as 0.035 at z$>$3. Therefore, the proved accuracy of the photometric redshifts allows us to push our analysis out to z$\lesssim$6, even if the number of sources at z$>$4 is relatively small (84 sources).

3 SED-fitting decomposition of 3 GHz sources

In this section, we fit the multiwavelength SEDs of our radio sources to disentangle the AGN emission from that related to the host-galaxy. It is well known that radio-selected samples contain distinct galaxy populations (e.g. Condon 1984; Windhorst et al. 1985; Gruppioni et al. 1999) in terms of star formation and AGN properties. Therefore, fitting the multiwavelength SEDs may provide meaningful results only if AGN and galaxy templates are both taken into account.

We used both the SED-fitting code magphys666The original magphys code is publicly available at this link: http://www.iap.fr/magphys/magphys/MAGPHYS.html (da Cunha et al. 2008), and the three-component SED-fitting code sed3fit by Berta et al. (2013), which accounts for an additional AGN component777The three-component SED-fitting code sed3fit can be retrieved from http://cosmos.astro.caltech.edu/page/other-tools. The aforementioned references provide for a full description of these SED-fitting codes. Here we briefly outline the main prescriptions that are relevant for our analysis.

The magphys code is designed to reproduce a variety of galaxy SEDs, from weakly star-forming to starbursting galaxies, over a wide redshift range888For extensive application of the magphys code in deriving physical properties of galaxies, see also Smith et al. (2012); Rowlands et al. (2014); Michałowski et al. (2014); Hayward & Smith (2015).. This code relies on the energy balance between the dust-absorbed stellar continuum and the reprocessed dust emission at infrared wavelengths. This recipe ensures that optical and infrared emission originating from star formation are linked in a self-consistent manner, but does not account for a possible AGN emission component. The three-component SED-fitting code presented by Berta et al. (2013) combines the emission from stars, dust heated by star formation, and a possible AGN-torus component from the library of Feltre et al. (2012, see also Fritz et al. 2006). This approach results in an effectively simultaneous three-component fit. For each best-fit parameter, the code provides a corresponding probability distribution function (PDF), which enables the user to obtain reliable confidence ranges for parameter estimates (see e.g. Calistro Rivera et al., in prep., for a similar SED-fitting technique).

We decomposed each observed SED by using the best available redshift (either spectroscopic or photometric, see Sect. 2.3) as input, and we derived integrated galaxy properties, such as SFR and $M_{\star}$, for each individual source. The SFR was derived from the total IR (rest 8-1000$\leavevmode\nobreak\ \mu$m) luminosity taken from the best-fit galaxy SED (i.e. corrected for a possible AGN emission), assuming a Kennicutt (1998) conversion factor scaled to a Chabrier (2003) IMF. We note that about 37% of our sample are not $\geq$3$\sigma$ detected in any Herschel bands. To obtain better constrained IR luminosities, we performed SED-fitting using the nominal PACS and SPIRE 3$\sigma$ upper limits, which are equal to 5.0 (100$\leavevmode\nobreak\ \mu$m), 10.2 (160$\leavevmode\nobreak\ \mu$m), 8.1 (250$\leavevmode\nobreak\ \mu$m), 10.7 (350$\leavevmode\nobreak\ \mu$m), and 15.4 (500$\leavevmode\nobreak\ \mu$m) mJy, including confusion noise (Lutz et al. 2011; Oliver et al. 2012). We modified the $\chi^{2}$ calculation to correctly account for those Herschel bands that have only upper limits, similar to the approach adopted by da Cunha et al. (2015).

As a sanity check, we verified that the IR luminosities based on our three-component fit are in good agreement with those calculated independently using a different set of IR templates (from Chary & Elbaz 2001; Dale & Helou 2002; Siebenmorgen & Krügel 2007; Polletta et al. 2007; Wuyts et al. 2008; Elbaz et al. 2011; Nordon et al. 2012, see Berta et al. 2013 for a comprehensive discussion). We briefly discuss the comparison with the Herschel fluxes derived through stacking in Sect. 4.4.3. The $M_{\star}$ is derived from the SED decomposition itself, which allows us to obtain robust estimates if the optical-NIR SED is dominated by the host-galaxy light (e.g. Bongiorno et al. 2012).

In order to quantify the relative incidence of a possible AGN component, we fitted each individual SED, both with the three-component approach and the magphys code. The fit obtained with the AGN is preferred if the reduced $\chi^{2}$ value of the best fit is significantly (at $\geq 99$% confidence level, on the basis of a Fisher test) smaller than that obtained from the fit without the AGN; see Delvecchio et al. (2014) for details. From our analysis, we found that 1 169 out of 7 729 radio sources (about 15%) show a $\geq 99$% significant AGN component in their best fit.

We extensively tested this technique against independent AGN indicators in the COSMOS field, such as MIR colours and X-rays (see Delvecchio et al. 2014). For instance, Lanzuisi et al. (2015) showed that the AGN radiative luminosities derived from SED decomposition were consistent (1$\sigma$=0.4 dex) with those calculated from X-ray spectra and assuming a set of bolometric corrections (e.g. Lusso et al. 2012). Moreover, the unprecedented accuracy of photometric redshifts and the photometric coverage exploited in this work, strengthened by our sizeable sample, further increase the reliability of our method. However, if the galaxy light outshines the AGN in the full optical-to-mm SED, this statistical technique becomes progressively less effective in identifying AGN.

4 Classification of 3 GHz radio sources

We combine SED-fitting decomposition (Sect. 3) with other multiwavelength AGN diagnostics to reach a more complete census of AGN in our sample. These additional indicators are taken from X-ray, MIR and radio data, which allow us to identify two main populations of radio-selected AGN in our sample: moderate-to-high radiative luminosity AGN (HLAGN) and low-to-moderate radiative luminosity AGN (MLAGN). Hereafter, we refer to these populations as HLAGN (X-ray, MIR, and SED-selected AGN) and MLAGN (radio-excess sources that are not HLAGN), as explained in the next sections.

This naming convention comes from the idea that the selection criteria based on SED-fitting, X-ray, and MIR data preferentially select higher luminosity AGN, where the term “luminosity” here refers to the AGN radiative luminosity (Lrad,AGN), which is a proxy of the SMBH accretion rate (BHAR; e.g. Alexander & Hickox 2012). This classification does not translate into a sharp threshold in the accretion efficiency (or Eddington ratio) between HLAGN and MLAGN, but rather reflects the reliability of the adopted diagnostics in identifying such AGN populations combined with the sensitivity of our survey at various wavelengths. Moreover, the tags “low to moderate” and “moderate to high” intentionally imply a potential overlap in Lrad,AGN between the two classes at various redshifts. However, at a given redshift, these differently selected AGN display significantly distinct distributions of AGN luminosity, as detailed in Sect. 4.4.1. Therefore, the present classification should be considered as observationally based, and aimed at dissecting our radio sources on the basis of AGN diagnostics that are known to be luminosity dependent. A detailed investigation on the distribution of radio-selected AGN as a function of their intrinsic Eddington ratio will be presented in a forthcoming paper (Delvecchio et al., in prep.).

In Sects. 4.1 and 4.2 we describe in more detail the multiwavelength diagnostics used to identify MLAGN and HLAGN, respectively, while in Sect. 4.4 we justify this naming convention by studying their Lrad,AGN distributions.

4.1 Moderate-to-high radiative luminosity AGN

As previously mentioned, the so-called HERG population identified in the local Universe consists of highly accreting SMBHs on the basis of the presence of high-excitation lines in their optical spectra (e.g. Smolčić et al. 2008), which implies radiatively efficient accretion. In order to detect potential HERG analogues in our sample, we combine SED decomposition (Sect. 3) with X-ray and MIR indicators. All these selection criteria are sensitive to an excess of emission likely arising from accretion onto the central SMBH rather than from star formation. As a consequence, despite the different biases intrinsic to each criterion, all of these criteria preferentially select higher radiative luminosity AGN. The selection criteria used to identify this AGN population are briefly summarised below. Smolcic et al. (2017a) provide a detailed description of the X-ray and MIR-based AGN indicators.

First, SED-fitting decomposition identifies 1 169 sources with $\geq$99% significant AGN component in their global SED, SED AGN hereafter (see Sect. 3).

Second, we used X-ray luminosities ($L_{\rm x}$) in the rest-frame [0.5–8] keV. The $L_{\rm x}$ estimates were calculated for the 903 X-ray detected sources, by assuming a fixed X-ray spectral slope $\Gamma=1.8$, and correcting for nuclear obscuration on the basis of the measured hardness ratio (e.g. Xue et al. (2010)). We identified 855 sources with X-ray luminosity $L_{\rm x}\geq$1042 erg s*-1* as X-ray AGN (e.g. Szokoly et al. 2004). We verified that the X-ray emission expected from recent $L_{\rm x}$–SFR relations (taken from Symeonidis et al. 2014) is always negligible for our X-ray AGN (about a few percent). On the basis of the aforementioned relation, about 30 X-ray detected sources with $L_{\rm x}<$1042 erg s*-1* show an X-ray excess and, therefore, could be considered low-luminosity X-ray AGN. However, we prefer to apply the same cut at $L_{\rm x}\geq$1042 erg s*-1* for all our sources to avoid potential contamination from outliers with respect to the $L_{\rm x}$– SFR relation.

Third, MIR colours can be very useful in identifying AGN, both unobscured and heavily obscured. Donley et al. (2012) proposed a conservative criterion to select AGN, on the basis of the MIR colour-colour diagram drawn from a combination of the four Spitzer-IRAC (3.6, 4.5, 5.8, and 8.0$\leavevmode\nobreak\ \mu$m) bands. We followed Eqs. (1) and (2) of their paper to identify AGN at z$<$2.7, while at higher redshift we applied the additional conditions stated in their Eqs. (3) and (4) to minimise the contamination from high-redshift starbursts without AGN. This method is highly reliable for bright AGN, but becomes incomplete at $L_{\rm x}<$1044 erg s*-1*. In total, 455 out of 7 729 radio sources (about 6%) satisfy the Donley et al. (2012) criterion, and therefore are classified as MIR AGN.

Hereafter, we will use the term “moderate-to-high radiative luminosity AGN” (HLAGN) to collectively refer to the union of X-ray, MIR, and SED-selected AGN identified in our sample, for a total of 1 604 objects (21% of the radio sample).

Figure 2 shows the percentages of HLAGN classified from each criterion: the percentage of AGN that fulfills all the criteria simultaneously is only about 14% of the full HLAGN population. This small overlap further suggests that different AGN diagnostics are sensitive to distinct AGN populations. This overlap increases with increasing X-ray luminosity, which is 7% for $10^{42}<L_{\rm x}<$1043 erg s*-1*, 25% for $10^{43}<L_{\rm x}<$1044 erg s*-1*, and 49% for $L_{\rm x}>$1044 erg s*-1*. These relatively small percentages are mainly driven by the incompleteness of the MIR classification, as the Donley et al. (2012) criterion is very conservative. We checked that the agreement between X-ray and SED-fitting diagnostics is as high as 21% for $10^{42}<L_{\rm x}<$1043 erg s*-1*, 52% for $10^{43}<L_{\rm x}<$1044 erg s*-1*, and 79% for $L_{\rm x}>$1044 erg s*-1*.

In Fig. 3 we illustrate some examples of best-fit SEDs, showing different levels of agreement between the AGN diagnostics described above. In all the panels, red circles indicate the (rest-frame) multiwavelength photometry, while downward pointing arrows set the 3$\sigma$ upper limits in the Herschel bands. Solid lines represent the best-fit templates of AGN (red), galaxy star formation (blue), and the sum of the two (black).

The left panel shows the SED of an unambiguous AGN, successfully identified from X-rays, MIR-colours, and SED-fitting decomposition. The central panel shows an AGN identified only from SED decomposition. Indeed, galaxies hosting heavily obscured AGN might be undetected in the X-rays, but also misclassified from MIR colours since the Donley et al. (2012) criterion is highly incomplete at $L_{\rm x}<10^{44}$ erg s*-1*. However, neither SED-fitting decomposition nor MIR colours can identify an AGN when the optical-IR SED is outshined by the host-galaxy light (right panel), although the X-ray luminosity suggests the presence of a moderately luminous X-ray AGN ($L_{\rm x}\sim$1043 erg s*-1*). We looked at the observed distribution of the X-ray to optical-UV index, defined as $\alpha_{\rm ox}$ = $-$Log[$L_{\rm 2\leavevmode\nobreak\ keV}$ / $L_{\rm 2500\leavevmode\nobreak\ \AA}$]/2.605, where $L_{\rm 2500\leavevmode\nobreak\ \AA}$ and $L_{\rm 2\leavevmode\nobreak\ keV}$ are the rest-frame monochromatic luminosities at 2500Å and 2 keV, respectively (e.g. Zamorani et al. 1981). We verified that the observed distribution of $\alpha_{\rm ox}$ for HLAGN identified solely from X-rays peaks at $\alpha_{\rm ox}\sim$1, unlike the average value $\alpha_{\rm ox}\sim$1.37 found for X-ray selected AGN in the COSMOS field (Lusso et al. 2010). The lower $\alpha_{\rm ox}$ suggests that HLAGN identified only from X-rays are optically fainter than the rest of X-ray AGN in the COSMOS field, as expected from their galaxy-dominated SEDs. This is also confirmed by the fact that in more than 80% of them, the optical-NIR photometry has been fitted without AGN templates when calculating the photometric redshifts (see Marchesi et al. 2016; M. Salvato et al. in prep.).

By using different and complementary tracers of highly accreting AGN, we can build a more representative (though not 100% complete) sample of HLAGN. Our analysis would certainly benefit from optical-NIR spectroscopy to identify AGN at lower intrinsic luminosities. Unfortunately, the spectral lines used to calculate the emission line ratios [O III]/H$\beta$ and [N II]/H$\alpha$) in the BPT diagram (from Baldwin et al. 1981) are detected only in a low percentage (about 5%) of our radio sample, mostly at z$<$0.5, as some optical lines (e.g. H$\alpha$) would be redshifted outside the observed spectral window at higher redshifts. For consistency, in this work we preferred to make use of AGN selection criteria that are applicable to the entire sample.

For our 1 604 HLAGN, we took the SFR and $M_{\star}$ estimates from the best-fit solution obtained with the three-component SED-fitting code by Berta et al. (2013). This approach allows us to account for a possible AGN contamination in galaxy parameter estimates and to study uncertainties and degeneracies, which is important when comparing galaxy properties between AGN and non-AGN hosts. However, we checked that the AGN contribution to the total (8-1000$\leavevmode\nobreak\ \mu$m) IR luminosity is very small (a few percent) for most of the HLAGN sample, as argued by previous studies of X-ray and IR-selected AGN (e.g. Mullaney et al. 2011; Santini et al. 2012; Rosario et al. 2012).

4.2 Low-to-moderate radiative luminosity AGN via radio excess

Radio sources that are not classified as HLAGN in Sect. 4.1 do not show evidence of AGN activity according to X-ray, MIR, or SED decomposition. However, this does not necessarily mean that SMBH accretion is not occuring at all, but rather that the AGN diagnostics described above might fail in detecting signatures of less efficient accretion episodes. As mentioned in the previous sections, radio observations are crucial for chasing such an elusive AGN population.

To identify lower radiative luminosity AGN in our sample, we first considered the 3 GHz selected sources that are not classified as HLAGN (i.e. 79%). For each of them, the $L_{\rm IR}$ obtained via SED fitting without AGN (da Cunha et al. 2008) has been converted to IR-based SFR (SFRIR hereafter) by assuming a Kennicutt (1998) scaling factor and a Chabrier (2003) IMF. To identify the possible AGN contribution in radio emission, we analysed the ratio between the 1.4 GHz radio luminosity $L_{\rm 1.4\,GHz}$ and the SFRIR for each source. Fig. 4 shows the distribution of their ratio (i.e. $L_{\rm 1.4\,GHz}$/SFRIR, in logarithmic scale) as a function of redshift (black points). Typical 1$\sigma$ uncertainties of the observed ratio are of the order of 0.15 dex. The histogram in the top right corner shows the distribution of our sources as a function of $L_{\rm 1.4\,GHz}$/SFRIR. The distribution is skewed towards high values of the ratio, suggesting that AGN activity might be contributing to the integrated radio emission. However, the average $L_{\rm 1.4\,GHz}$/SFRIR also increases with redshift, which partly explains the skewness of the observed distribution. To quantify these factors, we split our sample into seven redshift bins (0.01$<$z$<$0.3, 0.3$<$z$<$0.7, 0.7$<$z$<$1.2, 1.2$<$z$<$1.8, 1.8$<$z$<$2.5, 2.5$<$z$<$3.5, and 3.5$<$z$<$5.7) and fit each single distribution with a log-normal function that reproduces the peak and the negative part of the observed distribution well. The values of the peak (and dispersion) of the Gaussian function identified at each redshift bin are represented in Fig. 4 with blue filled circles (and relative errors). The position of the peak generally shifts to higher $L_{\rm 1.4\,GHz}$/SFRIR ratios with increasing redshift, which justifies our choice of a redshift-dependent approach. Moreover, some recent studies (e.g. Magnelli et al. 2015; Delhaize et al. 2017) have found a slight, but significant decrease of qIR (proportional to SFRIR/$L_{\rm 1.4\,GHz}$) with redshift through a careful treatment of non-detections via stacking or double-censored survival analysis.

We calculated the 3$\sigma$ deviation from the peak of the log-normal function at each redshift bin (see red open circles in Fig. 4). By fitting the open circles with a power-law function (red solid line in Fig. 4), we derived the analytical expression that describes a redshift-dependent threshold in radio excess as follows:

[TABLE]

From this expression, we set a threshold above which the radio emission shows a $>$3$\sigma$ excess compared to that expected from star formation. The aforementioned threshold identifies 1 333 low-to-moderate radiative luminosity AGN (hereafter MLAGN) via radio excess, corresponding to 17% of our radio sample. The percentage of MLAGN, which would have been identified through a single threshold at all redshifts, would be around 18% instead of the 17% found from Eq. 1; this would not affect our results.

The choice of a 3$\sigma$ radio excess imposed in each redshift bin ensures a negligible contamination from star-forming galaxies without radio excess (about 0.15%). On the other hand, our selection might miss a significant number of potential MLAGN in our sample, which is estimated to be around 1 000 sources (75% of the identified sample of MLAGN) based on the difference between the distribution below the threshold and the best-fitting Gaussian profile, in all redshift bins. A comparison with other definitions of radio excess found in the literature is presented in Sect. 4.3.

The Kennicutt (1998) $L_{\rm IR}$–SFRIR conversion assumes that the total IR luminosity arises entirely from optically thick, dust-obscured regions. While this assumption is reasonable in highly star-forming galaxies, such as those detected by Herschel (Wuyts et al. 2011; Magnelli et al. 2013), this is not true for weakly star-forming (or passive) systems, where a significant portion of the IR luminosity may originate from ($>$ few Gyr) old stellar populations (e.g. Groves et al. 2012). We verified that the unobscured SFR derived from the UV galaxy emission (e.g. Papovich et al. 2007) is around 5% of the obscured SFRIR, on average, therefore its contribution would not significantly affect our definition of radio excess presented above.

Moreover, our threshold is calibrated on a radio-selected sample, which is expected to bias the observed distribution towards higher values of $L_{\rm 1.4\,GHz}$/SFRIR compared to the true (i.e. unbiased) distribution. These arguments suggest that our definition of radio excess is likely to be fairly conservative and the number of MLAGN selected in this way should more properly be considered as a lower limit.

In summary, by combining multiwavelength AGN diagnostics, we managed to isolate two populations of AGN in our 3 GHz selected sample. We identified 1 604 HLAGN (21%) and 1 333 MLAGN (17%), which collectively make our sample of radio-selected AGN. The remainder of the sample (62%) is characterised in detail in Smolcic et al. (2017a). Moreover, we found that about 30% of HLAGN shows also a $\geq$3$\sigma$ radio excess. We checked that the relative numbers of HLAGN identified from each criterion and shown in Fig. 2 do not change between sources with and without significant radio excess. A more detailed analysis of the average galaxy and AGN properties for the aforementioned AGN classes (MLAGN, HLAGN, including the subsample with radio excess) is presented in Sect. 6.

4.3 Comparison with radio-based classifications

In this section we compare our source classification with other independent methods from the literature.

4.3.1 Comparison with VLBI sources

We cross-matched our 3 GHz selected sample of 7 729 sources with the Very Large Baseline Interferometry (VLBI) 1.4 GHz source catalogue from N. Herrera Ruiz et al. (in prep.). The authors targeted the radio sources selected by the 1.4 GHz VLA-COSMOS survey (Schinnerer et al. 2010) with VLBI at $\lesssim 0\aas@@fstack{\prime\prime}01$ angular resolution, reaching a 1$\sigma$ sensitivity of 10$\leavevmode\nobreak\ \mu$Jy beam*-1* in the central part of the field. They detected 468 radio sources at signal-to-noise ratio higher than 5.5. A total of 354 matches have been found within $0\aas@@fstack{\prime\prime}4$ (half-beam size of 3 GHz VLA observations). Given the high angular resolution, VLBI is sensitive to the radio emission on circum-nuclear scales (from $d\sim$20 pc at $z$=0.1 to $d\sim$80 pc at $z$=2), likely arising from an AGN.

Interestingly, we found that 91% of VLBI sources are classified as AGN by our method, where 55% are MLAGN and 36% are HLAGN. About 88% of the HLAGN also show a $>$3$\sigma$ radio excess compared to the SFRIR (Sect. 4.2). These notably high percentages of VLBI sources classified as AGN (both HLAGN and MLAGN) in our sample demonstrate the high reliability of our classification method.

4.3.2 Comparison with Bonzini et al. (2013)

In a similar study, Bonzini et al. (2013, B13 hereafter) carried out a panchromatic analysis of about 800 high-redshift (z$\leq$4) radio sources in the E-CDFS, selected with VLA at 1.4 GHz. They separated radio sources into radio-loud AGN (RL-AGN), radio-quiet AGN (RQ-AGN) and star-forming galaxies (SFGs) on the basis of the observed 24$\leavevmode\nobreak\ \mu$m-to-1.4 GHz flux ratio (also called q24,obs, see Donley et al. 2005). Briefly, B13 selected RL-AGN that are below the 2$\sigma$ deviation from the average q24,obs at a given redshift. In case q24,obs was an upper limit, the authors selected RL-AGN that are below the 1$\sigma$ deviation from the average value. Within and above the 2$\sigma$ deviation, they selected RQ-AGN that are not RL-AGN and at the same time fulfilling either X-ray or MIR diagnostics (similar to those discussed in Sect. 4.1); the rest of the sample was classified as SFGs. On top of these criteria, B13 applied further checks (see their Sect. 3.5.1) to improve the sample characterisation, which led the authors to reclassify 11 sources from SFGs to RQ or RL AGN in their sample.

Despite the larger area, our 3 GHz data in COSMOS are deeper than the E-CDFS data at 1.4 GHz, which allows us to compare our classification with B13 in the same flux density range. First, we scaled our 3 GHz flux density to the observed 1.4 GHz for each source (as discussed in Sect. 2.1). Secondly, we cut our sample at total 1.4 GHz flux $S_{\rm 1.4\leavevmode\nobreak\ GHz}>$ 37$\leavevmode\nobreak\ \mu$Jy as in B13, yielding 4 573 sources (around 59% of the full sample) and we computed q24,obs for all of them. Thirdly, we applied the same criteria of B13 to identify RQ-AGN, RL-AGN and SFGs in our sample (including their additional AGN diagnostics for consistency, see Sect. 4.4.2), and show the numbers in Table 1.

This comparison suggests that the HLAGN and MLAGN identified in this work fairly overlap with the RQ-AGN and RL-AGN classes, respectively, even if our classification tends to classify more objects as AGN (45%) than the B13 classification (32%). The RQ-AGN percentage found by B13 (24%) is higher than that listed in Table 1 (13%), and is likely driven by the higher incompleteness of X-ray observations in the COSMOS field towards moderately luminous X-ray AGN ($10^{42}<$ $L_{\rm x}<10^{44}$ erg s*-1*), especially at high redshifts (z$>$2). The aim of this comparison is not to invalidate either of the methods, but simply to clarify how different nomenclatures compare to each other. The main difference in the source classification lies in the different definitions of radio excess. In their work, B13 used a redshift-dependent threshold in q24,obs, which was calibrated on the M82 template. The dash-dotted line in Fig. 5 (top panel) shows the redshift-dependent threshold defined by B13. Black symbols indicate the distribution of our 3 GHz radio sources (7 729) as a function of q24,obs and redshift with our radio-excess sources (1 814 in total, being 1 333 in MLAGN and 481 in HLAGN) highlighted in green. Downward pointing arrows indicate (5$\sigma$) upper limits due to non-detection at 24$\leavevmode\nobreak\ \mu$m. The black dashed line sets the threshold q${}_{\rm 24,obs}<$ 0 adopted by Donley et al. (2005) to identify radio-excess sources. Del Moro et al. (2013) already pointed out that q24,obs is a reliable tracer of radio excess, although it is not complete. We confirm this statement, as most of the sources below the q24,obs threshold (by B13) also satisfy our radio-excess definition. However, the B13 criterion at z$>$1 becomes even more stringent than that proposed by Donley et al. (2005). This decreasing trend of q24,obs with redshift is driven by the shape of the M82 template SED, which is rather peculiar compared to the average SED of star-forming galaxies at 0$<$z$<$3, implying a steeper decline with redshift compared to what is observed in our sample. However, the percentage of z$>$2 sources in B13 is relatively small compared to our sample, which implies that the M82 template SED shape should not have a large effect on the source classification in the E-CDFS sample.

4.3.3 Comparison with Del Moro et al. (2013)

An alternative method to search for radio excess is by using the (rest-frame) FIR-to-radio flux ratio $q_{\rm FIR}$ (e.g. Sargent et al. 2010), where the FIR flux refers to the rest-frame range 42.5–122.5$\leavevmode\nobreak\ \mu$m, and the radio flux refers to the rest-frame 1.4 GHz. Recently, Del Moro et al. (2013) used $q_{\rm FIR}<$ 1.68 in star-forming galaxies within the GOODS-North field to identify sources with $>3\sigma$ radio excess.

As proposed by Del Moro et al. (2013, see their Figure 5), we calculate the qFIR for our 3 GHz sources and show in Fig. 5 (bottom panel) their distribution with redshift with respect to the threshold set by Del Moro et al. (2013). Del Moro et al. (2013) calibrated the qFIR threshold on a sample of sources detected at both 1.4 GHz and 24$\leavevmode\nobreak\ \mu$m, which are, therefore, on average more star-forming than our radio-selected galaxies in COSMOS. If limiting our sample to 24$\leavevmode\nobreak\ \mu$m detected sources, we estimate the percentage of radio-excess sources to be around 13%, which is in agreement with the percentage found by Del Moro et al. (2013). We checked that the average qFIR values of our 3 GHz sources also detected at 24$\leavevmode\nobreak\ \mu$m closely resemble those of Del Moro et al. (2013) in the common redshift range. This is expected, given that the radio-to-24$\leavevmode\nobreak\ \mu$m sensitivity limits are similar between the GOODS-North and COSMOS fields. For comparison with Del Moro et al. (2013), we calculated the $q_{\rm FIR}$ for our radio-excess sources and verified that 100% of these sources satisfy the condition $q_{\rm FIR}<$ 1.68, while around 72% of sources with $q_{\rm FIR}<$ 1.68 satisfy our definition of radio excess. This check further supports the reliability of our definition of radio excess.

4.4 Further tests of the source classification

The classification scheme presented in Sects. 4.1 and 4.2 was based on a few assumptions that we test in this section. In particular, we detail the motivation for our naming convention (Sect. 4.4.1), show how our classification would change if considering additional AGN diagnostics (Sect. 4.4.2), and compare the IR luminosities derived for Herschel-undetected sources against IR stacking (Sect. 4.4.3).

4.4.1 The choice of the naming convention

Our samples of HLAGN and MLAGN include AGN identified through different diagnostics. As shown in Fig. 2, only 14% of HLAGN meet simultaneously the three diagnostics (X-ray, MIR, and SED decomposition), suggesting that various AGN selection criteria might be sensitive to a broad range of AGN luminosities. We study the distribution of MLAGN and HLAGN as a function of AGN radiative luminosity (Lrad,AGN), as derived from SED-fitting decomposition, converted to BHAR following Eq. 1 from Alexander & Hickox (2012) and assuming a canonical mass-energy efficiency conversion of 10% (e.g. Marconi et al. 2004). Each estimate of Lrad,AGN is calculated from the corresponding best-fit AGN template (Sect. 3) obtained from SED-fitting decomposition. This parameter should be considered an upper limit if the AGN component is not significant from SED-decomposition, which is the case for all MLAGN (by definition) and for about 27% of HLAGN identified solely from X-ray or MIR criteria. In the latter case, the Lrad,AGN has been taken from the 95th percentile of the corresponding PDF obtained from the sed3fit code, which is equivalent to an upper limit at 90% confidence level.

Fig. 6 shows the Lrad,AGN distribution separately for MLAGN (red dashed line), HLAGN (blue solid line), and the subsample of HLAGN not identified from SED fitting (blue dashed line) in seven redshift bins. The vertical lines represent the median of the corresponding distribution at each redshift bin. The median values show a significant difference (around 1 dex) in Lrad,AGN between MLAGN and HLAGN at all redshifts. Despite the non-negligible overlap between the two distributions, we stress that the Lrad,AGN estimates for MLAGN consist of upper limits, which implies that the difference between the true distributions, at a given redshift, is even more significant.

This test suggests that, at each redshift, HLAGN are significantly more powerful than MLAGN and justifies the naming convention proposed in this work in a statistical sense.

4.4.2 Additional AGN diagnostics

As carried out in B13, we applied additional diagnostics to verify the robustness of our AGN identification method.

• •

Optical spectra: we searched for sources flagged as broad line AGN in the optical spectra taken from the zCOSMOS-Bright survey (Lilly et al. 2007, 2009); we found 15 in total, which were all pre-classified as HLAGN on the basis of X-ray and mid-IR criteria.

• •

VLBI sample: as mentioned in Sect. 4.3.1, 354 radio sources selected at 3 GHz were also identified in the VLBI sample (N. Herrera Ruiz et al. in prep.) available in the COSMOS field. Of these, 30 (8%) were not identified as either MLAGN or HLAGN, although they are likely to be AGN.

• •

Inverted radio spectra: we found a total of 11 radio sources not classified as radio-selected AGN, but that have inverted radio spectra ($\alpha>$0), where the spectral index $\alpha$ is set to the observed 1.4–3 GHz slope. This feature may indicate the presence of a compact radio core (e.g. Kellermann & Pauliny-Toth 1969).

• •

Hardness ratio: we found a total of seven X-ray sources that were not classified as AGN in our sample with a positive hardness ratio ($HR>$0), indicating the likely presence of obscured AGN (e.g. Brusa et al. 2010).

• •

Polycyclic aromatic hydrocarbon (PAH) features: sources lying in the IRAC colour-colour wedge discussed in B13 (see their Sect. 3.5.5) likely display PAH-dominated SEDs, which are typical of star-forming galaxies. We found 444 sources that fulfil this criterion with only five (26) that are classified as MLAGN (HLAGN) in our sample.

The above-mentioned criteria suggest that some of our sources might be misclassified. We found that 46 radio sources (30 from VLBI detection, 9 from inverted radio spectra, and seven from the hardness ratio) should be reclassified from non-AGN to AGN in our sample. However, these criteria are applicable to a very low percentage of our sample, meaning that by incorporating them we would likely introduce a bias against sources with no available diagnostics. Reclassifying these 46 sources (0.6% of our sample) would have no impact on our main results and conclusions. We also found 31 radio-selected AGN displaying PAH-dominated SEDs in the MIR; even though this criterion does not rule out a potential AGN contribution in other bands, the percentage of possibly misclassified AGN is minimal. For these reasons, we limited our AGN selection criteria to those presented in Sects. 4.1 and 4.2.

4.4.3 Testing radio excess with infrared stacking

As discussed in Sect. 3, the $L_{\rm IR}$ and SFR estimates are less robust for sources without Herschel detection, which constitute about one-third of our 3 GHz selected sample. For these sources, we fitted their SEDs using the nominal upper limits in the PACS and SPIRE bands (as described in Sect. 3). As a sanity check, we compared the Herschel fluxes used for SED fitting with those obtained via Herschel stacking. The difference between the two approaches enables us to test the robustness of our radio-excess definition.

We looked at the 79% of sources not classified as HLAGN, shown in Fig. 4, and considered only those without $\geq$3$\sigma$ Herschel detection in any PACS or SPIRE band, which amounts to 2 203 in total. As carried out in previous papers (e.g. Santini et al. 2012; Rosario et al. 2012; Bonzini et al. 2013), we split this sample in different redshift bins (see Table 4) and we stacked the PACS and SPIRE images at the optical-NIR position of each source in the same bin using a stacking tool from Béthermin et al. (2010)999This IDL routine is described in detail in Béthermin et al. (2010) and can be retrieved at https://www.ias.u-psud.fr/irgalaxies/downloads.php..

For PACS images, we performed a mean stacking on the residual maps at 100 and 160$\leavevmode\nobreak\ \mu$m, from which all detections were removed to avoid contamination by nearby brighter sources. Point spread function (PSF) photometry was performed on the final stacked images, using the uncertainty maps (Lutz et al. 2011) as weights. The stacked fluxes were corrected for aperture and correlated noise101010The full PACS documentation is available at www.mpe.mpg.de/resources/PEP/DR1_tarballs/.. The SPIRE images111111The SPIRE images for the COSMOS field were taken from the HeDAM database: http://hedam.lam.fr/HerMES/index/download. at 250, 350, and 500$\leavevmode\nobreak\ \mu$m were already given in units of surface brightness (Jy beam*-1*), hence we inferred the stacked flux directly from the value of the central pixel in the stacked image. To evaluate the uncertainty on the stacked fluxes, we performed a bootstrapping analysis (e.g. Shao et al. 2010; Santini et al. 2012; Rosario et al. 2012). Briefly, a set of $N$ sources, where $N$ is equal to the number of stacked sources per redshift bin, is randomly chosen 1 000 times, allowing repetition of the same source. To mitigate the possible contamination from a few brighter outliers, we set the stacked flux to the median of the distribution obtained from bootstrapping, while the error on the flux was drawn by the 16th and 84th percentiles of the same distribution. For PACS images, we also corrected the errors for high-pass filtering effects. We did not correct the stacked fluxes for possible blending in the optical-NIR images. In every redshift bin we obtained $>$2$\sigma$ detection in two to five Herschel bands. The final stacked fluxes and corresponding uncertainties are summarised in Table 4.

We re-fitted the 2 203 SEDs again with magphys, using the corresponding stacked fluxes as detections for each source. The resulting $M_{\star}$ estimates agree very well (no offset, 1$\sigma$ dispersion is about 0.15 dex) with those obtained in Sect. 3. The newly derived $L_{\rm IR}$ estimates are generally lower (by 50%) than those obtained in Sect. 3, but for 40% of the stacked sample we found slightly higher $L_{\rm IR}$. This is mainly because while in most cases the fluxes obtained from Herschel stacking (see Table 4) are lower than the flux values used for SED fitting, in a few bins they are instead slightly higher, especially at z$\gtrsim$1.5. However, we verified that using the $L_{\rm IR}$ estimates derived from stacking would only minimally affect the overall distribution of $L_{\rm 1.4\,GHz}$/SFRIR (Fig. 4). In particular, we checked that our classification would remain unchanged for 90% of the sources, while the remainder of the sample would move either from SFGs to MLAGN (5%) or vice versa (5%).

The purpose of this test was to quantify the impact of using different sets of Herschel fluxes on the source classification. The general agreement obtained between the two methods ensures the robustness of our classification. As a consequence, we decided to keep using the Herschel upper limits introduced in Sect. 3 through the rest of this work.

5 Catalogue description

The value-added catalogue presented in this section includes classification and selected physical properties used in this work for our 3 GHz radio sample with optical-NIR counterparts (7 729 sources in total). We also list the individual criteria used in this work to classify our radio sources (columns 14 to 17). This way, any user can easily retrieve our classification or adjust it to a different set of selection criteria. The catalogue will be made available through the COSMOS IPAC/IRSA database121212http://irsa.ipac.caltech.edu/Missions/cosmos.html. Here we describe its structure, following the same format of Table 5.

• •

(1) Identification number of the radio source (ID).

• •

(2) Right ascension (J2000) of the radio source.

• •

(3) Declination (J2000) of the radio source.

• •

(4) Best redshift available for the source.

• •

(5) Origin of the redshift: spectroscopic (“spec”) if available, photometric (“phot”) otherwise.

• •

(6) 3 GHz integrated radio flux density [$\mu$Jy].

• •

(7) 3 GHz (rest-frame) radio luminosity [log W Hz*-1*].

• •

(8) 1.4 GHz (rest-frame) radio luminosity [log W Hz*-1*], obtained as described in Sect. 2.1.

• •

(9) Star formation infrared (8-1000$\leavevmode\nobreak\ \mu$m rest-frame) luminosity derived from SED fitting [log $L_{\odot}$]. If the source is classified as HLAGN, this value represents the portion of the total infrared luminosity arising from star formation, while it corresponds to the total IR luminosity otherwise (see Sect. 3).

• •

(10) Flag for Herschel detection at $\geq$3$\sigma$, in at least one band: “true” if detected, “false” if only upper limits are available.

• •

(11) Stellar mass derived from SED-fitting decomposition [log $M_{\odot}$]. The value is drawn from the fit with AGN if the source is classified as HLAGN and otherwise from the fit without AGN. Calculated with a Chabrier (2003) IMF.

• •

(12) Star formation rate [$M_{\odot}$yr*-1*] obtained from the total infrared luminosity listed in column (9), assuming the Kennicutt (1998) conversion factor, and scaled to a Chabrier (2003) IMF.

• •

(13) Rest-frame [NUV-$r$] colour obtained from the best-fitting galaxy template and corrected for dust attenuation (AB magnitude).

• •

(14) X-ray-AGN: “1” if true, “0” otherwise.

• •

(15) MIR AGN: “1” if true, “0” otherwise.

• •

(16) SED-AGN: “1” if true, “0” otherwise.

• •

(17) Radio-excess: “1” if true, “0” otherwise.

• •

(18) Class: moderate-to-high radiative luminosity AGN (HLAGN), low-to-moderate radiative luminosity AGN (MLAGN), or neither of the two (empty space). A source is classified as HLAGN if (14)=1 $\vee$ (15)=1 $\vee$ (16)=1, while it is classified as MLAGN if (14,15,16)=(0,0,0) $\wedge$ (17)=1.

6 Results

In this section we present the average AGN and galaxy properties for our 3 GHz radio sources classified as AGN host galaxies. In particular, we focus our analysis on MLAGN (17%), HLAGN (21%), and also the subsample of HLAGN with radio excess (6%). We show, for these classes, the location in the SFR–$M_{\star}$ plane (Sect. 6.1) and the distributions of AGN and galaxy properties of their hosts (Sect. 6.2) at different cosmic epochs.

6.1 The SFR–$M_{\star}$ plane of radio-selected AGN

Fig. 7 shows the 2D density contours in SFR–$M_{\star}$ plane for our samples of MLAGN (red), HLAGN (blue), and for the subsample of HLAGN with radio excess (blue dashed contours). Black dots represent our joint sample of aforementioned AGN at different redshifts. The 2D density contours enclose the sources (from outer to inner contours) with density levels $>$35, $>$50, $>$68, $>$80, $>$90, and $>$95% of the maximum 2D density for a given class and redshift bin. The black dashed line marks the so-called main sequence (MS) of star-forming galaxies (taken from Whitaker et al. 2012), which is known to evolve positively with redshift (e.g. Noeske et al. 2007; Elbaz et al. 2011; Speagle et al. 2014; Schreiber et al. 2015).

The $M_{\star}$ and SFR estimates for each source were computed directly by the three-component fit, hence already correcting for a possible AGN contamination (Sect. 3). Typical 1$\sigma$ uncertainties on $M_{\star}$ and SFR are of the order of 0.1 dex, but for Herschel undetected sources the uncertainty in SFR is around 0.2 dex.

Table 2 summarises the number of sources shown in Fig. 7 for each class and redshift bin. These numbers show that the redshift distribution for both HLAGN and MLAGN peaks around z$\sim$1. The two AGN populations become comparable around z$\sim$1.5, while at higher and lower redshift the HLAGN generally outnumber MLAGN. The percentage of HLAGN with a $>$3$\sigma$ radio excess is roughly constant with redshift (around 25–35% of the HLAGN sample).

Our HLAGN and MLAGN appear to lie in different regions of the SFR–$M_{\star}$ plane, at various redshifts. At low redshift (z$<$0.3), the two AGN classes show rather distinct $M_{\star}$ distributions, where MLAGN are more than two times more massive than HLAGN (1011 $M_{\odot}$ versus 10${}^{10-10.5}\leavevmode\nobreak\ M_{\odot}$). This difference is consistent with that found by previous studies in the local Universe (e.g. Smolčić et al. 2009; Best & Heckman 2012). However, we notice that HLAGN with radio excess predominantly lie in the high-$M_{\star}$ tail of the HLAGN population and closely resemble the distribution of MLAGN at this redshift (z$<$0.3). Therefore, this subsample seems to show intermediate SFR and $M_{\star}$ distributions between the two AGN categories. At higher redshift (0.3$<$z$<$1.8), the two AGN populations show a larger overlap in the SFR–$M_{\star}$ plane than that observed at lower redshift. However, the bulk of HLAGN is generally located around the MS relation, while MLAGN preferentially lie in the lower part of the MS with typically lower SFRs (by a factor of 2–3) compared to HLAGN. The subsample of HLAGN with radio excess is mostly concentrated between the locations of the two main AGN classes. At even higher redshift (z$>$1.8) the overlap between the two distributions decreases, with HLAGN, largely located above the MS relation, having on average both higher M*⋆* and SFR than the MLAGN. A detailed description of the distributions of host-galaxy properties is given in Sect. 6.2. The highest redshift bin contains only a few tens of sources, which explains the noise seen in the 2D density contours.

This difference in SFR is consistent with the different percentages of Herschel-detected sources between the two AGN classes: 57% of HLAGN has a $>$3$\sigma$ detection in at least one Herschel band, while this percentage decreases to only 17% for MLAGN. These numbers slightly decrease with redshift for either classes, from 81% (27%) at z$<$0.3 to 50% (17%) in the highest redshift bin for HLAGN (MLAGN). These numbers further suggest that the host galaxies of HLAGN are typically star forming at all redshifts. Although this work does not aim to investigate the nature of the sources not classified as AGN, we checked that the subsample of sources that are neither HLAGN nor MLAGN (i.e. 62%) are mostly located on the MS relation at all redshifts, thus resembling the distribution of HLAGN in the SFR–$M_{\star}$ plane. This suggests that the remainder of our sample might consist of mostly star-forming galaxies.

A comprehensive study of the radio-AGN population in the SFR–M*⋆* plane has been presented by Bonzini et al. (2015). They found that most of RQ AGN show significant star formation in their hosts, and typically (75%) lie along the MS relation, likewise SFGs, at various redshifts. Moreover, Bonzini et al. (2015) found that the majority of RL AGN reside in less star-forming galaxies, which are predominantly located below the MS. Despite the different nomenclature and sample selection used by the authors (see Table 1 for a comparison), the qualitative agreement with their results is reassuring.

6.2 Physical properties of AGN hosts

In this section we investigated the distributions of galaxy and AGN properties for the populations of HLAGN and MLAGN. We applied a two-sample Kolmogorov-Smirnov (K-S) test to quantify the difference or similarity between two distributions. This statistical test allows us to determine if two input data sets could be drawn from a common parent distribution without any assumption about its shape. We run this test to evaluate the difference between HLAGN and MLAGN in terms of various galaxy properties. The probabilities P(K-S) obtained for each parameter and redshift bin are listed in Table 3, along with the corresponding $\sigma$ level. Smaller values of P(K-S) indicate a more significant difference between the two data sets.

6.2.1 Distribution of galaxy properties

Figs. 8, 9, and 10 show the distributions of $M_{\star}$, SFR, and rest-frame [NUV–$r$] colours, respectively, for the following classes: MLAGN (red), HLAGN (blue), and the subsample of HLAGN with radio excess (blue thicker distribution). Vertical lines show the median value of the corresponding distribution. The distributions are shown in seven redshift bins out to z$\lesssim$6 and are normalised to the highest maximum value of the two distributions.

As mentioned in Sect. 6.1, the $M_{\star}$ distributions of MLAGN at low redshift are skewed towards higher $M_{\star}$ compared to HLAGN, and the difference remains significant up to z$\sim$1 at $\gtrsim$99% level (see Table 3). At z$\sim$1.5 (4th redshift bin) the two distributions appear more similar. At higher redshifts (z$\sim$2), we observe a possible reversal of the $M_{\star}$ distributions with the HLAGN populating the high-$M_{\star}$ tail. At this redshift, the two-sample K-S test finds an almost 6$\sigma$ difference between the two distributions. However, we are not able to confirm or disclaim this trend at z$>$2.5, given the tentative significance (about 2$\sigma$) of the results obtained from the K-S test. A more detailed discussion and interpretation of these trends is presented in Sect. 7. As seen in the SFR–$M_{\star}$ plane (Fig. 7), the subsample of HLAGN with radio excess overlaps significantly with the distribution of MLAGN, showing intermediate $M_{\star}$ between the two AGN classes (except in the highest redshift bin).

In Fig. 9 we show the same plots for the SFR, obtained by integrating the best-fit galaxy template over the range (rest-frame) 8–1000$\leavevmode\nobreak\ \mu$m, and by assuming a Kennicutt (1998) scaling factor and a Chabrier (2003) IMF. As already seen in Fig. 7, we confirm that the HLAGN with radio excess populate the lower tail of the SFR distribution, overlapping significantly with MLAGN. The difference between HLAGN and MLAGN in SFR remains visible and $\gtrsim 99$% significant in the higher redshift bins, up to z$\sim$3, as well.

Rest-frame [NUV-$r$] colours were calculated from the best-fit galaxy template of each source and also corrected for dust attenuation. The distribution of [NUV-$r$] colours shown in Fig, 10 confirms that most of HLAGN have blue or green rest-frame optical colours ([NUV-$r$]$<$3.5; Ilbert et al. 2010) at all redshifts. On the other hand, MLAGN are more pronounced towards quiescent systems ([NUV-$r$]$>$3.5; Ilbert et al. 2010), at least up to z$\sim$1. This bimodality in the colour distributions becomes progressively less pronounced at higher redshift, also showing that the host-galaxies of MLAGN become, on average, more star forming with increasing redshift (see Sect. 4.4.3). It is recognised that the number density of quiescent galaxies (selected via optical colours) at $M_{\star}>$1010 $M_{\odot}$ decreases with increasing redshift (e.g Brammer et al. 2011, Ilbert et al. 2013). However, the difference in [NUV-$r$] between HLAGN and MLAGN remains highly significant up to z$\sim$3.5, while it disappears at 3.5$<$z$<$5.7. The subsample of HLAGN with radio excess shows intermediate colours between the rest of HLAGN and the population of MLAGN.

6.2.2 Distribution of AGN properties

Fig. 11 shows the distribution of both AGN radiative and mechanical power for our classes of AGN.

We calculated the AGN radiative power (Lrad,AGN) of HLAGN from the best-fit AGN template obtained with the three-component SED-fitting code sed3fit for sources both with and without radio excess. The typical uncertainties on Lrad,AGN are around 0.4 dex for sources with ($\geq$99%) significant AGN component (i.e. SED-AGN, see Sect. 4), while we took the Lrad,AGN from the 95th percentile of the corresponding PDF obtained from the sed3fit code for HLAGN not identified as such from SED decomposition (30% of HLAGN); this is equivalent to an upper limit at 90% confidence level on the AGN radiative luminosity. Fig. 11 shows the normalised distributions of Lrad,AGN, separately for HLAGN (blue) and for the subsample with radio excess (blue thicker). The distributions of Lrad,AGN cover a broad range ($>$3 dex) in each redshift bin, which is around 1042-45 erg s*-1* at z$<$0.3, and 1044-47 erg s*-1* at z$\sim$3. As mentioned in Sect. 4.4.1, the percentage of HLAGN identified from SED fitting does not depend on the presence of a radio excess. Therefore, applying the upper limits at 90% confidence level on Lrad,AGN does not affect the ratio of the Lrad,AGN distributions between HLAGN and their subsample with radio excess.

We calculated the rest-frame 1.4 GHz radio luminosity $L_{\rm 1.4\,GHz}$ for each source by scaling its radio flux from 3 GHz to 1.4 GHz and taking the observed 1.4–3 GHz spectral index $\alpha$, as explained in Sect. 2.1. The presence of a $>$3$\sigma$ radio excess suggests that a notable portion of the radio emission is not arising from star formation processes in the host, but possibly from the central SMBH. For this reason, each $L_{\rm 1.4\,GHz}$ measurement was scaled to the portion associated with AGN activity, based on the deviation of the observed $L_{\rm 1.4\,GHz}$ - to - SFRIR ratio from the peak of the Gaussian function (associated with star formation) at the corresponding redshift bin (blue points in Fig. 4).

We converted the AGN-related radio emission to AGN mechanical power (Lmech,AGN) of the radio jet, by assuming the redshift-independent relation by Willott et al. (1999), which is based on theoretical grounds and adopted in other studies (e.g. Merloni & Heinz 2008; La Franca et al. 2010, see Best & Heckman 2012 for a review). We used this relation expressed in terms of $L_{\rm 1.4\,GHz}$ (see Eq. 1 from Heckman & Best 2014). Willott et al. (1999) combined all of the uncertainties on this relation into a single factor, $f_{\rm W}$, which can range between 1 and 20. This scaling factor is still a matter of debate in the literature (Godfrey & Shabala 2016). Nonetheless, following the approach of numerous studies (e.g. Merloni & Heinz 2007; Smolčić et al. 2009; La Franca et al. 2010; Best & Heckman 2012; Pracy et al. 2016), we make the simplistic assumption that this relation holds at all radio luminosities $L_{\rm 1.4\,GHz}$ that are probed by our sample.

The normalised distributions of Lmech,AGN are shown in Fig. 11 for both MLAGN (orange) and HLAGN with radio excess (orange thicker distribution). The typical range in Lmech,AGN probed by the distributions is about 2 dex wide in each redshift bin, which is around 1041-43 erg s*-1* at z$<$0.3 and 1043-45 erg s*-1* at z$\sim$3. Fig. 11 shows the distributions of Lmech,AGN by taking $f_{\rm W}$=5, which is consistent with the relation derived by Daly et al. (2012). The vertical lines indicate the median value of the corresponding distribution. However, we calculated the range within which the median value could shift, by changing fW between fW=1 and fW=20, which is shown by the orange horizontal lines around the median.

Our analysis seems to suggest that the overall AGN properties observed for HLAGN with radio excess are similar to the rest of HLAGN in terms of radiative power and are also consistent with MLAGN in terms of mechanical power. The subsample of HLAGN with radio excess (6% of our parent 3 GHz radio sample) is particularly interesting because it enables a direct comparison between radiative and mechanical AGN power for the same sources. Despite the uncertainties on the relation proposed by Willott et al. (1999), we show that the AGN mechanical power in HLAGN with radio excess is typically lower than (or at most marginally comparable to) the AGN radiative power, depending on fW, although in some cases Lmech,AGN can exceed Lrad,AGN (see Heckman & Best 2014). While the AGN power of HLAGN occurs predominantly in radiative form, MLAGN display a substantial mechanical AGN luminosity component. These properties may suggest that HLAGN and MLAGN samples qualitatively resemble radio AGN types often referred to as radiative mode (or HERG) and jet mode (or LERG), respectively. In addition, we note that MLAGN have significantly lower Lrad,AGN than HLAGN with radio excess, despite both classes showing a relatively high radio loudness. As a consequence, a simple RL–RQ separation would not allow such direct insight into the fundamental properties of AGN.

7 Discussion

A radio-based selection allows us to study a mixture of galaxy populations that are powered by either star formation, AGN activity, or both. It is also crucial to exploit multiwavelength ancillary data to reach a more comprehensive perspective of the nature of our sources. In this work, we made use of this approach to derive integrated AGN and galaxy properties and compare them between different AGN classes, and over a wide range of radio luminosity and redshift. In this section, we discuss and interpret our findings in the context of current AGN and galaxy evolutionary scenarios.

7.1 Radio emission in HLAGN and MLAGN

The origin of radio emission in the sub-mJy radio population is still a matter of debate. Recent studies based on interferometric radio observations of sub-mJy radio sources have the potential to shed light on this issue. For example, Herrera Ruiz et al. (2016) analysed in detail the interferometric images of three RQ-AGN obtained with VLBI in the COSMOS field, which are part of the sample described in Sect. 4.3.1. The comparison between VLBI and VLA fluxes suggested that 50–75% of the radio emission in these sources is arising from non-thermal AGN activity. We note that these sources would have been classified as HLAGN with radio excess according to our method. Similar conclusions were reached by Chi et al. (2013) and Maini et al. (2016) on different samples of RQ-AGN, supporting the idea that some radio sources could be predominantly powered by AGN activity.

At low redshift (z$<$0.3), independent hints on the origin of radio emission in the sub-mJy radio population were provided by Kimball et al. (2011), who constructed the 6  GHz radio luminosity function from a sample of QSO host galaxies at 0.2$<$z$<$0.3. They concluded that radio emission in sources with 6  GHz luminosity $L_{\rm 6\leavevmode\nobreak\ GHz}>$1022.5 W Hz*-1* was AGN related, while in fainter sources it was mainly driven by star formation. However, the parent sample analysed by the authors consists of optically identified QSOs from the SDSS, which are systematically more powerful, relative to our 3 GHz sample of AGN in COSMOS, due to our smaller comoving volume covered at z$<$0.3. Support for radio emission that is powered by AGN activity was also provided by White et al. (2015), who studied a sample of RQ-QSOs at 1.4 GHz flux $S_{\rm 1.4\leavevmode\nobreak\ GHz}<$1 mJy. The authors stress, however, that their analysis may be biased towards the brightest optically identified QSOs.

A complementary view on this topic benefits from deeper radio surveys, which can push this analysis to higher redshifts and to intrinsically fainter radio sources. For example, Bonzini et al. (2015) and Padovani et al. (2015) investigated the origin of radio emission in RQ and RL AGN in the E-CDFS down to 37$\leavevmode\nobreak\ \mu$Jy (5$\sigma$). They found a mixture of AGN and SFGs contributing to the sub-mJy radio population, where RQ AGN is predominantly powered by star formation. We checked this by exploiting a larger sample of 3 GHz selected sources with optical-NIR counterparts, counting in total 1 604 HLAGN and 1 333 MLAGN (Sect. 4) out to z$\lesssim$6. The analysis presented in Sects. 4.1 and 4.2 suggests that roughly 70% of the HLAGN does not show a $\geq$3$\sigma$ radio excess, which might suggest that radio and infrared emission in HLAGN are commonly (to a certain amount) powered by star formation in their hosts, as proposed by previous studies (e.g. Morić et al. 2010; Baldi et al. 2013; Padovani et al. 2015). However, the radio excess detected for the remaining 30% is a potential signature of radio-selected AGN activity, possibly linked to jet-mode (or radio-mode) feedback, often referred to in the literature (e.g. Hardcastle et al. 2007; Best & Heckman 2012; Heckman & Best 2014). On the other hand, radio emission in our sample of MLAGN is predominantly arising from non-thermal radiation likely ascribed to AGN activity, rather than star formation in their hosts. These results agree with the conclusions presented by Padovani et al. (2015) and Bonzini et al. (2015) for a sample of RQ-AGN and RL-AGN, supporting the composite nature of the sub-mJy radio source population (e.g. Smolčić et al. 2008; Padovani et al. 2011; Baldi et al. 2014).

7.2 Radio AGN in the context of galaxy evolution

We attempt to interpret the nature of our HLAGN and MLAGN populations in the framework of AGN-galaxy evolution. As suggested by previous authors (e.g. Hardcastle et al. 2007; Smolčić 2009; Best & Heckman 2012; Padovani et al. 2015), in the local Universe the HERG and LERG classes show a clear dichotomy in terms of AGN and host-galaxy properties. These findings have been interpreted within a self-consistent evolutionary scenario, where HERG and LERG trace earlier and later stages, respectively, of galaxies’ life cycle (see Heckman & Best 2014 for a comprehensive review).

At higher redshift, Merloni & Heinz (2008) proposed a model to reproduce the kinetic and radiative luminosity function of AGN in which the highly efficient accretion onto the SMBH can produce both kinetic and radiative feedback (e.g. Veilleux et al. 2013), which are consistent with the AGN properties observed for our HLAGN with and without radio excess, respectively. Nevertheless, the power from weakly accreting SMBHs ($\lambda_{\rm Edd}\leq$ 10*-2*, also named “advection-dominated accretion flow”, ADAF; e.g. Blandford & Znajek 1977) is mainly in the form of kinetic feedback (Bower et al. 2006; Fanidakis et al. 2011), linking to the properties of our MLAGN population.

Semi-analytic models predict different accretion modes between highly and weakly accreting AGN. On the one hand, highly accreting AGN have been usually connected to a fast gas accretion mode in galaxy halos in which the free fall times are usually longer than the cooling times. On the other hand, weakly accreting AGN are in the regime of slow gas accretion, where cooling time is much larger than the free fall time (e.g. Fanidakis et al. 2011, 2012).

From an observational point of view, we found that galaxies hosting MLAGN are more massive, redder, and less star forming compared to HLAGN, at least up to z$\sim$1. In particular, the most massive galaxies ($M_{\star}\sim$1011 $M_{\odot}$) at these redshifts typically host MLAGN, while the $M_{\star}$ distributions of HLAGN and MLAGN become comparable at z$\sim$1.5 and display a reversal at z$\sim$2. This trend is unlikely to be driven by the incompleteness in M*⋆, as the optical-NIR selected sample in the COSMOS field is $>$80% complete at M${}_{\star}>$109.7M☉* out to z$\sim$4 (see Davidzon et al. 2017). We stress that this $M_{\star}$ behaviour is observed in our radio-selected sample, while it might not be the same for differently selected samples of AGN. A more comprehensive analysis combining panchromatic samples of X-ray, MIR, and radio-selected AGN would be crucial to test the widespread validity of this finding. This possible hint of “downsizing” (e.g. Cowie et al. 1996) links back to the known anti-hierarchical growth of galaxies over cosmic time with the most massive systems evolving earlier and faster than their lower mass counterparts (see also Bundy et al. 2006; Fontanot et al. 2009). In particular, this M*⋆* behaviour is expected if the most massive galaxies trigger higher radiative luminosity AGN activity earlier than less massive galaxies and then fade to lower radiative luminosity AGN at lower redshifts. The same qualitative argument is proposed in the evolution of AGN with the number density of powerful AGN ($L_{\rm x}>$1044 erg s*-1*) peaking earlier in cosmic time compared to lower luminosity AGN (e.g. Barger & Cowie 2005; Hasinger et al. 2005; Silverman et al. 2008; Ueda et al. 2014).

Different studies of AGN host galaxies have argued that AGN accretion preferentially occurs in gas-rich galaxies (Vito et al. 2014), and that the percentage of galaxies hosting X-ray AGN increases with infrared luminosity $L_{\rm IR}$ (e.g. Bongiorno et al. 2012; Santini et al. 2012). This is consistent with the increasing gas fraction observed in MS galaxies from low to high redshift (e.g. Daddi et al. 2010; Saintonge et al. 2012; Tacconi et al. 2013), and explained via the Schmidt–Kennicutt relation (Schmidt 1959; Kennicutt 1998). The increasing gas fraction from low to high redshift might explain the higher occurence of higher radiative luminosity AGN in massive galaxies ($M_{\star}\sim$1011 $M_{\odot}$) at higher redshift (z$\sim$2). Indeed, we showed that galaxies hosting HLAGN are mostly on the MS relation (see also Rosario et al. 2012), which implies a large availability of cold gas supplies, and possibly a more efficient fuelling mechanism of the central SMBH (i.e. with higher accretion rates), compared to the physical processes taking place in MLAGN.

According to this scenario, less star-forming galaxies are less likely to host an active SMBH. Interestingly, we found that most of MLAGN reside in weakly star-forming galaxies, which are typically located a factor of 2–3 below the MS. A plausible interpretation is that the difference in cold gas reservoirs leads HLAGN and MLAGN to be mainly powered by different accretion mechanisms. This raises the question of what triggers AGN activity in these two AGN populations. Shedding light on this issue requires a thorough investigation of the Eddington ratio distributions between these two AGN classes, which will be presented in a future work (I. Delvecchio et al., in prep.).

It is worth noting that HLAGN with radio excess show intermediate $M_{\star}$, [NUV-$r$] and SFR distributions between MLAGN and the rest of HLAGN, especially at z$<$1. Under the assumption that, in a stable phase, the AGN feedback occurs predominantly in either radiative or mechanical form, the population of HLAGN with radio excess might coincide with a transitional phase of AGN feedback.

According to semi-analytic models (e.g. Croton et al. 2006; Marulli et al. 2008, Hopkins et al. 2008), AGN feedback is one of the possible means to track the AGN host galaxies from the blue cloud of star-forming systems to the red sequence of passive galaxies, passing through a transition (often referred to as “green valley”), where the star formation is weaker but not yet stopped. According to this possible scenario, MLAGN and HLAGN with radio excess might represent intrinsically the same galaxies but that are observed at different stages of their AGN duty cycle, in which the energy produced via accretion onto the SMBH is emitted in either radiative or mechanical forms. The lower level of star formation in MLAGN might be a consequence of AGN-driven feedback, where the radio emission powered by the AGN could limit or hamper the galaxy star formation. This scenario is supported by studies of radio-selected AGN, where jet-induced feedback can strongly impact the molecular gas supplies of the host- galaxy (e. g. Feruglio et al. 2010; Morganti et al. 2013; Combes et al. 2013). In this context, it is possible that the population of HLAGN with radio excess probes a particular stage of the radio-mode feedback phase, where the molecular gas in the host galaxy is not yet depleted, and no evident impact in the integrated properties of the galaxy should be detectable during this transition (see Figs. 8, 9, and 10). This scenario is also supported by recent spectroscopic observations of powerful outflows detected in X-ray-MIR selected AGN, some of which show a significant radio excess (e.g. Perna et al. 2015; Lonsdale et al. 2015; Brusa et al. 2016). For these reasons, we stress that our sample of HLAGN with radio excess could be ideal to investigate the impact of AGN feedback, both radiative and mechanical, in a statistical sense and in a wide redshift range.

8 Conclusions

This work presents a multiwavelength analysis of radio-selected AGN host-galaxy properties out to z$\lesssim$6. Our sample consists of about 7 700 radio sources selected at 3 GHz in the COSMOS field, and cross-matched with optical-NIR counterparts. The exquisite photometry and redshifts available enabled us to use multiwavelength diagnostics to identify two main AGN populations in our sample: HLAGN (21%, out of which 30% also shows a $>$3$\sigma$ radio-excess) and MLAGN (17%). We analysed the average properties of their host galaxies at different cosmic epochs and summarise our main conclusions as follows:

We tested our source classification method against independent criteria used in recent radio-based studies (e.g. Del Moro et al. 2013; Bonzini et al. 2013; Padovani et al. 2015; N. Herrera Ruiz et al., in prep.), finding a good agreement and demonstrating the robustness of our method. 2. 2.

We provided a value-added catalogue containing the classification and the main physical properties discussed in this work for each radio source ($M_{\star}$, SFR, [NUV-$r$] colours, $L_{\rm 3\leavevmode\nobreak\ GHz}$ and $L_{\rm IR,SF}$). 3. 3.

Our HLAGN and MLAGN lie in different regions of the SFR–$M_{\star}$ plane, where the former are, on average, less massive and more star forming than the latter at various redshifts. We analysed in detail the observed distributions of various galaxy properties, finding significantly higher SFR and bluer [NUV-$r$] colours in HLAGN compared to MLAGN at all redshifts. Nevertheless, the $M_{\star}$ distribution is mainly populated by MLAGN at the highest $M_{\star}$ values ($M_{\star}$¿1011 $M_{\odot}$) at z$<$1, while the two AGN classes equally contribute to the highest $M_{\star}$ at z$\sim$1.5, and display a 6$\sigma$ reversal in the $M_{\star}$ behaviour at z$\sim$2. 4. 4.

Our results are consistent with radio emission predominantly arising from star formation in around 70% of HLAGN, while the remaining 30% shows a $\geq$3$\sigma$ radio excess that is likely attributable to AGN activity. The fractional AGN contribution to the radio emission in MLAGN is expected to be around 80–90%.

Overall, the differences in galaxy properties seen between these two AGN classes suggest that HLAGN and MLAGN samples trace two distinct galaxy populations in a wide range of redshifts. This might reflect the presence of two different driving mechanisms of AGN activity, which is possibly linked to the different availability of cold gas supplies in their hosts. In this scenario, the subsample of HLAGN with radio excess might coincide with a transitional phase during the AGN duty cycle, in which AGN activity occurs in both radiative and mechanical forms.

Acknowledgements.

The authors are grateful to the anonymous referee for his/her careful reading and useful comments, which improved the content of this manuscript. ID, VS, JD, MN, and OM acknowledge the European Union’s Seventh Framework programme under grant agreement 337595 (ERC Starting Grant, “CoSMass”). NB acknowledges the European Union’s Seventh Framework programme under grant agreement 333654 (CIG, ‘AGN feedback’). CL is funded by a Discovery Early Career Researcher Award (DE150100618). Parts of this research were conducted by the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020. MB and PC acknowledge supports from the PRIN-INAF 2014. MA acknowledges partial support from FONDECYT through grant 1140099. AK acknowledges support by the Collaborative Research Council 956, sub-project A1, funded by the Deutsche Forschungsgemeinschaft (DFG). This work was partially supported by NASA Chandra grant number GO3-14150C and GO3-14150B (FC, SM). ID is grateful to B. Magnelli for useful suggestions and to Takamitsu Miyaji for his help with CSTACK.

Appendix A Results from Herschel stacking

In Table 4, we list the fluxes and corresponding 1$\sigma$ uncertainties obtained by means of Herschel stacking. The results reported below refer to the sample of 2 203 radio sources without ($>$3$\sigma$) detection in any Herschel band, and not classified as HLAGN. See Sect. 4.4.3 for a more detailed description.

Appendix B Value-added 3 GHz radio catalogue

For guidance, we show 20 lines of the value-added catalogue used in this work in Table 5, following the same format introduced in Sect. 5.