LOFAR window on SF galaxies and AGN

The LOFAR window on star-forming galaxies and AGN – curved radio SEDs and IR-radio correlation at


We present a study of the low-frequency radio properties of star forming (SF) galaxies and active galactic nuclei (AGN) up to redshift . The new spectral window probed by the Low Frequency Array (LOFAR) allows us to reconstruct the radio continuum emission from 150 MHz to 1.4 GHz to an unprecedented depth for a radio-selected sample of galaxies in of the LOFAR Boötes field. Using the extensive multi-wavelength dataset available in Boötes and detailed modelling of the FIR to UV spectral energy distribution (SED), we are able to separate the star-formation (N=758) and the AGN (N=784) dominated populations. We study the shape of the radio SEDs and their evolution across cosmic time and find significant differences in the spectral curvature between the SF galaxy and AGN populations. While the radio spectra of SF galaxies exhibit a weak but statistically significant flattening, AGN SEDs show a clear trend to become steeper towards lower frequencies. No evolution of the spectral curvature as a function of redshift is found for SF galaxies or AGN. We investigate the redshift evolution of the infrared-radio correlation (IRC) for SF galaxies and find that the ratio of total infrared to 1.4 GHz radio luminosities decreases with increasing redshift: . Similarly, shows a redshift evolution following . Calibration of the 150 radio luminosity as a star formation rate tracer suggests that a single power-law extrapolation from is not an accurate approximation at all redshifts.

galaxies: active, evolution, photometry, starburst, infrared: galaxies, radio continuum: galaxies

1 Introduction

Radio selected samples of galaxies primarily consist of two populations: star-forming galaxies (hereafter SF galaxies) and active galactic nuclei (AGN). As both star formation activity and black-hole growth in AGN are processes closely related to the overall mass growth of galaxies (Shapley, 2011; Best & Heckman, 2012), radio-emitting galaxies therefore represent a unique laboratory for investigating the epoch of peak galaxy assembly, . While studies of these galaxies from the far-infrared (FIR) to the ultraviolet (UV) have contributed greatly to our understanding of galaxy evolution, statistically significant samples of radio-selected galaxies across cosmic time have only recently become available for exploration thanks to a new generation of radio facilities. At low radio frequencies the unprecedented sensitivity and resolution of the Low Frequency Array (LOFAR, van Haarlem et al., 2013) opens a new window for galaxy evolution studies.

Low frequency radio emission in SF galaxies and AGN have different physical origins, although both are thought to be dominated by synchrotron emission. In SF galaxies the synchrotron emission is powered by high-energy electrons and positrons (cosmic rays, CRs), accelerated in supernova remnants (SNRs), that emit when interacting with the diffuse magnetic field of the galaxy (Condon, 1992). Due to the short lifetime of the massive stars producing Type II and Type Ib supernovae, the synchrotron emission in SF galaxies is closely related to recent star formation, so that its emission (e.g. at 1.4 GHz) is widely used as a star formation tracer (Condon, 1992; Bell, 2003; Schmitt et al., 2006; Murphy et al., 2011).

In radio-selected AGN, the synchrotron radiation is ultimately powered by the central accreting black hole. However, the observed emission is believed to be emitted from regions which differ depending on the nature of the AGN (see Smolcic, 2016, for a review of AGN radio classifications). Fundamental physical differences have been seen between radio AGN classified as high-excitation and low-excitation radio galaxies (HERGs and LERGs, Hardcastle et al., 2007; Best & Heckman, 2012), where HERGs appear to be the dominant population at high radio luminosities (), while LERGs appear to dominate below this limit. Most HERGs are observed to consist of three large-scale structures: jets, hotspots and lobes (consistent with Fanaroff & Riley (1974) class II, FRII galaxies). High-energy CRs first emit while being accelerated in the relativistic jet, which transports them from the central AGN to shock regions called hotspots (e.g. Meisenheimer et al., 1989). Finally, the CRs expand away from the hotspot center, forming the large lobes of synchrotron emission (e.g. Krause et al., 2012).

The synchrotron emission in all of these processes is commonly described by a power-law function , where is the radio flux density at a given frequency, , and the corresponding power-law slope. Nevertheless, there are processes which can alter the shape of the spectra. While in SF galaxies the synchrotron emission is typically observed to follow a slope of (e.g. Gioia et al., 1982; Condon, 1992; Yun et al., 2001), intrinsic changes in the CR energy distributions or environmental and ISM processes (free-free absorption, ionization losses, and synchrotron self-absorption) can significantly alter the spectral shape of this emission (Lacki, 2013, and references therein), intrinsic changes in the CR energy distributions or environmental and ISM process . In AGN, both spectral ageing (e.g. Harwood et al., 2015) and the relative brightness of the different components (e.g. core, jets, hotspots and lobes) play an important role in shaping the integrated radio SEDs of radio AGN (e.g. Hardcastle, 2009).

Multi-frequency spectral studies of the integrated radio SED of AGN are not widespread in the literature (though there are a few exceptions, Laing & Peacock, 1980; Ker et al., 2012; Singh et al., 2013; Kharb et al., 2016; Mahony et al., 2016). Recent literature has focused on morphological studies through spectral index maps (e.g. Harwood et al., 2015; Vardoulaki et al., 2015). A luminosity dependence of the spectral shape of the integrated radio spectrum was studied e.g. by Laing & Peacock (1980), where they suggest radio AGN with have spectra that steepen at low frequencies, while brighter AGN show the opposite trend. Similarly, Whittam et al. (2016) observed spectral flattening for radio-selected galaxies at the high-frequency end (15.7 GHz) and suggest this may be due to the cores of Fanaroff & Riley (1974) class I sources (FRI) becoming dominant at these high frequencies. Deep multi-frequency radio observations of representative samples of galaxies are needed to study the radio SED and understand the physical processes shaping it.

Observational studies of the spectral properties of radio SF galaxies have so far focused on the local universe as sensitivity usually prevents the detection of statistically significant samples of galaxies at higher redshifts. An exception to this is the radio spectral slope study presented by Ibar et al. (2009, 2010) for a sample of submillimetre galaxies (SMGs). They found no redshift evolution of the spectral indices for SMGs, although due to the sparse coverage of the radio SED they could not rule out the presence of curvature. For local SF galaxies, spectral flattening towards low radio frequencies has been observed (), with thermal absorption or intrinsic synchrotron curvature as plausible explanations for this (Israel & Mahoney, 1990; Clemens et al., 2010; Marvil et al., 2015). Theoretical models explaining an alternative picture to the simple power-law shape which includes spectral curvature have also been developed (e.g. Lacki, 2013).

In SF galaxies the radio emission is typically calibrated to trace star formation based on a tight empirical correlation with the IR radiation (IR-radio correlation (IRC), de Jong et al., 1985; Helou et al., 1985). In the local universe the IRC has been observed to be roughly linear across more than three orders of magnitude in FIR-luminosity () (Yun et al., 2001; Magnelli et al., 2015), for different galaxy classes (from dwarf to ultra-luminous IR galaxies, ULIRGs), and in star-forming regions within galaxies (e.g. Dumas et al., 2011; Tabatabaei et al., 2013). The basic understanding of the IRC relies on both the cold dust IR emission and the radio emission being tracers of recent star formation (e.g. calorimeter model and conspiracy model, Voelk, 1989; Lacki et al., 2010, respectively). The possibility of a redshift evolution of the IRC has motivated an extensive debate from the observational point of view: while several studies claim the existence of significant redshift evolution(Seymour et al., 2008; Ivison et al., 2010b; Magnelli et al., 2015; Delhaize et al., 2017), a few studies suggest that such an evolution is a product of selection biases (e.g. Appleton et al., 2004; Ibar et al., 2008; Jarvis et al., 2010; Sargent et al., 2010; Bourne et al., 2011) or a dust temperature dependence (e.g. Smith et al., 2014). In theoretical studies different trends as a function of redshift have also been discussed (e.g. Schober et al., 2016; Schleicher & Beck, 2013; Lacki et al., 2010).

In this work we take advantage of the unique sensitivity and resolution of LOFAR (van Haarlem et al., 2013; Röttgering et al., 2011), combining deep 150 MHz radio observations (Williams et al., 2016) and the wealth of ancillary data available in the Boötes field to address some of the key outstanding questions regarding the radio emission of galaxies, namely, what are the star formation and AGN contributions to the radio continuum at low frequencies?; does the low-frequency radio emission of SF galaxies exhibit a correlation with the IR luminosity as tight as that observed at 1.4 GHz?; and if so, how does the IRC evolve with redshift at both 1.4 GHz and 150 MHz? Answering these questions is crucial to test the reliability of photometric-redshift estimations based on the radio–IR SEDs (e.g. Yun & Carilli, 2002; da Cunha et al., 2015). Finally, this study will allow us to investigate the low-frequency radio emission as a SFR tracer of galaxy populations at higher redshifts.

This paper is structured as follows. In section 2 we present the multi-frequency radio and FIR to UV data used for this study, while section 3 describes the selection strategy for our sample. In section 4 we discuss the SED analysis and classification into SF galaxies and AGN populations. The analysis of the spectral slope and curvature for the two populations is described in section 5. Next, the IRC in our data is investigated in section 6 and section 7 discusses the low frequency radio emission as a SFR diagnostic. Finally, section 8 summarizes our findings. Throughout the paper we adopt a concordance flat -CDM cosmology with , , and (Komatsu et al., 2009; Planck Collaboration et al., 2014) and all quoted magnitudes assume the AB system (Oke & Gunn, 1983) unless otherwise specified.

2 Survey data

The National Optical Astronomy Observatory (NOAO) Deep Wide-field survey (NDWFS; Jannuzi & Dey, 1999) targeted the sky seen towards the constellation of Boötes as one of its deep extragalactic fields. Originally the NDWFS covered in the optical and near-infrared B R, I and K bands. Since then, the Boötes Field has been surveyed across the electromagnetic spectrum, including deep X-ray (Murray et al., 2005), mid-infrared (Jannuzi et al., 2010) and far-infrared (Oliver et al., 2012) photometric observations as well as extensive spectroscopic surveys (e.g. Kochanek et al., 2012). The field therefore represents one of the richest multi-wavelength datasets among the wide deep extragalactic surveys and is complementary to the deep radio observations provided by LOFAR.

2.1 LOFAR 150 MHz

The first source catalogue used for our sample selection is based on 150 MHz radio observations of the Boötes field (Williams et al., 2016) taken with the LOFAR High Band Antennae (HBA). The calibration and imaging were achieved with the ‘Facet’ calibration scheme presented by van Weeren et al. (2016), which corrects for direction dependent effects (DDEs) caused by the ionosphere and imperfect knowledge of LOFAR station beam shapes. The resulting image, as presented by Williams et al. (2016), covers 19 deg with an rms noise of Jy beam in the central region of the field and a resolution of . These values represent up to more than one order of magnitude improvement compared to the images existing at this wavelength (Intema et al., 2011; Williams et al., 2013). The source extraction for the catalogue was done using the Python Blob Detection and Source Measurement software, (PyBDSM: Mohan & Rafferty, 2015), which performs Gaussian fitting to decompose radio interferometry images, grouping the Gaussians together into individual sources where appropriate. The LOFAR 150 MHz radio source catalogue contains 6276 sources detected with a peak flux density threshold of 5 within the coverage region shown in Fig. 1. The LOFAR 150 MHz luminosity distributions for our total sample classified into SF galaxies and AGN are shown in the lower panel of Figure 2.

