Galaxy Evolution in the Radio Band: The Role of Star-forming Galaxies and Active Galactic Nuclei

The first ingredient is constituted by the global SFR function ${dN}/d\,\mathrm{log}{\dot{M}}_{\star }$, namely, the number density of galaxies per logarithmic bin of SFR $[\mathrm{log}{\dot{M}}_{\star },\,\mathrm{log}{\dot{M}}_{\star }+d\,\mathrm{log}{\dot{M}}_{\star }]$ at given redshift z. This has been accurately determined by Mancuso et al. (2016a, 2016b) by exploiting the most recent determinations of the evolving galaxy luminosity functions from far-IR and UV data.

In a nutshell, UV data have been dust-corrected according to the local empirical relation between the UV slope β_UV and the IR-to-UV luminosity ratio IRX (see Meurer et al. 1999; Calzetti et al. 2000), which is also routinely exploited for high-redshift galaxies (see Bouwens et al. 2009, 2015, 2016a, 2016b). For SFGs with intrinsic SFR ${\dot{M}}_{\star }\gtrsim 30\,{M}_{\odot }$ yr⁻¹ the UV data, even when dust-corrected via the UV slope–IRX relationship, strongly underestimate the intrinsic SFR, which is instead better probed by far-IR observations. This is because high SFRs occur primarily within heavily dust-enshrouded molecular clouds, while the UV slope mainly reflects the emission from stars obscured by the diffuse, cirrus dust component (Silva et al. 1998; Efstathiou et al. 2000; Efstathiou & Rowan-Robinson 2003; Coppin et al. 2015; Reddy et al. 2015; Mancuso et al. 2016a). On the other hand, at low SFR ${\dot{M}}_{\star }\,\lesssim 10\,{M}_{\odot }$ yr⁻¹ the dust-corrected UV data efficiently probe the intrinsic SFR. Moreover, in late-type galaxies at z ≲ 1 the far-IR emission itself can be contributed by the cirrus component, heated by the general radiation field from evolved stellar populations. To correct for such an effect, which otherwise may cause the SFR inferred from far-IR data to be appreciably overestimated, we have adopted the prescription by Clemens et al. (2013). In Figure 1 we report the overall data compilation from far-IR and dust-corrected UV observations. The luminosity L and SFR ${\dot{M}}_{\star }$ scale have been related using $\mathrm{log}{\dot{M}}_{\star }/{M}_{\odot }\,{\mathrm{yr}}^{-1}\approx -9.8+\mathrm{log}L/{L}_{\odot }$, a good approximation for both far-IR and (intrinsic) UV luminosities.

Zoom In Zoom Out Reset image size

Then we determined a smooth, analytic representation of the SFR function in terms of the standard Schechter shape

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{d\,\mathrm{log}{\dot{M}}_{\star }}({\dot{M}}_{\star },z)={ \mathcal N }(z){\left[\displaystyle \frac{{\dot{M}}_{\star }}{{\dot{M}}_{\star ,c}(z)}\right]}^{1-\alpha (z)}\,{e}^{-{\dot{M}}_{\star }/{\dot{M}}_{\star ,c}(z)},\end{eqnarray} \tag{ 1 }$$

characterized at any given redshift z by three parameters, namely, the normalization ${ \mathcal N }$, the characteristic SFR ${\dot{M}}_{\star ,c}$, and the faint-end slope α. We determine the values of the three Schechter parameters over the range z ∼ 0–10 in unitary redshift bins by performing an educated fit to the data. Specifically, UV data are fitted for SFRs ${\dot{M}}_{\star }\lesssim 30\,{M}_{\odot }$ yr⁻¹ since in this range dust corrections based on the β_UV are reliable, while far-IR data are fitted for SFRs ${\dot{M}}_{\star }\gtrsim {10}^{2}\,{M}_{\odot }$ yr⁻¹ since in this range dust emission is largely dominated by molecular clouds and reflects the ongoing SFR. To obtain a smooth yet accurate representation of the SFR functions at any redshift, we find it necessary to (minimally) describe the redshift evolution for each parameter p(z) of the Schechter function as a third-order polynomial in log-redshift $p(z)\,={p}_{0}+{p}_{1}\,\xi +{p}_{2}\,{\xi }^{2}+{p}_{3}\,{\xi }^{3}$, with $\xi =\mathrm{log}(1+z)$. The values of the parameters $\{{p}_{i}\}$ are reported in Table 1. The resulting SFR functions for representative redshifts z ≈ 0, 1, 3, and 6 are illustrated in Figure 1.

Table 1.  SFR Function Parameters

Parameter	p0	p1	p2	p3
$\mathrm{log}{ \mathcal N }(z)$	−2.13	−8.90	18.07	−9.58
$\mathrm{log}{\dot{M}}_{\star ,c}(z)$	0.72	8.56	−10.08	2.54
$\alpha (z)$	1.12	3.73	−7.80	5.15

Note. We have adopted the Meurer–Calzetti law to compute dust correction for UV data and the Clemens et al. (2013) prescriptions to subtract cirrus emission from low-redshift (z ≲ 1) far-IR data.

Download table as:  ASCIITypeset image

In Mancuso et al. (2016a, 2016b) and Lapi et al. (2017) we have validated the global SFR functions against independent data sets, including galaxy number counts at significant submillimeter/far-IR wavelengths, redshift distributions of gravitationally lensed galaxies, galaxy stellar mass function via the continuity equation, main sequence of SFGs, cosmological evolution of the average SFR and gamma-ray burst rates, and high-redshift observables including the history of cosmic reionization.

The second ingredient is constituted by deterministic evolutionary tracks for the history of star formation and BH accretion in an individual galaxy, gauged on a wealth of multiwavelength observations and inspired theoretically by the in situ coevolution scenario. This envisages star formation and BH accretion in galaxies to be essentially in situ, time-coordinated processes (e.g., Lapi et al. 2006, 2011, 2014; see also Lilly et al. 2013; Aversa et al. 2015; Mancuso et al. 2016a, 2016b), triggered by the early collapse of the host dark matter halos, but subsequently controlled by self-regulated baryonic physics and in particular by energy/momentum feedbacks from supernovae and AGNs.

