© ESO, 2015

1. Introduction

Radio-loud active galactic nuclei (AGNs) are a key component driving galaxy evolution; the feedback their expanding radio jets provide is essential for preventing large-scale cluster cooling flows and halting the growth of massive elliptical galaxies (e.g. Fabian et al. 2006; Best et al. 2006, 2007; Croton et al. 2006; Bower et al. 2006). To understand the timescales upon which these processes occur, it is important to first understand the evolution of the radio luminosity function (RLF) to high-redshift. An early measurement of this came from Dunlop & Peacock (1990), who found a higher value for the space density (by two to three orders of magnitude) of both flat and steep-spectrum radio AGN than is seen locally. They also saw the first indication of an expected higher redshift density decline at z ~ 2.5, corresponding to the build-up of these objects in the early Universe. However, their work, and that of subsequent studies (e.g. Shaver et al. 1996; Jarvis et al. 2001; Waddington et al. 2001), lacked the necessary depth and volume needed to unambiguously measure this high-redshift behaviour.

This situation improved with the development of the Combined EIS-NVSS Survey of Radio Sources (CENSORS; Best et al. 2003), a survey designed to maximise the coverage of steep-spectrum radio sources close to the high-redshift break in the RLF. Rigby et al. (2011, hereafter R11) used CENSORS, combined with additional radio source samples, source counts, and determinations of the local RLF to investigate the space density evolution of the P_{1.4 GHz} > 10²⁵ W/Hz steep-spectrum population via grid-based modelling with no prior assumptions included about the high-redshift behaviour. This robustly identified the post-peak space density decline in the RLF, and found that this turnover appears to be luminosity-dependent; at lower radio powers (P_{1.4 GHz} = 10²⁵⁻²⁶ W/Hz) the space densities peak at z ≳ 1, but the peak moves to higher redshift for the more luminous objects (z ≳ 3 for P_{1.4 GHz} > 10²⁷ W/Hz). A luminosity dependence in the position of the steep-spectrum RLF peak can be interpreted as a sign of cosmic downsizing, where the most massive black holes form at earlier epochs than their less massive counterparts. This has also been seen for other AGN populations, selected at other radio, optical, far-infrared, and X-ray wavelengths (e.g. De Zotti et al. 2010; Hasinger et al. 2005; Richards et al. 2005; McAlpine et al. 2013; Delvecchio et al. 2014), and also reproduced in simulations of black hole growth (e.g. Fanidakis et al. 2012; Hirschmann et al. 2012, 2014).

Steep-spectrum radio sources can be split into two distinct populations: typically powerful objects with standard accretion of cold gas and likely to be merger driven (cold-mode); and typically weaker sources with radiatively inefficient accretion directly from their hot gas haloes (jet-mode or hot-mode). Locally, the latter dominate below P_{1.4 GHz} ~ 10²⁶ W/Hz, but both classes are present at all powers (Best & Heckman 2012). The major limitation of R11 was the lack of constraint at lower radio powers (P ≲ 10²⁵ W/Hz) out to z ≳ 2. The R11 power range was thus only able to probe the RLF regimes where the cold-mode sources are in the majority, and therefore we could not draw any conclusions about the relative contributions of the two classes across cosmic time. This paper builds and expands on this initial foundation by including radio powers an order of magnitude fainter than R11, allowing the RLF evolution to be measured across 24 ≤ log P_{1.4 GHz} ≤ 28, to high redshift. This improved coverage means that the apparent luminosity-dependence of the peak space density can be investigated for the first time from the brightest sources down to the regime where the dominant steep-spectrum population is shifting between the two types.

The layout of the paper is as follows. Section 2 describes the modelling method used to determine the best-fitting steep spectrum RLFs. Section 3 presents the SXDF radio source sample. The results of incorporating this into the modelling process are given in Sect. 4. Section 5 discusses the luminosity-dependence of the position of the RLF peak. Finally, the results are summarised and interpreted in Sect. 6. Throughout this paper values for the cosmological parameters of H₀ = 70 km s^-1 Mpc^-1, Ω_m = 0.3, and Ω_Λ = 0.7 are used and the radio spectral index α is defined as S_ν ∝ ν^−α.