2.2 Radio photometry for the Boötes field

The construction of radio SEDs for the sources in the Boötes field is one of the main aims of this paper. Conveniently, the Boötes field has been covered by several previous radio surveys. Among these, the high sensitivity and the resolution of the VLA-P Survey at 325 MHz, GMRT observations at 608 MHz and the deep WSRT catalogue at 1.4 GHz, make them the best data sets to complement the LOFAR data in our study.

The VLA-P catalogue (Coppejans et al., 2015) used in this investigation, was drawn from a 324.5 MHz image of a radius of 2.05 in the NOAO Boötes field using Karl G. Jansky Very Large Array (VLA) P-band observations. The image covered a single pointing and has a resolution of 5.65.1 arcsec with a central noise of 0.2 mJy beam increasing to 0.8 mJy beam at the edge of the image. The source extraction for the catalogue was done using the PyBDSM software, (PyBDSM: Mohan & Rafferty, 2015). The source detection threshold for the construction of the catalogue for this study is 3 . This is a reliable detection limit since our study includes only sources from the catalogue which have LOFAR detected counterparts. As discussed by Coppejans et al. (2015), the VLA-P catalogue was matched to WENSS (Rengelink et al., 1997) to check the absolute flux density scale and the primary beam correction.

The GMRT image at 608 MHz is a mosaic constructed from four pointings covering 1.95 deg of the Boötes field (project code 28_064). The mosaic has a resolution of arcsec and a noise level of Jy beam. Primary flux density calibration was done with 3C286 using the wideband low-frequency flux density standard of Scaife & Heald (2012). The source extraction for the catalogue was also performed using the PyBDSM software, (PyBDSM: Mohan & Rafferty, 2015).

The 1.4 GHz data (de Vries et al., 2002) are drawn from the deep (1612 hr) Westerbork Synthesis Radio Telescope (WSRT) observations of the approximately 6.68 deg Boötes Deep Field. The image covers 42 discrete pointings and has a limiting sensitivity of Jy beam and a resolution of 1327 arcsec. The source extraction in the public catalogue was done using automated routines described in detail by Rengelink et al. (1997), which consist basically of Gaussian fitting to islands of detected brightness in the radio map. The source detection threshold for the inclusion of the flux densities in the catalogue is 5 times the local rms, resulting in the full catalogue containing 3172 sources. Since the Boötes field has been covered by previous radio surveys at this wavelength such as the NRAO VLA Sky Survey (NVSS; Condon et al., 1998) and the Faint Images of The Radio Sky at 20 cm Survey (FIRST; Becker et al., 1995), these data were used to calibrate the survey flux densities and positions.

Assuming a spectral index of -0.7, the LOFAR measurements offer sensitivity values comparable to the 1.4 GHz WSRT map presented by de Vries et al. (2002) (rms: 28 Jy beam) and resolution values comparable to the 365 MHz VLA-P map published by Coppejans et al. (2015) (), which are the best radio data available at the respective frequencies to date.

Band Detections [%]
Selection LOFAR 150 MHz 100%
I 806 nm 100%

NUV 300 nm 93%
u 365 nm 99%
optical B 445 nm 100%
R 658 nm 100%
NIR z 900 nm 99%
Y 1020 nm 99%
J 1222 nm 100%
H 1630 nm 100%
K 2190 nm 100%
MIR IRAC1 3.6 m 100%
IRAC2 4.5 m 100%
IRAC3 5.8 m 100%
IRAC4 8.0 m 100%
MIPS24 24.0 m 100%

SPIRE 250 m 83%
SPIRE 350m 79%
SPIRE 500 m 67%
RADIO WSRT 1.4 GHz 67% (88%3)4
GMRT 610 MHz 23%5
VLA-P 325 MHz 63% (88%)
Table 1: Multiwavelength coverage of the LOFAR+I-band selected sample in the Boötes field. The column of ‘Detections’ specifies which fraction of the final sample has counterparts in the respective bands.
Figure 1: Spatial coverage of the optical and radio photometry including the I-band selected photometry presented by Brown et al. (2007) and the four deep catalogues at 150, 325, 610, 1400 MHz, available for the Boötes field. The dashed area is the area chosen for our study.

2.3 Optical and Infrared Photometry

The second source-catalogue used for our sample selection is the combined I-band-selected psf-matched photometry catalogue presented by Brown et al. (2007). The photometric bands included in the catalogue are presented in Table 1 covering a wide range of wavelengths, spanning from 0.15 to 24 m.

This catalogue includes the original NDWFS observations in the optical and near- infrared; B R, I and K bands. With an absolute positional uncertainty of arcsec, the I-band images reach depths of 24.9 AB magnitude ( within 2 arcsec diameter aperture).

Mid-infrared (MIR) counterparts are drawn from the Spitzer Deep Wide-field Survey (SDWFS; Ashby et al., 2009), using the InfraRed Array Camera (IRAC) instrument on the Spitzer Space Telescope, providing images at 3.6, 4.5, 5.8, and 8.0 m. Infrared 24 m photometry from Spitzer was provided by the Multiband Imaging Photometer (MIPS) AGN and Galaxy Evolution Survey (MAGES; Jannuzi et al. (2010)), while deep J, H, K photometry was drawn from the NOAO Extremely Wide-field Infrared Imager (NEWFIRM; Autry et al., 2003) survey. A fraction of the total Boötes field is covered by the z-Boötes survey (Cool, 2007), providing photometry in the z-band for 7.62 deg of the field, with some gaps in the coverage due to the 10 arcmin gaps in the instrument CCDs. Y and U bands photometry is provided using images from the Large Binocular Camera (LBC) mounted on the Large Binocular Telescope (LBT). Finally, NUV photometry (1800-2750 Å) has been included from the GALEX/GR6 surveys (Martin et al., 2003; Bianchi et al., 2014).

As described in detail by Brown et al. (2007), this psf-matched catalogue wass constructed by regridding and smoothing the individual survey images corresponding to the u-, B-,R-, I-, z-, Y-, J-, H- and K-bands to a common scale so that the sources’ point spread functions (PSF) are Moffat profiles. Fluxes are extracted from these images for all the sources with I-band detections using SExtractor (Bertin & Arnouts, 1996), while for the remaining bands (IRAC and MIPS bands) aperture fluxes were extracted. Regions surrounding very extended galaxies and saturated stars were excluded.

In total, the combined catalogue provides photometry for k sources, including galaxies and AGN with I-band magnitudes .

Figure 2: Redshift and LOFAR 150 MHz luminosity distributions for the SF galaxies (red) and AGN (blue) populations. The classification of the total sample into SF galaxies and AGN is explained in detail in section 4.2.

2.4 Far-infrared Photometry: HerMes DR3

Far-infrared photometry is a key ingredient for this study, as it allows us to constrain total infrared luminosities and trace SFRs of the sources. To increase the infrared spectral coverage of the final sample, we added SPIRE data at 250, 350 and 500 m from the Herschel Multi-tiered Extragalactic Survey (HerMes; Oliver et al., 2012). One of the greatest challenges in measuring fluxes at such long wavelength is the high confusion noise, a result of the large SPIRE beam sizes. Since one effect of confusion is that it increases the positional uncertainty of sources (e.g., Hogg, 2001), cross-identification with other wavelengths becomes very challenging. To deal with this issue, we use the third data release (DR3) cross-identification catalogues, which are selected with positional priors at 24 m and obtained with the technique described by Roseboom et al. (2010). This technique consists basically of performing cross-identifications in map-space so as to minimise source blending effects by using a combination of linear inversion and model selection techniques. In this way they produce reliable cross-identification catalogues based on Spitzer MIPS 24µm source positions, giving significantly greater accuracy in the flux density compared to other traditional source recovery methods and so recovering a much larger fraction of faint SPIRE sources.

The FIR photometry was included in the multi-wavelength data by cross-correlating the source selection catalogue described above with the Hermes DR3 catalogue, using the task TSKYMATCH2 from the STILTS software package (TOPCAT implementation, Taylor, 2006). Using the optical position of the sources in the Brown et al. (2007) catalogue, TSKYMATCH2 returned the SPIRE fluxes associated to all sources with MIPS 24m to optical separations of r arcsec (less than half the FWHM of the MIPS 24m beam).

The percentage of sources with counterparts at different wavelengths are listed in Table 1.

2.5 Redshifts

The spectroscopic AGN and Galaxy Evolution Survey (AGES) has covered 7.7 deg of the Boötes Field providing spectroscopic redshifts for 23 745 galaxies and AGN. However, given that the majority of sources at in the AGES catalogue are QSOs, it is crucial to estimate robust photometric redshifts, as this study focuses on both AGN and star forming galaxies. Spectroscopic redshifts are available for around 45 and 35 per cent of our total selection of SF galaxies and AGN respectively (the selection and classification is described in sections 3 and 4) . Photometric redshifts are estimated as described below for the remaining fraction of the sources.

Figure 3: Test for completeness of the LOFAR selection. The solid lines depict the redshift evolution of the expected radio fluxes for different SF galaxies and are colour-coded by SFR. The 150 MHz luminosities where calculated based on the IR-radio correlation, assuming and a spectral index of for the extrapolation from 1.4 GHz to 150 MHz. The dashed line corresponds to the 150 MHz flux evolution as a function of redshift for AGN of . The grey shaded area displays the flux limit implied in our radio selection, which correspond to 5

Photometric redshifts are provided by the optimised catalogue described by Duncan et al. (in prep.), produced using the Boötes photometry of Brown et al. (2007), presented above. The innovative method used for this redshift catalogue consists of combining three different -estimation methods, by using the EAZY photometric redshift software (Brammer et al., 2008) customized with three different template sets: one set of stellar only templates (EAZY default library; Brammer et al., 2008) and two sets including AGN and QSO contributions (SWIRE; Polletta et al., 2007) and (Atlas of Galaxy SEDs; Brown et al., 2014). This methodology has been motivated by the fact that comparisons of -estimation techniques have shown that for a suite of multiple z-estimates, the ensemble average estimates offer significant statistical improvements in redshift accuracy compared to the estimates of any single set of predictions (Dahlen et al., 2013; Carrasco Kind & Brunner, 2014). These three individual estimates were then combined using a Hierarchical Bayesian combination method (Dahlen et al., 2013), as an alternative to a straight addition of the probability distributions of the three estimates. The main advantage of this method is that it determines the consensus probability P for each object, given the possibility that the individual measured probability distributions may be wrong. These results were also optimised using zero-point offsets calculated from the spectroscopic redshift sample. The redshift distributions for the total sample classified into SF galaxies and AGN are shown in the upper panel of Figure 2.