In a nutshell, during the early stages of a galaxy's evolution, the competition between gas condensation and energy/momentum feedback from supernovae and stellar winds regulates the SFR. In low-mass galaxies the SFR is small, ${\dot{M}}_{\star }\lesssim {\rm{a}}$ few tens M_⊙ yr⁻¹, and it slowly decreases over long timescales of several gigayears because of progressive gas consumption. On the other hand, in high-mass galaxies huge gas reservoirs can sustain violent, almost constant SFR ${\dot{M}}_{\star }\gtrsim {10}^{2}\,{M}_{\odot }$ yr⁻¹, while the ambient medium is quickly enriched with metals and dust; the galaxy behaves as a bright submillimeter/far-IR source. After a time τ_b ∼ some 10⁸ yr the SFR is abruptly quenched by the energy/momentum feedback from the central supermassive BH, and the environment is cleaned out; thereafter the stellar populations evolve passively and the galaxy becomes a red and dead early-type.

From the point of view of the central BH, during the early stages plenty of gas is available from the surroundings, so that considerable accretion rates sustain mildly super-Eddington emission with Eddington ratios $\lambda \equiv L/{L}_{\mathrm{Edd}}\gtrsim 1;$ radiation trapping and relativistic effects enforce radiatively inefficient, slim-disk conditions (see Begelman 1979; Li 2012; Madau et al. 2014). During these early stages, the BH bolometric luminosity is substantially smaller than that of the host SFG, but it increases exponentially. After a time τ_b ≳  a few 10⁸ yr, the nuclear power progressively increases to values similar to or even exceeding that from star formation in the host galaxy. As mentioned above, strong energy/momentum feedback from the BH removes interstellar gas and dust while quenching star formation; the system behaves as an optical quasar. Residual gas present in the central regions of the galaxy can be accreted onto the BH at progressively lower, sub-Eddington accretion rates. When the Eddington ratio falls below a critical value around λ ≲ 0.3 (see McClintock et al. 2006), the disk becomes thin, yielding the standard spectral energy distributions (SEDs) observed in type 1 AGNs. Eventually, the BH activity ceases because of gas exhaustion in the nuclear region. At low redshift z ≲ 1, especially within a rich environment, gravitational interaction or even a galaxy merger can temporarily rekindle a starburst and the BH activity. A schematic evolution of the SFR and BH accretion rate as a function of the galaxy age is reported in Figure 2.

Zoom In Zoom Out Reset image size

On this basis, we have computed the relative time spent by the AGN in a given logarithmic bin of bolometric luminosity L_AGN, i.e., the AGN duty cycle, as

$$\begin{eqnarray}&&\displaystyle \frac{d\delta }{d\,\mathrm{log}{L}_{\mathrm{AGN}}}({L}_{\mathrm{AGN}},z| {\dot{M}}_{\star })\approx \displaystyle \frac{{\tau }_{\mathrm{ef}}+{\tau }_{\mathrm{AGN}}}{{\tau }_{{\rm{b}}}}\,\mathrm{ln}10;\end{eqnarray} \tag{ 2 }$$

here τ_ef is the e-folding time (depending on the Eddington ratio λ and radiative efficiency) during the early AGN phase, τ_AGN is the characteristic time of the declining AGN phase, and τ_b is the duration of the star formation period before the AGN quenching. In Mancuso et al. (2016a, 2016b) such parameters have been set by comparison with observations, including the Eddington ratio distributions at different redshifts, the fraction of host SFGs in optically/X-ray-selected quasars, the fraction of AGN hosts with given stellar mass as a function of the Eddington ratio, the BH mass function via the continuity equation, the main sequence of SFGs, and the AGN coevolution plane (i.e., bolometric AGN luminosity vs. SFR or stellar mass). Note that the AGN duty cycle depends on the average SFR through the above parameters, since at the end of the evolution of the galaxy the central BH-to-stellar mass ratio M_BH/M_⋆ must take on the locally observed values ≈10⁻³ (e.g., Kormendy & Ho 2013; McConnell & Ma 2013; Shankar et al. 2016). We defer the reader to the Mancuso et al. (2016b) paper for a detailed description of the above parameter values and their dependence on the average SFR, which we adopt in full here.

In terms of the duty cycle, we can now map the SFR functions into the AGN bolometric luminosity functions as

$$\begin{eqnarray}\displaystyle \frac{{dN}}{d\,\mathrm{log}{L}_{\mathrm{AGN}}}({L}_{\mathrm{AGN}},z) & = & \displaystyle \int d\,\mathrm{log}{\dot{M}}_{\star }\,\displaystyle \frac{{dN}}{d\,\mathrm{log}{\dot{M}}_{\star }}\,\displaystyle \frac{d\delta }{d\,\mathrm{log}{L}_{\mathrm{AGN}}}\\ & & \times ({L}_{\mathrm{AGN}},z| {\dot{M}}_{\star }).\end{eqnarray} \tag{ 3 }$$

The outcome for representative redshifts z ≈ 0, 1, 3, and 6 is illustrated in Figure 3 and compared with a data compilation from optical and hard X-ray observations. The data have been converted to bolometric luminosity using the Hopkins et al. (2007) corrections, while the corresponding number densities have been corrected for obscured (also Compton-thick) AGNs after Ueda et al. (2014). The pleasing agreement between our determination and the data confirms that the AGN duty cycle is correctly determined.

Zoom In Zoom Out Reset image size

