The XXL survey LII : The evolution of radio AGN luminosity function determined via parametric methods from GMRT, ATCA, VLA and Cambridge interferometer observations

A set of three wide shallow surveys were the $7$C, $6$CE and $3$CRR fields, observed with the Cambridge Low-Frequency Synthesis Telescope and the Cambridge Interferometer. See Willott et al. (2001) and references therein for details about the surveys, which we briefly summarize below.

The $7$C field was observed at frequencies of $151\ \mathrm{MHz}$, detecting sources above $0.5\ \mathrm{Jy}$ (Willott et al. 2001). The complete area of the survey, which consists of three distinct regions: $7$C-I, $7$C-II and $7$C-III, equalls $72.22\ \mathrm{deg}^{2}$ (i.e. $0.022\ \mathrm{sr}$). The spectral indices were determined using multifrequency radio data for $7$C-I, $7$C-II, and $38\ \mathrm{MHz}$ $8C$ data for $7$C-III (Lacy et al. 1999), resulting altogether in a mean spectral index of $\alpha\approx-0.64$. The redshift information was derived from follow-up optical and near-infrared observations (Willott et al. 2001). Most of the redshift were determined spectroscopically ($\approx 85\%$) while the remaining sources have photometric redshifts. The number of sources in the $7$C catalogue equals $128$.

The $6$CE survey at $151\ \mathrm{MHz}$ covered an area of $\approx 340\ \mathrm{deg}^{2}$ ($0.103\ \mathrm{sr}$) capturing sources with flux density $2\ \mathrm{Jy}<S_{151\ \mathrm{MHz}}<3.93\ \mathrm{Jy}$ (Rawlings et al. 2001). The catalogue contained $59$ sources, of which all but three have spectroscopically determined redshifts. The mean spectral index, obtained by a polynomial fit to the multi-frequency data (Rawlings et al. 2001) equaled $\alpha\approx-0.51$. For more details on the catalogue see Rawlings et al. (2001).

The $3$CRR catalogue, from observations at $178\ \mathrm{MHz}$ spans an area of $\approx 13,900\ \mathrm{deg}^{2}$ ($4.23\ \mathrm{sr}$), with a detection limit of $10.9\ \mathrm{Jy}$. The redshift information is present for all $173$ sources in the sample as describd in Willott et al. (1999). The spectral index, calaculated at rest frame $151\ \mathrm{MHz}$ had a mean of $\alpha\approx-0.67$.

The $1.4\ \mathrm{GHz}$ rest-frame radio luminosities for the $7$C, $6$CE and $3$CRR used here were computed from flux, redshift and spectral index values given in the corresponding catalogs, using a newer cosmology defined in Sec. 1 ($H_{0}=70\ \mathrm{kms^{-1}Mpc^{-1}}$, $\Omega_{m}=0.3$, and $\Omega_{\Lambda}=0.7$).

The observations of the XXL-North field were performed at $610\ \mathrm{MHz}$ with the Giant Metrewave Radio Telescope (GMRT). The observations were divided into two distinct parts: the inner part of the field spanned an area of $11.9\ \mathrm{deg}^{2}$ (the XMM-Large Scale Structure, XMM-LSS field). The data were taken from an earlier study by Tasse et al. (2007), and then re-reduced as discussed in Šlaus et al. (2020) (hereafter XXL Paper XLI). These observations reach a mean rms of $200\ \mathrm{\mu Jy\ beam^{-1}}$. The remaining $18.5\ \mathrm{deg}^{2}$, were observed by XXL Paper XXIX. They have a mean rms of $45\ \mathrm{\mu Jy\ beam^{-1}}$. The number of observed sources in both parts of the field was $5,434$, using a signal-to-noise ratio of $S/N\geq 7$. The spectral indices for both parts of the field were obtained by matching the catalogue with another radio catalogue, the NRAO Very Large Array Sky Survey (NVSS) at $1400\ \mathrm{MHz}$ (Condon et al. 1998) as described in XXL Paper XXIX. The mean spectral index equaled $-0.65$ for the inner part of the field, and $-0.75$ for the outer, due to the difference in survey depth. The $157$ multi-component sources were created from components as described in XXL Paper XLI. More details on the observations and the corresponding catalog can be found in XXL Paper XXIX.