3 Sample Selection

The selection strategy for the sample of galaxies and AGN in the Boötes Field involved several steps. First, sources detected both in LOFAR 150 MHz and I-band were selected to build the primary source catalogue, which is additionally complemented by multi-wavelength ancillary data. Second, we constrained this catalogue to sources with positions within the region of overlap among the 150 MHz, 325 MHz and 1.4 GHz radio maps (see Section 2.2 and Fig. 1). Finally, we classified these sources into SF galaxies and AGN-dominated sources using the AGNfitter algorithm (Calistro Rivera et al., 2016) and restrict the posterior analysis on the sub-samples separately.

3.1 LOFAR + I-band selection

The identification of optical counterparts to the LOFAR detections from the Brown et al. (2007) catalogue is described in detail by Williams et al., in prep. We use the likelihood ratio (LR) method (Richter, 1975) to quantify the probability of an I-band detected galaxy being the true counterpart of the radio emission observed in the LOFAR map. The fully developed method used is described by Tasse et al. (2008). From the 3317 LOFAR-detected sources, which lie within the 9 deg boundary of the Brown et al. (2007) catalogue, 2326 (70 per cent) are found to have unique optical counterparts.

Fig. 3 shows the implications of our radio selection for the nature of our sample. Lines correspond to the expected radio fluxes at 150 MHz, based on the IR-radio correlation, assuming and a spectral index of for the extrapolation from 1.4 GHz to 150 MHz. We choose to use a value of based on the observed median value of found for our sample (see Fig. 13). The flux cut from our radio selection is represented by the grey shaded area. Due to our selection, at our sample appears to be complete for SF galaxies of M/yr (yellow line). At redshifts only galaxies with /yr are represented in a complete way according to the LOFAR selection. We expect to be complete for AGN of at all redshifts. We thus conclude that at redshifts our sample is limited to highly SF galaxies (starbursts) and radio-selected AGN.

Figure 4: Example SED of dusty SF galaxy from our sample at z=1.2. The multi-wavelength coverage of our total sample is shown here, from the low-frequency radio till the UV regime. The radio data points represent the radio SED covered in this study: 150, 325, 610 MHz and 1.4 GHz. The FIR-UV photometry is decomposed into physical components of the host galaxy and the AGN using the code AGNfitter. The integrated luminosities used in our classification scheme were computed by integrating the SED within the areas represented here as the grey-shaded rectangles. The left rectangle covers the integration limits for the main IR components: the galactic cold dust emission (green line) and the hot dust emission from the AGN torus (purple line). The right rectangle is the integration area for the optical main components: the stellar emission (orange line) and the accretion disk emission (blue line).The thin red line represents the total fitted SED. The dashed red and grey lines correspond to the synchroton and thermal contributions to the radio emission, respectively.

3.2 Selection for radio SED construction

The radio SEDs constructed from the data described in Section 2.2 are built up mostly with three data points: LOFAR 150 MHz, VLA-P 365 MHz and WSRT 1400 MHz (Fig. 1), complemented with a partial coverage of GMRT observations at 610 MHz. From the 2326 LOFAR sources with optical counterparts, 1955 sources lie within the region of overlap among the 150 MHz, 325 MHz and 1.4 GHz radio maps required for our study. To study the shape of the radio continuum we need reliable spectral index measurements along the total frequency coverage. As spectral index measurements are prone to be easily affected by systematic offsets between different bands, special care should be exercised in having an homogeneous flux scales and in dealing with non-detections. We dedicate this section to this purpose, review all systematics intrinsic to each data set and corroborate the homogeneity of the flux scaling.

Flux scales

The combination or comparison of radio maps at different frequencies requires them to undergo a standard calibration scaling. For this study we adjust the fluxes of the four different bands to the scale that is most accurate at low radio frequencies, that of Scaife & Heald (2012). Since the main uncertainties of using a standard calibration scale are at frequencies lower that 1 GHz we focus our investigation on the LOFAR and VLA-P data sets. However, we note that the scaling of the high-frequency range covered by the GMRT and WSRT datasets was verified to also be consistent with the scale by Scaife & Heald (2012).

To investigate the reliability of LOFAR fluxes at 150 MHz, Williams et al. (2016) have compared high signal-to-noise sources to three other data sets with known calibrations in the Boötes Field: NVSS (at 1.4 GHz), WENSS (Rengelink et al., 1997) (at 365 MHz) and VLSSr (Lane et al., 2014) (at 74 MHz), where the WENSS fluxes had been scaled a priori by a factor of 0.9 to be consistent with the Scaife & Heald (2012) scale. The inference of the correction factor for the LOFAR fluxes consists of calculating the spectral indices between the lower (74 MHz) and higher (325, 1400 MHz) frequencies and predicting the LOFAR flux density at the central frequency (150 MHz). The mean flux density ratio between the predicted and uncorrected LOFAR flux was found to be . Finally, as described by Williams et al. (2016), this factor was used to adjust the LOFAR flux densities to the Scaife & Heald (2012) flux scale.

The reliability of the VLA-P flux-scale was tested using two independent approaches which yielded consistent results. First, we investigated the calibrator source (3C216) that lies in the Boötes field. We calculated the ratio the 325 MHz flux measured during the calibration of the map (which assumes the flux scale of Perley&Butler 2010) and the 325 MHz flux predicted by Scaife & Heald (2012). A correction factor of 0.91 was needed to adjust the calibrated fluxes to the standard scale of Scaife & Heald (2012). In a second approach, we compared high signal-to-noise VLA-P sources with their counterpart fluxes from the WSRT, GMRT and the LOFAR catalogues. The spectral indices between the lower (150 MHz) and higher frequencies (610 MHz, 1.4 GHz) were calculated and we could predict the VLA-P flux density at the central frequency (325 MHz). The mean flux density ratio between the predicted and uncorrected VLA-P flux was , in agreement with the data-independent method presented above.


Although the sensitivity of the LOFAR image is equivalent to the WSRT sensitivity under the assumption of a spectral index of , 33 per cent of the LOFAR selected sources are not detected in WSRT. This issue is also relevant for the VLA-P data, which are less sensitive () than the equivalent LOFAR image at this frequency for the same assumed (). We find that 37 per cent of the LOFAR selected sources are not detected in VLA-P.

To estimate the fluxes of the sources undetected in WSRT and VLA-P we use forced photometry technique. This process consists of extracting aperture fluxes from the radio maps at the location of known LOFAR selected sources. To calibrate this process we first applied the forced photometry on catalogued sources. We found that aperture fluxes were sensitive to artifacts of the imaging process (e.g. negative fluxes) and to reduce this effect, we used apertures smaller than the WSRT and VLA-P beam sizes. Finally, we used our results on catalogued sources to derive a correction factor, which was applied on the undetected sources. The flux uncertainties on the forced photometry were calculated by repeating the same flux extraction process as above using the corresponding rms-maps.

Although the aperture fluxes have been determined to be reliable, the results on radio continuum that are based on forced photometry (see Section 5), are carefully differentiated from those of SEDs fully sampled with detected sources, using forced photometry fluxes just as upper limits. Moreover, for more specific tests in the following sections we apply conservative flux cuts to minimize the fraction of non-detections in our data.

4 Analysis of AGN and SF galaxies with SED-fitting

The FIR to UV SEDs of the sources in our sample were decomposed into different physical contributions through fitting the multi-wavelength photometry using the SED-fitting algorithm AGNfitter. Some examples of the fitting output are included in Appendix B. An advantage of using AGNfitter is that it is based on a Markov Chain Monte Carlo (MCMC) technique and infers the probability density functions (PDFs) of the physical parameters. This way it provides a robust calculation of their uncertainties and improving the recognition of correlations and degeneracies among them.

The total active galaxy model in AGNfitter consists of the superposition of the host galaxy emission and the nuclear AGN emission. The host galaxy emission is modelled as a combination of a stellar component and the reprocessed emission of cold/warm dust in starburst regions. At nuclear scales, the AGN emission is modelled as a combination of an accretion disk component (Big Blue Bump, BBB) and a hot dust ’torus’ component.

A large set of relevant physical parameters for the galaxy (SFR, ) and AGN (N-torus, ) are calculated, some of which are detailed below. However, for full details we refer the reader to Calistro Rivera et al. (2016). From the 1955 LOFAR-I-band selected sources within the region of overlap between the multi-frequency radio data (section 3.2), 1542 ( per cent) can be classified as 758 SF galaxies and 784 AGN.

Figure 5: Classification into AGN and SF galaxies using ratios of integrated luminosities of galactic versus AGN components. Integrated luminosities were computed from component templates fitted with AGNfitter. The red-shaded area shows the region populated by sources which are considered SF galaxies, while the blue-shaded areas are populated by sources classified as AGN. The data points are shown in black and their error bars are shown in red or blue, depending if their classified as SF galaxies or AGN, respectively.

4.1 Stellar masses and star-formation rates

Stellar masses and SFRs () were estimated by fitting the observed data to Bruzual & Charlot (2003) templates, assuming a Chabrier (2003) initial mass function (IMF) and an exponentially declining star formation history modulated by the time scale parameter . The distribution of stellar masses derived for our sample has a median value of , where the error bars correspond to the scatter given by the 16th and 84th percentiles. The templates were corrected for dust absorption assuming the Calzetti et al. (1994) reddening law. Additional SFRs from the IR-emission were calculated using integrated infrared luminosities of cold dust templates () (Dale & Helou, 2002; Chary & Elbaz, 2001). To derive the SFRs from the fitted cold dust emission we use the calibrations presented by Murphy et al. (2011), which are an updated version from those by Kennicutt (1998). These estimates are described in detail in Calistro Rivera et al. (2016).

4.2 Classification into AGN and SF galaxies

Figure 6: Comparison of the radio-independent classification by AGNfitter to the IR-radio correlation. IR luminosities are inferred as the luminosities of the cold dust emission component in the sources’ SEDs integrated within the range 8-1000 m, while radio luminosities are k-corrected observations at 1.4 GHz, assuming a spectral index of . The solid line represents the canonical value of the correlation , while the shaded area show a scatter of . The red and blue data points represent the values for galaxies classified as SF galaxies and AGN respectively.

The classification strategy consisted of comparing the disentangled contributions of the AGN and the host galaxy in both the optical-UV and the infrared regimes. It is important to note that this classification may still be contaminated by low excitation radio AGN (LERGs), which lack AGN signatures in their non-radio SED. In section 4.3 we will describe how we corrected for these potential mis-classifications.

In the optical-UV, the AGN contribution arises from the BBB luminosity , which is then compared to the galaxy stellar emission , integrating both inside the same frequency range 0.1 m. In the MIR, the AGN luminosity contribution is produced mainly by the torus and is compared to the cold dust emission in starburst regions of the host galaxy , as well integrating both components within the same frequency range 1 m. A representation of the SED decomposition and the integration ranges is shown in Fig. 4.

