© ESO, 2017

1. Introduction

Understanding massive galaxy formation is one of the major quests of modern astrophysics, and optimally studied via synergy between theory (magnetohydrodynamic and semi-analytic simulations) and (panchromatic) observations. Our current understanding of massive galaxy formation requires energetic outflows from radio-luminous active galactic nuclei (AGN) to suppress stellar growth in the most massive galaxies (e.g., Benson et al. 2003; Bower et al. 2006; Croton et al. 2006, 2016). This is mostly based on the advance of semi-analytic models in the last decades.

Semi-analytic models are used to study the properties of galaxies formed in cosmological models within dark matter halos by describing ongoing physical processes, such as gas cooling, star formation, merger rates, and various types of feedback. One of their major past challenges was related to substantial overprediction of the quantity of the most massive galaxies in the Universe (see, e.g., Steinmetz 1997). This was solved by introducing a process, dubbed radio-mode AGN feedback (e.g., Croton et al. 2006).

By now radio-mode AGN feedback has become a standard and key ingredient in semi-analytic models that enables reproduction of the observed galaxy properties well (Croton et al. 2006, 2016; Bower et al. 2006; Sijacki et al. 2007). The feedback is assumed to be related to radio AGN outflows as the main source that heats the halo gas surrounding a massive galaxy, thereby quenching its star formation and limiting growth, and thus avoids creating overly massive galaxies. Introducing such a feedback process was motivated by radio AGN that were observed to i) have massive galaxy hosts (e.g., Best et al. 2005; Kauffmann et al. 2008; Smolčić et al. 2009) and ii) interact with, i.e., heat the intra-cluster medium on large (group/cluster) scales, potentially solving the so-called cooling-flow problem in galaxy clusters (see Fabian 2012, for a review).

Initial implementations of radio-mode feedback were aimed to solve the cooling-flow problem and generalized based on simple phenomenological descriptions (Croton et al. 2006) or a sharp cutoff in cooling (Bower et al. 2006) that were both applied beyond a critical halo mass threshold. For example, in the Croton et al. (2006) model the source of heating was related to low-luminosity radio activity caused by hot gas accretion onto the central supermassive black hole (SMBH), once a static hot halo had formed around the host galaxy. The halo virial mass threshold beyond which a static hot halo forms was taken to be ≳2.5 × 10¹¹M_⊙ (see their Fig. 2).

A more complex cycle between gas cooling and radio-mode AGN heating was recently implemented within the updated Croton et al. (2006) model, i.e., the Semi-Analytic Galaxy Evolution (SAGE) model (Croton et al. 2016). In this model cooling and heating of the halo gas has been more directly coupled and the phenomenological treatment of radio-mode feedback has been put on more physical grounds by assuming accretion onto central SMBHs hosted by massive haloes of spherical, Bondi-Hoyle type (Bondi 1952), and has been scaled by a radio-mode efficiency parameter that modulates the strength of the accretion and subsequently the radio-mode feedback.

From an observational point of view, radio-mode AGN feedback, and particularly its cosmological relevance for massive galaxy formation, as suggested by the semi-analytic models, is difficult and challenging to test. While observations of large, spatially resolved radio galaxies inducing cavities (i.e., buoyantly rising bubbles) in the hot X-ray emitting intra-cluster medium (ICM) are possible for a small number of well-studied, nearby systems (e.g., Bîrzan et al. 2004, 2008; O’Sullivan et al. 2011), at higher redshifts the observational data is limited as most of the systems are unresolved in both the radio and X-ray bands. This necessitates the use of scaling relations to convert the monochromatic radio luminosity to kinetic luminosity inferred either on the basis of well-studied nearby systems (Merloni & Heinz 2007; Bîrzan et al. 2008; Cavagnolo et al. 2010; O’Sullivan et al. 2011; Godfrey & Shabala 2016) and assumed to hold at high redshifts (Smolčić et al. 2009; La Franca et al. 2010; Pracy et al. 2016), or on theoretical grounds, including various assumptions drawn from well-studied and resolved radio galaxies (Willott et al. 2001; Daly et al. 2012; Godfrey & Shabala 2016), and again assumed to hold for populations of radio AGN that have been observed in deep radio surveys over a large redshift range (Best et al. 2006, 2014).

A commonly used approach to estimate the cosmic evolution of radio-mode feedback, i.e., that of the volume averaged kinetic luminosity density of radio AGN, is based on constraining the cosmic evolution of radio AGN luminosity functions, convolved with a scaling relation between the monochromatic and kinetic luminosities (Merloni & Heinz 2008; Cattaneo & Best 2009; Smolčić et al. 2009; 2015; Best et al. 2014; Pracy et al. 2016). Hence, the cosmic evolution of radio AGN has direct implications on constraining the redshift evolution of radio-mode AGN feedback.

Past studies have shown that radio AGN evolve differentially, consistent with cosmic downsizing; low radio-luminosity sources evolve less strongly than high radio-luminosity sources and the number density peak occurs at higher redshift for higher luminosity radio AGN (Dunlop & Peacock 1990; Willott et al. 2001; Waddington et al. 2001; Rigby et al. 2011). Consistent with the evolution of optically and X-ray selected quasars (e.g., Schmidt et al. 1995; Silverman et al. 2008; Brusa et al. 2009), high-luminosity radio AGN (L_{1.4 GHz}> 2 × 10²⁶ W Hz^-1) show a strong positive density evolution with redshift out to z ~ 2, beyond which their comoving volume density starts decreasing (Dunlop & Peacock 1990; Willott et al. 2001). Intermediate- and low-luminosity radio AGN (L_{1.4 GHz}~ 1.6 × 10²⁴ − 3 × 10²⁶ W Hz^-1) have been shown to evolve slower with a comoving volume density turnover occurring at a lower redshift (z ~ 1 − 1.5; Waddington et al. 2001; Clewley & Jarvis 2004; Sadler et al. 2007; Donoso et al. 2009; Smolčić et al. 2009).

The goal of recent studies has been to constrain the cosmic evolution of the low-luminosity radio AGN reaching toward the epoch of reionization (z ≲ 6; e.g., McAlpine et al. 2013; Smolčić et al. 2015; Padovani et al. 2015). This is becoming possible only now via deep panchromatic surveys, including radio data obtained with the recently upgraded radio facilities, such as the Karl. G. Jansky Very Large Array (VLA). In this context, the VLA-COSMOS 3 GHz Large Project (Smolčić et al. 2017b) is to date the deepest radio continuum survey over a relatively large, 2 square degree field. Combined with the rich COSMOS multiwavelength data (Scoville et al. 2007), the Large Project thus provides a unique data set to study the evolution of the faintest observable radio AGN since z ~ 5 to date and put the results in context of radio-mode AGN feedback as an assumed cosmologically relevant process for massive galaxy formation. This is the aim of the work presented here.

In Sect. 2 we describe the VLA-COSMOS 1.4 and 3 GHz projects and the 3 GHz selected radio AGN sample used here. In Sect. 3 we derive the 1.4 GHz rest-frame radio luminosity functions for our AGN and model their evolution out to z ~ 5. In Sect. 4 we consider the obtained results in the context of radio-mode AGN feedback, comparing them with the SAGE semi-analytic model, and in Sect. 5 we discuss various unknowns remaining in both, the observational and analytic approaches. We summarize in Sect. 6. Throughout the paper we adopt H₀ = 70, Ω_M = 0.3,Ω_Λ = 0.7. We define the radio spectral index, α, via S_ν ∝ ν^α, where S_ν is the flux density at frequency ν. We use a Chabrier (2003) initial mass function.

2. Data and radio AGN sample