2. Modelling process

The modelling method in R11 makes no assumptions about the shape and evolution of the RLF. It takes as input five radio galaxy redshift distributions, measurements of the radio source counts over 0.05 mJy ≤ S_{1.4 GHz} ≤ 94 Jy, and determinations of the RLF. The radio samples are selected such that they provide a good coverage of the radio luminosity – redshift plane out to z ~ 4 (Fig. 1).

These input data are used to constrain co-moving space densities, ρ, determined at various points on a P-z grid of radio luminosities and redshifts from which the cosmic evolution of the radio galaxy population can be read directly. This grid is made up of three components, steep-spectrum, flat-spectrum, and star-forming sources, that represent the different radio source populations present. Since the steep-spectrum component is the dominant population at the luminosities and redshifts considered in R11, it is the only one allowed to vary in the modelling process (see Sect. 3 for further discussion). The constant star-forming grid is created by evolving the local star-forming galaxy luminosity function of Sadler et al. (2002) by P/ (1 + z)^2.5 to z = 2 and by 0.06P at z> 2; the constant flat-spectrum grid is taken as the median of the set of evolutionary models presented in Dunlop & Peacock (1990) for this source population. The starting point for the steep-spectrum grid is formed by evolving the local AGN RLF (Sadler et al. 2002) by (1 + z)³ in density.

The grids cover a range of 19.25 ≤ log P_{1.4 GHz} ≤ 29.25 in radio luminosity, equally spread in steps of Δlog P_{1.4 GHz} = 0.5, and 0.1 ≤ z ≤ 6.0 in redshift, evaluated at z = 0.1,0.25,0.5,1.0,2.0,3.0,4.0,6.0. The space density at any intermediate P-z grid point can simply be interpolated from its four neighbours. Figure 1 illustrates the distribution of grid points, together with the corresponding constraints offered by the input data. Unconstrained densities are excluded from the fitting process.

The best-fitting space density grid is determined using the amoeba algorithm for downhill simplex minimisation (Nelder & Mead 1965), run in a multi-stage loop with varying scaling and tolerance parameters for the initial steps. The goodness of fit for a particular model grid is the combined likelihood of the predicted source counts, redshift samples, and local luminosity functions compared to the real data. The uncertainties in the model grid are taken as the marginalized error calculated using the inverse Hessian matrix.

Full details of both the modelling process and the input datasets used can be found in R11.

3. New input data: the Subaru/XMM-Newton Deep Field radio sample

The faintest radio galaxy sample used in R11 was the AGN subsample of the VLA-COSMOS survey (Smolčić et al. 2008), with a 1.4 GHz flux density limit of 100 μJy. This dataset was only complete to z ≤ 1.3, which – since the next faintest sample had a 2 mJy flux density limit – severely limited the high-redshift constraint in the P-z grid at low radio powers. This restricted the RLF evolution analysis to log P_{1.4 GHz} ≥ 25. To improve this constraint a radio survey of comparable depth but better redshift coverage, carried out in the Subaru/XMM-Newton Deep Field (SXDF; Simpson et al. 2006, 2012), is now included in the modelling. This covers 0.81 sq. deg. and contains a complete sample of 505 sources with S_{1.4 GHz} ≥ 100 μJy. Robust redshifts, both spectroscopic (51%) and photometric, exist for 99% of the sample. The photometric redshifts for the remaining seven sources are contaminated by foreground objects and are therefore excluded here.

At this flux density level ≳50% of the sources are expected to be star-forming galaxies (e.g. Norris et al. 2013). Since this paper is concerned with the behaviour of the steep-spectrum population (which comprises the majority of the remainder), the catalogue is restricted to radio-loud AGNs only, as described below. In principle, this step is not necessary as the evolution of the star-forming RLF in the modelling process should be accounted for by the treatment applied to the star-forming grid. However, uncertainties in the evolution parameters used for this, combined with the non-negligible contribution of this population in this survey, means that it is prudent to perform this step to limit the effect of any errors.