As presented in Fig. 5, in order to classify the total sample into SF galaxies and AGN we define the following scheme:

If the ratio of AGN torus to galaxy cold dust luminosities is smaller than one (), we consider the source’s emission in the IR is dominated by star formation and we classify it as a SF-galaxy. The optical/UV contributions do not change the classification significantly, since the results for the 94 per cent of the sample which satisfy , also satisfy . In the cases with (only around 6 per cent of the sample), the direct emission of the AGN accretion disk appears to dominate over the stellar emission. However, since their FIR emission is dominated by the dust component in star-forming regions, these sources may represent the fraction of highly obscured star forming galaxies (e.g. Casey et al., 2014) and are thus classified as star-forming galaxies.

If the ratio of AGN torus to galaxy cold dust luminosities is greater than one () the source is considered an AGN according to our scheme. In the optical, the classification as AGN is less stringent and the source may host an AGN both in the case were , where the accretion disk direct emission dominates the optical and if , where the source may present obscuration at nuclear scales and is considered an obscured AGN.

To test this classification, we take advantage of the availability of the radio data and use the infrared-radio correlation (IRC; e.g. Yun et al., 2001), which is a property observed almost exclusively by SF galaxies and radio-quiet AGN. Fig. 6 shows an excellent agreement between the radio independent classification method based on SED-fitting compared to the IRC results. This illustrates the capability of our method in classifying the sources into AGN and galaxies. Prior to our correction for contamination by LERGs our samples consists of 810 SF galaxies and 732 AGN, both populating the redshift range from 0.05 to 2.5.

4.3 Contamination by Low Excitation Radio Galaxies

A consequence of classifying AGN and galaxies based on their SEDs from the FIR-UV is that the population of LERGs are prone to be misclassified as galaxies, due to the lack of AGN signature in their non-radio SED. This mis-classification can be partially corrected by the assumption that LERGs would populate the area in the IRC typically covered by AGN, despite being classified as SF galaxies. In Fig. 6 LERG candidates are the red data points which are outliers from the IRC area (, as estimated below in section 6).

We applied a conservative correction for LERGs contamination and reclassified the LERG candidates into the AGN class by choosing the values below the 2.5 sigma region of the IRC value to be the division line (, as observed for our sample). This correction finds a total of 52 LERGs in our sample, which represents around 7 per cent of the total AGN population in our data. To test our correction we investigate other properties of the 52 sources identified as LERGs. We find that although LERGs have IR luminosities dominated by SF processes by definition , the dominance is clearly weak () compared to the values for SF galaxies (, were the errorbars represent the scatter given by the 14th and 86th percentiles). Since these sources present a tendency towards the values of the HERGs classified as AGN (), these results support the reliability of our correction. For a detailed discussion of the properties of mass-selected LERGs and HERGs populations within Boötes see Williams et al. (in prep).

We note that this correction applies for sources where AGN and star-formation activity are mutually exclusive processes, but this is not generally expected since AGN activity may occur in strongly star-forming galaxies (Stevens et al., 2003; Rees et al., 2016). Moreover, since the IRC is a central topic of this work, the classification of our sources needs to be as independent of it as possible to avoid biases in the results. Through the chosen limits, we consider our correction is conservative enough for our purposes.

Our final sample consists of 758 SF galaxies and 784 AGN, which are well covered in the redshift range , while the redshift range is slightly affected by incompleteness due to the flux limit of the data.

4.4 Comparison of the classification strategy with literature values

A classification based solely on the SED-fitting output yields that the relative source fractions of SF galaxies and AGN is 49 and 51 per cent respectively, given our LOFAR-I-band selected sample. A direct comparison to other observational and theoretical studies on this fraction is a complex task, since it is highly dependent on the sensitivity of the surveys (e.g. Jarvis & Rawlings, 2004; Appleton et al., 2004; Ibar et al., 2008; Simpson et al., 2012). For instance, a similar LOFAR-selected sample presented by Hardcastle et al. (2016) shows a number density split of SF galaxies to AGN of approximately 30/70. The lower SF galaxies fraction found can be explained due to the shallower radio data used in that study but also due to their different selection criteria (based on the FIR-radio correlation only). However, by using their selection criteria in our sample we still recover our original fraction of 50/50. Similarly, a study of radio sources selected at 3 GHz in the COSMOS field (Delvecchio et al., 2017) finds a SF galaxies/AGN ratio of based on an equivalent SED-fitting classification.

Investigations of the SF galaxy fraction in 1.4 GHz radio surveys (e.g. Seymour et al., 2008; Bonzini et al., 2013; Smolčić et al., 2008) have found that while SF galaxies show predominance at low flux densities, the AGN population start dominating at 100 Jy. Assuming a spectral index value of this turning point corresponds to Jy at 150 MHz, which is close to the detection limit of our sources. Based on these assumptions, the observed 49/51 ratio in our study is thus in agreement with the ratio expected for the flux regime around the dominance turning point.

5 The radio SEDs of SF galaxies and AGN

We characterize the radio SED for galaxies and AGN by calculating the distributions of different spectral indices, which we define as:


For the total sample we calculate the most general spectral index using only the two extreme points of the total frequency range studied, from 150 MHz to 1.4 GHz.

More detailed spectral properties such as curvature are studied by calculating adjacent frequency pairs, and in Section 5.2. To reliably measure spectral indices over narrow frequency ranges requires higher signal to noise data. We therefore use a subsample with LOFAR fluxes above 2 mJy in order to minimize the effect of non-detections on these results. Finally, the spectral study is refined for a fraction of the total sample which has GMRT spatial coverage at 610 MHz in Section 5.3. For this sub-sample ( per cent of the total) the distributions of three different spectral indices are calculated for adjacent frequency pairs and .

5.1 Spectral index – all sources

Figure 7: Spectral index for SF galaxies (upper panel) and AGN-dominated galaxies (lower panel). Values corresponding to the 16th, 50th and 84th percentiles are calculated for the total sample. Detected and undetected sources in 1.4 GHz are depicted in red and orange respectively for SF galaxies, and dark blue and sky blue for AGN. The solid line displays the fitted Gaussian with mean and intrinsic scatter .

The distribution of the spectral index in Fig. 7 includes the complete LOFAR-I-band selected sample. Note that we use forced photometry in the cases of non-detections in 325 MHz or 1.4 GHz. We obtain the median and scatter values for the distributions of SF galaxies and AGN: and respectively, derived according to the 16th, 50th and 84th percentiles of the total samples.

For a proper interpretation of the width of the distribution, whether it corresponds to measurement errors or intrinsic spread of the spectral index distribution, we fit a Gaussian distribution to the total samples. Using an MCMC algorithm (emcee: Foreman-Mackey et al., 2013) we infer the parameters and (mean value and intrinsic spread) that characterise the intrinsic distribution of the parameters, independently of the measurement errors. In order to avoid that extreme outliers have a strong effect on the fit we restrict the fit on 99.7 percent of the population distributed around the median value. We choose this very conservative cut in order not to alter the shape of the distributions, which are characterised by longer tails than a normal distribution. The logarithm of the likelihood function obtained through the fitting of these parameters, assuming the measurement errors are Gaussian, is given by


where (spread of the intrinsic distribution) and (measurement errors) are coupled. Strictly speaking, the distribution described by Eq. 2 is not a Gaussian, but a weighted sum of Gaussians with varying widths. To visualize the shape of the intrinsic distribution, Gaussian distributions with the and intrinsic inferred above are over-plotted on the histograms of Fig. 7.6

Figure 8: Spectral indices and for SF galaxies (upper panel) and AGN-dominated galaxies (lower panel). Detected and undetected sources in 1.4 GHz are depicted in red and orange respectively for SF galaxies, and dark blue and sky blue for AGN. Median and scatter values corresponding to the 16th, 50th and 84th percentiles are calculated for all detected sources in VLA-P and WSRT ( per cent of the total). The solid line displays the fitted Gaussian with mean and intrinsic scatter . A schematic representation of the observed curvature is also shown in the central panel for SF galaxies (red line) and AGN (blue line).

The distribution of the SF galaxy sample is best fitted with a mean value and intrinsic scatter of , while the spectral index distribution for AGN presents and . Although the median values of the SF-galaxy population are slightly higher than the AGN population, and the fitted mean values show a similar trend, both populations are consistent within the error bars. To test the difference between both distributions for significance we used the non-parametric two-sample Kolmogorov-Smirnov (KS) test and found that the distributions are statistically different at a confidence level greater than 99.9 per cent for all sources and even greater when considering only detections. In conclusion we find that using the total spectral index to characterize the radio SEDs of galaxies, a statistically significant difference between starburst-dominated and the AGN-dominated sample is found, where the total sample of SF galaxies show a slightly steeper spectra than AGN with a difference in their median values of . The resulting values agree with previous calculations of the low-frequency spectral slope in the literature (starbursts: Condon, 1992; Ibar et al., 2008; Ivison et al., 2010b; Marvil et al., 2015), (AGN: Singh et al., 2013).

Figure 9: Spectral indices , and for SF galaxies (upper panel) and AGN-dominated galaxies (lower panel). These distributions corresponds to the fraction of our sample with GMRT coverage (see Fig. 1), which makes up per cent of the total sample. Detected and undetected sources in 1.4 GHz are depicted in red and orange respectively for SF galaxies, and dark blue and sky blue for AGN. Values corresponding to the 16th, 50th and 84th percentiles are calculated for all detected sources in VLA-P and WSRT ( per cent of the total). A schematic representation of the observed curvature is also shown in the central panel for SF galaxies (red line) and AGN (blue line).

5.2 Spectral indices and

Higher order spectral properties such as curvature and higher order features can be studied by adding the central P-band observations at 325 MHz to the data set.

Biased curvature observations due to non-detections

Before including the central data point, we need to model our expectations according to the characteristics of the available data. Following the distribution of presented above, a single slope spectrum around and for galaxies and AGN would imply observing similar values for the spectral indices and , when adding the central data point at 325 MHz. However, measurement uncertainties can easily alter these observations and produce biased conclusions on curvature. An important bias is the detection level of the central band, which can severely compromise the interpretation of steep spectra undetected in this band, shifting the median to flatter spectra towards lower frequencies. We simulated the behaviour of our sample of LOFAR selected sources assuming that their spectral indices follow the ones observed in the histograms of Fig. 7, and proceed to mimic similar systematics to those of our sample to investigate the expected values after the addition of the 325 MHz data point. The distributions resulting from the simulations show a clear offset in these values compared to the initial distribution; demonstrating the biasing effect of using upper limits or excluding undetected sources without a proper treatment.

We would like to point out that this is important for the interpretation of results on radio continuum studies constructed from data of diverse sensitivity similar to ours (see e.g.: Marvil et al., 2015; Clemens et al., 2010). To overcome these biases we proceed with a cut in LOFAR fluxes and select sources above 2 mJy, reducing the percentage of undetected sources in VLA-P from 41 per cent to 14 per cent for SF galaxies and 27 per cent to 10 per cent for AGN, leaving us still with a large enough sample for a statistical study with 189 SF galaxies and 421 AGN.

