The LOFAR window on star-forming galaxies and AGN - curved radio SEDs and IR-radio correlation at $0 < z < 2.5$

TL;DR

This study uses LOFAR data to analyze the radio spectral energy distributions and infrared-radio correlation of star-forming galaxies and AGN up to redshift 2.5, revealing spectral curvature differences and redshift evolution in the IR-radio ratio.

Contribution

It provides the first detailed analysis of low-frequency radio SEDs and IR-radio correlation evolution for a large, multi-wavelength selected galaxy sample across cosmic time.

Findings

SF galaxies show spectral flattening at low frequencies.

AGN SEDs become steeper at lower frequencies.

The IR-radio ratio decreases with redshift for both populations.

Abstract

We present a study of the low-frequency radio properties of star forming (SF) galaxies and active galactic nuclei (AGN) up to redshift $z=2.5$. 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 $1542$ galaxies in $\sim 7~ \rm{deg}^2$ of the LOFAR Bo\"otes field. Using the extensive multi-wavelength dataset available in Bo\"otes 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: $ q_{\rm 1.4GHz} = (2.45 \pm 0.04) \times (1+z)^{-0.15 \pm 0.03} $. Similarly, $q_{\rm 150MHz}$ shows a redshift evolution following $ q_{\rm 150GHz} = (1.72 \pm 0.04) \times (1+z)^{-0.22 \pm 0.05}$. Calibration of the 150 MHz radio luminosity as a star formation rate tracer suggests that a single power-law extrapolation from $q_{\rm 1.4GHz}$ is not an accurate approximation at all redshifts.

Figures

Tables

Table 1. 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.

	Band	$\nu \left[Hz\right],\lambda \left[m\right]$	Detections [%]
Selection	LOFAR	150 MHz	100%
	I	806 nm	100%
UV	NUV	300 nm	93%
	u	365 nm	99%
optical	BW	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 $\mu$m	100%
	IRAC2	4.5 $\mu$m	100%
	IRAC3	5.8 $\mu$m	100%
	IRAC4	8.0 $\mu$m	100%
	MIPS24	24.0 $\mu$m	100%
FIR	SPIRE	250 $\mu$m	83%
	SPIRE	350$\mu$m	79%
	SPIRE	500 $\mu$m	67%
RADIO	WSRT	1.4 GHz	67% (88%111After the flux cut applied for the curvature analysis (section 5.2.2))222Forced photometry included for non-detections
	GMRT	610 MHz	23%333Due to spatial partial coverage, see Fig. 1
	VLA-P	325 MHz	63% (88%${}^{\text{<a><span>1</span></a>}}$)${}^{\text{<a><span>2</span></a>}}$

Equations

$${\alpha }_{{\nu }_1}^{{\nu }_2}=\frac{\log (S_{{\nu }_1}/S_{{\nu }_2})}{\log ({\nu }_2/{\nu }_1)},$$

$$ln\mathcal{L}=C-\frac{1}{2}\sum \left(ln({\sigma }^2+e_i^2)+\frac{(x_i-\mu {)}^2}{{\sigma }^2+e_i^2}\right),$$

$$q_{irc}=\log \left(\frac{L_{\mathrm{IR}}/(3.75\times 10^{12}\mathrm{Hz})}{L_{\mathrm{radio}}}\right),$$

$$\mathcal{L}=\sum \left({\sigma }_{\mathrm{tot}}+\frac{(q_{\mathrm{obs}}-(C\times (1+z{)}^{\gamma }){)}^2}{{\sigma }_{\mathrm{tot}}^2}\right),$$

$${\sigma }_{\mathrm{tot}}=\sqrt{{\sigma }_{\mathrm{int}}^2+\Delta q_{obs}^2}$$

$$q_{1.4}(z)=(2.45\pm 0.04)\times (1+z{)}^{-(0.15\pm 0.03)}$$

$$R^2=\frac{{\sum }_i(q_{\mathrm{obs}-\mathrm{i}}-q_{1.4}(z_i){)}^2}{{\sum }_i(q_{\mathrm{obs}-\mathrm{i}}-⟨q_{\mathrm{obs}}⟩{)}^2}$$

$$q_{150,exp}=q_{1.4}+\log {\left(\frac{1400}{150}\right)}^{\alpha }$$