The AGN/star-formation separation is done by calculating the ratio of mid-infrared, S_{24 μm}, to radio flux density via the k-corrected q parameter: q = log (S_{24 μm}/S_{1.4 GHz}). Following Donley et al. (2005), AGNs are defined as sources with q< 0, indicating excess radio emission over that from star-formation. Applying this cut results in a final sample of 179 sources. The q parameter adopted is higher than the q = −0.23 value used previously to select radio-loud objects (e.g. Ibar et al. 2008). However, inspection of Fig. 6 in Simpson et al. (2012) shows that the higher cut is more appropriate for this sample. In practice, the difference in the modelling results from changing q is negligible within the errors.

The AGN cut also removes radio quiet quasars from the sample. These also begin to be a larger fraction of the radio population at S_{1.4 GHz} < 1 mJy, and display radio emission that may arise mainly from star-formation instead of black hole accretion, although this is still a topic of active debate (Smolcic et al. 2015, and references therein). Their contribution to the overall population is not explicitly included here, but rather assumed to be covered by the star-forming grid.

The redshift distribution of the sample is shown in Fig. 2, together with the 12th order polynomial used to represent it in the modelling process. This step is necessary as the low number of sources at high-redshift makes a direct calculation of the χ² difficult (see R11 for further details).

The VLA-COSMOS sample previously used in the modelling process is now replaced here by the SXDF sample described above.

4. Model results

4.1. Dataset comparison

The independent redshift distribution generated from the best-fit steep-spectrum P-z grid is shown in Fig. 2; this agrees well with the real dataset at z> 0.3. Below this the SXDF might miss resolved sources (Simpson et al. 2012), which may explain the apparent overprediction of the model at low redshift. The figure also shows the small relative contribution to the total population from the un-fitted flat-spectrum grid.

The total number of AGN (i.e. flat- plus steep-spectrum sources) given by the model is 199 ± 31 (188 ± 29 at z ≥ 0.3), which compares well to the actual figure of 179 (178 at z ≥ 0.3). In addition, the predicted number of star-forming galaxies (377) broadly agrees with the 326 observed in the SXDF sample and suggests that the treatment of this population in the modelling is sufficient and does not give rise to additional uncertainty in the grid results.

The agreement between the model and the other input datasets is also good; however, these comparisons are not explicitly shown here as they are essentially unchanged from the results presented in R11.

4.2. Model RLFs

The predicted variation of space density with redshift for eight radio power bins (24 ≤ log P ≤ 28) is shown in Fig. 3. We note that the ρ for a particular P-z point is only shown if it is constrained by both the source counts and the redshift distributions (Fig. 1).

The main benefit of including the SXDF sample in the modelling comes at the lowest powers (log P< 25.5), where it allows the RLF behaviour to be studied at z> 1. Inspection of Fig. 3 shows that decreases in space density are also present here, though these tend to occur at lower redshifts than for higher power ranges; this is discussed further in Sect. 5 and is consistent with the recent results of Padovani et al. (2015).

As part of the modelling process the original P-z grid is converted into a grid of flux density and redshift to allow easier comparison with the source counts and redshift distributions. To make this S-z grid a spectral index, α, is needed; following R11, this is assumed to vary with redshift as α = 0.83 + 0.4log (1 + z) (Ubachukwu et al. 1996), with an additional uncertainty of α ± 0.2 at each redshift to account for any inter-bin variations. Figure 3 shows the range of best-fit results returned when this assumption is changed to either a constant but very steep spectrum (α = 1.5), a constant value of 0.8 (the mean value at low redshift), or a stronger increase with redshift (α = 0.8 + 0.25z). The general effect of using the steeper α values is to increase the overall density values and weaken, though not remove, the z ≳ 1 space density turnovers as sources move to bins in the S-z grid corresponding to higher radio powers. The assumption of the shallower value of α = 0.8 has the opposite effect. The small spread in results seen at 26 < log P< 26.5 following this α variation suggests that this is the radio power bin that is most strongly constrained in the modelling process.

5. Luminosity-dependence of the high-redshift turnover

Inspection of Fig. 3 suggests that, as in R11, the position of the space density peak is luminosity dependent. This can be clearly seen if the model steep-spectrum RLFs are shown together, as in Fig. 4. The addition of the SXDF sample as a constraint to the RLF modelling allows this behaviour to be probed to lower radio powers than previously possible.