Results on flux-selected sources with S

Histograms of the low and high frequency spectral indices for SF galaxies and AGN are presented in Fig. 8. We quote the median values and scatter based on the 16th and 84th percentiles of the distributions and fit a Gaussian to the distribution of sources with VLA-P detections (red and dark blue areas for galaxies and AGN respectively).

The upper panel of Fig. 8 shows the median and scatter values of the spectral index distributions for SF galaxies: and . A simple comparison of the median values of the two different spectral index distributions (similarly, using mean values from the fitted Gaussian curves) shows in general a slight flattening towards lower frequencies. Using the non-parametric two-sample KS-test, we find that the difference between the low- and high frequency distributions of SF galaxies is highly significant at a level greater than 99.99 per cent (p-value). Our statistical study thus concludes that the spectral continuum for our SF galaxy sample is consistent with a curved spectrum, showing a slight flattening towards lower frequencies.

Similarly, the lower panels of Fig. 8 present the median and scatter values for the population of AGN: and . In contrast to the SF galaxies population, the medians of the distributions and fitted mean values in the lower panels of Fig. 8 show a steepening towards lower frequencies with a difference in the mean values of the Gaussian fits of . Also here, we test this observation for statistical significance using the two-sample KS-test. We find that the sample of galaxies hosting AGN present curvature in their radio SED with a confidence level greater than 99.99 per cent (p-value). Our study thus concludes that radio SEDs of AGN in our sample show a statistically significant steepening in their radio continuum going to lower frequencies.

These results will be discussed in a physical context in section 5.6. We would like to remark that the significant differences found between the spectral curvature of SF galaxies and AGN in Fig. 8 confirm that the observed curvature is a real spectral feature and not due to calibration issues in our sample.

5.3 Spectral indices , and

We take advantage of the partial availability of GMRT data at 610 MHz and construct low-frequency radio SEDs with four data points for the fraction of sources that lie inside the GMRT coverage (green region in Fig. 1). The GMRT subsample constitutes 24 per cent of the total sample of LOFAR-I-band selected sources with LOFAR fluxes mJy, which implies a total of 198 sources. Only 11 per cent have VLA-P non-detections and are replaced by forced-photometry. Taking into account detected sources this sub-sample consists of 41 SF galaxies and 101 AGN.

As expected from a solely spatial cut, the spectral index distributions present similar median values to the total field, with general spectral indices for SF galaxies of , and for AGN of . Adding one more data point to the radio SEDs, we calculate the three spectral indices , and following Eq. 1. Constructing the total SED, as sketched in the central panel of Fig. 9, we are able to recover spectral index ratios which suggest a consistent behaviour (within the larger statistical errors) with our previous finding in Fig. 8.

Figure 10: Evolution of spectral index as a function of redshift for SF galaxies (red) and AGN (blue). The data points are the mean values of , grouped in bins of equal density. The error bars represent bootstrap errors.

The upper and lower panels of Fig. 9 show three spectral index distributions for low, medium and high frequency pairs for both SF galaxies and AGN, respectively. While the distributions’ median values for SF galaxies ( and ) suggest a slight tendency to have flatter spectra towards low frequencies, the distributions’ medians for AGN ( and ) are consistent with a steepening of the continuum going towards lower frequencies. A quick examination to these characteristic values suggest these results are similar to the results of section 5.2.2, finding curvature as sketched in the middle panel of Fig. 8.

To test this systematic curvature for statistical significance, a KS-test on the difference of the spectral index distributions was performed. Comparing all three spectral index distributions at adjacent frequency pairs we find that the curvature is not statistically significant, since the KS-test could not reject the null-hypothesis due to the small statistics of the sub-sample with GMRT coverage (40 SF galaxies and 102 AGN). However, a KS-test on the distributions of the spectral indices at the extremes of the frequency coverage, and , shows that a difference between these distributions is statistical significant for the AGN population at a level of 98 per cent, consistent with curvature and confirming the results of the previous section. For the SF-galaxy population this is not the case. This may be due to the small sample used for the study in this section (40 sources), which weakens the statistics.

5.4 Redshift evolution of the spectral curvature

In Figs. 10 and 11 we study the redshift evolution of the total spectral index and the curvature parameter , respectively. For this study we use the samples with detections in both VLAP and WSRT (90 per cent of the sample with S ), which correspond to the characteristic values estimated from the distributions in Fig. 8.

Investigating the change in total slope as a function of redshift is important to validate the approximation of a single power law as the k-correction for radio luminosities. This assumption holds as long as the SED slope remains constant at the rest-frame frequencies for low redshift (150 MHz-1.4 GHz at ) and high redshift (450 MHz-4.2 GHz at ) sources. Fig. 10 shows that evolution of the slope as a function of redshift cannot be clearly recognized for either populations but shows a large scatter around the canonical value of . Specifically for AGN, some data points in Fig. 10 may suggest an unclear trend as a function of redshift. However, using a linear fitting to the data points we find a redshift evolution very close to 0 (slope= ). We also tested the robustness of this finding using different binning schemes and the result remains consistent. This result is expected since the rest-frame frequency range covered by our sample should be completely dominated by synchrotron emission and the strong flattening from thermal emission should not be visible below frequencies around 10 GHz (e.g. Gioia et al., 1982).

Figure 11: Evolution of the spectral curvature parameter () as a function of redshift for SF galaxies (red) and AGN (blue). The data points are the mean values of the curvature parameter, binned in groups of equal density in redshift. The error bars represent bootstrap errors.

A non-evolving mean spectral slope disfavours a scenario where the synchrotron spectral shape would be significantly affected by large scale properties of the galaxy that present a strong redshift evolution, such as integrated luminosities or SFRs. On the contrary, the absence of such an evolution suggests that the synchrotron emission would be determined by other rather local properties affecting the CR injection mechanisms that are not coupled to an evolution in redshift. Examples of such processes are CR acceleration mechanisms and energy losses confined to sub-kpc structure around individual supernova remnants, galactic magnetic fields, and the immediate surrounding interstellar medium. We suspect these local properties establish the shape of the radio SED, while differences in global properties may primarily affect the normalisation of the spectrum. This result is consistent with the findings of Magnelli et al. (2015) and Ivison et al. (2010a), who found no significant evolution of the spectral index with redshift across for SF galaxies.

Fig. 11 shows the redshift evolution of the spectral curvature, parametrized as the difference between spectral indices and , binned in six different redshift bins of equal sample sizes. The most remarkable observation is that a contrast between the curvature of both populations is clearly present at all redshifts. While SF galaxies show continuously positive values of , AGN present negative values, consistent with Figs. 8 and 9.

We find no clear evolution of the curvature with redshift for both SF galaxies and AGN populations, as shown in Fig. 11. We quantify this by fitting straight lines of different slopes for both binned populations. We find that the slope for the SF galaxies population is consistent with zero (), while the slope for AGN is found to be . We test these results for significance by calculating the coefficient of determination statistics , which describes the proportion of the variance in the curvature parameter that is predictable from the redshift evolution given by the fit. From the results of the coefficient of determination we find that the fits for the SF galaxies and AGN populations account for less than 13 and 40 per cent of the scatter respectively. Given the very small slope values, the relatively large errors and the low values of , we conclude that we find no significant redshift evolution of the spectral curvature.

Finally, we explore the redshift dependence of the total sample if undetected sources were also included. Although non-detections produce no significant differences to the total characterization of the curvature due to their low fraction (10 percent), we find that they imprint a weak redshift dependence on the results. This effect occurs since luminosities of undetected fluxes are calculated as the product of a narrow distribution of fluxes close to the noise limit and the squared luminosity distance which imprints a strong redshift dependence. This observation is especially important for samples of large non-detection fractions, as we will discuss in section 6.

5.5 Spectral curvature dependence on other properties

To investigate the origin of the spectral curvature for galaxies and AGN, we studied its dependence on the different physical properties inferred by SED-fitting. For the SF galaxy sample, no dependence on IR luminosity was found and consequently on . Similarly we investigated the evolution of spectral curvature as a function of specific star formation rate sSFR, estimated by the SED fitting routine, and found no significant dependence.

For the AGN sample, a significant evolution was found as a function of the integrated luminosity of the torus component , which represents a proxy for the total AGN power. As shown in Fig. 12, this suggests that more powerful AGN seem to present weaker curvature, i.e. flatter spectra. We quantify this by fitting straight lines of different slopes to the z-evolution of the binned AGN population presented in Fig. 12. We estimate the slope to be , where the fit accounts for 93 per cent of the scatter according to the coefficient of determination of this fit. Since the spectral curvature of AGN does not present a redshift dependence as significant as this luminosity dependence (Fig. 11), our results suggest that AGN-luminosity is the driver of this spectral curvature trend. We discuss these observations in section 5.6.

5.6 Comparison with previous spectral index studies

In this section radio SEDs from 150 MHz to 1.4 GHz have been constructed for the total population of LOFAR-selected sources with fluxes above 2 mJy. Our results suggest the presence of oppositely directed curvature in the radio continuum of star-forming compared to AGN-dominated galaxies. While SF galaxies present a slight flattening towards lower frequencies suggesting positive curvature, AGN radio SEDs reveal a systematic steepening towards lower frequencies, implying negative curvature. Moreover, no redshift evolution of the spectral curvature was found for both the SF galaxy and AGN populations. Comparing our results to previous findings in the literature is not a trivial task due to the unparalleled depth of our observations, particularly at low frequencies, which allows us to study the radio SED properties of sources to higher redshifts and lower luminosities than previous studies.

SF galaxies

Figure 12: Curvature observed in AGN plotted against the torus luminosity , integrated in the wavelength range 1-40 m. The data points are the mean values of the curvature parameter, binned in groups of equal density in . The error bars represent bootstrap errors.

Several studies of the radio SEDs of nearby SF galaxies agree with our results that positive curvature is not an unusual spectral property at low frequencies (e.g. Marvil et al., 2015; Murphy et al., 2013; Williams & Bower, 2010; Clemens et al., 2010; Condon, 1992; Israel & Mahoney, 1990), although few differences are to be noted. For instance, the average spectrum of Marvil et al. (2015) nearby sources is curved following a change of per logarithmic frequency decade towards lower frequencies, which is a stronger flattening than that found in this study (). One possible explanation is that Marvil et al. (2015) include a large fraction of non detections at the central wavelengths ( per cent for two of the four frequencies sampled), which is an important source of biased curvature if the results of our simulation in Section 5.2.1 applies here as well. Similarly, Clemens et al. (2010) studied 20 nearby LIRGS and ULIRGS from 244 to 8.4 GHz, finding complex curvature in their source as well as claiming not only flattening but a strong turn over at around 1 GHz, which when extrapolated towards lower frequencies, implies a positive spectral index at frequencies GHz. With the availability of LOFAR data at 150 MHz we can confirm this is not the case for most of the SF galaxies present in our field but the curvature is rather small, presenting a slightly flatter slope but still of negative sign.