We anticipate that to properly address the radio emission in RQ systems, it will be convenient to compute the luminosity function of SFGs hosting an AGN with X-ray emission above a given threshold L_X,min. This quantity is given by

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{d\,\mathrm{log}{\dot{M}}_{\star }}({\dot{M}}_{\star },z| \,\gt \,{L}_{{\rm{X}},\min })=\displaystyle \frac{{dN}}{d\,\mathrm{log}{\dot{M}}_{\star }}({\dot{M}}_{\star },z)\\ &&\quad \times {\displaystyle \int }_{\gt {L}_{{\rm{X}},\min }}d\,\mathrm{log}{L}_{\mathrm{AGN}}\,\displaystyle \frac{d\delta }{d\,\mathrm{log}{L}_{\mathrm{AGN}}}({L}_{\mathrm{AGN}},z| {\dot{M}}_{\star });\end{eqnarray} \tag{ 4 }$$

the threshold ${L}_{{\rm{X}},\min }\approx {10}^{42}$ erg s⁻¹ will be employed, since this is the value commonly adopted by observers (e.g., Bonzini et al. 2013; Padovani et al. 2015) to clearly discern the nuclear X-ray emission from that associated with star formation ${L}_{{\rm{X}},\mathrm{SFR}}\approx 7\times {10}^{41}$ erg ${{\rm{s}}}^{-1}\,({\dot{M}}_{\star }/{10}^{2}\,{M}_{\odot }\,{\mathrm{yr}}^{-1})$ (see, e.g., Vattakunnel et al. 2012). The resulting SFR functions are illustrated as dashed lines in Figure 1. Most of the AGNs with significant X-ray powers are hosted at z ≳ 1 by galaxies with SFRs ${\dot{M}}_{\star }\gtrsim {10}^{2}\,{M}_{\odot }$ yr⁻¹; however, their number density is smaller than that of the global star-forming population by a factor of 10⁻¹, which reflects the nearly constant behavior of the SFR versus the exponential growth of the BH accretion rate during most of the galaxy lifetime, quantified by the AGN duty cycle.

The radio emission associated with star formation comprises two components that are well known to correlate with the SFR (see Condon 1992; Bressan et al. 2002; Murphy et al. 2011): free–free emission (fully dominating at frequencies ν ≳ 30 GHz) emerging directly from H ii regions containing massive, ionizing stars; and synchrotron emission resulting from relativistic electrons accelerated by supernova remnants.

As to the free–free emission, we use the classic expression (see Murphy et al. 2011; Mancuso et al. 2015)

$$\begin{eqnarray}{L}_{\mathrm{ff}} & \approx & 3.75\times {10}^{26}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}\,\displaystyle \frac{{\dot{M}}_{\star }}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}\\ & & \times {\left(\displaystyle \frac{T}{{10}^{4}{\rm{K}}}\right)}^{0.3}\,g(\nu ,T){e}^{-h\nu /{kT}},\end{eqnarray} \tag{ 5 }$$

where $g(\nu ,T)$ is the Gaunt factor

$$\begin{eqnarray}&&g(\nu ,T)\\ &&=\mathrm{ln}\left\{\exp \left[5.960-\displaystyle \frac{\sqrt{3}}{\pi }\,\mathrm{ln}\left({Z}_{i}\,\displaystyle \frac{\nu }{\mathrm{GHz}}{\left(\displaystyle \frac{T}{{10}^{4}{\rm{K}}}\right)}^{-1.5}\right)\right]+e\right\},\end{eqnarray} \tag{ 6 }$$

approximated according to Draine (2011), and the quantity ${e}^{-h\nu /{kT}}$ tentatively renders electron energy losses. This equation reproduces the Murphy et al. (2012) calibration at 33 GHz for a pure hydrogen plasma (Z_i = 1) with temperature T ≈ 10⁴ K; we adopt these values in the following.

As to the synchrotron emission, the calibration with the SFR is a bit more controversial since it involves complex and poorly understood processes such as the production rate of relativistic electrons, the fraction of them that can escape from the galaxy, and the magnetic field strength. We use the calibration proposed by Murphy et al. (2011, 2012) and then adopted in the widely cited review by Kennicutt & Evans (2012),

$$\begin{eqnarray}{L}_{\mathrm{sync}} & \approx & 1.9\times {10}^{28}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}\,\displaystyle \frac{{\dot{M}}_{\star }}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}{\left(\displaystyle \frac{\nu }{\mathrm{GHz}}\right)}^{-{\alpha }_{\mathrm{sync}}}\\ & & \times {\left[1+{\left(\displaystyle \frac{\nu }{20\mathrm{GHz}}\right)}^{0.5}\right]}^{-1}\,F[{\tau }_{\mathrm{sync}}(\nu )],\end{eqnarray} \tag{ 7 }$$

where ${\alpha }_{\mathrm{sync}}\approx 0.75$ is the spectral index (e.g., Condon 1992), the term in square brackets renders spectral aging effects (see Banday & Wolfendale 1991), and the function $F(x)=(1-{e}^{-x})/x$ takes into account synchrotron self-absorption in terms of the optical depth (e.g., Kellermann 1966; Tingay & de Kool 2003)

$$\begin{eqnarray}&&{\tau }_{\mathrm{sync}}\approx {(\nu /{\nu }_{\mathrm{self}})}^{-{\alpha }_{\mathrm{sync}}-5/2}\end{eqnarray} \tag{ 8 }$$

that is thought to become relevant at frequencies $\nu \lesssim {\nu }_{\mathrm{self}}\,\approx 200\,\mathrm{MHz}$.