In order to obtain the source redshifts, the catalogue was cross-matched with a multiwavelength catalogue from Fotopoulou et al. (2016), hereafter XXL Paper VI, using only the subset of the catalog that has identifications in the Spitzer Infrared Array Camera (IRAC) Channel 1 band at $3.6\ \mathrm{\mu m}$ (PI M. Bremer, limiting magnitude of 21.5 AB), to obtain uniform depth. The redshifts of the IRAC-detected sources were determined photometrically using the full multi-wavelength data (Fotopoulou, in prep.). Details on the redshift accuracy can be found in XXL Paper XLI. The IRAC survey covered roughly $80\%$ of the radio field. After further removing the noisy edges of the radio map, the area of the inner part of the field equaled $6.3\ \mathrm{deg}^{2}$ and the area of the outer $14.2\ \mathrm{deg}^{2}$.

The $25\ \mathrm{deg}^{2}$ of the XXL-South field were observed with the Australia Telescope Compact Array (ATCA), at $2.1\ \mathrm{GHz}$ (XXL Paper XVIII). The observations reached a depth of $\approx 41\ \mathrm{\mu Jy\ beam^{-1}}$. The details of the observations are described in XXL Paper XVIII. The catalogue consists of $6,239$ single component sources and an aditional $48$ sources composed of multiple components. The spectral indices were determined by matching the $2.1\ \mathrm{GHz}$ catalogue with the Sydney University Molonglo Sky Survey (SUMSS) at $843\ \mathrm{MHz}$ (Bock et al. 1999) reaching sources with peak flux density of $6\ \mathrm{mJy}$. After taking into account the bias that arises from a high detection limit of the SUMSS survey, the median spectral index was estimated at $\alpha\approx-0.75$.

The catalogue was furthermore matched with a multiwavelength catalogue using a likelihood technique as described in Ciliegi et al. (2018) (hereafter XXL Paper XXVI). The multiwavelength catalogue contained data from near-infrared and optical up to X-ray data (for details see XXL Paper VI). The cross-correlation process resulted in $4,770$ optical/NIR counterparts, with $414$ of them also detected in the X-ray band (XXL Paper XXVI). Since $12$ of these sources were classified as stars, they were removed from the sample, resulting in a catalogue of $4,758$ sources (Butler et al. 2018b, hereafter XXL Paper XXXI). There are $1,110$ spectroscopic redshifts and $3,648$ photometric redshifts listed in the catalogue (XXL Paper XXXI, XXL Paper VI). The details concerning the accuracy of the photometric redshifts and the overall redshift distribution of the sample can be found in XXL Paper XXXI. The median spectral index of the matched catalogue is flatter and equals $-0.45$ (XXL Paper XXXI).

The deepest radio survey used in this works is the VLA-COSMOS $3\ \mathrm{GHz}$ Large Project (Smolčić et al. 2017a). The area of the observed field also covered by multiwavelength data, is $2\ \mathrm{deg}^{2}$ and the detection limit at $5\sigma$ equaled $11.5\ \mathrm{\mu Jy\ beam^{-1}}$. The full catalogue contained $10,830$ sources of which $67$ are multi-component. The spectral indices were derived from cross-correlation with the $1.4\ \mathrm{GHz}$ Joint catalogue by Schinnerer et al. (2010), using survival analysis to account for the bias of different detection limits, as described in Smolčić et al. (2017a). The mean spectral index was estimated at $\alpha=-0.73$.

The radio catalogue was further matched with a multiwavelength catalogue as described in Smolčić et al. (2017). This resulted in $\approx 93\%$ of the sources obtaining a counterpart ($8035/8696$ in the unmasked part of the field). For $7778$ of these sources there exists a redshift estimate. Of these, $2740$ are spectroscopic ($\approx 34\%$), while the other $5123$ are photometric. For details on the redshift estimation see Delvecchio et al. (2017) and Smolčić et al. (2017).