Spectral flattening in SF galaxies hints towards physical conditions which, integrated over the total galaxy, appear to be a common property in our SF galaxies sample and are expected to remain constant at all redshifts. These conditions can be intrinsic to the production of the cosmic ray electrons (CRE) responsible for the synchrotron emission component (Basu et al., 2012) or environmental, connected to the extreme interstellar medium (ISM) in star-forming regions into which the CREs are injected (Lacki, 2013). Consistent with the second hypothesis, Murphy et al. (2013) find increasing flattening of radio spectral index of galaxies with increasing specific star formation rate (sSFR). We find such a dependence neither for sSFR nor for SFR. Similarly, Basu et al. (2015) find that the flattening of the spectral index of nearby SF galaxies is unlikely to be caused by thermal free-free absorption but instead suggest a dependence on gas surface density. Due to the lack of gas surface density tracers in our sample, we can not test for this.


Other attempts have been made in the construction of radio SED for radio AGN, but these have been limited to small samples and/or single spectral index studies (e.g. Laing & Peacock, 1980; Singh et al., 2013; Mauch et al., 2013; Kharb et al., 2016; Hardcastle et al., 2016). Understanding the physical processes which drive integrated radio spectra of AGN is not an easy task due to the multiple physical conditions contributing to its total emission, as multi-component structures and orientation effects, and hence studies have focused on spectral index maps ( maps) of spatially resolved sources (Vardoulaki et al., 2015; Harwood et al., 2015). For instance, in the case of FRI and FRII sources, the spectral shape can be dramatically different between jet, lobes to hot spots due to spectral ageing (Harwood et al., 2015). While hot-spots are sources of constant CRE supply and present flat spectral indices (Laing & Peacock, 1980), lobes consist of earlier-ejected material which has lost energy, producing considerably steeper radio spectral indices. The overall spectrum thus depends on the history and environment of the radio source.

An important question to understand the integrated radio spectra is which of the physical components dominate the total emission, at different frequency ranges, and across different orders of magnitude in luminosity? A basic consensus among the observations is that radio AGN present generally flatter spectra in the GHz regime than SF galaxies (Murphy et al., 2013), similar to our finding in Section 5 and, mainly for sources at and frequencies , in Fig. 8. As hot-spots can be the dominant source in total flux density in luminous sources (Jenkins & McEllin, 1977), this could suggest that these frequencies are dominated by hotspot emission. Consistent with our observation in Fig. 12, Laing & Peacock (1980) finds a luminosity dependence of the integrated spectral shapes for a sample of radio sources selected at 178 and 2700 MHz, arguing that more powerful sources have flatter spectra than weaker sources, while weaker sources have spectra which steepen at low frequencies (but see also Rees et al., 2016, for different results). While their study is limited to a much smaller sample and luminosity range due to their shallower data, we are able to confirm the luminosity dependence, ruling out a redshift dependence. The steepening at low frequencies observed for our weaker AGN sources could be explained by steep-spectrum components dominating in this frequency regime, while flat-spectrum components become relatively more important at higher frequencies, causing the spectral index to flatten. These observations are in good agreement with the results presented by Whittam et al. (2016). They report spectral flattening for radio-selected galaxies at the high-frequency end (15.7 GHz) and suggest this may be due to the cores of Fanaroff & Riley (1974) class I sources (FRI) becoming dominant at these high frequencies.

However, considering high flux uncertainties is especially relevant for AGN due to the different combinations of resolution and source extraction methods used, as they are extended sources. In contrast, this is not a relevant issue for SF galaxies, which are mostly compact sources. These uncertainties might weaken the strength of our conclusions on AGN and better radio data and larger samples are needed to confirm this finding.

6 The infrared-radio correlation (IRC)

A remarkable property intrinsic to the radio SED of SF galaxies is the tight correlation with infrared luminosity, which has been observed at a few radio bands across five orders of magnitude in luminosity (Yun et al., 2001). In this section we will characterize the IRC for the 1.4 GHz and the LOFAR data at 150 MHz and investigate its evolution with redshift for our sample. For this study we include all detected and undetected sources in both the 1.4 GHz WSRT and the Herschel-SPIRE data. As explained in detail in section 3.2.2, undetected fluxes are estimated using forced photometry measurements on the positions of LOFAR-detected sources. In this section we discard AGN and consider only sources selected as SF galaxies from the total LOFAR-I-band sample, which comprises 810 galaxies.

The IRC is parametrized by the widely used value , which is defined as the ratio


where is the total rest-frame infrared luminosity (in ) of the cold dust SED template integrated along the rest-frame frequency range m (equivalent to Helou et al., 1985; Bell, 2003; Ivison et al., 2010b; Sargent et al., 2010; Bourne et al., 2011; Murphy et al., 2011). L is the luminosity at the radio frequency to be studied, in this case 1.4 GHz or 150 MHz in . The factor is the frequency corresponding to m, used in the definition to make a dimensionless quantity.

Figure 13: Histograms representing the distributions of for the subsample of SF galaxies only. The values are calculated

Radio rest-frame luminosities are computed using k-corrections inferred assuming , where is the mean value observed for the SF galaxies distribution (see Fig. 7). We use a fixed value for this study, since the available data does not sample the radio SED curvature accurately enough for a useful k-correction.

In contrast to most previous studies, our unique infrared coverage and the large dust models library allow us to robustly compute the total infrared luminosities. This approach has several advantages over using bolometric corrections based on monochromatic luminosities and single models as discussed by Smith et al. (2014). Moreover our MCMC approach in SED fitting allows us to infer robust uncertainties for the luminosities and consequently for q.

Figure 14: q-value for the IRC corresponding to radio luminosities at 1.4 GHz plotted against redshift. Black error bars depict the observed values within the 2 region of the -distribution in Fig. 13, while red lines correspond to the fitted q-values inferred by the equation in the legend, taking into account uncertainties on the parameters calculated through MCMC sampling. The pink shaded area correspond to the fitted intrinsic scatter of the correlation.

6.1 The IRC at 1.4 GHz and its redshift evolution

While 1.4 GHz luminosities of SF galaxies clearly follow the IRC (Fig. 6), in the upper panel of Fig. 13 is scattered in a distribution with a median value of , where the error bars represent the error on the median. Fig. 13 reveals a small number of outlier galaxies exhibiting a radio excess, which is the consequence of contamination as part of our IR-based classification. Note that these were carefully excluded for the redshift-evolution study in the next section. This distributions agrees well with the literature within the errors, although it is slightly lower than most of the local values observed for local SF galaxies (e.g. using FIR emission: Yun et al. (2001), Helou et al. (1985), using total IR emission as in our case: Bell (2003) ). In this section we investigate which parameter drives this offset and is able to explain the scatter observed.

In Fig. 14 we test the evolution of (1.4 GHz) as a function of redshift. To quantify a redshift evolution we fit the function: , where is a constant and the variable driving the evolution. We include only 86 per cent of sources, which lie within the region of the distribution in Fig. 13, so that the fit is not driven by outliers. To account for measurement uncertainties in both variables and intrinsic scatter, the linear regression to infer the parameters is based on the likelihood function




The fitting results for q in Fig. 14 suggest an evolution with redshift given by


with an intrinsic scatter of .

As expected, after the fit for redshift evolution the value at is higher than the median value of the total sample by . To infer the proportion of the variance in the distribution of Fig. 13, which is explained by the -evolution, we infer the coefficient of determination , defined in this case as,


where q are the observed q-values and q are values calculated following Eq. 6. We calculate that the redshift evolution accounts for per cent of the variance in the total distribution, while most of the intrinsic scatter is still not explicable.

A dependence of the q evolution on the evolution around 1.4 GHz is unlikely, since this would imply a clear evolution of in Fig. 10, which is not observed. Moreover, since the spectral curvature discussed in our study is observed only at lower frequencies this should not affect the k-correction for the 1.4 GHz luminosities. To investigate whether this redshift dependency is rather a luminosity dependence of the we separated the total sample into luminosity bins of similar sample sizes and were able to confirm a redshift evolution inside most of the luminosity bins, discarding such a degeneracy (see plots in Appendix A for details). Finally, we test whether the z-evolution of is related to a changing ratio of AGN and SF contribution to the IR emission (). We find no correlation between these quantities indicating that the redshift evolution is not due to AGN contamination.

We discuss possible explanations to the offset observed between the mean of our sample () and the local mean IRC found by (Bell, 2003, )). An important factor to take into account in IRC studies is its sensitivity to selection biases (Sargent et al., 2010). A common selection bias in radio-selected samples arises from samples being incomplete at IR data for sources with faint radio fluxes (Ibar et al., 2008; Sargent et al., 2010). This might not affect our sample, since the FIR emission expected for the 150 MHz fluxes in our sample (based on the IRC) is considerably above the flux limit of the XID SPIRE fluxes (Roseboom et al., 2010). Non-detections comprise less than 17 per cent and are compensated using forced photometry. Moreover, as discussed above, no direct dependence of the z-evolution of on radio luminosity could be found (Appendix). We note that by using a radio-selected sample we might be losing the FIR-dominated sources which drop out at the radio at the high redshift end of our sample and could shift the distribution to lower values. However, Casey et al. (2012) remarked that the XID SPIRE fluxes are complete in relation to a pure SPIRE selected sample and that the fraction of sources that could drop out in the radio at is negligible. Since per cent of the sources in our sample have redshifts , the possible bias against faint radio sources is not significant.

The mean value of our total distribution is lower () than the local mean IRC presented by Bell (2003). However, from the fitted redshift-evolution in Equation 6 the predicted value corresponding to the local universe would be which is well within the error on the mean quoted by (Bell, 2003, )). The tension between our observations and those by Bell (2003) is therefore reduced. The offset between the mean values can be primarily explained by the redshift evolution and the different redshift ranges probed in the two studies.

Figure 15: q-value for the IRC corresponding to radio luminosities at 150 MHz plotted against redshift. Black error bars depict the observed values within the 2 region of the -distribution in Fig. 13, while red lines correspond to the fitted q-values inferred by the equation in the legend, taking into account uncertainties on the parameters calculated through MCMC sampling. The pink shaded area correspond to the fitted intrinsic scatter of the correlation.

The picture of decreasing values as a function of redshift is being supported by a number of studies based on FIR and radio data sets of variate properties. Although some theoretical studies predict the opposite trend (Schober et al., 2016) and few observational studies claim to find no observable trend (Pannella et al., 2015), a long list of recent Herschel-FIR and Spitzer-IR based studies have observed a redshift evolution of (e.g., Ivison et al., 2010b; Bourne et al., 2011; Casey et al., 2012; Jarvis et al., 2010; Smith et al., 2014; Magnelli et al., 2015; Delhaize et al., 2017). However, some of these studies have chosen to favour a no-trend scenario due to the modest evolution compared to the large size of their measurement errors or other observed dependences. It is reassuring that our results are consistent with studies that take into account the treatment of biases due to selection as presented by Ivison et al. (2010b), which compare different selections to test their results and Magnelli et al. (2015), which is based on a mass-selected galaxy sample. They find an evolution that goes as and respectively. Similarly, Delhaize et al. (2017) find a redshift evolution of using highly sensitive 3GHz VLA observations and Herschel-FIR data and accounting for radio and IR non-detections through survival analysis.