Given that for the Chabrier IMF the far-IR luminosity is given by ${L}_{\mathrm{FIR}}\,[{\rm{W}}]\approx 3\times {10}^{36}\,{\dot{M}}_{\star }\,[{M}_{\odot }\,{\mathrm{yr}}^{-1}]$ in terms of the SFR, the above calibrations Equations (5) and (7) yield a far-IR versus 1.4 GHz correlation parameter q_FIR ≡ log(L_FIR/3.75 × 10¹² W) − log(L_1.4GHz/W Hz⁻¹) ≈ 2.77; this is slightly higher than the classic value q_FIR ≈ 2.35 (see Yun et al. 2001) but in excellent agreement with the recent determinations by Novak et al. (2017) and Delhaize et al. (2017). Note that at ν ≈ 1.4 GHz a synchrotron-to-free–free luminosity ratio L_synch/L_ff ≈ 5.4 is found, somewhat lower than the classic value ≈8 quoted by Condon (1992) from his analysis of the M82 SED, but in good agreement with more recent data and models (see Bressan et al. 2002; Murphy et al. 2011, 2012; Obi et al. 2017).

Following Mancuso et al. (2015), we also take into account the lower efficiency in producing synchrotron emission by galaxies with small SFRs ${\dot{M}}_{\star }\lesssim {\rm{a}}$ few M_⊙ yr⁻¹ (see Bell 2003), on correcting Equation (7) as

$$\begin{eqnarray}&&{L}_{\mathrm{sync},\mathrm{corr}}=\displaystyle \frac{{L}_{\mathrm{sync}}}{1+{({L}_{0,\mathrm{sync}}/{L}_{\mathrm{synch}})}^{\zeta }},\end{eqnarray} \tag{ 9 }$$

with ζ ≈ 2 and ${L}_{0,\mathrm{sync}}\approx 3\times {10}^{28}$ erg s⁻¹ Hz⁻¹. We anticipate that this correction is necessary to reproduce the local 1.4 GHz luminosity function at small radio powers ${L}_{1.4\mathrm{GHz}}\lesssim {L}_{0,\mathrm{sync}}$.

Note that other phenomena may contribute to change the synchrotron luminosity in specific frequency and redshift range, and that in the lack of a consensus physical understanding and detailed modeling we decide not to include in our fiducial approach. First, at low frequencies in a dense medium where relativistic and thermal electrons spatially coexist, the synchrotron emission can be absorbed; this would be described by an additional multiplicative factor ${e}^{-{\tau }_{\mathrm{ff}}}$ in Equation (7), where

$$\begin{eqnarray}&&{\tau }_{\mathrm{ff}}\approx {\left(\displaystyle \frac{T}{{10}^{4}{\rm{K}}}\right)}^{-1.35}{\left(\displaystyle \frac{\nu }{1.4\mathrm{GHz}}\right)}^{-2.1}\,\displaystyle \frac{\mathrm{EM}}{6\times {10}^{6}\,\mathrm{pc}\,{\mathrm{cm}}^{-6}}\end{eqnarray} \tag{ 10 }$$

is the free–free absorption optical depth (Condon 1992; Bressan et al. 2002) in terms of the emission measure EM of the plasma. However, on average synchrotron emission in SFGs is produced on spatial scales much larger than that of the thermal electrons. Moreover, in most of the local galaxy population the physical conditions would assume free–free absorption to become relevant only at very low frequencies ν ≲ 100 MHz, although there are some controversial cases related to starburst cores where absorption has been reported to be effective even at ν ≲ 1 GHz (e.g., Vega et al. 2008; Schober et al. 2017). Note also that in high-redshift SFGs the expected increase in the average density of the medium is easily offset by the z-dependence induced in the rest frame ${\tau }_{\mathrm{ff}}\propto {\nu }^{-2.1}{(1+z)}^{-2.1}$.

Second, at high redshift the relativistic electrons producing synchrotron can lose energy owing to a number of processes, most noticeably inverse Compton scattering off cosmic microwave background (CMB) photons (e.g., Murphy 2009; Lacki & Thompson 2010; Schober et al. 2017), whose energy density grows with redshift as ${(1+z)}^{4}$. However, the amplitude of the breakdown and the redshifts at which it happens vary in connection with the assumed properties of high-redshift galaxies. In particular, Bonato et al. (2017) have shown that such effects are constrained to be minor at least out to redshift z ≲ 3. In fact, Smith et al. (2014) have measured a direct dependence of the radio-to-monochromatic-far-IR luminosity ratio L_{1.4 GHz}/L_{250 μm} with dust temperature, which could lead to some balancing of the inverse Compton losses. We stress that at high rest-frame frequencies the free–free emission (not affected by the above) will dominate anyway over the synchrotron, so that the correction to the total radio power will be small.

Third, an open issue concerns the evolution with redshift of the normalization in the radio luminosity versus SFR relation, namely, the ${q}_{\mathrm{FIR}}$ parameter mentioned above. Some studies found it to be unchanged or to undergo only minor variations with redshift (e.g., Ibar et al. 2008; Bourne et al. 2011; Mao et al. 2011; Smith et al. 2014), while others have reported a significant, albeit weak, evolution (e.g., Seymour et al. 2009; Basu et al. 2015; Magnelli et al. 2015; Delhaize et al. 2017; Novak et al. 2017). This can be rendered on multiplying Equation (7) by a factor ${10}^{{q}_{0}[1-{(1+z)}^{-{q}_{1}}]}$; Magnelli et al. (2015) find q₀ = 2.35 and q₁ = 0.12, Novak et al. (2017) report q₀ = 2.77 and q₁ = 0.14, and Delhaize et al. (2017) suggest q₀ = 2.9 and q₁ = 0.19. On the other hand, Bonato et al. (2017) have shown that the evolution by Magnelli et al. (2015) is marginally consistent with the 1.4 GHz radio luminosity functions and deep counts down to μJy levels; in Section 4 we will revise the issue in view of the most recent data on the redshift-dependent radio luminosity function (Novak et al. 2017) and counts (Prandoni et al. 2017) of SFGs.