Another study used in this work is the Combined EIS–NVSS Survey Of Radio Sources (CENSORS) by Brookes et al. (2008), containing $150$ radio sources. The catalogue comes from the spectroscopic observations of the $1.4\ \mathrm{GHz}$ sources observed by Best et al. (2003a). The area of observations equalled $6\ \mathrm{deg}^{2}$, and the detection limit reached $7.2\ \mathrm{mJy}$. The catalogue contained spectroscopic redshifts for $71\%$ of the sample, while the remaining redshifts were determined via the K-z relation, as described by Brookes et al. (2008). The spectral indices were lacking and were set to a constant value of $-0.8$, following the assumptions made in Brookes et al. (2008). Four sources from the sample whose radio emission comes from star formations were classified by Brookes et al. (2008). The rest of the sample consists of AGN sources.

The Best, Röttgering $\&$ Lehnert (BRL) sample is defined from the Molonglo Reference Catalogue, with the spectroscopic data compiled and observed by Best et al. (1999), Best et al. (2000), Best et al. (2003a) and Best et al. (2003b). The radio observations were performed at $408\ \mathrm{MHz}$, with a detection limit of $5\ \mathrm{Jy}$. The area of observations was bounded by declination $\delta\in[-30^{o},10^{o}]$, excluding the Galactic plane $|b|>10^{o}$

The sample by Wall & Peacock (1985) covers an area of $\approx 32,000\ \mathrm{deg}^{2}$ ($9.81\ \mathrm{sr}$). The observations were performed via ANRAO/Parkes, NRAO/Greenbank, and MPIfR/Bonn. The frequency of observations equalled $2.7\ \mathrm{GHz}$, with a detection limit of $2\ \mathrm{Jy}$ at this frequency. The catalogue consists of $233$ sources, with $171$ of them ($73\ \%$) having measured redshifts. The remaining redshifts were estimated from $V$-band magnitudes, as described by Wall & Peacock (1985). The spectral indices of the sources were determined from additional $5\ \mathrm{GHz}$ data. Given the high detection limit of the survey, we assumed it consists purely of AGNs.

The Combined NVSS-FIRST Galaxies (CoNFIG) sample (Gendre & Wall 2008) comes from observations at $1.4\ \mathrm{GHz}$, selecting NRAO Very Large Array (VLA) Sky Survey (NVSS) sources from the north field of the Faint Images of the Radio Sky at Twenty centimetres (FIRST) survey, spanning $\approx 5000\ \mathrm{deg}^{2}$ ($1.5\ \mathrm{sr}$). The catalogue contains $274$ sources, selected above a detection limit of $1.3\ \mathrm{Jy}$. Redshifts are determined for $89\%$ of the sample. The redshifts are mosltly spectroscopic ($230$ sources), and some determined via the $R-z$ relation ($14$ sources), as described by Gendre & Wall (2008). The spectral indices present in the catalogue were determined from observations at lower frequencies, up to $178\ \mathrm{MHz}$.

A crucial step in the process of Bayesian parameter estimation is determining the likelihood function. We followed here Marshall et al. (1983) (see also Christlein et al. 2009 for a more detailed derivation). By dividing the complete luminosity-redshift space into infinitesimal cells $\mathrm{d}z\mathrm{d}L$ and assuming that each cell is small enough to contain up to one source, invoking the Poisson distribution, the probability of observing $N$ sources of the complete sample is:

where $\lambda$ is the expected number of sources per bin. The first product goes over the complete sample of $N$ sources, while the second one takes into account that all the remaining cells must remain empty. The expected density of sources in a given luminosity bin $\mathrm{d}L$ is given by the luminosity function $\phi(z,L)$:

By taking the customary logarithm of the probability and rearranging the sums we obtain:

where $\phi_{i}$ and $V_{i}$ are associated with a particular source of the catalogue. The first sum goes over the observed sample. The second sum, which was turned into an integral, covers the whole available $(z,L)$ space. The limits of the integral therefore follow the detection limit of the survey, i.e. it numbers all the cells where in principle a source could be observed. It therefore also follows that the integral equals the total predicted number of sources above the detection limit of the survey (Christlein et al. 2009). The log-Likelihood, $\ln\mathcal{L}$, is defined as (e.g. Marshall et al. 1983):

The expression can be further simplified by noting that not all terms depend on the luminosity function parameters. The first term of the last equation can be divided into:

The second term of this relation does not depend on the luminosity function parameters and, as such, provides the $\ln L$ relation with a constant value not important in the minimisation process. It can therefore be omitted. We have finally:

This expression is the one found commonly in the literature (e.g. Kelly et al. 2008, Yuan et al. 2020). Furthermore, this expression can be generalised naturally to multiple fields $j$ with different detection limits and observational areas as:

where the first sum covers all the sources from all the composite fields and each integral in the second sum reaches the depth of the corresponding field as denoted by the lower limit. If the incompleteness of the survey near the detection limit is significant, it can be included as a separate completeness function, as described in the next subsection.

The completeness corrections of each survey can be introduced naturally by using a smooth detection limit which is a function of flux, instead of an abrupt cutoff. The corrections for each survey were taken from their respective papers, as described below.

For the XXL-North field the correction is given in XXL Paper XXIX. This correction corresponds to the one arising from noise near the detection limit. In XXL Paper XLI we also introduced another correction arising from the losses during the matching of radio data with the multiwavelength catalogue which were not negligible. This correction is a function of redshift.

For the COSMOS field, the correction can be found in Smolčić et al. (2017a) in their Figure 16 or Table 02. Finally, as seen from XXL Paper XVIII and Willott et al. (2001) the other catalogues can be considered complete. Therefore for these catalogues no corrections were included.

A feature of the Bayesian formalism discussed above is the ability to compare the fit between different models. A direct comparison is obtained by calculating the odds ratio (Thrane & Talbot 2019):

where we have re-introduced the evidence from Sect. 4, $E=p(D|M,I)$ (Liddle 2007). If, furthermore, there is no model preferred by the priors, the odds ratio reduces to a ratio of evidences called the Bayes factor:

This ratio is also commonly expressed as a difference in log-scale. It should be noted that this method naturally selects models with a lower number of parameters, or in other words it incorporates Occam’s razor (Liddle 2007).

Another, more approximative, method is the Akaike information criterion (AIC; Akaike 1974, Liddle 2007, XXL Paper VI). The definition of the AIC value is:

where $k$ denotes the number of parameters. The model with the lower value of AIC is the one that corresponds to a better fit. It can be seen that this method also penalises a larger number of parameters. In other words, the Occam’s razor is included naturally within the comparison.

Lastly we have the Bayesian information criterion (BIC; Schwarz 1978, Liddle 2007, XXL Paper VI), similar to AIC, and defined as:

where $N$ is the number of data points. This value is the numerical approximation for the Bayes factor. From the above expression it is also immediately clear that for a sufficiently large $N$, the BIC penalty for a large number of parameters is stronger than that of AIC.

The shape of the local luminosity function is usually described by a power-law with an exponential cut-off (Saunders et al. 1990, Sadler et al. 2002, Smolčić et al. 2009):

where the base of the logarithm is $10$. Here $L^{*}$ is the break luminosity, $\Phi^{*}$ the normalisation and $\sigma$ the high-luminosity slope. Another possible choice would have been the double power law used often for radio and X-ray AGN samples (Dunlop & Peacock 1990, Mauch & Sadler 2007, Smolčić et al. 2017b, XXL Paper VI). We discuss our choice in more detail in the discussion. Lastly, another possibility is the bimodal model from Willott et al. (2001) which has a different form for the low and high luminosity end of the sample. This model is discussed separately in Sect. 6.3 below.