$$q_{150}(z)=(1.72\pm 0.04)\times (1+z{)}^{-(0.22\pm 0.05)}$$

$$\left(\frac{\mathrm{SF}{\mathrm{R}}_{\mathrm{IR}}}{{\mathrm{M}}_{\odot }\mathrm{y}{\mathrm{r}}^{-1}}\right)=3.88\times 10^{-44}\left(\frac{{\mathrm{L}}_{\mathrm{IR}}}{{\mathrm{erg}\mathrm{s}}^{-1}}\right).$$

$$\left(\frac{\mathrm{SF}{\mathrm{R}}_{150\mathrm{MHz}}(\mathrm{z})}{M_{\odot }y{r}^{-1}}\right)=1.455\times 10^{-24}\times 10^{{\mathrm{q}}_{150}(\mathrm{z})}\times \left(\frac{{\mathrm{L}}_{150}}{{\mathrm{erg}\mathrm{s}}^{-1}}\right),$$

Full text

††pagerange: The LOFAR window on star-forming galaxies and AGN – curved radio SEDs and IR-radio correlation at $0<z<2.5$–9††pubyear: 2002

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, $1<z<3$ . 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 ($L_{1.4\rm GHz}>10^{26}\,\rm W\,Hz^{-1}$), 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 $S(\nu)\propto\nu^{\alpha}$, where $S(\nu)$ is the radio flux density at a given frequency, $\nu$, and $\alpha$ 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 $\alpha\sim-0.7$ (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 $L_{1.4\,\rm GHz}<10^{25}\rm\,W\,Hz^{-1}$ 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 ($\nu<1\rm\,GHz$), 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 ($10^{9}<L_{\mathrm{IR,(8-1000}\mu m)}[\rm L_{\odot}]<10^{12.5}$) (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 $\Lambda$-CDM cosmology with $H_{0}=70\,\rm km\,s^{-1}\,Mpc^{-1}$, $\Omega_{\mathrm{m}}=0.3$, and $\Omega_{\Lambda}=0.7$ (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 $\rm 9~{}deg^{2}$ in the optical and near-infrared B${}_{\rm W},$ 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 deg2 with an rms noise of $\sigma_{\rm rms}\sim 120-150~{}\mu$Jy beam*-1* in the central region of the field and a resolution of $5.6\times 7.4{\rm~{}arcsec}$. 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$\sigma_{\rm rms}$ 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.6$\times$5.1 arcsec with a central noise of 0.2 mJy beam*-1* increasing to 0.8 mJy beam*-1* 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 $\sigma_{\rm rms}$. 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 deg2 of the Boötes field (project code 28_064). The mosaic has a resolution of $5\times 5$ arcsec and a noise level of $\sigma_{\rm rms}\sim 40-70~{}\mu$Jy beam*-1*. 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 (16$\times$12 hr) Westerbork Synthesis Radio Telescope (WSRT) observations of the approximately 6.68 deg2 Boötes Deep Field. The image covers 42 discrete pointings and has a limiting sensitivity of $\sigma_{\rm rms}\sim 28\,\mu$Jy beam*-1* and a resolution of 13$\times$27 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 $\sim$-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 $\mu$Jy beam*-1*) and resolution values comparable to the 365 MHz VLA-P map published by Coppejans et al. (2015) ($5.6\times 5.1{\rm~{}arcsec}$), which are the best radio data available at the respective frequencies to date.

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 $\mu$m.

This catalogue includes the original NDWFS observations in the optical and near- infrared; B${}_{\rm W},\,$ R, I and K bands. With an absolute positional uncertainty of $<0.1$ arcsec, the I-band images reach depths of 24.9 AB magnitude ($5\sigma_{\rm rms}$ 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 $\mu$m. Infrared 24 $\mu$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, Ks 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 deg2 of the field, with some gaps in the coverage due to the 10 arcmin gaps in the instrument CCDs. Y and Uspec 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-, Bw-,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 $\sim 830$k sources, including galaxies and AGN with I-band magnitudes $I<24$.

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 $\mu$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 $\mu$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 24$~{}\mu$m to optical separations of r${}_{24-opt}<3$ arcsec (less than half the FWHM of the MIPS 24$~{}\mu$m 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 deg2 of the Boötes Field providing spectroscopic redshifts for 23 745 galaxies and AGN. However, given that the majority of sources at $z>1$ 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.

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 $z_{\rm phot}$-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 $z_{\rm phot}$-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 $z_{\rm phot}$ 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 $z_{\rm phot}$ estimates. The main advantage of this method is that it determines the consensus probability P${}_{z_{\rm phot}}$ 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 deg2 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 $q_{1.4\,\rm GHz}=2.3$ and a spectral index of $\alpha=-0.7$ for the extrapolation from 1.4 GHz to 150 MHz. We choose to use a value of $q_{1.4\,\rm GHz}=2.3$ based on the observed median value of $q_{1.4}$ 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 $z<1.3$ our sample appears to be complete for SF galaxies of $SFR>300$ M*⊙*/yr (yellow line). At redshifts $1.5<z<2.3$ only galaxies with $SFR>1000~{}M_{\odot}$/yr are represented in a complete way according to the LOFAR selection. We expect to be complete for AGN of $\rm L_{150\rm MHz}>10^{25}W\,Hz^{-1}$ at all redshifts. We thus conclude that at redshifts $z>1.5$ our sample is limited to highly SF galaxies (starbursts) and radio-selected AGN.

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.

3.2.1 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 $1.01\pm 0.1$. 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 $0.91\pm 0.02$, in agreement with the data-independent method presented above.

3.2.2 Non-detections

Although the sensitivity of the LOFAR image is equivalent to the WSRT sensitivity under the assumption of a spectral index of $\alpha\sim 0.7$, 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 ($\sigma_{\rm rms}\sim 0.2~{}{\rm~{}mJy}~{}\rm beam^{-1}$) than the equivalent LOFAR image at this frequency for the same assumed $\alpha$ ($\sigma_{\rm rms}\sim 0.074~{}\rm mJy~{}beam^{-1}$). 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, $M_{*}$) and AGN (NH-torus, $L_{\rm bol}$) 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 ($\sim 79$ per cent) can be classified as 758 SF galaxies and 784 AGN.

4.1 Stellar masses and star-formation rates

Stellar masses and SFRs ($\textup{SFR}_{\rm opt}$) 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 $\tau$. The distribution of stellar masses derived for our sample has a median value of $3.01^{+5.01}_{-2.50}\times 10^{10}\rm M_{\odot}$, 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 ($\textup{SFR}_{\rm IR}$) (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

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 $L_{\rm BB}$, which is then compared to the galaxy stellar emission $L_{\rm GA}$, integrating both inside the same frequency range 0.1 $<\lambda<10\,\mu$m. In the MIR, the AGN luminosity contribution is produced mainly by the torus $L_{\rm TO}$ and is compared to the cold dust emission in starburst regions of the host galaxy $L_{\rm SB}$, as well integrating both components within the same frequency range 1 $<\lambda<40\,\mu$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 ($\log(L_{\rm TO}/L_{\rm SB})<0$), 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 $\log(L_{\rm TO}/L_{\rm SB})<0$, also satisfy $\log(L_{\rm BB}/L_{\rm GA})<0$. In the cases with $\log(L_{\rm BB}/L_{\rm GA})>0$ (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 ($\log(L_{\rm TO}/L_{\rm SB})>0$) 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 $\log(L_{\rm BB}/L_{\rm GA})>0$, where the accretion disk direct emission dominates the optical and if $\log(L_{\rm BB}/L_{\rm GA})<0$, 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 ($q\sim 2.3\pm 0.5$, 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 ($q=2.3-2.5\sigma=0.8$, 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 $\log(L_{\rm TO}/L_{\rm SB})<0$, the dominance is clearly weak ($\log(L_{\rm TO}/L_{\rm SB})\sim-0.50\pm 0.51$) compared to the values for SF galaxies ($\log(L_{\rm TO}/L_{\rm SB})\sim-1.08\pm 1.04$, 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 ($\sim 1.4\pm 1.2$), 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 $z=[0.05,1.7]$, while the redshift range $z=[1.7,2.5]$ 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 $\sim$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 $\sim 60/40$ 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 $\mu$Jy. Assuming a spectral index value of $\sim-0.7$ this turning point corresponds to $\sim 600\,\mu$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:

[TABLE]

For the total sample we calculate the most general spectral index $\alpha^{1400}_{150}$ 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, $\alpha^{1400}_{325}$ and $\alpha^{325}_{150}$ 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 ($25$ per cent of the total) the distributions of three different spectral indices are calculated for adjacent frequency pairs $\alpha^{1400}_{610},\alpha^{610}_{325}$ and $\alpha^{325}_{150}$.

5.1 Spectral index $\alpha^{1400}_{150}$ – all sources

The distribution of the spectral index $\alpha^{1400}_{150}$ 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: $\alpha^{1400}_{150}=-0.78^{+0.24}_{-0.24}$ and $\alpha^{1400}_{150}=-0.73^{+0.33}_{-0.30}$ 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 $\mu$ and $\sigma$ (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

[TABLE]

where $\sigma$ (spread of the intrinsic distribution) and $e_{i}$ (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 $\mu$ and intrinsic $\sigma$ inferred above are over-plotted on the histograms of Fig. 7.444The 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.

The distribution of the SF galaxy sample is best fitted with a mean value $\mu=-0.729\pm 0.010$ and intrinsic scatter of $\sigma=0.215\pm 0.007$, while the spectral index distribution for AGN presents $\mu=-0.664\pm 0.0112$ and $\sigma=0.294\pm 0.008$. 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 $\alpha^{1400}_{150}$ 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 $\delta\alpha^{1400}_{150}=-0.048\pm 0.011$. 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).

5.2 Spectral indices $\alpha^{325}_{150}$ and $\alpha^{1400}_{325}$

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.

5.2.1 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 $\alpha^{1400}_{150}$ presented above, a single slope spectrum around $-0.73$ and $-0.67$ for galaxies and AGN would imply observing similar values for the spectral indices $\alpha^{1400}_{325}$ and $\alpha^{325}_{150}$, 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.

5.2.2 Results on flux-selected sources with S${}_{150}>2.0{\rm~{}mJy}$

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: $\alpha^{325}_{150}=-0.63^{+0.57}_{-0.49}$ and $\alpha^{1400}_{325}=-0.74^{+0.27}_{-0.41}$. 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$<10^{-11}$). 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: $\alpha^{325}_{150}=-0.80^{+0.55}_{-0.69}$ and $\alpha^{1400}_{325}=-0.56^{+0.42}_{-0.35}$. 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 $\Delta\alpha\sim-0.35\pm 0.05$. 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$<10^{-11}$). 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 $\alpha^{325}_{150}$, $\alpha^{610}_{325}$ and $\alpha^{1400}_{610}$

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 $S_{150}>2$ 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 $\alpha^{1400}_{150}=-0.80\pm 0.27$, and for AGN of $\alpha^{1400}_{150}=-0.70\pm 0.31$. Adding one more data point to the radio SEDs, we calculate the three spectral indices $\alpha^{325}_{150}$, $\alpha^{610}_{325}$ and $\alpha^{1400}_{610}$ 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.

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 ($\alpha^{1400}_{610}=-0.82,\alpha^{610}_{325}=-0.89$ and $\alpha^{325}_{150}=-0.65$) suggest a slight tendency to have flatter spectra towards low frequencies, the distributions’ medians for AGN ($\alpha^{1400}_{610}=-0.59,\alpha^{610}_{325}=-0.63$ and $\alpha^{325}_{150}=-0.67$) 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, $\alpha^{1400}_{610}$ and $\alpha^{325}_{150}$, 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 $\alpha^{1400}_{150}$ and the curvature parameter $\alpha^{325}_{150}-\alpha^{1400}_{325}$, respectively. For this study we use the samples with detections in both VLAP and WSRT ($\sim$90 per cent of the sample with S${}_{150}>2.0{\rm~{}mJy}$ ), 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 $z\sim 0$) and high redshift (450 MHz-4.2 GHz at $z=2$) 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 $\alpha^{1400}_{150}=-0.7$. 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= $-0.05\pm 0.03$). 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).

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 $\alpha^{\rm 1.4GHz}_{\rm 610MHz}$ with redshift across $0<z<3$ for SF galaxies.

Fig. 11 shows the redshift evolution of the spectral curvature, parametrized as the difference between spectral indices $\alpha^{325}_{150}$ and $\alpha^{1400}_{325}$, 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 $\alpha^{325}_{150}-\alpha^{1400}_{325}$, 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 ($-0.055\pm 0.103$), while the slope for AGN is found to be $0.147\pm 0.096$. We test these results for significance by calculating the coefficient of determination statistics $R^{2}$, 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 $R^{2}$, 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 $\textup{SFR}_{\rm IR}$. 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 $L_{\rm TO}$, 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 $0.168\pm 0.019$, 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.

5.6.1 SF galaxies

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 $\Delta\alpha=0.2$ per logarithmic frequency decade towards lower frequencies, which is a stronger flattening than that found in this study ($\Delta\alpha=0.1$). One possible explanation is that Marvil et al. (2015) include a large fraction of non detections at the central wavelengths ($>60$ 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 $<1$ 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.

5.6.2 AGN

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 ($\alpha$ 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 $z>0.7$ and frequencies $\nu>325{\rm~{}MHz}$, 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 $\sim$810 galaxies.

The IRC is parametrized by the widely used value $q_{\rm irc}$, which is defined as the ratio

[TABLE]

where $L_{\rm IR}$ is the total rest-frame infrared luminosity (in $\rm~{}erg~{}s^{-1}$) of the cold dust SED template integrated along the rest-frame frequency range $8-1000\,\mu$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). Lradio is the luminosity at the radio frequency to be studied, in this case 1.4 GHz or 150 MHz in $\rm erg~{}s^{-1}~{}Hz^{-1}$. The factor $3.75\times 10^{12}\rm~{}Hz$ is the frequency corresponding to $80~{}\mu$m, used in the definition to make $q_{irc}$ a dimensionless quantity.

Radio rest-frame luminosities are computed using k-corrections inferred assuming $S_{\nu}\propto\nu^{\alpha}$, where $\alpha$ is the mean value observed for the SF galaxies distribution $\alpha=-0.73$ (see Fig. 7). We use a fixed $\alpha$ 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 qirc.

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 $q_{1.4}$ is scattered in a distribution with a median value of $q_{1.4}=2.28^{+0.017}_{-0.021}$, 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 $q_{1.4}$ 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: $q_{1.4}=2.34\pm 0.26$ Yun et al. (2001), $q_{1.4}=2.3$ Helou et al. (1985), using total IR emission as in our case: $q_{1.4}=2.64\pm 0.26$ 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 $q_{1.4}$ (1.4 GHz) as a function of redshift. To quantify a redshift evolution we fit the function: $q_{1.4}=C\times(1+z)^{\gamma}$, where $C$ is a constant and $\gamma$ the variable driving the evolution. We include only 86 per cent of sources, which lie within the $2\sigma$ 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

[TABLE]

where

[TABLE]

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

[TABLE]

with an intrinsic scatter of $\sigma_{\rm int}=0.53$.

As expected, after the fit for redshift evolution the value $q_{1.4}$ at $z=0$ is higher than the median value of the total sample by $\delta=0.11$. To infer the proportion of the variance in the distribution of Fig. 13, which is explained by the $z$-evolution, we infer the coefficient of determination $R^{2}$, defined in this case as,

[TABLE]

where qobs-i are the observed q-values and q${}_{\rm 1.4}(z)$ are values calculated following Eq. 6. We calculate that the redshift evolution accounts for $\sim 6$ per cent of the variance in the total distribution, while most of the intrinsic scatter is still not explicable.

A dependence of the q${}_{\rm 1.4}-z$ evolution on the $\alpha-z$ evolution around 1.4 GHz is unlikely, since this would imply a clear evolution of $\alpha_{\rm tot}$ 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 $q_{1.4}$ 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 $q_{1.4}$ is related to a changing ratio of AGN and SF contribution to the IR emission ($L_{\rm TO}/L_{\rm SB}$). 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 $q_{1.4}$ of our sample ($q_{1.4}\sim 2.28\pm 0.02$) and the local mean IRC found by (Bell, 2003, $q_{1.4}\sim 2.64\pm 0.26$)). 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 $q_{1.4}$ 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 $q_{1.4}$ 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 $z<2$ is negligible. Since $>90$ per cent of the sources in our sample have redshifts $z<2$, the possible bias against faint radio sources is not significant.

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

The picture of decreasing $q_{1.4}$ 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 $q_{1.4}$ (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 $q_{1.4}\propto(1+z)^{-0.12}$ and $q_{1.4}\propto(1+z)^{-0.15}$ respectively. Similarly, Delhaize et al. (2017) find a redshift evolution of $q_{1.4}\propto(1+z)^{-0.19}$ 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 $q_{150}$ in the lower panel of Fig. 13. The distribution of $q_{150}$ presents a median value and error on the median of $q_{150}=1.544^{+0.015}_{-0.023}$.

It is customary to calculate the value of $q_{irc}$ for different radio frequencies interpolated from the well calibrated value $q_{1.4}$ 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

[TABLE]

Assuming a unique power law value from 1400 to 150 MHz of $\alpha=-0.73$, we obtain an expected q-value of $q_{150,exp}=1.59$. This value is well in agreement with the distribution observed in the lower panel Fig. 13 ($q_{150,obs}=1.54$) given the errors.

Also here we consider sources with $q$-values within the 2$\sigma$ region of the distribution in Fig. 13 and investigate the evolution of $q_{150}$ with redshift. We find the relation

[TABLE]

with an intrinsic scatter of $\sigma_{\rm int}=0.53$. Proceeding with a similar analysis as for $q_{1.4}$, we calculate the $R^{2}$ values for this regression and obtain that the redshift evolution accounts for $\sim 8.4$ 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 $q_{150}$ is that its redshift dependence, $q_{150}\propto(1+z)^{-0.22\pm 0.05}$, is slightly stronger ($\delta\gamma=0.07$) than that observed for $q_{1.4}$, with $q_{1.4}\propto(1+z)^{-0.15\pm 0.04}$. We explore the connection between the stronger redshift dependence of $q_{150}$ and our results on the spectral curvature described in section 5.4.

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 $q_{1.4}\propto(1+z)^{-0.22\pm 0.04}$, which is a similar result to the one obtained from $q_{150}$. The study comprising exclusively sources undetected in 1.4 GHz returns an evolution of $q_{1.4}\propto(1+z)^{-0.07\pm 0.04}$. From these results we conclude that the weaker redshift evolution of $q_{1.4}$ compared to $q_{150}$ 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 $q_{1.4}$ 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 ($8-1000\,\mu\rm m$, see section 4 for details ) to compute the SFRs throughout the paper following:

[TABLE]

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:

[TABLE]

where the factor $(1.455\times 10^{-24}\times 10^{\rm q_{150}(z)})$ is equivalent to $2.450\times 10^{-24}$, $2.289\times 10^{-24}$, $2.149\times 10^{-24}$ and $1.965\times 10^{-24}$ for redshifts $z$ around 0.1, 0.5, 1.0 and 2.0, respectively.

In Fig. 16 we show how the redshift evolution of the $q_{150}$-value suggests a more accurate SFR estimation compared to assuming a unique $q_{150}$ value. We over-plot the observed median SFRIR 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 $q_{150}$-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 $\sim 1500$ galaxies at redshifts of $0.05<z<2.5$. 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 $\alpha^{1400}_{150}=-0.73$. In comparison, AGN were found to exhibit a flatter spectral slope, $\alpha^{1400}_{150}=-0.67$. Although the difference is small compared to the width of the observed $\alpha^{1400}_{150}$ 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 $q_{1.4}(z)=2.45\pm 0.04\times(1+z)^{-0.15\pm 0.03}$, where $q_{1.4}$ is the ratio between the total infrared luminosity ($L_{\rm{IR}},8-1000\,\mu m$) and the 1.4 GHz radio luminosity ($L_{\textup{1.4 GHz}}$). The observed evolution accounts only for $\sim 6$ per cent of the variance of the distribution. The large observed intrinsic scatter of $\sigma_{int}=0.53$ 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 $q_{150}(z)=1.72\pm 0.04\times(1+z)^{-0.22\pm 0.05}$. A comparison of the observed $q_{150}$ values and those extrapolated from $q_{1.4}$, as a function of redshift, give similar results but with some scatter around $z\sim 1$.

• •

After converting $L_{\rm IR}$ to $\rm SFR_{\rm IR}$ 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.

Appendix B SED-fitting examples

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