The average total radio power for a given SFR is the sum of the contribution from synchrotron and free–free emission ${\bar{L}}_{\nu }({\dot{M}}_{\star })={L}_{\mathrm{synch},\mathrm{corr}}+{L}_{\mathrm{ff}}$. In Figure 4 we show these quantities as a function of the SFR for three representative frequencies ν ≈ 0.15, 1.4, and 10 GHz, and as a function of the frequency for three different values of the SFR. It is seen that the free–free emission increasingly dominates over synchrotron in moving toward high frequencies ν ≳ 10 GHz owing to its flatter spectrum, and at small SFRs ${\dot{M}}_{\star }\lesssim {\rm{a}}$ few M_⊙ yr⁻¹ owing to the inefficiency of synchrotron emission after Equation (9). At low frequencies ν ≲ 200 MHz the synchrotron emission is also suppressed because of the self-absorption process.

Zoom In Zoom Out Reset image size

We consider, consistently with observations (see Condon 1992; Bressan et al. 2002), a Gaussian scatter of ${\sigma }_{\mathrm{log}L}\approx 0.2$ dex around the average ${\bar{L}}_{\nu }({\dot{M}}_{\star })$ relationship. The radio luminosity function of SFGs is then obtained as

$$\begin{eqnarray}\displaystyle \frac{{{dN}}_{\mathrm{SF}}}{d\,\mathrm{log}{L}_{\nu }}({L}_{\nu },z) & = & \displaystyle \frac{1}{\sqrt{2\pi }\,{\sigma }_{\mathrm{log}L}}\,\displaystyle \int d\,\mathrm{log}{\dot{M}}_{\star }\displaystyle \frac{{dN}}{d\,\mathrm{log}{\dot{M}}_{\star }}\\ & & \times ({\dot{M}}_{\star },z){e}^{-{[\mathrm{log}{L}_{\nu }-\mathrm{log}{\bar{L}}_{\nu }({\dot{M}}_{\star })]}^{2}/2{\sigma }_{\mathrm{log}L}^{2}}.\end{eqnarray} \tag{ 11 }$$

According to the framework discussed in Section 2.2 and illustrated in Figure 2, during the early phase of a massive galaxy's evolution, the SFR is sustained at high, nearly constant values, while the BH mass is small but increases rapidly. After a few e-folding times, the X-ray power from the nucleus overwhelms that associated with star formation, so that the AGN is clearly detectable in X-rays at luminosities L_X ≳ 10⁴² erg s⁻¹. This is the case for many SFGs selected in the far-IR band and then followed up in X-rays (e.g., Mullaney et al. 2012; Delvecchio et al. 2015; Rodighiero et al. 2015). However, the situation in the radio band is likely very different.

Since the BH accretes at high rates from large gas reservoirs, slim-disk conditions develop, featuring rather low radiative efficiency because of photon trapping and relativistic effects (see Begelman 1979; Li 2012; Madau et al. 2014); the accretion rates can be substantially super-Eddington ${\dot{M}}_{\mathrm{BH}}\gg {L}_{\mathrm{Edd}}/{c}^{2}$, but the emitted luminosity is only moderately above Eddington with $\lambda \equiv L/{L}_{\mathrm{Edd}}\,\gtrsim$ a few. The accretion is nearly spherical and chaotic; hence, the spins of the BH and of the disk stay small, and rotational energy cannot be easily funnelled into jets to power radio emission via the Blandford & Znajek (1977) or the Blandford & Payne (1982) mechanisms (e.g., Meier 2002; Jester 2005; Fanidakis et al. 2011). In addition, in this phase the BH is growing but still too small to originate large-scale AGN outflows and winds (indeed, star formation is ongoing in the host), so that even nuclear radio emission from AGN-driven shocks is not expected. In fact, the absence of appreciable nuclear radio emission in strongly star-forming, gas-rich SFGs at high redshift appears to be confirmed by recent observations (see Ma et al. 2016; Heywood et al. 2017).

All in all, in this phase the AGN may be detectable in X-rays, but it is almost RS, so that any radio emission from the system should be mainly ascribed to the SFR in the host SFG. The number density of SFGs hosting an RS AGN detectable in X-rays can be easily estimated on the basis of Equation (4) as

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{\mathrm{RQ} \mbox{-} \mathrm{SF}}}{d\,\mathrm{log}{L}_{\nu }}({L}_{\nu },z)=\displaystyle \frac{1}{\sqrt{2\pi }\,{\sigma }_{\mathrm{log}L}}\,\displaystyle \int d\,\mathrm{log}{\dot{M}}_{\star }\displaystyle \frac{{dN}}{d\,\mathrm{log}{\dot{M}}_{\star }}\\ &&\quad \times ({\dot{M}}_{\star },z| \gt {L}_{X,\min }){e}^{-{[\mathrm{log}{L}_{\nu }-\mathrm{log}{\bar{L}}_{\nu }({\dot{M}}_{\star })]}^{2}/2{\sigma }_{\mathrm{log}L}^{2}},\end{eqnarray} \tag{ 12 }$$

where ${L}_{X,\min }\approx {10}^{42}$ erg s⁻¹ is usually considered by observers as the selection threshold for the presence of the AGN. We checked that adopting as a threshold the X-ray emission L_X,SFR ≈ 7 × 10⁴¹ erg ${{\rm{s}}}^{-1}\,({\dot{M}}_{\star }/{10}^{2}\,{M}_{\odot }\,{\mathrm{yr}}^{-1})$ associated with star formation after the calibration by Vattakunnel et al. (2012) does not change appreciably the outcome.

During the late evolution of a massive galaxy with age exceeding some 10⁸ yr, the BH has grown to large masses and originates outflows that can quench the star formation in the host. Meanwhile, the BH accretion rates decline to sub-Eddington levels, and the accretion disk becomes thin and radiatively efficient (see Shakura & Sunyaev 1973). The BH and the accretion disk spin up rapidly, and rotational energy can be easily funnelled into jets (e.g., Blandford & Znajek 1977; Blandford & Payne 1982).