The position of the peak redshift in each radio power range is determined by fitting a polynomial (generally of order 2, but of order 3 where necessary) to the model space density distributions in the P-z grid. This overcomes the ambiguity arising from the discrete nature of the P-z grid, and the wide redshift bins used. Figure 5 shows the results of this for the best-fitting steep-spectrum grid presented in Sect. 4.2, together with the spread arising from repeating this process for the grids determined using altered spectral indices. The small variation in values at log P = 26.25 supports the previous identification of 26 < log P< 26.5 as the most constrained radio power bin in the model grid.

The luminosity-dependence of the peak redshift, z_peak, is clear with the most powerful radio sources peaking at earlier times than their weaker counterparts. However, the extension to lower radio powers reveals an apparent flattening of the relationship: below log P ≲ 26.5 the results are consistent with a constant value of z_peak.

5.1. Comparison with quasar behaviour

The luminosity-dependent behaviour of the RLF peak redshift seen here for steep-spectrum radio sources has also been observed in other AGN populations (e.g. De Zotti et al. 2010; Hasinger et al. 2005; Richards et al. 2005; Delvecchio et al. 2014). A direct comparison of the z_peak changes is possible for flat-spectrum quasars, using the luminosity function (QLF) models of Hopkins et al. (2007); these are derived from a large set of observed QLFs, from the mid-infrared to X-ray.

The optical B-band (0.44 μm) is chosen for extracting the quasar behaviour, for ease of comparison to the luminosities in the radio waveband. It should be noted, however, that the variation in peak redshift calculated using the bolometric QLF instead is qualitatively similar. The Hopkins et al. model is queried for a luminosity range over which it is reliable (20 < log P_B−band < 27), and the z_peak value determined in each case. The subsequent result is shown in Fig. 6, together with that already found for the steep-spectrum sources. The ratio of the two curves is also shown as a function of peak redshift.

The shape of the z_peak−log P_{B − band} curve derived for the quasar population is broadly similar to that of the radio galaxies, with indications of flattening at 24 < log P_ν < 26, before continuing to decrease at lower powers. The log (P_{1.4 GHz}/P_{B − band}) ratio begins to increase from a constant value of ~1.5 at z> 3, within the uncertainties; it is interesting to note that log (P_{1.4 GHz}/P_{B − band}) ≥ 1.5 is the historical definition of a radio-loud quasar (Kellermann et al. 1989).

5.2. The relation to cosmic downsizing

The obvious physical interpretation of the luminosity-dependent space density peak seen here is cosmic downsizing, in which the most massive black holes have formed by z ~ 4, with the less massive forming at later times. This apparently contradictory anti-hierarchical growth can be explained if the dominant mechanism of AGN fuelling changes with cosmic time from cold gas accretion via major mergers (cold-mode) to radiatively inefficient accretion directly from hot gas haloes (jet-mode or hot-mode). Simulations of AGN evolution support this picture; they show that hierarchical black hole growth is able to produce the observed downsizing trend if inefficient accretors come to dominate at low redshift (Fanidakis et al. 2012; Hirschmann et al. 2014).

Locally, hot-mode sources make up the majority of the radio-galaxy population at log P_{1.4 GHz} ≲ 26 (Best & Heckman 2012), though both types are generally found over the range of radio powers studied. In comparison, by z ~ 0.7 the space densities of the two are similar at low radio powers, with the cold-mode sources becoming more dominant at log P_{1.4 GHz} ≳ 26 (Best et al. 2014). The flattening seen here in the z_peak vs. log P_{1.4 GHz} relation occurs at corresponding radio power to this transition, lending further support to the assertion that the RLF evolution seen here directly represents the interplay between the two populations. It also suggests that the weaker, likely hot-mode, sources undergo little evolution in space density; this finding is broadly replicated in Best et al. (2014), which directly catalogued sources into the two classes.