6.2 The IRC at 150 MHz and its redshift evolution

One main focus of this study is to investigate for the first time the IRC for radio emission at low frequencies, which is crucial for LOFAR observations and relevant as well for future SKA data. Similarly to the previous section, we calculate the value of the IRC parameter for the 150 MHz data and present the distribution of in the lower panel of Fig. 13. The distribution of presents a median value and error on the median of .

It is customary to calculate the value of for different radio frequencies interpolated from the well calibrated value at 1.4 GHz (e.g., Yun et al., 2001; Murphy et al., 2011) and under the assumption of a single power law radio emission. Such a value would be calculated as


Assuming a unique power law value from 1400 to 150 MHz of , we obtain an expected q-value of . This value is well in agreement with the distribution observed in the lower panel Fig. 13 () given the errors.

Also here we consider sources with -values within the 2 region of the distribution in Fig. 13 and investigate the evolution of with redshift. We find the relation


with an intrinsic scatter of . Proceeding with a similar analysis as for , we calculate the values for this regression and obtain that the redshift evolution accounts for per cent of the variance in the total distribution, while most of the intrinsic scatter is still not accounted for.

One remarkable property of the redshift evolution of is that its redshift dependence, , is slightly stronger () than that observed for , with . We explore the connection between the stronger redshift dependence of and our results on the spectral curvature described in section 5.4.

Figure 16: SFR inferred from the total infrared luminosity as a function of radio luminosity at 150 MHz. Scatter points are the observed median SFR values of all SF-sources, colour-coded by redshift bin. Redshifts were binned in groups of equal density, where median values of each bins are labeled in the legend. The transparent lines following the same colour-coding correspond to the calculated SFR values using Eq.11.The dotted line shows the SFR values, when a non-evolving q-value is assumed.

In Fig. 11 we have shown that the spectral curvature is not a function of redshift when considering only the detected sources in 1.4 GHz (90 per cent of total). We added that the inclusion of the remaining fraction of undetected sources (10 per cent) would imprint a weak redshift dependence on the results (section 5.4). However, in that case we clarified that this finding produced no significant differences to the total characterization of the curvature. In contrast to section 5.4, for the present IRC study we use the total LOFAR-I-band selected sample for which 33 percent of sources have non-detections at 1.4 GHz (Table 1). These non-detections will have a stronger impact here than in the spectral curvature study.

We quantify this impact by investigating the redshift evolution of subsamples of 1.4 GHz detections and non-detections exclusively. For the subsample of detected sources we find an evolution of , which is a similar result to the one obtained from . The study comprising exclusively sources undetected in 1.4 GHz returns an evolution of . From these results we conclude that the weaker redshift evolution of compared to is an effect of non-detections in 1.4 GHz rather than a result of any evolution in the spectral curvature (see Fig. 11, which illustrates the lack of redshift evolution). We thus conclude that the redshift evolution in both IRC parameter distributions are consistent to a large extent. These observations suggest that the strength of the redshift dependence found for can be considered a lower limit and that a stronger evolution is possible.

7 The low frequency radio luminosity as a SFR tracer

Calibrating non-thermal radio continuum emission as a SFR tracer is done mostly based on the IRC due to the reliability of the IR emission as an extinction-free SFR diagnostic (see also Brown et al. (in prep.) and Gürkan et al. (in prep.) for other studies on this at low redshifts). Since this is especially important for our sample due to the dusty environments observed in highly SF galaxies, we use the total IR data (, see section 4 for details ) to compute the SFRs throughout the paper following:


as presented as well by Murphy et al. (2011) and Bell (2003). Combining Equations 10, 3 and 9 one can express the 150 MHz luminosity as a SFR as:


where the factor is equivalent to , , and for redshifts around 0.1, 0.5, 1.0 and 2.0, respectively.

In Fig. 16 we show how the redshift evolution of the -value suggests a more accurate SFR estimation compared to assuming a unique value. We over-plot the observed median SFR values of all SF-sources (scatter points) to the inferred SFR values using Eq.11 (transparent lines), both colour-coded by redshift bins. In contrast we plot the SFR values, as assuming a unique -value at all redshifts as a dotted black line. As can be seen in Fig. 16, values derived from Eq.11 give a better description of the data.

8 Summary

In this paper we have investigated two properties of radio selected SF galaxies and AGN, namely their radio spectral energy distributions and the IR-radio correlation (IR), and how these evolve across the bulk of cosmic history. Our study is based on a radio and I-band selected sample from the Boötes field and comprises galaxies at redshifts of . Taking advantage of the rich multi-wavelength coverage available for the Boötes field, we constructed panchromatic SEDs for our sample, incorporating UV/optical to far-infrared observations. The sources were then classified into star-formation and AGN-dominated galaxies using a multi-wavelength SED-fitting method (AGNfitter code, Calistro Rivera et al., 2016).

We have characterized the statistical radio spectral index and spectral curvature properties for our total samples of SF galaxies and AGN using deep radio data across a wide range in radio frequency (150, 325, 610 MHz and 1.4 GHz). Additionally, intrinsic physical properties of the galaxy and AGN were investigated by calibrating spectral properties with parameters inferred from multi-wavelength SED-fitting. We arrive at the following conclusions:

  • SF galaxies exhibit an average radio spectral slope of . In comparison, AGN were found to exhibit a flatter spectral slope, . Although the difference is small compared to the width of the observed distributions we find it to be statistically significant.

  • Observations of radio SEDs indicate a difference in the curvature of the radio continuum for SF galaxies compared to AGN-dominated galaxies. SF galaxies exhibit a slightly curved radio spectrum which flattens at frequencies 325 MHz. Conversely, AGN-dominated radio galaxies display a curved radio spectrum which shows a systematic steepening towards low frequencies below 1.4 GHz. In both cases we find these results to be significant based two-sample KS-tests.

  • No evolution of the curvature with redshift is observed for both the SF galaxies and AGN samples.

  • In SF galaxies we found no dependence of the curvature on star-formation rate or specific star-formation rate.

  • A correlation between the radio spectral curvature and torus luminosity is observed for AGN, with low luminosity AGN exhibiting a higher degree of steepening than high luminosity sources.

In addition to studying the radio SEDs, the extensive mid- to far-IR data available for our sample allowed us to also study in detail the infrared radio correlation (IRC) for the star-forming galaxies:

  • An evolution with redshift of the IRC for 1.4 GHz radio luminosities is observed, following , where is the ratio between the total infrared luminosity () and the 1.4 GHz radio luminosity (). The observed evolution accounts only for per cent of the variance of the distribution. The large observed intrinsic scatter of is still predominantly unexplained. We discuss that the inclusion of radio non-detections in this sample may bias the result towards a weaker redshift dependence and thus conclude that this evolution can be considered a lower limit for the IRC at 1.4 GHz.

  • Similarly, an evolution of the IRC with redshift for 150 MHz radio luminosities is found following . A comparison of the observed values and those extrapolated from , as a function of redshift, give similar results but with some scatter around .

  • After converting to to calibrate the low frequency radio luminosity as a measure of SFR, we show that the observed IRC redshift evolution is potentially an important factor for the future use of radio-luminosity to estimate un-biased star-formation rates.

9 Acknowledgements

The authors thank the anonymous referee for the very careful reading of the manuscript and the many insightful comments and suggestions that improved the paper. We thank Jacqueline Hodge, Frank Israel and Michael Brown for useful discussions and suggestions. GCR, KJD and HJR acknowledge support from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007- 2013) /ERC Advanced Grant NEW-CLUSTERS-321271. WLW and MJH acknowledge support from the UK Science and Technology Facilities Council [ST/M001008/1]. PNB is grateful for support from the UK STFC via grant ST/M001229/1. MJJ acknowledges support from UK Science and Technology Facilities Council and the South African SKA Project EKM acknowledges support from the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020. GJW Gratefully acknowledges support from The Leverhulme Trust. This research made use of ASTROPY, a community-developed core Python package for astronomy (Astropy Collaboration et al. 2013) hosted at http://www.astropy.org/.

Appendix A Luminosity vs redshift dependence of the IRC

Figures 17 and 18 show that the redshift evolution is found also within subsamples of similar luminosity, suggesting the trend observed in Figs. 14 and 15 is not the consequence of a dependence on luminosity.

Figure 17: q-value for the IRC corresponding to the radio at 1.4 GHz plotted against redshift for different luminosity bins. Error bars correspond to bootstrap errors.
Figure 18: q-value for the IRC corresponding to the radio at 150 MHz plotted against redshift for different luminosity bins. Error bars correspond to bootstrap errors.

Appendix B SED-fitting examples

Examples of SED fitting of galaxies in our sample using the AGNfitter code.