However, the jets driven by thin-disk accretion are rather ineffective in producing radio emission with respect to the advection-dominated flows powering low-z, steep-spectrum RL objects⁷ (e.g., Meier 2002; Jester 2005; Fanidakis et al. 2011). This appear to be confirmed by the observed anticorrelation between the ratio of the jet power and the accretion luminosity versus the Eddington ratio (see Fernandes et al. 2011; Punsly & Zhang 2011; Sikora et al. 2013; Rusinek et al. 2017). In most of the instances, the net result of a thin-disk accretion will be RQ AGNs with weak, small-scale jets. Note, however, that when the BH mass is very large, thin-disk conditions may also originate a flat-spectrum RL AGN, as can be the case for some blazars observed out to high redshift (see Ghisellini et al. 2013); on the other hand, these sources constitute a minority (≲10%) both of the radiatively efficient AGNs and of the overall RL AGN population (see Section 3.4).

Alternatively to the small-scale jet origin, the radio emission of RQ AGNs may be traced back to shock fronts associated with the AGN-driven outflow (e.g., Zakamska & Greene 2014; Nims et al. 2015), or to winds that originated from the outermost portion of the thin accretion disk (see Blundell et al. 2001; King et al. 2013), or to electron acceleration via magnetic reconnection in the thin-disk corona (see Laor & Behar 2008; Raginski & Laor 2016). Whatever mechanism operates, in this phase of a massive galaxy's evolution the radio emission from the system should be mainly ascribed to the nuclear activity typical of an RQ AGN.

The number density of RQ AGNs can be estimated as follows. We start from the AGN bolometric luminosity function of Equation (3), convert the bolometric power in X-rays via the Hopkins et al. (2007) correction, and then derive the AGN radio power by using the relation between rest-frame X-ray and 1.4 GHz radio luminosity observed for a sample of (mainly) RQ AGNs by Panessa et al. (2015; see also Brinkmann et al. 2000),