For our analysis we employ the VLA-COSMOS 3 GHz Large Project (Smolčić et al. 2017b) and the VLA-COSMOS 1.4 GHz Large and Deep Projects (Schinnerer et al. 2004, 2007, 2010). The VLA-COSMOS 3 GHz Large Project provides radio data at 10 cm wavelength within the COSMOS 2 square degree field down to an average (1σ) sensitivity of 2.3 μJy/beam over a 0.75″ resolution element. The source catalog lists 10 830 sources with signal-to-noise ratios (S/N) ≥ 5 (see Smolčić et al. 2017b, for details). Within the VLA-COSMOS 1.4 GHz Large Project the 2 square degree COSMOS field has been observed with the VLA down to a uniform root mean square (rms) of ~10–15 μJy/beam at a resolution of 1.5″. The VLA-COSMOS Deep Project has added further VLA 1.4 GHz observations to the inner square degree yielding an rms in this area of ~7 (12) μJy/beam over a 1.5″ (2.5″) resolution element. The joint catalog lists 2864 sources with S/N ≥ 5 at 1.5″ and/or 2.5″ resolution in the Large and/or Deep Projects. The number of sources detected in both the 3 and 1.4 GHz catalogs is 2530 (see Smolčić et al. 2017b).

Multiwavelength counterparts of the VLA-COSMOS 3 GHz Large Project sources have been identified by Smolčić et al. (2017a). For the cross-correlation they used the i) most up-to-date COSMOS photometric redshift catalog (COSMOS2015 hereafter; see Laigle et al. 2016 for details); ii) i-band selected catalog (Capak et al. 2007); and iii) 3.6 μm selected catalog (Sanders et al. 2007). The latter two were used to supplement potentially missed counterparts due to a nondetection in the COSMOS2015 catalog. To assure the most reliable photometric redshifts and a clean selection function here we used the most secure counterparts identified in the COSMOS2015 catalog within an area of 1.77 square degrees, uncontaminated by bright stars or saturated objects (see also Delvecchio et al. 2017; Delhaize et al. 2017; and Novak et al. 2017). Within this area 92% (8,035/8696) of the 3 GHz radio sources were associated with (COSMOS2015, i-band or IRAC) counterparts; 96% (7729) of these have counterparts drawn from the COSMOS2015 catalog and, thus, are the most precise photometric redshifts available (see Laigle et al. 2016; see also Fig. 4 in Smolčić et al. 2017a).

Using the full COSMOS (ultraviolet to millimeter) multiwavelength data set a three-component SED fitting procedure was performed for each galaxy taking into account the energy balance between the ultraviolet (UV) and infrared (IR) emissions of the galaxy and the contribution from a potential AGN component (see Delvecchio et al. 2017, for details; see also Smolčić et al. 2017a). These fits yielded an estimate of the total IR luminosity that was associated with the best-fit star-forming galaxy template (i.e., the AGN contribution was subtracted, if present), which was then converted to a star formation rate via the Kennicutt (1998) conversion, assuming a Chabrier (2003) initial mass function. To identify galaxies whose radio emission predominantly arises from the AGN, galaxies with radio excess were identified via a redshift dependent (3σ) threshold in the distribution of the ratio of the logarithms of radio luminosity (which is directly proportional to star-formation rate; e.g., Condon 1992) over the IR-derived star formation rate (see Fig. 5 in Smolčić et al. 2017a).

We define our radio AGN sample as the sample of 3 GHz sources with COSMOS2015 counterparts, with photometric or spectroscopic redshifts, detected in the near-infrared (NIR), which exhibit radio emission in (>3σ) excess of that expected from their hosts’ (IR-based) star formation rates. The sample contains 1814 sources and the chosen criterion assures that at least 80% of their radio emission is due to the AGN component. About 50% (899) of these are also detected at 1.4 GHz, allowing us to derive their spectral indices directly. In Fig. 1 we show the redshift distribution of our radio AGN sample, the 1.4 GHz rest-frame radio luminosity and stellar mass as a function of redshift. The sample reaches out to a redshift of z ~ 5 and comprises 1.4 GHz radio luminosities in the range of 10²² − 10²⁷ W Hz^-1 and typical stellar masses within ~3 × (10¹⁰ − 10¹¹)M_⊙. About 3% of the sources have estimated stellar masses that are lower than ~3 × 10⁹M_⊙. Although radio selection and X-ray selection preferentially select high-mass galaxies (e.g., Smolčić 2009; Smolčić et al. 2009), this fraction is very similar to the fraction of low-mass galaxies in samples of X-ray selected AGN (see, for example, Fig. 1 in Bongiorno et al. 2016, where the host galaxy mass function of the XMM-selected AGN in COSMOS is analyzed). The classification of these low-mass AGN as radio-excess sources seems to be statistically robust.

Fig. 1 Redshift distribution (top panel), 1.4 GHz rest-frame radio luminosity (middle panel), and stellar mass (bottom panel) as a function of redshift for our radio-excess AGN sample.

3. Radio AGN luminosity functions and their cosmic evolution

In this section we derive the AGN luminosity functions out to z ~ 5 (Sect. 3.1) and model the cosmic evolution of radio AGN in the COSMOS field (Sect. 3.2). We further present the cosmic evolution of their number and luminosity densities (Sect. 3.3) and a comparison with results from the literature (Sect. 3.4).

3.1. Radio AGN luminosity functions out to z ~ 5

We computed the volume densities (for a given redshift range) as described in detail by Novak et al. (2017). We followed the 1 /V_max procedure and corrected for a combined set of radio incompletenesses, including radio detection, noise, and resolution biases, as well as for the incompleteness of the counterpart catalog due to radio sources without assigned NIR counterparts (the latter being overall less than 10%). We refer to Sect. 3.1. in Novak et al. (2017) for a detailed description of the procedure (see also their Fig. 2).

The rest-frame 1.4 GHz luminosities¹ were computed using the observed-frame 3 GHz flux densities. If the source was also detected at 1.4 GHz the inferred spectral index was used in the computation. If the source was undetected at 1.4 GHz we assumed a spectral index of α = − 0.7, which corresponds to the median value derived for all radio sources detected at 3 GHz and also taking the limits (i.e., nondetections at 1.4 GHz) into account via survival analysis (see also Smolčić et al. 2017b). Spectral indices used in the further analysis are shown in Fig. 2. A trend toward steeper spectra can be seen with increasing redshift. Without deep observations at other radio frequencies we cannot state whether this trend is real (e.g., owing to spectral steepening of the radio spectrum at higher rest-frame frequencies, as sampled for our high-redshift sources) or a bias due to our flux limited observations.

Fig. 2 Spectral indices of our radio-excess AGN sample as a function of redshift (gray points). Median values with interquartile ranges for sources with measured 1.4–3 GHz spectral indices in different redshift bins are shown with blue circles. Mean values for all sources (i.e., also undetected at 1.4 GHz for which we manually set α = − 0.7) in same redshift bins are shown with red squares.

The 1.4 GHz radio luminosity functions (LFs) for our radio AGN, separated into nine redshift ranges out to z ~ 5, are shown in Fig. 3, and tabulated in Table 1. Redshift bins were chosen to be large enough to contain a statistically significant number of galaxies and to mitigate possible photometric uncertainties (i.e., sources falling into the wrong bin). We report LFs using the median luminosity in each luminosity bin, while the error bars show the width of the luminosity bin.