This assertion can be tested by predicting the dual population space density evolution, taking as a starting point the well-defined local RLFs for the two classes (Best & Heckman 2012). For the evolution of the hot-mode population, we follow Best et al. (2014) in assuming that their evolution broadly follows that of the space density of massive quiescent galaxies, which constitute their host galaxies. Best et al. (2014) showed that the space density of these galaxies remains relatively flat out to z ~ 0.8, and then rapidly decreases as (1 + z)^-6.5. We modify this to allow for an increase in the space density as (1 + z) out to the peak redshift at z = 0.8, as this is a better match to the evolution of the 10²⁴ <P_{1.4 GHz} < 10²⁵ RLFs in the model grid.

For the cold-mode population, we determine the luminosity-dependent peak redshift by comparison with the optical quasars. Figure 6 suggests that there is a roughly constant ratio between the luminosities of radio galaxies and quasars at high values of z_peak where the cold-mode sources dominate. Assuming this is true for this population, z_peak for a particular radio luminosity is therefore taken as the peak of the corresponding B-band quasar luminosity. For this population, we allow a more rapid pre-peak increase in space density of (1 + z)³, which matches the evolution seen at the highest radio luminosities in the model grid. The post-peak decrease is again taken as (1 + z)^-6.5.

The predicted RLFs for the two populations are then combined, sampled using the same luminosity and redshift bins as in the model grid, and the peak redshift extracted.

The variation in peak redshift with radio luminosity arising from this simple model is shown in Fig. 5. It does well at reproducing the overall behaviour derived from the best-fit model grid. At the faintest (log P_{1.4 GHz} < 24.5) luminosities the predicted RLFs are dominated by the hot-mode sources, whilst the cold-mode sources dominate at log P_{1.4 GHz} > 26.5. In between, the RLFs are a combination of the two classes. This is reflected in the changing shape of the z_peak − log P_{1.4 GHz} curve: in this scenario the plateau at 24.5 < log P_{1.4 GHz} < 26.5 arises from the transition from the cold-mode to the hot-mode populations.

6. Summary and conclusions

The results presented in this paper have demonstrated the benefits of extending the P–z plane coverage to lower radio powers when measuring the evolution of the steep-spectrum RLF. The addition of the updated SXDF sample to the constraining datasets used in the modelling process shows that the high-redshift turnover in the RLF persists to log P_{1.4 GHz} > 24 – an order of magnitude fainter than could be measured previously.

The redshift at which the peak space density occurs is found to be luminosity-dependent at log P_{1.4 GHz} ≳ 26, with the most powerful sources peaking at earlier times than their weaker counterparts. This mirrors what was seen when considering quasar optical luminosity functions, which is to be expected as the radio population is dominated by cold-mode (i.e. quasar-like) sources at these luminosities. The presence of a mean ratio for these powers (log [ P_{1.4 GHz}/P_{B − band} ] ~ 1.5), which remains constant with z_peak, again suggests that the radio and optical observations are seeing the same objects. This picture changes at lower radio luminosities (log P_{1.4 GHz} ≲ 26). Here the value of the peak redshift remains constant within the uncertainties, whereas it continues to fall for the optical population. This implies an increase in the P_{1.4 GHz}/P_{B − band} ratio and suggests that the dominant radio population is changing from cold-mode to hot-mode.

Simple models for the space density evolution of the hot-mode and cold-mode populations are able to reproduce the observed luminosity dependence in the peak redshift, reinforcing the conclusion that this behaviour arises from the transition between the two populations. To properly understand the relative contributions of the hot-mode and cold-mode sources to the RLF evolution, their behaviour should be minimised separately in the fitting. This is easy to do with this grid-based modelling method, but the input data needed are currently lacking across the required redshift range. Future large, deep radio surveys in regions with extensive multiwavelength coverage will make this possible, whilst also allowing the coverage of the P-z plane to extend to radio powers of log P_{1.4 GHz} < 24.

Acknowledgments

E.E.R. acknowledges financial support from NWO (grant number: NWO-TOP LOFAR 614.001.006). E.E.R. and J.A. thank the Leiden/ESA astrophysics program for summer students (LEAPS) which supported J.A. in Leiden. The authors would like to thank Chris Simpson for sharing the q parameters in the SXDF catalogue.

References

All Tables

All Figures