$$\begin{eqnarray}{L}_{\nu ,\mathrm{AGN}} & \approx & 5.7\times {10}^{-22}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}{\left(\displaystyle \frac{{L}_{X,\mathrm{AGN}}}{\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{1.17}\\ & & \times {\left(\displaystyle \frac{\nu }{1.4\mathrm{GHz}}\right)}^{-{\alpha }_{\mathrm{AGN}}},\end{eqnarray} \tag{ 13 }$$

where α_AGN ≈ 0.7 is the spectral index for an optically thin synchrotron emission. Consistently with observations (e.g., Brinkmann et al. 2000; Panessa et al. 2015), we consider a scatter ${\sigma }_{\mathrm{log}{L}_{X}}\approx 0.4$ around the resulting average relationship ${\bar{L}}_{\nu ,\mathrm{AGN}}({L}_{\mathrm{AGN}})$. The related statistics is given by

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{\mathrm{RQ} \mbox{-} \mathrm{AGN}}}{d\,\mathrm{log}{L}_{\nu }}({L}_{\nu },z)=\displaystyle \frac{1}{\sqrt{2\pi }\,{\sigma }_{\mathrm{log}{L}_{X}}}\,\displaystyle \int d\,\mathrm{log}{L}_{\mathrm{AGN}}\,\displaystyle \frac{{dN}}{d\,\mathrm{log}{L}_{\mathrm{AGN}}}\\ &&\quad \times ({L}_{\mathrm{AGN}},z){e}^{-{[\mathrm{log}{L}_{\nu }-\mathrm{log}{\bar{L}}_{\nu ,\mathrm{AGN}}({L}_{\mathrm{AGN}})]}^{2}/2{\sigma }_{\mathrm{log}{L}_{X}}^{2}}.\end{eqnarray} \tag{ 14 }$$

RL AGNs are galaxies with clear signs of intense AGN activity in the radio band. In fact, RL systems constitute a mixed bag of objects with rather different properties in terms of accretion levels, excitation line emission, and variability timescales (see review by Padovani 2016; Tadhunter 2016).

Steep-spectrum RL AGNs are characterized by optically thin synchrotron radio emission from (large-scale) relativistic jets; they are typically associated with low-redshift (z ≲ 1) activity of very massive BHs at the center of massive early-type galaxies. These sources feature low Eddington ratios λ ≲ 10⁻² likely enforced by advection-dominated accretion flows (e.g., Narayan & Yi 1994; Meier 2002; Fanidakis et al. 2011). The associated large-scale jets have been shown to affect the thermodynamics of the surrounding intracluster medium (see reviews by McNamara & Nulsen 2007; Cavaliere & Lapi 2013) and even to lift huge molecular gas reservoirs, possibly promoting star formation (see Russell et al. 2017). Flat-spectrum RL AGNs are instead characterized by more compact radio emission, typically associated with radiatively efficient, thin-disk conditions (see review by Heckman & Best 2014; see also Massardi et al. 2016).

RL AGNs are well known to dominate the bright portion of the radio counts above 0.5 mJy at 1.4 GHz. As such, they are not our main interest in this paper, which is mainly focused on the faint radio counts to be explored via next-generation surveys. However, for comparison with the total radio luminosity function and counts observed to date, we empirically include them in our analysis.

RL AGNs constitute a small fraction of the overall AGN population, reaching at most ≲10% for powerful optically or X-ray-selected quasars (see Williams & Rottgering 2015); for many of them the fueling mechanism is highly stochastic and orientation effects are relevant. Thus, a description similar to that we pursued for RQ AGNs, based on the AGN bolometric luminosity function and on the radio versus X-ray luminosity correlation (see Sambruna et al. 1999; Fan & Bai 2016), is not viable. Therefore, we return to the empirical description of the cosmological evolution for RL objects by Massardi et al. (2010), which has been extensively tested against a wealth of data on luminosity function and redshift distributions at least out to redshift z ≲ 3. For the reader's convenience we provide a brief account of the Massardi et al. (2010) description here.

These authors consider two flat-spectrum populations with different evolutionary properties, namely, flat-spectrum radio quasars and BL Lacs, and a single steep-spectrum population; for sources of each population a simple power-law spectrum is adopted, ${S}_{\nu }\propto {\nu }^{-\alpha }$ with ${\alpha }_{\mathrm{FSRQ}}={\alpha }_{\mathrm{BLLac}}=0.1$ and α_SS = 0.8. The comoving luminosity function at a given redshift is described by a double power law

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{d\,\mathrm{log}{L}_{\nu }}({L}_{\nu },z)=\displaystyle \frac{{n}_{0}}{{[{L}_{\nu }/{L}_{c}(z)]}^{a}+{[{L}_{\nu }/{L}_{c}(z)]}^{b}},\end{eqnarray} \tag{ 15 }$$

and its redshift evolution is rendered in terms of a pure luminosity evolution

$$\begin{eqnarray}\mathrm{log}{L}_{c}(z) & = & \mathrm{log}{L}_{c}(0)+2\,{k}_{\mathrm{ev}}\,z\,{z}_{\mathrm{top}}\\ & & \times [1-{(z/{z}_{\mathrm{top}})}^{{m}_{\mathrm{ev}}}/(1+{m}_{\mathrm{ev}})],\end{eqnarray} \tag{ 16 }$$

where

$$\begin{eqnarray}&&{z}_{\mathrm{top}}={z}_{\mathrm{top},0}+\displaystyle \frac{\delta {z}_{\mathrm{top}}}{1+{L}_{c}(0)/{L}_{\nu }}\end{eqnarray} \tag{ 17 }$$

is the redshift at which L_c(z) reaches its maximum.

This empirical rendition is characterized by eight parameters: n₀, a, b, ${L}_{c}(0)$, k_ev, m_ev, z_top,0, and $\delta {z}_{\mathrm{top}}$, with different values for each of the three populations considered (flat-spectrum radio quasars, BL Lacs, and steep-spectrum radio quasars). The parameters have been determined by Massardi et al. (2010) by fitting the luminosity function and redshift distributions from various surveys; we defer the reader to the Massardi et al. paper for a full description of this procedure and for the resulting parameter values at ν = 1.4 GHz (see in particular their Table 1), which we adopt here.

It is worth mentioning that at frequencies ν ≲ 1 GHz the RL counts turn out to be largely dominated by steep-spectrum sources; for ν ≳ few GHz this is still true for ${S}_{\nu }\lesssim 1$ Jy, while at higher fluxes flat-spectrum sources starts to contribute appreciably.

We compute the differential number counts in the radio band by integrating over redshift the luminosity functions above

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{d\mathrm{log}{S}_{\nu }\,d{\rm{\Omega }}}({S}_{\nu })=\int {dz}\,\displaystyle \frac{{dV}}{{dz}\,d{\rm{\Omega }}}\,\displaystyle \frac{{dN}}{d\,\mathrm{log}{L}_{\nu }}({L}_{\nu (1+z)},z),\end{eqnarray} \tag{ 18 }$$

where the flux is given by

$$\begin{eqnarray}&&{S}_{\nu }=\displaystyle \frac{{L}_{\nu (1+z)}\,(1+z)}{4\pi \,{D}_{L}^{2}(z)}\end{eqnarray} \tag{ 19 }$$

in terms of the cosmological volume per unit solid angle ${dV}/{dz}\,d{\rm{\Omega }}$ and of the luminosity distance D_L(z). The redshift distribution is the integrand of the previous expression, in turn integrated over the luminosities above the one corresponding to a lower flux limit ${S}_{\nu ,\mathrm{lim}}$ via Equation (19).

We take into account strong galaxy–galaxy lensing of SFGs by using the amplification distribution ${dp}/d\mu$ derived in Lapi et al. (2012). The lensed differential counts are obtained as follows:

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{\mathrm{lens}}}{d\mathrm{log}{S}_{\nu }\,d{\rm{\Omega }}}({S}_{\nu })\\ &&\quad =\displaystyle \int {{dz}}_{s}\,\displaystyle \frac{1}{\langle \mu \rangle }\,{\displaystyle \int }^{{\mu }_{\max }}d\mu \,\displaystyle \frac{{dp}}{d\mu }\,\displaystyle \frac{{dN}}{d\mathrm{log}{S}_{\nu }\,d{\rm{\Omega }}}({S}_{\nu }/\mu ),\end{eqnarray} \tag{ 20 }$$

where the maximum amplification μ_max ≈ 25, as appropriate for extended sources of a few kiloparsecs, is adopted; the factor $\langle \mu \rangle$ at the denominator can be approximated to 1 in the case of large-area surveys, as considered here.

It will be of fundamental importance to test our expectations of Section 3 concerning the radio emission from SFGs and radio AGNs in different stages of a massive galaxy's evolution, by looking at large samples of radio sources with multiband coverage (e.g., X-ray, far-IR, and radio). We present below three specific examples.

The first is focused on the Eddington ratio distributions of the supermassive BHs hosted in SFGs, RQ AGNs, and RL AGNs. In Figure 14 we illustrate the schematic evolution with galactic age of the Eddington ratio λ, based on the star formation and BH accretion histories presented in Figure 2, and the corresponding λ-distributions. We expect the BHs hosted in SFGs to show a quiet narrow Eddington ratio distribution centered around (mildly super-) Eddington values λ ≳ 1. The same holds for RS AGNs, with a distribution slightly offset toward smaller values of λ, which after the peak of BH activity starts decreasing. RQ AGNs are expected to feature a much broader distribution skewed toward λ ≲ 0.3, the value allowing the development of a thin-disk accretion suitable for nuclear radio emission (see Section 3). Finally, low-z steep-spectrum RL AGNs, featuring low accretion rates, should feature a distribution shifted toward very small values λ ≲ 0.01; however, given the lack of a detailed physical understanding of their radio emission processes, we do not attempt definite predictions for this class.

Zoom In Zoom Out Reset image size

The redshift dependence of the distributions for the different classes from z ∼ 0 to 2 is mild and dictated by the redshift evolution of the Eddington ratio in the early stages of the evolution, which is required by independent data sets on the AGN luminosity functions (see discussion in Section 2). It is remarkable that our predictions for RQ AGNs are in agreement with the observational determinations, though within large uncertainties, by Panessa et al. (2015) and Padovani et al. (2015). This adds further validation to our overall picture for the radio emission from SFGs and RQ AGNs. A further, important but challenging test would consist in observationally determining the λ-distribution for pure SFGs hosting an RS AGN (X-ray detected) and confronting the outcome with our predictions for this population.

The second example is focused on the locus occupied by SFGs and AGNs on the main-sequence diagrams: SFR versus stellar mass, SFR versus X-ray luminosity, and ratio of X-ray luminosity to SFR versus stellar mass. In Figure 15 we place in such diagrams at z ∼ 2 the data with radio information by Padovani et al. (2015), highlighting the population of SFGs, RQ AGNs, and RL AGNs. Other data sets referring to mass/far-IR-selected galaxies (Rodighiero et al. 2015), X-ray-selected AGNs (Stanley et al. 2015), mid-IR-selected AGNs (Xu et al. 2015), and optically selected quasars (Netzer et al. 2016) are also reported for completeness, although we note the caveat that the detection thresholds in SFR and X-ray luminosity are slightly different among them and with respect to the Padovani et al. sample.

Zoom In Zoom Out Reset image size

In addition, we show three typical evolutionary tracks based on the star formation and BH accretion histories presented in Figure 2, which correspond to values of SFR ${\dot{M}}_{\star }\approx 30$, 300, and 1000 M_⊙ yr⁻¹ at the time when the AGN activity peaks. The shaded area shows the average relationships computed as in Mancuso et al. (2016b), taking into account the number density of galaxies and AGNs and the relative time spent by individual objects in different portions of the evolutionary tracks. We expect galaxies without signs of nuclear activity to be rather young objects, featuring stellar masses appreciably smaller than implied by the average relationship at given SFR (corresponding to ages ≲  few 10⁸ yr), and X-ray luminosities L_X ≲ 10⁴² erg s⁻¹ mostly dominated by star formation. On the other hand, we expect RQ AGNs to be more evolved objects with stellar masses lying closer to the average relationship at given SFR (corresponding to ages ≲10⁹ yr). They are also expected to host an X-ray-detectable AGN with L_X ≳ 10⁴² erg s⁻¹, RS when the star formation is still sustained, and progressively radio-active when star formation is in the way of getting quenched. Finally, low-z steep-spectrum RL AGNs are hosted mainly by galaxies in passive evolution, so that their star formation activity is often undetected at all.

Our expectations on SFGs and RQ AGNs are indeed consistent with the current multiwavelength data with radio classification from Padovani et al. (2015). In some detail, most of the objects classified as SFGs have young ages ≲  few 10⁸ yr and lie to the left of the average main-sequence relationship at a given SFR; moreover, they have only upper limits on L_X and on L_X/SFR ratios. Contrariwise, most of the objects classified as RQ AGNs have ages of several 10⁸ yr, stellar masses consistent with the average main-sequence relationship, X-ray luminosities L_X ≳ 10⁴² erg s⁻¹, and L_X/SFR ratios appreciably higher than for SFGs.

At z ≲ 0.6 this picture has been partially confirmed over the large Herschel-ATLAS fields by Gurkan et al. (2015). However, at z ≳ 2 the limited area around 0.3 deg² of the Padovani et al. sample does not allow us to populate with considerable statistics the diagrams at stellar masses ≳10¹¹ M_⊙ and X-ray luminosities L_X ≳ 10⁴⁴ erg s⁻¹. We expect to find there an appreciable number of objects classified as RQ AGNs (but not of SFGs) with at least partially quenched SFR. It would be interesting to check such a trend with multiband data of comparable quality on larger areas, as is the case for the Prandoni et al. (2017) WSRT observations over 6.6 deg² in the framework of the LH Project, and as will become routinely possible with the advent of SKA and its precursors.

The third example is focused on disentangling the relative contribution from the active nucleus and from large-scale star formation to the radio emission of individual RQ AGNs. In Figure 16 we illustrate the locus occupied by RQ AGNs in a diagram where the radio luminosity 1.4 GHz from the nucleus is plotted against that from star formation. First of all, we plot the data for RQ quasars at z ∼ 1 by White et al. (2017); these authors followed up an optically selected quasar sample with radio (FIRST) and far-IR (Herschel) observations to probe the relative contribution from the nucleus and from star formation to the radio emission.

Zoom In Zoom Out Reset image size

As for our predictions, we show three typical evolutionary tracks based on the star formation and BH accretion histories presented in Figure 2, which correspond to values of the peak AGN bolometric luminosities L_AGN ≈ 3 × 10⁴⁵, 10⁴⁶, and 3 × 10⁴⁶ erg s⁻¹, approximately the same range sampled by White et al. (2017). During the early stages of a galaxy's evolution, star formation is nearly constant, while the AGN luminosity is exponentially increasing, to originate a vertical track in the diagram; after the AGN luminosity peaks, both the star formation and the AGN luminosity decrease, and the galaxy moves to the left part of the diagram with a roughly flat track (the detailed shape depends on τ_AGN; see Section 2.2). The shaded area shows the average relationship computed as in Mancuso et al. (2016b), taking into account the number density of AGNs with different luminosities and the relative time spent by individual objects in different portions of the evolutionary tracks.

Our expectations are in good agreement with the results from White et al. (2017) for RQ quasars, which tend to cluster close to the peaks of the individual evolutionary tracks, and lie within the average relationship (with its scatter) predicted by Mancuso et al. (2016b). All in all, for most of them the radio emission is found to be mainly contributed by the active nucleus. It would be interesting to test further our predictions (e.g., to independently constrain the timescale τ_AGN; see above) with larger samples spanning a wider luminosity range and attaining a higher sensitivity in the far-IR and radio bands.