Figure 19: SED-fitting examples: Photometric data points for each source are plotted as black circular markers with error bars. We pick 8 different realizations from the parameters posterior probability distributions and over-plot them in order to visualize the effect of the parameters’ uncertainties on the SEDS. The SED shapes of the physical components are presented as solid lines: the galactic cold dust emission (green), the hot dust emission from the AGN torus (purple), the stellar emission (orange) and the accretion disk emission (blue).The linear combination of these, the ’total SED’, is depicted as a red line.


  1. pagerange: The LOFAR window on star-forming galaxies and AGN – curved radio SEDs and IR-radio correlation at 9
  2. pubyear: 2002
  3. After the flux cut applied for the curvature analysis (section 5.2.2)
  4. Forced photometry included for non-detections
  5. Due to spatial partial coverage, see Fig. 1
  6. The slight apparent mismatch between the fit and the data (especially for the AGN sample) is explained by the presence of the small bulk of outlier sources at the high frequency tail of the spectral index distribution, by undetected sources (orange / sky blue bars) being weighted less than detections (red /dark blue bars) due to their larger errors (following Eq. 2) and by the fact that a Gaussian fit is only an approximate description of the distribution.


  1. Appleton P. N., et al., 2004, ApJS, 154, 147
  2. Ashby M. L. N., et al., 2009, ApJ, 701, 428
  3. Autry R. G., et al., 2003, in Iye M., Moorwood A. F. M., eds, Proc. SPIEVol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes. pp 525–539, doi:10.1117/12.460419
  4. Basu A., Roy S., Mitra D., 2012, ApJ, 756, 141
  5. Basu A., Beck R., Schmidt P., Roy S., 2015, MNRAS, 449, 3879
  6. Becker R. H., White R. L., Helfand D. J., 1995, ApJ, 450, 559
  7. Bell E. F., 2003, ApJ, 586, 794
  8. Bertin E., Arnouts S., 1996, A&AS, 117, 393
  9. Best P. N., Heckman T. M., 2012, MNRAS, 421, 1569
  10. Bianchi L., Conti A., Shiao B., 2014, Advances in Space Research, 53, 900
  11. Bonzini M., Padovani P., Mainieri V., Kellermann K. I., Miller N., Rosati P., Tozzi P., Vattakunnel S., 2013, MNRAS, 436, 3759
  12. Bourne N., Dunne L., Ivison R. J., Maddox S. J., Dickinson M., Frayer D. T., 2011, MNRAS, 410, 1155
  13. Brammer G. B., van Dokkum P. G., Coppi P., 2008, ApJ, 686, 1503
  14. Brown M. J. I., Dey A., Jannuzi B. T., Brand K., Benson A. J., Brodwin M., Croton D. J., Eisenhardt P. R., 2007, ApJ, 654, 858
  15. Brown M. J. I., et al., 2014, ApJS, 212, 18
  16. Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
  17. Calistro Rivera G., Lusso E., Hennawi J. F., Hogg D. W., 2016, ApJ, 833, 98
  18. Calzetti D., Kinney A. L., Storchi-Bergmann T., 1994, ApJ, 429, 582
  19. Carrasco Kind M., Brunner R. J., 2014, MNRAS, 442, 3380
  20. Casey C. M., et al., 2012, ApJ, 761, 140
  21. Casey C. M., Narayanan D., Cooray A., 2014, Phys. Rep., 541, 45
  22. Chabrier G., 2003, PASP, 115, 763
  23. Chary R., Elbaz D., 2001, ApJ, 556, 562
  24. Clemens M. S., Scaife A., Vega O., Bressan A., 2010, MNRAS, 405, 887
  25. Condon J. J., 1992, ARA&A, 30, 575
  26. Condon J. J., Cotton W. D., Greisen E. W., Yin Q. F., Perley R. A., Taylor G. B., Broderick J. J., 1998, AJ, 115, 1693
  27. Cool R. J., 2007, ApJS, 169, 21
  28. Coppejans R., Cseh D., Williams W. L., van Velzen S., Falcke H., 2015, MNRAS, 450, 1477
  29. Dahlen T., et al., 2013, The Astrophysical Journal, 775, 93
  30. Dale D. A., Helou G., 2002, ApJ, 576, 159
  31. Delhaize J., et al., 2017, preprint, (arXiv:1703.09723)
  32. Delvecchio I., et al., 2017, preprint, (arXiv:1703.09720)
  33. Dumas G., Schinnerer E., Tabatabaei F. S., Beck R., Velusamy T., Murphy E., 2011, AJ, 141, 41
  34. Fanaroff B. L., Riley J. M., 1974, MNRAS, 167, 31P
  35. Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  36. Gioia I. M., Gregorini L., Klein U., 1982, A&A, 116, 164
  37. Hardcastle M. J., 2009, in Saikia D. J., Green D. A., Gupta Y., Venturi T., eds, Astronomical Society of the Pacific Conference Series Vol. 407, The Low-Frequency Radio Universe. p. 121
  38. Hardcastle M. J., Evans D. A., Croston J. H., 2007, MNRAS, 376, 1849
  39. Hardcastle M. J., et al., 2016, preprint, (arXiv:1606.09437)
  40. Harwood J. J., Hardcastle M. J., Croston J. H., 2015, MNRAS, 454, 3403
  41. Helou G., Soifer B. T., Rowan-Robinson M., 1985, ApJ, 298, L7
  42. Hogg D. W., 2001, AJ, 121, 1207
  43. Ibar E., et al., 2008, MNRAS, 386, 953
  44. Ibar E., Ivison R. J., Biggs A. D., Lal D. V., Best P. N., Green D. A., 2009, MNRAS, 397, 281
  45. Ibar E., Ivison R. J., Best P. N., Coppin K., Pope A., Smail I., Dunlop J. S., 2010, MNRAS, 401, L53
  46. Intema H. T., van Weeren R. J., Röttgering H. J. A., Lal D. V., 2011, A&A, 535, A38
  47. Israel F. P., Mahoney M. J., 1990, ApJ, 352, 30
  48. Ivison R. J., et al., 2010a, MNRAS, 402, 245
  49. Ivison R. J., Magnelli B., Ibar E. e. a., 2010b, A&A, 518, L31
  50. Jannuzi B. T., Dey A., 1999, in Weymann R., Storrie-Lombardi L., Sawicki M., Brunner R., eds, Astronomical Society of the Pacific Conference Series Vol. 191, Photometric Redshifts and the Detection of High Redshift Galaxies. p. 111
  51. Jannuzi B., et al., 2010, in American Astronomical Society Meeting Abstracts #215. p. 513
  52. Jarvis M. J., Rawlings S., 2004, New A Rev., 48, 1173
  53. Jarvis M. J., et al., 2010, MNRAS, 409, 92
  54. Jenkins C. J., McEllin M., 1977, MNRAS, 180, 219
  55. Kennicutt Jr. R. C., 1998, ARA&A, 36, 189
  56. Ker L. M., Best P. N., Rigby E. E., Röttgering H. J. A., Gendre M. A., 2012, MNRAS, 420, 2644
  57. Kharb P., Srivastava S., Singh V., Gallimore J. F., Ishwara-Chandra C. H., Ananda H., 2016, MNRAS, 459, 1310
  58. Kochanek C. S., et al., 2012, ApJS, 200, 8
  59. Komatsu E., et al., 2009, ApJS, 180, 330
  60. Krause M., Alexander P., Riley J., Hopton D., 2012, MNRAS, 427, 3196
  61. Lacki B. C., 2013, MNRAS, 431, 3003
  62. Lacki B. C., Thompson T. A., Quataert E., 2010, ApJ, 717, 1
  63. Laing R. A., Peacock J. A., 1980, MNRAS, 190, 903
  64. Lane W. M., Cotton W. D., van Velzen S., Clarke T. E., Kassim N. E., Helmboldt J. F., Lazio T. J. W., Cohen A. S., 2014, MNRAS, 440, 327
  65. Magnelli B., et al., 2015, A&A, 573, A45
  66. Mahony E. K., et al., 2016, MNRAS, 463, 2997
  67. Martin C., et al., 2003, in Blades J. C., Siegmund O. H. W., eds, Proc. SPIEVol. 4854, Future EUV/UV and Visible Space Astrophysics Missions and Instrumentation.. pp 336–350, doi:10.1117/12.460034
  68. Marvil J., Owen F., Eilek J., 2015, AJ, 149, 32
  69. Mauch T., Klöckner H.-R., Rawlings S., Jarvis M., Hardcastle M. J., Obreschkow D., Saikia D. J., Thompson M. A., 2013, MNRAS, 435, 650
  70. Meisenheimer K., Roser H.-J., Hiltner P. R., Yates M. G., Longair M. S., Chini R., Perley R. A., 1989, A&A, 219, 63
  71. Mohan N., Rafferty D., 2015, PyBDSM: Python Blob Detection and Source Measurement, Astrophysics Source Code Library (ascl:1502.007)
  72. Murphy E. J., et al., 2011, The Astrophysical Journal, 737, 67
  73. Murphy E. J., Stierwalt S., Armus L., Condon J. J., Evans A. S., 2013, ApJ, 768, 2
  74. Murray S. S., et al., 2005, ApJS, 161, 1
  75. Oke J. B., Gunn J. E., 1983, ApJ, 266, 713
  76. Oliver S. J., et al., 2012, MNRAS, 424, 1614
  77. Pannella M., et al., 2015, ApJ, 807, 141
  78. Planck Collaboration et al., 2014, A&A, 571, A16
  79. Polletta M., et al., 2007, ApJ, 663, 81
  80. Rees G. A., et al., 2016, MNRAS, 455, 2731
  81. Rengelink R. B., Tang Y., de Bruyn A. G., Miley G. K., Bremer M. N., Roettgering H. J. A., Bremer M. A. R., 1997, A&AS, 124
  82. Richter G. A., 1975, Astronomische Nachrichten, 296, 65
  83. Roseboom I. G., et al., 2010, MNRAS, 409, 48
  84. Röttgering H., et al., 2011, Journal of Astrophysics and Astronomy, 32, 557
  85. Sargent M. T., et al., 2010, ApJ, 714, L190
  86. Scaife A. M. M., Heald G. H., 2012, MNRAS, 423, L30
  87. Schleicher D. R. G., Beck R., 2013, A&A, 556, A142
  88. Schmitt H. R., Calzetti D., Armus L., Giavalisco M., Heckman T. M., Kennicutt Jr. R. C., Leitherer C., Meurer G. R., 2006, ApJS, 164, 52
  89. Schober J., Schleicher D. R. G., Klessen R. S., 2016, preprint, (arXiv:1603.02693)
  90. Seymour N., et al., 2008, MNRAS, 386, 1695
  91. Shapley A. E., 2011, ARA&A, 49, 525
  92. Simpson C., et al., 2012, Monthly Notices of the Royal Astronomical Society, 421, 3060
  93. Singh V., Shastri P., Ishwara-Chandra C. H., Athreya R., 2013, A&A, 554, A85
  94. Smith D. J. B., et al., 2014, MNRAS, 445, 2232
  95. Smolcic V., 2016, preprint, (arXiv:1603.05687)
  96. Smolčić V., et al., 2008, ApJS, 177, 14
  97. Stevens J. A., et al., 2003, Nature, 425, 264
  98. Tabatabaei F. S., et al., 2013, A&A, 552, A19
  99. Tasse C., Le Borgne D., Röttgering H., Best P. N., Pierre M., Rocca-Volmerange B., 2008, A&A, 490, 879
  100. Taylor M. B., 2006, in Gabriel C., Arviset C., Ponz D., Enrique S., eds, Astronomical Society of the Pacific Conference Series Vol. 351, Astronomical Data Analysis Software and Systems XV. p. 666
  101. Vardoulaki E., et al., 2015, A&A, 574, A4
  102. Voelk H. J., 1989, A&A, 218, 67
  103. Whittam I. H., Riley J. M., Green D. A., Davies M. L., Franzen T. M. O., Rumsey C., Schammel M. P., Waldram E. M., 2016, MNRAS, 457, 1496
  104. Williams P. K. G., Bower G. C., 2010, ApJ, 710, 1462
  105. Williams W. L., Intema H. T., Röttgering H. J. A., 2013, A&A, 549, A55
  106. Williams W. L., et al., 2016, MNRAS, 460, 2385
  107. Yun M. S., Carilli C. L., 2002, ApJ, 568, 88
  108. Yun M. S., Reddy N. A., Condon J. J., 2001, ApJ, 554, 803
  109. da Cunha E., et al., 2015, ApJ, 806, 110
  110. de Jong T., Klein U., Wielebinski R., Wunderlich E., 1985, A&A, 147, L6
  111. de Vries W. H., Morganti R., Röttgering H. J. A., Vermeulen R., van Breugel W., Rengelink R., Jarvis M. J., 2002, AJ, 123, 1784
  112. van Haarlem M. P., et al., 2013, A&A, 556, A2
  113. van Weeren R. J., et al., 2016, ApJS, 223, 2
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description