Fig. 3 Radio AGN 1.4 GHz rest-frame luminosity functions in nine redshift bins out to z ~ 5 for radio-excess AGN in the COSMOS field. The volume densities derived from our radio sample with the Vmax method are shown by the filled black dots. A 1σ confidence interval of the combined luminosity and density evolution in individual redshift bins is shown as the gray shaded area. Pure luminosity evolution and pure density evolution 2-parameter models are shown with blue and red lines, respectively (see text for details). The dotted curve in each panel indicates the local luminosity function (Eq. (2)). The vertical line in each panel shows the 5σ sensitivity limit (assuming a spectral index of α = − 0.7) at the high-redshift end of each redshift bin. Various results from the literature, as indicated in the legend, are also shown (each scaled to the median redshift in each bin using the evolution model as reported in the corresponding study).

3.2. Evolution of the radio AGN luminosity function

The cosmic evolution of an astrophysical population is usually parameterized by density and luminosity evolution of its local luminosity function as $\mathrm{\Phi }\mathrm{\left(}\mathit{L,z}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{{\mathit{\alpha }}_{\mathit{D}}}\mathrm{\times }{\mathrm{\Phi }}_{\mathrm{0}}\left[\frac{\mathit{L}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{{\mathit{\alpha }}_{\mathit{L}}}}\right]\mathit{,}$(1)where α_D and α_L are the characteristic density and luminosity evolution parameters, respectively, L is the luminosity, Φ(L,z) is the luminosity function at redshift z, and Φ₀ is the local luminosity function. The analytic form for the local radio AGN LF, adopted here, and also shown in Fig. 3, is taken from Mauch & Sadler (2007) and parametrized with two power laws ${\mathrm{\Phi }}_{\mathrm{0}}\mathrm{\left(}\mathit{L}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\Phi }}^{\mathrm{\ast }}}{\mathrm{(}{\mathit{L}}^{\mathrm{\ast }}\mathit{/}\mathit{L}{\mathrm{)}}^{\mathit{\alpha }}\mathrm{+}\mathrm{(}{\mathit{L}}^{\mathrm{\ast }}\mathit{/}\mathit{L}{\mathrm{)}}^{\mathit{\beta }}}\mathit{,}$(2)where the parameters are the normalization ${\mathrm{\Phi }}^{\mathrm{\ast }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{0.4}}{\mathrm{10}}^{-5.5}$ Mpc^-3 dex^-1 (scaled to the base of dlog L), the knee position L^∗ = 10^24.59 W Hz^-1, and the bright and faint end slopes α = − 1.27, and β = − 0.49, respectively. Mauch & Sadler (2007) derived their AGN LF using 2661 detections in the 6dFGS–NVSS field with a median redshift of med(z) = 0.073 and a span of six decades in luminosities. With such a sample they were able to constrain well both the faint and bright end of the local AGN LF. For consistency with other studies in the literature and to provide a broad overview, we here take the LF of the radio AGN population as derived by Mauch & Sadler (2007), and model its evolution as it is usually done (using Eq. (1), see also below). Also, as discussed in Sect. 1 a two-population model may be more appropriate for modeling the evolution of radio AGN. This is discussed further in Sect. 5.

In Fig. 4 we show the best-fit pure density (PDE; α_L = 0) and pure luminosity (PLE; α_D = 0) evolutions for each redshift bin (also tabulated in Table 2), which can be considered as the two extreme cases of evolution. For a conservative approach the outlying, lowest luminosity bins at z< 1.3 were ignored in the fitting process.