The evolution of the aforementioned local LFs is given by both density and luminosity evolution as (e.g. Smolčić et al. 2017b):

where $\alpha_{D}$ and $\alpha_{L}$ are the parameters quantifying density and luminosity evolution respectively. In this work we refer to this model as Sadler$+02$. We also tested pure density evolution (PDE) and pure luminosity evolution (PLE) where only one evolution parameter is different from zero. Another parametrization of luminosity function evolution can be found in Novak et al. (2017):

where $\beta$ parameters quantify the change of evolution with redshift. We refer to this model as Novak$+18$. In order to account for the difference in evolution between the high and low luminosity end of the sample, we also investigated the luminosity dependent density evolution (LDDE; Schmidt & Green 1983, Ueda et al. 2003) following XXL Paper VI:

where:

where $L_{a}$ is the luminosity where the evolution changes according to relation (19), $z_{c}$ redshift after which the evolution changes and $p_{1,2}$ the parameters of evolution. Lastly, as introduced by Massardi et al. (2010) and Bonato et al. (2017), the evolution of sources can be described by a luminosity dependent luminosity evolution model (LDLE). Where the break luminosity $L^{*}$ evolves as:

and:

Here $k_{evo}$ and $m_{ev}$ are free parameters of evolution, while $z_{top,0}$ and $\delta z_{top}$ determine the redshift where the evolution changes.

A special case for both the local luminosity function and its evolution is the bimodal model taken from Willott et al. (2001). We refer to it in this work as Willott$+01$. The shape and evolution of the sample have a different analytical form for the high and low luminosity end of the sample. Following Smolčić et al. (2009), in this work we used model ”C” from Willott et al. (2001), being the most flexible one, defined as:

where $\Phi_{l}$ is the low luminosity end of the function:

and $\Phi_{h}$ is the high luminosity end:

Here $L^{*}$ denotes the break luminosity, $\Phi_{0}$ the normalisation and $\alpha$ the slope of the luminosity functions. Parameter $z_{0}$ is the redshift at which the evolution changes. These parameters exist separately for the high and low luminosity end of the sample as denoted by the extra indices $h,l$. Parameters $k_{l}$, $z_{h1}$ and $z_{h2}$ quantify the evolution.

Before using the observed data discussed in Sect. 2, we tested the methodology by using simulated data. For this purpose we created custom Python codes which created catalogues of mock-observations starting from an assumed LF. We then tried to re-create the assumed LF by modeling the LFs using our methodology on simulated mock-observations. We tested both the parametric method described in Sect. 4 and the method of maximum volumes described in Sect. 5. The test were performed for a wide range of different fields and their combinations, starting with a wide range of assumed LFs. The results were always in good agreement with the starting luminosity function, within the range defined by the uncertainties defined via $90\%$ quantiles for a sample of LFs drawn randomly from the posterior. This provided us with confidence that the methodology used in this work is sound.

As an example, we describe here the process performed on a Schechter LF model, with a superposition of PDE and PLE evolution (named Sadler$+02$ within this work). The area of the field was set to $40.46\ \mathrm{deg}^{2}$, and the detection limit to $50\ \mathrm{\mu Jy}$. This resulted in a simulated catalogue of $6378$ mock sources above the detection limit, created by randomly selecting sources via the assumed LF. The starting parameters of the LF are given in Tab. 3. The scatter in redshifts was set here to be negligible, but a finite uncertainty in redshifts, via hierarchical bayesian interference, was also tested. The parameter modeling was performed on this simulated data set, using the same codes later used on observational data. The retrieved parameters are shown also in Tab. 3. The detection limit in this example was a step-function i.e. the completeness corrections were not present. We also assumed a mean spectral index of $-0.7$. The codes were also tested for non-negligible completeness corrections, by introducing completeness correction as a separate function during the integration of log-Likelihood. The methodology was tested for all the models described in Sect. 6 and on different areas and depths of mock-catalogues. The parameters of the LFs were always retrieved successfully.