To fit a simple, continuous model to the data, we follow Novak et al. (2017), and add a redshift dependent term to the α_L, and α_D parameters in Eq. (1), and model the evolution of the local luminosity function using the following form: $\mathrm{\Phi }\mathrm{\left(}\mathit{L,z,}{\mathit{\alpha }}_{\mathit{L}}\mathit{,}{\mathit{\beta }}_{\mathit{L}}\mathit{,}{\mathit{\alpha }}_{\mathit{D}}\mathit{,}{\mathit{\beta }}_{\mathit{D}}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{{\mathit{\alpha }}_{\mathit{D}}\mathrm{+}\mathit{z}\mathrm{·}{\mathit{\beta }}_{\mathit{D}}}\mathrm{\times }{\mathrm{\Phi }}_{\mathrm{0}}\left[\frac{\mathit{L}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{{\mathit{\alpha }}_{\mathit{L}}\mathrm{+}\mathit{z}\mathrm{·}{\mathit{\beta }}_{\mathit{L}}}}\right]\mathit{,}$(3)where α_L, β_L, α_D, and β_D are the four free parameters. For pure luminosity evolution (α_D = β_D = 0) the χ² minimization procedure yields best-fit parameter values of α_L = 2.88 ± 0.82, β_L = − 0.84 ± 0.34, while for pure density evolution (α_L = β_L = 0) the best-fit parameters are α_D = 2.00 ± 0.18, β_D = − 0.60 ± 0.14. Fitting for all four parameters simultaneously yields a strong degeneracy between the parameters and mainly differs from the two two-parameter fits at the high-luminosity end in high-redshift bins (z> 2.1). As discussed in more detail in Sect. 5, the volume densities in these particular bins are the most sensitive to the assumption of a simple synchrotron power law for the K correction to rest-frame 1.4 GHz luminosity from the observed 3 GHz flux density. Hence, as the typical AGN spectrum at radio frequencies in deep surveys, such as the VLA-COSMOS 3 GHz survey, to date is not constrained, and further, given the degeneracies inherent to the four-parameter fit, we hereafter adopt the simple continuous two-parameter models as a representation of the evolution of the VLA-COSMOS 3 GHz (radio-excess) AGN. Also, simple, pure luminosity (density) evolution models are commonly used in the literature (e.g., Sadler et al. 2007; Smolčić et al. 2009; Donoso et al. 2009; McAlpine et al. 2013; Padovani et al. 2015).

Fig. 4 Best-fit parameters for the evolution of the local LF with redshift. The best-fit pure luminosity (density) evolution in a given redshift bin is shown by the black points in the top (bottom) panel. The blue (red) line in the top (bottom) panel shows the simple, continuous 2-parameter pure luminosity (density) evolution model (see text for details). Gray shaded areas in both panels correspond to the 68% confidence interval for a combined luminosity and density evolution. The large uncertainty in the combined fit is due to parameter degeneracy.

Fig. 5 Redshift evolution of the number density (top panel), and luminosity density (bottom panel) for our radio AGN (for models as indicated in the panels; see text for details). The number density was integrated with a lower boundary of L1.4 GHz= 1022 W Hz-1. Decreasing the boundary would systematically shift the number densities to higher values (for example, by a factor of 3–5 if L1.4 GHz= 1021 W Hz-1 were used). The luminosity density integral has a much weaker dependence on the integration limits (if it is broad enough) as the bulk of the integral is constrained close to (L∗,Φ∗).

3.3. Cosmic evolution of the number and radio luminosity densities

In Fig. 5 we show the redshift evolution of the number and luminosity densities using our two-parameter evolution models and the best-fit pure luminosity and density evolutions in each redshift bin. The number density in a given redshift bin was obtained by integrating the corresponding luminosity function over the logarithm of luminosity, ^∫Φ(L_{1.4 GHz}) d(log L_{1.4 GHz}), and the luminosity density by integrating the product of luminosity and volume density over the logarithm of luminosity, i.e., ${}^{\mathrm{\int }}\mathit{L}_{\mathrm{1.4}\mathrm{GHz}}\mathrm{\times }\mathrm{\Phi }\mathrm{\left(}{\mathit{L}}_{\mathrm{1.4}\mathrm{GHz}}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathrm{\left(}\log {\mathit{L}}_{\mathrm{1.4}\mathrm{GHz}}\mathrm{\right)}$.

The number density of our radio AGN increases by a factor of ~2–3 from z = 0 to z ~ 1.5, beyond which it decreases, reaching a number density equivalent to that in the local Universe at z ~ 3.5 and further decreasing beyond this redshift. The luminosity density shows a similar behavior. It rises by a factor of ~2–4 from z = 0 to z ~ 1.5, beyond which it decreases, at z ~ 5 reaching a value about an order of magnitude lower than the luminosity density derived at z = 0.

3.4. Comparison with the literature

Based on the 2SLAQ Luminous Red Galaxy Survey containing ~400 galaxies in a volume-limited sample at 0.4 <z< 0.7,Sadler et al. (2007) have found that their radio AGN (L_{1.4 GHz}≈ 10²⁴ − 10²⁷ W Hz^-1) undergo significant evolution since z ~ 0.7, parametrized with a pure luminosity evolution parameter α_L = 2.0 ± 0.3. Donoso et al. (2009) found fully consistent results, but with considerably smaller error bars as they used a sample of over 14 000 radio AGN. These derivations are in very good agreement with the luminosity function derived here (see 0.4 <z< 0.7 bin in Fig. 3) and with the pure luminosity evolution we find for the VLA-COSMOS radio-excess AGN, detected at 3 GHz (see top panel of Fig. 4).

McAlpine et al. (2013) studied the evolution of faint radio AGN in the VIDEO-XMM3 field out to z ~ 2.5. By fitting a combined evolution of star-forming and AGN galaxies, they found a slightly weaker evolution than that inferred by Sadler et al. (2007, see Table 4 in McAlpine et al. 2013. The authors argue that this is because the evolution of the VIDEO-XMM3 radio AGN is primarily driven by the higher redshift range (0.9 <z< 2.5), not constrained by the 2SLAQ survey. The pure luminosity evolution inferred by McAlpine et al. (2013, α_L = 1.2 ± 0.2 is consistent with the average α_L derived here in the equivalent (0.9 <z< 2.5) redshift range (for example, for z = 1.75, α_L + z·β_L = 1.4 ± 1.0 for our pure luminosity two-parameter fit; see top panel of Fig. 4). In each redshift range the VIDEO-XMM3-based volume densitites of their AGN (shown in Fig. 3), which were identified during the photometric-redshift estimation process via template fitting (see their Sect. 6.2.1.), are slightly below those derived here, particularly at the high-luminosity end in each redshift bin (see Fig. 3). This underestimate of the number density of AGN at the high-luminosity end, and especially at z> 1, however, has already been observed and reported by McAlpine et al. (2013, see their Fig. 8).

Based on the VLA-COSMOS 1.4 GHz survey Smolčić et al. (2009) have derived luminosity functions for their rest-frame color selected AGN out to z = 1.3, which is also shown in Fig. 3. Overall, there is good agreement between the derivations based on the (shallower) 1.4 GHz survey, and the 3 GHz survey used here. Smolčić et al. (2009) inferred a pure luminosity (density) evolution out to z = 1.3 for their AGN of α_L = 0.8 ± 0.1 (α_D = 1.1 ± 0.1). While their pure density evolution is in agreement with that derived here, the pure luminosity evolution is here derived to be slightly higher; i.e., α_L(z = 1.3) = 1.6 ± 0.3, and α_L + z·β_L = 1.8 ± 0.9 for z = 1.3, and our pure luminosity two-parameter fit. The differences are due to a combination of factors, such as i) the ~3 times increased sensitivity of the VLA-COSMOS 3 GHz, compared to the 1.4 GHz Large projects; ii) the different local luminosity functions used in the two studies; while the Mauch & Sadler (2007) local luminosity function is used here, the Sadler et al. (2002) local luminosity function was used by Smolčić et al. (2009); and iii) differently selected AGN. Smolčić et al. (2009) selected the radio AGN by requiring host galaxies with red rest-frame colors (Smolčić et al. 2008) that mimic the spectroscopically based identification (Baldwin et al. 1981) commonly used in the local Universe, while we use the radio-excess criterion to identify our AGN, which includes both host galaxies with red rest-frame colors, but also AGN selected via X-ray, and mid-infrared (MIR) criteria within blue host galaxies (see Sect. 5, and Fig. 7 below; see also Delvecchio et al. 2017).

Bonzini et al. (2013) identified a sample of (z ≤ 4) radio-loud AGN in the E-CDFS field by requiring that at a given redshift, they lie below the 2σ deviation from the average observed 24 μm-to-1.4 GHz flux ratio obtained using the M82 galaxy template. This criterion is similar to that used here to select radio-excess AGN (see Sect. 4.3.2. in Delvecchio et al. 2017), and thus in Fig. 3 we also compare the luminosity functions derived here with those derived by Padovani et al. (2015) for the radio-loud AGN in the E-CDFS field. Overall, the E-CDFS-based luminosity functions agree well with those derived here out to z = 1.3. Beyond this redshift a higher degree of discrepancy is observed. This is likely related to small number statistics given the size of the E-CDFS survey (0.3 square degrees) resulting in most of the E-CDFS volume density values at z> 1.3 constrained by bins containing fewer than three sources (cf. Table 5 in Padovani et al. 2015 and Table 1 here).

4. Radio-mode feedback considerations

In this section we derive the cosmic evolution of the radio AGN kinetic luminosity density (Sect. 4.1), and compare this value with the SAGE semi-analytic model (Sect. 4.2).

4.1. Cosmic evolution of the radio AGN kinetic luminosity density

The AGN studied here have, by construction, an observed excess of radio emission relative to that expected from star formation processes within their host galaxies. Thus, the origin of the observed 3 GHz radio emission in these sources is related to the AGN itself, i.e., to the energy released through mass accretion onto the central SMBH. This energy can be efficiently radiated away and/or channeled in kinetic form via collimated jets expanding through and beyond the host galaxy and is observable at radio wavelengths. The observed radio emission in our sources (albeit mostly unresolved or barely resolved at 0.75″ resolution) can thus be predominantly attributed to the synchrotron emission of relativistic particles within the jet structures of the sources (i.e., core, jets, and lobes). However, the observed radio emission amounts to only a fraction of the total kinetic luminosity² of the jets, while a much larger fraction is stored as the internal energy of the jets and lobes and is lost to the environment via the work carried out by the expanding radio jets (e.g., Willott et al. 1999). The latter is of particular interest as the energy (per unit time) deposited into the surroundings and dissipated can be directly linked to the heating of the surrounding medium generated by the jets of such AGN (see McNamara & Nulsen 2007, 2012, for reviews).

A simple scaling relation between the monochromatic radio and kinetic luminosities has been long sought after. In Appendix A we give a detailed overview of such scaling relations, commonly used in the literature (Willott et al. 1999; Bîrzan et al. 2004; 2008; Merloni & Heinz 2007; Cavagnolo et al. 2010; O’Sullivan et al. 2011; Daly et al. 2012; Godfrey & Shabala 2016). We hereafter adopt the relation derived by Willott et al. (1999) at 151 MHz rest-frame frequency and converted to 1.4 GHz rest-frame luminosity (Heckman & Best 2014), given as $\log {\mathit{L}}_{\mathrm{kin}}\mathrm{\left(}{\mathit{L}}_{{\mathrm{L}}_{\mathrm{1.4}\mathrm{GHz}}}\mathrm{\right)}\mathrm{=}\mathrm{0.86}\mathrm{·}\log {\mathit{L}}_{{\mathrm{L}}_{\mathrm{1.4}\mathrm{GHz}}}\mathrm{+}\mathrm{14.08}\mathrm{+}\mathrm{1.5}\log {\mathit{f}}_{\mathrm{W}}\mathit{,}$(4)where L_kin is the kinetic luminosity in units of W, L_{L1.4 GHz} is the 1.4 GHz rest-frame luminosity in units of W/Hz, and f_W is an uncertainty parameter (see below).

Willott et al. (1999) estimated the kinetic energy using minimum energy arguments, i.e., computing the minimum energy stored in the lobes to produce the observed synchrotron luminosity, considering the age of the source and efficiency of conversion of the kinetic luminosity into the internal energy of the observed synchrotron emission. These authors have folded all uncertainty factors in the calculation, such as departure from minimum energy conditions, uncertainty in the energy in nonradiating particles, and the composition of the jet, into one parameter, f_W, estimated to lie in the range of f_W ≈ 1 − 20. For f_W = 15, the normalization is close to that of the jet kinetic luminosities computed through X-ray observations of galaxy clusters with cavities induced in the hot, X-ray emitting intracluster gas by the radio jets and lobes (Bîrzan et al. 2004; 2008; Merloni & Heinz 2007; Cavagnolo et al. 2010; O’Sullivan et al. 2011), and for f_W = 4 the normalization matches the relation derived for powerful radio galaxies using equations of strong shock physics (see Daly et al. 2012, and references therein). We, however, also consider the full range of the uncertainty parameter (f_W = 1 − 20; see below).

The kinetic luminosity density at a given redshift is computed as the integral of the kinetic luminosity density over (the logarithm of the) monochromatic 1.4 GHz radio luminosity, ${\mathrm{\Omega }}_{\mathrm{kin}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{\int }{\mathit{L}}_{\mathrm{kin}}\mathrm{\left(}{\mathit{L}}_{{\mathrm{L}}_{\mathrm{1.4}\mathrm{GHz}}}\mathrm{\right)}\mathrm{\times }\mathrm{\Phi }\mathrm{\left(}{\mathit{L}}_{{\mathrm{L}}_{\mathrm{1.4}\mathrm{GHz}}}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathrm{(}\log {\mathit{L}}_{\mathrm{1.4}\mathrm{GHz}}\mathrm{)}\mathit{.}$(5)In Fig. 6 we show the cosmic evolution of the kinetic luminosity density derived using our two-parameter continuous evolution model (Eq. (3)) and the scaling relation between monochromatic radio luminosity and kinetic luminosity of Willott et al. (1999, Eq. (4 with uncertainty parameters f_W = 4 and 15, and the range encompassed by the extreme values of f_W (= 1,20). The various f_W values only systematically shift the derived volume-averaged kinetic luminosity density as a function of redshift. Furthermore, the derivation for f_W = 1 can be considered as a robust lower limit, as any other scaling relation available in the literature would have resulted in systematically higher values (see Appendix A for details). Under the assumptions made, the kinetic luminosity density rises by about a factor of three from z = 0 to z ~ 1.5 and decreases thereafter by close to two orders of magnitude by z ~ 5. This holds for both of our (pure luminosity and pure density) two-parameter models, as also illustrated in Fig. 6.

Fig. 6 Cosmic evolution of the kinetic luminosity density for our radio AGN. The kinetic luminosity density was derived using Eq. (5), and the Willott et al. (1999) relation between monochromatic (1.4 GHz rest-frame) radio luminosity and kinetic luminosity (Eq. (4)) with fW = 4 (pure luminosity evolution, PLE, model: blue full curve and blue shaded area showing the 1σ uncertainty), fW = 15 (purple full curve, and purple-hatched area showing the 1σ uncertainty), and fW = 1, and 20, with folded 1σ fitting errors (lower and upper gray dashed curves, respectively). Also shown is the kinetic luminosity density assumed in the SAGE semi-analytic cosmological model (green dash-dotted curve; Croton et al. 2016). For comparison, the evolution of the kinetic luminosity density for our radio AGN for the pure density evolution, PDE, model using only fW = 4 (Eq. (A.1)) is also shown with the red full curve and red hatched area. To avoid overcrowding the panel, the PDE result for fW = 15 is not shown, but it would be only systematically higher compared to the fW = 4 PDE model result and coincident with that for PLE with the corresponding fW = 15 value. See also Fig. A.2 for results based on other scaling relations.

4.2. Comparison with the SAGE semi-analytic model

In Fig. 6 we compare the redshift evolution of the COSMOS radio AGN kinetic luminosity density with that from the SAGE model (Croton et al. 2016). The SAGE model is a significant update of the semi-analytic model by Croton et al. (2006), including implementation of a more complex cycle between gas cooling and the radio-mode AGN heating. It is assumed that hot gas accretes onto the central black hole following the Bondi-Hoyle accretion (Bondi 1952), but including a, so-called, radio-mode efficiency parameter (κ_R; see Eq. (16) in Croton et al. 2016). This parameter is used to modulate the strength of black hole accretion (ṁ), assumed to be related to the luminosity of the black hole in radio mode in the standard way, as L = ηṁc², where η = 0.1 is the standard efficiency of mass to energy conversion, and c is the speed of light. In the SAGE model this luminosity is taken as the source of heating that offsets the energy losses of the cooling gas, such that the heating rate is simply the ratio between the luminosity and the specific (per unit mass) energy of gas in the hot halo (see Eq. (17) and Sect. 9.1. in Croton et al. 2016). The accretion luminosity assumed in the SAGE model can be taken as an equivalent to the, here derived, jet kinetic luminosity assuming that the bulk of the accretion energy is channeled in kinetic (rather than radiative) form. Hence, in Fig. 6 we can compare the redshift evolution of the COSMOS radio AGN kinetic luminosity density with the volume-averaged accretion luminosity responsible for radio-mode feedback in the SAGE model, which is equivalent to that obtained by Croton et al. (2006, their Fig. 3), but scaled by the radio-mode efficiency parameter (κ_R = 0.08; see Sect. 9.1.1. in Croton et al. 2016).

As seen in Fig. 6 at z ≳ 1, we find a remarkably similar slope of the redshift evolution of the COSMOS radio AGN kinetic luminosity density and the radio-mode accretion luminosity used in the SAGE model with the absolute values on the same order for f_W = 4 (the uncertainty parameter in the radio luminosity – kinetic luminosity relation; Eq. (4)). At z< 1 we find a steeper evolution than that in the SAGE model, for which our z = 0 value is a factor of 3–4 lower (f_W = 4). For f_W = 15 the kinetic luminosity density is at every redshift systematically higher than the radio-mode accretion luminosity in the model, while for the robust lower limit value (f_W = 1) the observationally based kinetic luminosity density is at each redshift about an order of magnitude lower than that in the SAGE model. Overall, for the most likely values of the normalization of the radio luminosity – kinetic luminosity relation (f_W = 4,15; see previous section, and Appendix A), the redshift evolution of the COSMOS radio AGN kinetic luminosity density is either on the same scale as that in the model or is systematically higher by about an order of magnitude or even higher if other scaling relations are used (see Fig. A.2). This would suggest that, in both cases, the energy deposited by faint radio AGN into their environment may be sufficient to offset the cooling energy losses, as postulated in the model. However, there is still a non-negligible number of simplifications and unknowns inherent in both the semi-analytic models, and the observational results, as discussed in the next section.

5. Discussing the unknowns

Observational studies of radio-selected AGN find two radio-luminous AGN populations with distinct host galaxy and AGN properties (Smolčić 2009; Smolčić et al. 2009, 2015; Smolčić 2016; Hardcastle et al. 2007; Buttiglione et al. 2010; Heckman & Best 2014; Bonzini et al. 2015; Padovani et al. 2015; Padovani 2016; Tadhunter 2016). One population is consistent with the standard, unified model AGN picture, in which the accretion occurs in a radiatively efficient manner at high Eddington rates (1 − 10% ≲ λ_Edd ≲ 100%). Evidence has been presented that this class is fueled by the cold intragalactic medium (IGM) phase and that it is not too likely to launch collimated jets (observable at radio wavelengths). The second population, however, exhibits radiatively inefficient accretion related to low Eddington ratios (λ_Edd ≲ 1 − 10%) and may be fueled by the hot phase of the IGM. It has been shown that this population is highly efficient in collimated jet production. The observed difference in Eddington ratio between the two populations can be linked to the switch between the standard accretion flow model, i.e., radiatively efficient, geometrically thin (but optically thick) disk accretion flow (Shakura & Sunyaev 1973), and a radiatively inefficient, geometrically thick (but optically thin) accretion flow (Esin 1997; Narayan et al. 1998), occurring at accretion rates below a certain Eddington ratio (1–10%; Rees et al. 1982; Narayan & Yi 1994; Meier 2002; Fanidakis et al. 2011). Further differences have been found demonstrating that the population of radiatively efficient accretors dominates the radio-AGN number densitites at the bright end (e.g., L_{1.4 GHz}≳ 10²⁶ W Hz^-1 for z< 0.3; Pracy et al. 2016) and that this population evolves more rapidly with cosmic time (e.g., Willott et al. 2001; Pracy et al. 2016).

The radio AGN studied here have been selected based on an excess of their radio emission relative to that expected from the (IR-based) star formation rates in their host galaxies, thus assuring that ≳80% of the radio emission arises from the AGN core, jet, and/or lobe component. In Fig. 7 we plot the absolute and fractional contributions of X-ray-, MIR-selected AGN, and the remaining AGN (not identified via X-ray or MIR emission), those hosted by red, quiescent galaxies, and those hosted by galaxies with blue or green rest-frame colors, implying substantial star formation activity in the hosts (see Smolčić et al. 2017a, for details). The X-ray and MIR regimes provide an efficient approach to identify radiatively efficient AGN, while red, quiescent host galaxies of radio AGN are shown to contain AGN with systematically lower radiative AGN luminosities (see, e.g., Fig. 7 in Smolčić et al. 2017a); further, this selection can be used to trace radiatively inefficient AGN (at least in the local Universe; e.g., Smolčić 2009).

As seen in Fig. 7, at z ≲ 1 our radio-excess AGN are composed of similar fractions of i) red, quiescent galaxies and ii) X-ray, MIR, and those hosted by star-forming galaxies, but not identified via X-ray or MIR emission. Beyond z = 1 the fraction of red, quiescent galaxies decreases steeply to a minimal fraction by z ~ 2. Hence, we can conclude that AGN with observational signatures of radiatively efficient accretion flows comprise a non-negligible fraction of our radio-excess AGN (30–40% X-ray- and MIR-selected AGN). If only a fraction of the radio-excess AGN hosted by blue or green, star-forming galaxies also contained radiatively efficient AGN, which could be the case given the expected cold gas supply in such, actively star-forming galaxies (e.g., Vito et al. 2014), the fraction of radiatively efficient AGN in our radio-excess sample would rise even further. This implies that our radio-excess AGN are likely a mix of radiatively efficient and inefficient black hole accretion flows, shown in other studies to evolve differently with cosmic time (e.g., Willott et al. 2001; Best et al. 2014; Pracy et al. 2016).

Fig. 7 Absolute (top panel) and fractional (bottom panel) contributions of various subpopulations (as indicated in the bottom panel) to the full radio-excess AGN sample (black curve) as a function of redshift.

The volume densities derived here for our radio AGN are, at the high-luminosity end in every redshift range, somewhat higher than our best-fit pure luminosity/density two-parameter models (see Fig. 3). This could potentially be attributed to the rise given the faster evolution of the high-luminosity radio AGN that are consistent with the radiatively efficient accretion flow population. However, these particular, high-luminosity bins at each redshift are the most sensitive to the assumption of a simple synchrotron power law for the conversion from the observed 3 GHz flux density to rest-frame 1.4 GHz luminosity. For example, assuming a spectral index of α = − 0.7 for all our AGN instead of using the observed value, if available, which steepens toward high redshift (see Fig. 2), this would have no effect on the low-luminosity volume densities in the high-redshift bins, but the high-luminosity bin values would decrease to be consistent with the pure luminosity (density) two-parameter models.

If the typical radio AGN spectrum were to steepen toward the high-frequency end owing to synchrotron energy losses (e.g., Miley 1980), as sampled by the observed frequency of 3 GHz at z> 2 (corresponding to rest-frame frequencies of >9 GHz) then, for example, a broken power law, rather than a simple, single power-law assumption would be more appropriate for the K correction. On the other hand, in this case, the simple, single power-law assumption of α = − 0.7 could result in a more correct 1.4 GHz rest-frame luminosity value than the directly observed higher spectral index values.

Although our radio-excess AGN are likely a mix of radiatively efficient and inefficient black hole accretion flows it could be assumed that, if heating from radio outflows is indeed occurring in these systems, this heating operates in a similar fashion in the various AGN within our sample, as postulated in the SAGE model. In the SAGE model it is taken that cooling flows develop in any halo with mass ≳2.5 × 10¹¹M_⊙, and further, that Bondi-Hoyle accretion of the hot (subdominant) component of the cooling gas, which fills the space between the (dominant) cold cloud component, is responsible for the luminosity of the SMBH in radio mode. Taking the stellar mass – halo mass relation of Behroozi et al. (2010), the halo mass threshold beyond which radio-mode feedback occurs can be converted to a stellar mass of ≳5 × 10⁹M_⊙. The stellar masses of our radio AGN satisfy this criterion (see Fig. 1), and can, hence, be taken as a representative population of that assumed in the SAGE model to be responsible for radio-mode feedback. In that respect, the similarity between the here derived redshift evolution of the kinetic luminosity density and that assumed in the SAGE semi-analytic model (at least at z> 1) is astonishing. As a further step a more complex treatment of radio AGN in semi-analytic models, taking into account the two dominant, thin- and thick-disk, accretion flows and the spins of SMBHs, would be desirable (e.g., Fanidakis et al. 2011). However, from an observational perspective likely the largest source of uncertainty in studies of radio-mode feedback is the uncertainty of the monochromatic radio luminosity to kinetic luminosity conversion. The relation used here contains an uncertainty factor on the order of 20 (f_W in Eq. (4); Willott et al. 1999) and various relations commonly used in the literature are summarized in Appendix A. The kinetic luminosity has been demonstrated to be not only a function of monochromatic radio luminosity, but also of other parameters such as the synchrotron electron age of the AGN and its environment (Bîrzan et al. 2008; Hardcastle & Krause 2014; Kapinska et al. 2015). A better constraint of the relation as a function of redshift and other physical parameters is still awaited. In this context, deep radio surveys, observed in a large enough range of radio frequency, allowing us to compute synchrotron ages, could provide an important step toward better understanding of the relevance of radio-mode AGN feedback in massive galaxy formation.

6. Summary and conclusions

We have used a sample of 1814 radio AGN selected from the VLA-COSMOS 3 GHz Large Project, and with NIR-detected (COSMOS2015 catalog) counterparts with photometric or spectroscopic redshifts. These were selected to have radio emission in >3σ excess of that expected from the IR-based star formation rates of their hosts. Such a criterion assures that at least 80% of their radio emission is due to the AGN component. This sample of AGN reaches out to a redshift of z< 6, comprises 1.4 GHz rest-frame radio luminosities in the range of 10²² − 10²⁷ W Hz^-1, and the typical stellar masses of the AGN are within ~3 × (10¹⁰ − 10¹¹)M_⊙.

Using the 1 /V_max method we derived the luminosity functions of the radio AGN out to z ~ 5. We further constrained the evolution of this population with continuous models of pure density and pure luminosity evolutions finding best-fit parameters of Φ^∗ ∝ (1 + z)^{(2.00 ± 0.18) − (0.60 ± 0.14)z} and L^∗ ∝ (1 + z)^{(2.88 ± 0.82) − (0.84 ± 0.34)z}, respectively. We find a turnover in number and luminosity densities of the population at z ≈ 1.5.

The 1.4 GHz luminosity was converted to kinetic luminosity via an analytically motivated relation. Taking into account the full range of uncertainty, we derived the cosmic evolution of the kinetic luminosity density provided by our AGN, which we compare to the radio-mode AGN feedback assumed in the SAGE model. We find that the kinetic luminosity exerted by our radio AGN may be high enough at each cosmic epoch since z ~ 5 to balance the radiative cooling of the hot gas, as assumed in the model. However, although our findings support the idea of radio-mode AGN feedback as a cosmologically relevant process in massive galaxy formation, many simplifications in both the observational and semi-analytic approaches still remain and need to be resolved before robust conclusions can be reached.

Acknowledgments

This research was funded by the European Union’s Seventh Frame-work programs under grant agreements 333654 (CIG, “AGN feedback”). V.S., M.N., I.D., and L.C. acknowledge support from the European Union’s Seventh Frame-work program under grant agreement 337595 (ERC Starting Grant, “CoSMass”). M.B. acknowledges support from the PRIN-INAF 2014. E.S. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 694343).

References

Appendix A: Scaling relations between monochromatic radio luminosity and kinetic luminosity

Here we summarize the scaling relations between monochromatic radio luminosity and kinetic jet luminosity commonly used in the literature and illustrated in Fig. A.1.

A commonly used technique to estimate the kinetic luminosity is via X-ray observations of galaxy clusters hosting radio galaxies, whose jets and lobes, which succeeded in propagating beyond their host galaxies and into the intracluster gas, induce cavities in the hot X-ray-emitting gas and inflate buoyantly rising bubbles (see, e.g., Bîrzan et al. 2004, 2008; Allen et al. 2006; Cavagnolo et al. 2010; O’Sullivan et al. 2011). As the total energy within a bubble is the sum of its internal energy and the work carried out by its inflation, the kinetic luminosity is given by the ratio of this total energy and the age of the bubble.

Bîrzan et al. (2008) derived kinetic luminosities for a sample of 24 radio galaxies (L_{1.4 GHz}< 10²⁸ W/Hz, z< 0.4) in clusters or groups and found a rather flat (A = 0.35) slope of the correlation (see Eq. (A.1) below), albeit with a large scatter of σ ~ 0.84 dex ³. This sample was augmented by Cavagnolo et al. (2010) and O’Sullivan et al. (2011) to better constrain the faint end (L_{1.4 GHz}~ 10²⁰ − 10²⁴ W Hz^-1), which resulted in a revised steeper (A ≈ 0.6 − 0.9) slope of the relation with a scatter of σ ~ 0.7 dex. Using similar samples of radio AGN in galaxy clusters (L_{1.4 GHz}≈ 10²⁰ − 10²⁶ W/Hz; Allen et al. 2006; Rafferty et al. 2006) Merloni & Heinz (2007) derived a relationship between the kinetic and (5 GHz) radio core luminosities, corrected for Doppler boosting, finding a slope of A = 0.81 (with a scatter of σ ≈ 0.5 dex). The high-luminosity end of this relation was supplemented using a sample of powerful FR II galaxies (z = 0.056 − 1.8; L_{1.4 GHz}> 10²⁷ W Hz^-1) by Daly et al. (2012). They find a slope of A = 0.84 and no evidence for a redshift evolution of this relation.

Recently the relations based on the data used in Bîrzan et al. (2004), Cavagnolo et al. (2010), and O’Sullivan et al. (2011) were reassessed by Godfrey & Shabala (2016) who have found only a weak correlation with slope A = 0.27 between the two quantities when distance to the systems in the sample (i.e., Malmquist bias) is taken into account (see Eq. (A.1) and their Table 2)⁴. However, contrary to the inferred flatter slope of the distant-dependent relation, based on theoretical grounds Godfrey & Shabala (2016) have predicted a much steeper slope (A ≳ 0.5), which is consistent with those inferred by all of the abovementioned studies. These authors have attributed this discrepancy to either additional model parameters being correlated with the radio jet luminosity or a systematic bias in the X-ray cavity luminosity calculation that induces a flattening in the observed slope.

Further evidence for a steeper slope (A ≈ 0.8 − 0.9) of the scaling relation arises from the analytic derivation by Willott et al. (1999). Willott et al. (1999) estimated the kinetic energy using minimum energy arguments, i.e., computing the minimum energy stored in the lobes to produce the observed synchrotron luminosity, considering the age of the source and the efficiency of conversion of the kinetic luminosity into the internal energy of the observed synchrotron emission. These authors inferred a slope of the relation of A = 0.86. Willott et al. (1999) used a reference frequency of 151 MHz, and they folded all uncertainty factors in the calculation (such as departure from minimum energy conditions, uncertainty in the energy in nonradiating particles, and the composition of the jet) into one parameter, f_W, estimated to lie in the range of f_W ≈ 1 − 20. We here adopted the conversion to a reference frequency of 1.4 GHz as presented by Heckman & Best (2014) (their Eq. (1)), noting that conversions using single power-law spectrum assumptions add additional uncertainties (see also below).

The various scaling relations, shown in Fig. A.1, can be summarized as $\log {\mathit{L}}_{\mathrm{kin}}\mathrm{=}\mathit{A}\mathrm{·}\log \frac{{\mathit{L}}_{\mathrm{1.4}\mathrm{GHz}}}{{\mathrm{10}}^{\mathrm{24}}\mathrm{W}\mathit{/}\mathrm{Hz}}\mathrm{+}\mathit{B,}$(A.1)where L_kin is the kinetic luminosity in units of W, L_{1.4 GHz} is the monochromatic radio luminosity at rest-frame 1.4 GHz in units of W/Hz, and A and B are constants as determined by the various authors as follows: $\begin{array}{ccc}\mathrm{\left(}\mathit{A,B}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \\ & \\ & \\ & \\ & \\ & \end{array}& & \end{array}$where f_W is the factor defined by Willott et al. (1999) encompassing all the possible systematic erros (f_W ≈ 1 − 20) and d_L is the luminosity distance expressed in units of Mpc.

Merloni & Heinz (2007) originally based their scaling relation using 5 GHz core luminosities in units of W (i.e., ${\mathit{L}}_{\mathit{\nu }}^{\mathrm{core}}\mathrm{=}\mathit{\nu }\mathrm{·}{\mathit{L}}_{\mathrm{nu}}^{\mathrm{core}}$, where ν = 5 GHz, and ${\mathit{L}}_{\mathrm{nu}}^{\mathrm{core}}$ is in units of W/Hz). We scaled this relation in such a way that the total rest-frame 1.4 GHz luminosity of the source is taken as input. This was carried out by converting between rest-frame 5 and 1.4 GHz luminosity using a spectral index of α = − 0.5 (Kimball & Ivezić 2008) and assuming that the core luminosity on average corresponds to half of the total luminosity of the AGN as supported by our VLA-COSMOS radio AGN also detected by the VLBA-COSMOS survey at 0.005″ angular resolution (Herrera Ruiz et al., in prep.). We caution, however, that such a conversion may not be applicable at the high-luminosity end (for sources such as, e.g., Cygnus A, not sampled well within the 2 square degreee COSMOS field), and that the conversion, thus, may be luminosity dependent, which is not accounted for here. Similarly, the original expression derived by Daly et al. (2012) is based on total 178 MHz luminosities converted here to monochromatic 1.4 GHz luminosities using a spectral index of α = − 0.7.

Fig. A.1 Various scaling relations between monochromatic 1.4 GHz luminosity (L1.4 GHz) and kinetic luminosity (Lkin) from the literature, as indicated in the top left of the panel and described in detail in the text (see Eq. (A.1)). Also shown, to guide the eye, are the data used by Bîrzan et al. (2004) and O’Sullivan et al. (2011) (L1.4 GHz≈ 1020 − 1028 W Hz-1; filled symbols). The analytically derived Willott et al. (1999) relation shown is that for fW = 15 (see text for details) in agreement with the relations based on a combination of X-ray cluster/group and radio data (symbols). For fW = 4, it would agree with the Daly et al. (2012) relation inferred for powerful radio galaxies (L1.4 GHz> 1027 W Hz-1).

In Fig. A.2 we show the kinetic luminosity density as a function of redshift for our pure luminosity or density evolution two-parameter models and using the various scaling relations. Given the similar (A ≈ 0.8) slope of most of the relations the resulting kinetic luminosity density approximately only systematically shifts along the y-axis. The Willott et al. (1999) relation taking f_W = 1 represents a robust lower limit, and, when considering only the A = 0.6 − 0.9 slope relations, the highest values are obtained using the Merloni & Heinz (2007) relation.

The result based on the Godfrey & Shabala (2016) relation is not shown in Fig. A.2 as it is dependent on the lower integration boundary given the shallow slope of the relation (A = 0.27). Furthermore, it is distance dependent and constrained only based on a sample out to z ≤ 0.23. Thus, extrapolating this relation beyond z ~ 0.2 and up to z ~ 5 would yield extremely high kinetic luminosity values. For example, for L_{1.4 GHz}= 10²⁵ W Hz^-1 the extrapolated Godfrey & Shabala (2016) relation would yield L_kin = (0.8,1.6,2.5,3.7) × 10⁴⁰ W for z = 2,3,4,5, respectively, while the maximum kinetic luminosity for the same input parameters based on the other relations would be L_kin = 1.6 × 10³⁸ W (Merloni & Heinz 2007 relation), thus about two order of magnitude lower. However, out to z ~ 0.1 the Godfrey & Shabala (2016) parametrization yields results that are consistent with the range of those obtained using the other relations considered, and at z> 0.3 this parametrization results in kinetic luminosity densities that are increasingly higher by more than an order of magnitude compared to the results based on the other relations.

Fig. A.2 Top (bottom) panel: kinetic luminosity density as a function of redshift for our pure luminosity (density) evolution 2-parameter model with the various scaling relations commonly applied in the literature, as indicated in each panel. For comparison the radio model related SMBH accretion luminosity from the SAGE semi-analytic model (Croton et al. 2016) is also shown.

All Tables

All Figures

	Fig. 1 Redshift distribution (top panel), 1.4 GHz rest-frame radio luminosity (middle panel), and stellar mass (bottom panel) as a function of redshift for our radio-excess AGN sample.
In the text

	Fig. 2 Spectral indices of our radio-excess AGN sample as a function of redshift (gray points). Median values with interquartile ranges for sources with measured 1.4–3 GHz spectral indices in different redshift bins are shown with blue circles. Mean values for all sources (i.e., also undetected at 1.4 GHz for which we manually set α = − 0.7) in same redshift bins are shown with red squares.
In the text

Fig. 3 Radio AGN 1.4 GHz rest-frame luminosity functions in nine redshift bins out to z ~ 5 for radio-excess AGN in the COSMOS field. The volume densities derived from our radio sample with the Vmax method are shown by the filled black dots. A 1σ confidence interval of the combined luminosity and density evolution in individual redshift bins is shown as the gray shaded area. Pure luminosity evolution and pure density evolution 2-parameter models are shown with blue and red lines, respectively (see text for details). The dotted curve in each panel indicates the local luminosity function (Eq. (2)). The vertical line in each panel shows the 5σ sensitivity limit (assuming a spectral index of α = − 0.7) at the high-redshift end of each redshift bin. Various results from the literature, as indicated in the legend, are also shown (each scaled to the median redshift in each bin using the evolution model as reported in the corresponding study).

In the text

Fig. 4 Best-fit parameters for the evolution of the local LF with redshift. The best-fit pure luminosity (density) evolution in a given redshift bin is shown by the black points in the top (bottom) panel. The blue (red) line in the top (bottom) panel shows the simple, continuous 2-parameter pure luminosity (density) evolution model (see text for details). Gray shaded areas in both panels correspond to the 68% confidence interval for a combined luminosity and density evolution. The large uncertainty in the combined fit is due to parameter degeneracy.

In the text

Fig. 5 Redshift evolution of the number density (top panel), and luminosity density (bottom panel) for our radio AGN (for models as indicated in the panels; see text for details). The number density was integrated with a lower boundary of L1.4 GHz= 1022 W Hz-1. Decreasing the boundary would systematically shift the number densities to higher values (for example, by a factor of 3–5 if L1.4 GHz= 1021 W Hz-1 were used). The luminosity density integral has a much weaker dependence on the integration limits (if it is broad enough) as the bulk of the integral is constrained close to (L∗,Φ∗).

In the text

Fig. 6 Cosmic evolution of the kinetic luminosity density for our radio AGN. The kinetic luminosity density was derived using Eq. (5), and the Willott et al. (1999) relation between monochromatic (1.4 GHz rest-frame) radio luminosity and kinetic luminosity (Eq. (4)) with fW = 4 (pure luminosity evolution, PLE, model: blue full curve and blue shaded area showing the 1σ uncertainty), fW = 15 (purple full curve, and purple-hatched area showing the 1σ uncertainty), and fW = 1, and 20, with folded 1σ fitting errors (lower and upper gray dashed curves, respectively). Also shown is the kinetic luminosity density assumed in the SAGE semi-analytic cosmological model (green dash-dotted curve; Croton et al. 2016). For comparison, the evolution of the kinetic luminosity density for our radio AGN for the pure density evolution, PDE, model using only fW = 4 (Eq. (A.1)) is also shown with the red full curve and red hatched area. To avoid overcrowding the panel, the PDE result for fW = 15 is not shown, but it would be only systematically higher compared to the fW = 4 PDE model result and coincident with that for PLE with the corresponding fW = 15 value. See also Fig. A.2 for results based on other scaling relations.

In the text

	Fig. 7 Absolute (top panel) and fractional (bottom panel) contributions of various subpopulations (as indicated in the bottom panel) to the full radio-excess AGN sample (black curve) as a function of redshift.
In the text

Fig. A.1 Various scaling relations between monochromatic 1.4 GHz luminosity (L1.4 GHz) and kinetic luminosity (Lkin) from the literature, as indicated in the top left of the panel and described in detail in the text (see Eq. (A.1)). Also shown, to guide the eye, are the data used by Bîrzan et al. (2004) and O’Sullivan et al. (2011) (L1.4 GHz≈ 1020 − 1028 W Hz-1; filled symbols). The analytically derived Willott et al. (1999) relation shown is that for fW = 15 (see text for details) in agreement with the relations based on a combination of X-ray cluster/group and radio data (symbols). For fW = 4, it would agree with the Daly et al. (2012) relation inferred for powerful radio galaxies (L1.4 GHz> 1027 W Hz-1).

In the text

	Fig. A.2 Top (bottom) panel: kinetic luminosity density as a function of redshift for our pure luminosity (density) evolution 2-parameter model with the various scaling relations commonly applied in the literature, as indicated in each panel. For comparison the radio model related SMBH accretion luminosity from the SAGE semi-analytic model (Croton et al. 2016) is also shown.
In the text