LOFAR-Bo\"otes: Properties of high- and low-excitation radio galaxies at $0.5

This paper presents a study of the redshift evolution of radio-loud active galactic nuclei (AGN) as a function of the properties of their galaxy hosts in the Bo\"otes field. To achieve this we match low-frequency radio sources from deep $150$-MHz LOFAR observations to an $I$-band-selected catalogue of galaxies, for which we have derived photometric redshifts, stellar masses and rest-frame colours. We present spectral energy distribution (SED) fitting to determine the mid-infrared AGN contribution for the radio sources and use this information to classify them as High- versus Low-Excitation Radio Galaxies (HERGs and LERGs) or Star-Forming galaxies. Based on these classifications we construct luminosity functions for the separate redshift ranges going out to $z = 2$. From the matched radio-optical catalogues, we select a sub-sample of $624$ high power ($P_{150\mathrm{\,MHz}}>10^{25}$ W Hz$^{-1}$) radio sources between $0.5 \leq z<2$. For this sample, we study the fraction of galaxies hosting HERGs and LERGs as a function of stellar mass and host galaxy colour. The fraction of HERGs increases with redshift, as does the fraction of sources in galaxies with lower stellar masses. We find that the fraction of galaxies that host LERGs is a strong function of stellar mass as it is in the local Universe. This, combined with the strong negative evolution of the LERG luminosity functions over this redshift range, is consistent with LERGs being fuelled by hot gas in quiescent galaxies.


INTRODUCTION
The evolution of radio-loud Active Galactic Nuclei (RL AGN) is closely entwined with that of their host galaxies and the central supermassive black holes that power them. The ability of the expanding radio lobes of RL AGN to do work on the surrounding intra-cluster medium provides a important 'feedback' mechanism by which a central black hole can regulate or extinguish star forma-E-mail: w.williams5@herts.ac.uk (WLW) tion in its parent galaxy (see e.g. Best et al. 2006Best et al. , 2007Bower et al. 2006;Croton et al. 2006;Fabian et al. 2006;Cattaneo et al. 2009). Over recent years RL AGN have come to be classified based on their Eddington-scaled accretion rates, with sources on either end of the scale exhibiting very different charactersitics (Best & Heckman 2012;Son et al. 2012;Russell et al. 2013;Mingo et al. 2014;Gürkan et al. 2014;Fernandes et al. 2015).
RL AGN with high Eddington-scaled accretion rates experience radiatively efficient accretion of cold gas via an accretion disc (e.g. Shakura & Sunyaev 1973) and therefore appear as 'quasars' (Silk & Rees 1998), with emission across the electromagnetic spectrum (e.g. Barthel 1989;Antonucci 1993;Urry & Padovani 1995). In the literature these are variously refered to as 'cold mode' or 'radiative mode' or 'high-excitation' sources because they are characterised by strong optical emission lines. They are typically hosted by lower mass, bluer galaxies in less dense environments (e.g. Tasse et al. 2008b;Janssen et al. 2012). While the most powerful radio sources tend to be high excitation radio galaxies (HERGs), they are in fact found at all radio powers (Best & Heckman 2012). Due to their strong evolution with redshift this mode is likely important in cutting off star formation at high redshifts and thus setting up the tight black hole vs bulge mass relation that is observed locally (Magorrian et al. 1998). At the low, or radiatively inefficient, end of the Eddington-scaled accretion rate spectrum radio galaxies are found to have no or weak emission lines (Hine & Longair 1979;Laing et al. 1994;Jackson & Rawlings 1997) and are thought to be fuelled by hot gas accreting directly onto the supermassive black hole (Hardcastle et al. 2007), e.g. via advection dominated accretion flows (ADAFs, e.g. Narayan & Yi 1995). Typically hosted by higher mass, redder galaxies and occurring in more dense environments (Best et al. 2005a), these sources have no mid-infrared emission or optical obscuration from dust (Whysong & Antonucci 2004;Ogle et al. 2006), they have no accretion-related X-ray emission Evans et al. 2006) and their radio powers tend to be low. Forming the bulk of the population in the local Universe, these low excitation radio galaxies (LERGs) are otherwise refered to as 'hot mode', 'radio mode' or 'jet mode' in the literature. LERGs have a direct link between the black hole and its hot gas fuel supply and can maintain elliptical galaxies at lower redshifts as 'old, red and dead' (e.g. Best et al. 2006) and can prevent strong cooling flows in galaxy clusters (e.g. Fabian et al. 2006). For a comprehensive review on the current understanding of the HERG/LERG populations see Heckman & Best (2014) and McNamara & Nulsen (2012) and references therein.
It is well known that, within the local universe (z < ∼ 0.3), the RL fraction, i.e. the fraction of galaxies hosting a RL AGN, is strongly dependent on the stellar mass of the host galaxies ( f RL ∝ M 2.5 * , Best et al. 2005b;Tasse et al. 2008b;Janssen et al. 2012;Simpson et al. 2013), increasing to > 30 per cent at stellar masses above 5 × 10 11 M for radio luminosities > 10 23 W Hz −1 . However this mass-dependence of the entire population is driven by that of LERGs which dominate the RL AGN population at these redshifts (Best et al. 2006). The RL fraction for HERGs has a much shallower mass-dependence, f RL ∝ M 1.5 * (Janssen et al. 2012). Furthermore, Janssen et al. (2012) have shown that the fraction of RL AGN for the two classes have different dependencies not only on the stellar mass of the host galaxies, but also on properties such as colour and star formation rate (SFR): red (passive) galaxies are a factor of a few times more likely to host lower power LERGs than blue (star-forming) galaxies of the same stellar mass; blue galaxies show a higher probability of hosting HERGs at all radio luminosities than red galaxies; and for blue galaxies, the likelihood of hosting either radio AGN type is a strong positive function of the SFR. It is clear that the presence of cold, star-forming gas in a galaxy clearly enhances the probability of its central BH becoming a RL AGN. This means that some LERG activity, especially at high radio luminosities, is not solely related to hot halo gas accretion and is consistent with it being produced at low accretion rates by either hot or cold gas (Heckman & Best 2014). A key open question is how the radio galaxy populations and RL fraction for each depends on host galaxy masses and colours at higher redshifts. As a first step, in studying the RL fraction at z ≈ 1 − 2, Williams & Röttgering (2015) found more than an order of magnitude increase in the fraction of lower mass galaxies (M * < 10 10.75 M ) which host RL AGN with radio powers P 1.4 GHz > 10 24 W Hz −1 compared to the local Universe.
Optical spectra are the key discriminator between HERGs and LERGs. Based on SDSS sepctroscopy, Best & Heckman (2012) built the largest sample of HERGs/LERGs in the local Universe, but this is harder to do at higher redshifts. Best et al. (2014, hereafter B14) provided the first sample of intermediate redshift (z < 1) objects that are spectroscopically classified as HERGs and LERGs. Since then, Pracy et al. (2016) have classified a much larger sample, but still probing only out to a redshift of about one. To build large high-redshift samples requires a method independent from spectroscopy for the separation of HERGs and LERGs. Quasarselection techniques based on Mid-infrared (MIR) colours (Stern et al. 2005;Donley et al. 2012;Stern et al. 2012) fail to select all high excitation sources and selections based on X-ray emission alone (e.g. Hickox et al. 2009) miss obscured and weaker sources. In this paper we classify a sample of RL AGN as HERGs and LERGs on the basis of their broad-band spectral energy distributions (SEDs), and study the RL fractions, radio luminosity functions and colour and mass dependencies for the two classes of RL AGN at intermediate redshifts of 0.5 ≤ z < 2. Preliminary results were presented in . This paper is structured as follows: the LOFAR 150-MHz radio data is described in Section 2 and the multi-wavelength datasets and catalogues we use are described in Section 3. In Section 4 we use SED fitting to determine photometric redshifts and galaxy parameters for the sample of optical galaxies. Section 5 describes our method for identifying optical counterparts to the LOFAR radio sources. In Section 6 we describe further SED fitting to classify sources from this RL AGN sample as HERGs and LERGs. Section 7 describes the selection of a well-defined subsample of RL AGN and presents an analysis of the properties of the RL AGN, including the RL fraction and luminosity functions of HERGs and LERGs. Throughout this paper we use AB magnitudes and a concordance cosmology with Ω M = 0.3, Ω Λ = 0.7, and H 0 = 70 km s −1 Mpc −1 . The spectral index, α, is defined as S ν ∝ ν α , where S is the source flux density and ν is the observing frequency. We assume a spectral index of −0.7 unless otherwise stated.

RADIO DATA
The low-frequency radio data are described by Williams et al. (2016), but we provide a brief summary here. The 8 hr observation was taken with the LOw Frequency ARray (LOFAR; van Haarlem et al. 2013) using the High Band Antennae (HBA) and covering the frequency range 130-169 MHz, with a central frequency of ≈ 150 MHz. Particular care was taken in the calibration and imaging to correct for direction-dependent effects (DDEs) caused by the ionosphere and imperfect knowledge of the LOFAR station beam shapes. This DDE calibration and imaging was achieved with the 'Facet' calibration scheme presented by van Weeren et al. (2016). The resulting image covers 19 deg 2 , with an rms noise of ≈ 120 − 150 µJy beam −1 . Assuming a spectral index of −0.7, the sensitivity of this map is comparable to the 28 µJy beam −1 rms of the WSRT 1.4 GHz-image made by de Vries et al. (2002). However, LOFAR's superior resolution of 5.6 × 7.4 arcsec (compared to 13 × 27 arcsec at 1.4 GHz), combined with its positional accuracy of < 1 arcsec, makes it significantly better for the optical identi- Figure 1. Coverage diagram for the Boötes field. The black circle shows the LOFAR 150-MHz coverage. The blue polygon shows the main I-selected psf-matched catalogue region, which is covered completely by both the ND-WFS (B W RIK) and SDWFS (3.6, 4.5, 5.8, and 8.0 µm). It covers a total of 9.2 deg 2 , when regions contaminated by bright stars are excluded. The red squares show the zBoötes coverage, which has some gaps. The orange circles show the GALEX NUV coverage. There is a small area not covered by the NEWFIRM survey (J, H, and K s , shown in cyan) and the LBT/LBC survey (U spec , and Y , shown in light green). fication of the radio sources. The LOFAR 150-MHz radio source catalogue contains 6 276 sources detected with a peak flux density threshold of 5σ , where σ is the local rms noise. The radio coverage is shown as a circle in Fig. 1.

MULTI-WAVELENGTH DATA
The Boötes field is among the widest of the famous deep extragalactic fields and was first observed as one of two fields within the National Optical Astronomy Observatory (NOAO) Deep Wide Field Survey (NDWFS; Jannuzi et al. 1999). Since then it has been surveyed across the electromagnetic spectrum. We describe here the surveys and datasets that are used in this work.

Combined Photometry Catalogue
The primary catalogue that we make use of is the combined I-bandselected psf-matched photometry catalogue presented by Brown et al. (2007Brown et al. ( , 2008. This catalogue includes 15 bands spanning 0.14-24 µm and combines several different surveys. These include the original optical (B W , R, and I) and the NIR (K) survey, surveys with the Spitzer Space Telescope at 3.6, 4.5, 5.8, and 8.0 µm (SDWFS; Ashby et al. 2009), and 24 µm (MAGES; Jannuzi et al. 2010), NUV and FUV surveys from GALEX, and deeper J, H, K s , and z band surveys. Brown et al. (2007) have constructed a combined psf-matched catalogue by regridding and smoothing the individual released sur-vey images to a common scale so that the stellar point-spread function (PSF) is a Moffat profile with a full width at half-maximum (FWHM) of 1.35 arcsec and β = 2.5 for the B w -, R-, I-, Y -, H, K-, and K s -bands and with a FWHM of 1.6 arcsec for the u-, zand J-bands. PSF fluxes are extracted from these images for all the sources in the I-band using SExtractor (Bertin & Arnouts 1996). For the remaining bands, aperture fluxes were extracted. Regions surrounding very extended galaxies and saturated stars were excluded. The final sample area is 9.2 deg 2 . The geometry of the Boötes field is shown in Fig. 1.

Additional Multi-wavelength Coverage
Boötes is part of the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012), which includes photometry using the Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010) instrument at 250µm, 350µm, and 500µm. Within HerMES, Boötes has 'level 5' coverage of 3.25 deg 2 to 5σ noise levels of 13.8, 11.3, and 16.4 mJy and 'level 6' coverage of 10.57 deg 2 to 5σ noise levels of 25.8, 21.2, and 30.8 mJy. In this paper, we use the maps (Levenson et al. 2010) from the fourth data release (DR4).
The AGN and Galaxy Evolution Survey (AGES; Kochanek et al. 2012) has provided redshifts for 23 745 galaxies and AGN across 7.7 deg 2 of the Boötes field. The AGES spectra were obtained for random sparse samples of normal galaxies brighter than m I < 20 mag (significantly deeper than SDSS). Additional samples of AGN, selected in the radio, X-ray, IRAC mid-IR, and MIPS 24 µm, were targeted to fainter limiting magnitudes (m I < 22.5 mag for point sources). The survey used the Hectospec instrument (Fabricant et al. 2005) on the MMT to obtain 3700-9200Å spectroscopy at a spectral resolution of 6Å (R ≈ 1000 Kochanek et al. 2012;Cool et al. 2012). The median redshift of the galaxies in the survey is z = 0.3, with 90 per cent of the redshifts in the range 0.085 z 0.66. However, the spectroscopic redshift completeness for the matched LOFAR sources believed to be at z > 1 is less than 50 per cent. For this reason we derive photometric redshifts, described in the following section. AGES also provides photometric redshifts, calculated using the LRT code by Assef et al. (2010) that fits a combination of an early-type, late-type, star forming, and (obscured) AGN to the observed broadband SEDs. The photometry they used is a subset of that used in this work.

SED FITTING
For the 888, 956 optical sources in the Brown et al. (2007) psf-matched photometry catalogue with m I ≤ 24 mag and FLAG DEEP = 1, and for which we have either spectroscopic or photometric redshifts, we fit their spectral energy distributions (SED) to determine galaxy parameters, including stellar mass, star formation rates and colours. Prior to any fitting, the photometry catalogue was filtered to remove catastrophic outliers, i.e. flux densities lower (higher) than 2.5 mag (1 mag) compared to the two adjacent filters were flagged (and not used in later fitting). These cutoffs were chosen to be sufficiently extreme not to flag any reasonable spectral emission or absorption features and by comparing to two adjacent filters bona fide spectral breaks are not flagged. About 1 − 2 per cent of the photometry points were flagged in this way.

Photometric Redshifts
Photometric redshifts are provided by the hybrid photometric redshift method presented by Duncan et al. (in prep) and , based on the Brown et al. (2007) photometry catalogue. The redshifts are derived by combining template-based estimates with additional Gaussian process estimates (Almosallam et al. 2016a,b) trained for subsets of the sample population, specifically infrared-, X-ray-, and optically-selected AGN as well as the remaining galaxy population. The three different template-based estimations were calculated following the methodology presented by , using the EAZY software (Brammer et al. 2008) and three different template sets: one set of stellar-only templates, the EAZY default library (Brammer et al. 2008), and two sets including both stellar and AGN/QSO contributions, the XMM-COSMOS templates (Salvato et al. 2008) and the Atlas of Galaxy SEDs (Brown et al. 2014).
The multiple individual z phot estimates were then combined using a Hierarchical Bayesian method (Dahlen et al. 2013), as an alternative to a straight addition of the probability distributions of the z phot estimates. The main advantage of this method is that it determines the consensus probability P(z 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 and the posterior redshift predictions calibrated such that they accurately represent the uncertainties in the estimates.

Comparison with AGES Redshifts
While the quality of the photometric redshifts is analysed in detail by Duncan et al. (in prep), we provide here a brief overview because the quality of the photometric redshifts is fundamental for the subsequent analysis. In Fig. 2 we show a comparison between the Duncan et al. (in prep) z phot and z spec for the sources with good AGES spectroscopic redshifts (with a signal-to-noise > 5). In general, the photometric redshifts compare well to the spectroscopic redshifts, although we note that this comparison is primarily from galaxies at z spec < 1.0. Galaxies that are > 3σ outliers from the one-to-one relation based on their redshift errors from the consensus z phot estimates are considered catastrophic outliers, the fraction of which is 1.2 per cent. As a measure of the accuracy of the photometric redshifts, we consider two quantities, computed after excluding the catastrophic outliers. The first goodness measure is the standard dispersion, σ z /(1 + z), defined by The second is the normalized median absolute deviation, or NMAD, of the residuals, defined as NMAD(∆z) = 1.48 × Median(∆z), where ∆z = (z phot − z spec )/(1 + z spec ). We measure σ z /(1 +z) = 0.11 and NMAD = 0.028. It is well known that photometric redshifts are poorly determined for AGN (e.g. Brodwin et al. 2006;Rowan-Robinson et al. 2008;Assef et al. 2010), and should preferably be fit using different methods (e.g. Salvato et al. 2009Salvato et al. , 2011. We compare the z phot and z spec for normal galaxies and AGN separately in Fig. 2. For this we use the sources flagged as AGN by Assef et al. (2010), which is based on their having a significant contribution by an AGN SED template. Excluding the galaxies selected as AGN in AGES, we find that the photometric redshifts are more accuarate for normal galaxies, with σ z /(1 + z) = 0.045 and NMAD = 0.030. Considering only the AGES AGN, we find σ z /(1+z) = 0.17 and NMAD = 0.065. This is comparable with the redshifts determined by Assef et al. (2010), with σ z /(1 + z) = 0.04 for normal galaxies and σ z /(1 + z) = 0.18 for point-like AGN. For comparison, the most accurate photometric redshifts available in the literature typically have σ z/(1 + z) 0.01 (e.g. Ilbert et al. 2009;Muzzin et al. 2013), but using 30 bands of broad, intermediate and narrow width.

Stellar Masses, Star Formation Rates and Rest-frame Colours
Stellar population parameters are determined by fitting galaxy SEDs using the FAST code (Kriek et al. 2009), based on the Bruzual & Charlot (2003) models. We assume solar metallicity, a Chabrier (2003) initial mass function (IMF), and a Calzetti et al. (2000) dust extinction law. The template SEDs are constructed in the standard way (see e.g. Muzzin et al. 2013), assuming exponentially declining star formation histories (SFHs) of the form SFR ∝ exp(t/τ), where t is the time since the onset of star formation and τ is the e-folding star formation timescale in units of Gyr. All galaxies are fitted assuming their redshift is the z spec from AGES or, where none is available, the consensus z phot estimate. In all, four parameters are determined per galaxy: τ , t, A V (the V band extinction), and a normalization. The stellar mass (M * ) is then determined from mass-to-light ratio of the best-fit SED multiplied by the best-fit normalization of the SED.
Rest-frame colours are derived using INTERREST (Taylor et al. 2009) with the consensus z phot estimates. We determine colours for the 0.1 u and 0.1 r bands, defined as the AB magnitudes in the SDSS u and r bands at z = 0.1. These colours allow straightforward comparison to SDSS results (e.g. Blanton et al. 2003c,b;Kauffmann et al. 2003;Blanton et al. 2003a).

OPTICAL IDENTIFICATION OF RADIO SOURCES
In this section we describe the identification of optical counterparts, from the I-band-selected optical catalogue described in Section 3, matched to the LOFAR radio sources, described in Section 2. We use a statistical technique to determine the probability that an Iband optical source is the true host of a particular radio source. Prior to this, we inspect the radio-optical images (radio contours of each radio source overlaid upon the corresponding I-band image) and classify their radio morphologies into the different categories described below.

Visual Classification
In order to identify the host galaxies of radio sources, the true location of the host galaxy with respect to the radio source should be known. Following Best et al. (2003) and Tasse et al. (2008a) we determine a strong subjective prior on this location for each source by visually inspecting all the radio-optical images and dividing them into the following classes based on the radio morphology: Class 1: For these sources the radio emission is expected to be coincident with the optical emission (although the optical emission may be below the detection limit). This occurs in sources such as starburst galaxies, compact core-dominated radio sources or radio sources where the radio core can be clearly identified. In these cases, the errors on the radio and optical positions can be used in a statistical way to identify the optical counterpart of each radio source. We consider all relatively small (usually single-component) Figure 2. Photometric redshifts from the Boötes I-selected catalogue vs. spectroscopic redshifts from the AGES catalogue. Left Only galaxies not indicated as AGN in the AGES catalogue are plotted. Right Sources indicated to be AGN by the AGES SED fitting. The solid and dotted blue curves show respectively the median and rms dispersion of δ z = (z phot − z spec )/(1 + z spec ), within 11 logarithmic-spaced bins across the spectroscopic redshift range. sources in this category, even if they are slightly resolved and take into account the larger uncertainties on the radio positions in the likelihood ratio analysis in the next section. We note that some radio sources appear resolved because of some bandwidth-and timesmearing in the LOFAR image (see Williams et al. 2016). Class 2: In the case that no radio core is identified (such as for classical double lobe FRII (Fanaroff & Riley 1974) radio sources, only a weak prior can be considered for the optical host position. The position of the host and associated errors are estimated based on the flux-weighted centroid of the multiple Gaussian fitting components, as described in more detail by Best et al. (2003). For very large such sources the error regions become large and these are then considered as Class 3 sources below. Class 3: When the radio source is large or very asymmetric, the flux weighted radio centroid and associated errors can be very far from the real optical host. We use the combination of radio morphology and optical properties (such as an elongated lobe pointing to a bright optical object), to infer the position of the optical counterpart. These sources are matched (or left without an optical match where none is obvious) visually on a case-to-case basis and the statistical method described below cannot be used. Class 4: These are clearly resolved and diffuse radio sources whose morphology is not suggestive of jets. This includes 'radio halos' and 'relics', typically found in clusters. These sources have been excluded from further analysis. Class 5: When the radio source overlaps a bright saturated source, we have classified the source as Class 5. These sources likely have contaminated photometry and have been excluded in further analysis.

Likelihood Ratio
For the Class 1 and Class 2 sources we employ a statistical method to determine the optical counterparts to the radio source. We use the likelihood ratio (LR) method (Richter 1975) to determine the probability that an I-band optical source is the true counterpart of a particular radio source. The LR method has been further developed by Prestage & Peacock (1983); Benn (1983); Wolstencroft et al. (1986) and Sutherland & Saunders (1992). Here we use the methodology outlined by Tasse et al. (2008a). The probability that an optical I-band source is the true optical counterpart of a given radio source is determined from the LR (Sutherland & Saunders 1992;Tasse et al. 2008a), defined as: where m is the I-band magnitude of the optical candidate, θ (< m) is the pre-determined probability that a radio source has an observed optical counterpart with magnitude < m, and ρ(< m) is the surface number density of objects with magnitude < m. The parameter r is the uncertainty-normalised angular distance between the radio core and the optical host candidate, defined as r 2 = (∆α/σ α ) 2 + (∆δ /σ δ ) 2 , where ∆ is the positional difference, σ is the uncertainty, and α and δ are the right ascension and declination respectively. For each α and δ , the uncertainty is the quadratic sum of the uncertainty on the radio position, σ radio , and on the optical position, σ opt . We adopt an optical astrometry accuracy of σ opt ≈ 0.35 arcsec, independent of the magnitude m I . The accuracy of the radio position, σ radio , is different for every source and depends on the signal-to-noise ratio and the Gaussian fitting parameters (Williams et al. 2016). The probability P id (i) of the i-th candidate being a true identification is: where θ (m lim ) is the fraction of radio sources having detected optical counterparts at the limiting magnitude of the survey, i refers to the candidate under consideration and j runs over the set of all possible candidates. We estimate the association probability assuming that θ and ρ depend only on the object magnitude m, which is taken as the I-band magnitude of the optical candidate. For each radio source we calculate the density function ρ(m) within 2 arcmin of the radio source centroid, in order to account for the variation of the surface density with position, or clustering of optical sources.
To estimate the function θ (< m), we follow the methodology of Tasse et al. (2008a). This involves simulating random radio and optical catalogues with a known fraction of radio-optical matches and comparing the simulated radio-optical separation distribution to the real distribution. We consider discrete I-band magnitude cuts in the interval 13 < I < 24 with an increment ∆m I = 0.2. For each of these cuts an optical catalogue with uniformly-distributed positions is generated. A corresponding radio catalogue is generated assuming a given fraction, θ (< m), of radio sources have an optical counterpart (i.e. the same position as a source in the optical catalogue), while the remainder have uniformly-distributed positions. All the radio and optical positions are then shifted by Gaussian-distributed offsets in right ascension and declination with the standard deviations given by the respective positional uncertainties. The distribution of the angular distance between radio sources and their closest object in the optical catalogue is then computed and compared to the real distribution through a Kolmogorov-Smirnov test. The fraction, θ (< m), corresponding to the maximum Kolmogorov-Smirnov probability is retained. For each I-band magnitude cut, the test is repeated 10 times, to estimate an error on θ (< m).

Radio-Optical Match Results
Of the 6 267 sources in the LOFAR 150-MHz catalogue, 3 894 lie within the boundary of the optical catalogue and may therefore have potential optical counterparts. Based on the visual inspection of the radio and optical images, we separated eight sources that had been grouped by the original source detection algorithm (i.e. where PYBDSF grouped two Gaussians into one source, but the optical images suggest these are two Class 1 sources instead of one Class 2 source). The majority, 3, 403, of these sources were classified as Class 1 (87 per cent), 177 sources (4.5 per cent) were classified as Class 2, 4 sources (< 1 per cent) were classified as Class 3, and 24 sources (< 1 per cent) were classified as Class 4 (diffuse sources) or Class 5 (sources with bad optical photometry). Some examples of the Class 1 and 2 sources with LR-matched optical sources are shown in Fig. A1 and Fig. A2 respectively in Appendix A. In the following analysis for the Class 1 and 2 sources with LR-matched optical sources we select only the match with the maximum probability, where there is more than one possible optical identification, and only the sources with P id > 0.7. Of the total of 3 902 sources, we found at least one optical counterpart for 2 428 sources (76 per cent) of which 2, 835 have m I < 24 mag. Fig. 3 shows the redshift distribution of all the matched radiooptical sources. A small number (30) of sources, not shown, have photometric redshifts in the range 3 < z < 6 are not shown here. We show also the predictions from the SKA simulated skies (Wilman et al. 2008) constrained by our observed coverage area and rms distribution. The distributions of the simulated sources are in very good agreement with the observed distributions, at least up to z < 2. Above this redshift there are indications of incompleteness in our matched sample as sources fall below our optical detection threshold. The dotted lines in this plot show the distribution of sources with spectroscopic redshifts from AGES -the low completeness of which motivates the need for a complete sample with photometric redshifts. However, due to the AGES selection criteria (see Kochanek et al. 2012), most of the sources at z > 1 with spectroscopic redshifts are AGN. Thus, the sources most likely to have poor photometric redshifts are more likely to have spectroscopic redshifts. Throughout the rest of this paper we adopt the z spec from AGES where possible, otherwise we use the consensus z phot estimates. The radio power versus redshift for these sources is shown in Fig. 4.

Contamination
In order to estimate the level of contamination by random matches, we generated 15 radio catalogues by randomising the positions of the sources in the real radio catalogue. We then cross-matched these 15 random radio catalogues with the optical sources in the same manner as described in the previous section. The distribution of optical identifications in stellar mass are plotted in Fig. 5. The contamination is high for sources with low stellar masses, likely driven by the higher surface density of faint optical galaxies with low stellar masses. The total contamination is ≈ 10 per cent for all the sources with stellar masses M * < 10 12 M . However, for sources with stellar masses M * < 10 9 M the contamination exceeds 90 per cent. We therefore do not consider optical identifications with objects with stellar masses below this value in later analysis.

PANCHROMATIC SED FITTING
While Gürkan et al. (2014) have suggested that for high radio luminosities, a single cut in 22 µm flux can be used to separate LERGs and HERGs, this may not be the case at lower luminosities (Janssen et al., in preparation). We therefore attempt to separate these sources based on further SED fitting to determine the relative contributions of the AGN and galaxy components at IR wavelengths. For this fitting we have included the FIR fluxes of these sources at 250 µm, 350 µm, and 500 µm by matching to the HerMES catalogue for the Boötes Field (Oliver et al. 2012). The Herschel fluxes were found by extracting the flux densities and errors directly from the DR4 maps at the optical source positions in the manner described by Hardcastle et al. (2013). The FIR fluxes are particularly important here to constrain the separation between the star-forming and AGN components. In order to decompose the SEDs of all the matched radio-optical sources we fit all the available multiwavelength photometry, including FIR, using the MCMC-based algorithm AGNFITTER (Calistro Rivera et al. 2016). Calistro Rivera et al. used AGNFITTER to separate star-forming galaxies and AGN. An example fitted SED is shown in Fig. 6. We note here that this fitting is dependent on the photometric redshifts and we have incorporated the full photometric redshift PDFs in the AGNFITTER analysis. For each source with a photometric redshift we produce 100 samples from the photometric redshift PDF, run AGNFITTER for each sample, and combine the PDFs for the individual AGNFITTER parameters. For sources with spectroscopic redshifts, we use those as a single sample.
The advantage of using AGNFITTER is that it infers the posterior probability density functions (PDFs) of the fitting parameters. This allows correlations and degeneracies among parameters to be recognised and allows for a robust calculation of the uncertainties for the inferred parameter values, and allows us to fold in the photometric redshift PDFs into the analysis. As described in Calistro Rivera et al. (2016) and following their nomenclature, the total model in AGNFITTER is the sum of the emission from the host galaxy and nuclear AGN. The host galaxy emission is represented by both stellar emission (GA) and the reprocessed light from  2000) dust extinction law covering a broad range in star-formation rates, including quiescent galaxies. The SB models used are the templates from Chary & Elbaz (2001) and Dale & Helou (2002), again covering a range in star-formation rates. A fit to a quiescent galaxy will yield a negligible SB component. The stellar templates come from the models of Bruzual & Charlot (2003) with a Chabrier (2003) initial mass function. The nuclear hot-dust torus models are taken from Silva et al. (2004).
A small fraction of sources (22, ≈ 1 per cent) have very poor fits, i.e. have AGNFITTER likelihood values < −100. These are excluded in further analysis. Some examples of the AGNFITTER SEDs with components in the three redshift intervals are shown in Appendix B in Fig. B1 and Fig. B2 for sources with good quality fits (quantified by likelihood values close to −1) and in Fig. B3 for sources with poor fits (quantified by likelihood values −20).
We compare the stellar masses and SFRs returned by AG-NFITTER to those we have derived using FAST (see Section 4.2). This comparison is shown in Fig. 7. We do not use SFRs in the subsequent analysis, but do use the fitted starburst components from AGNFITTER in classifying radio sources as star-forming. While the two codes are used to fit the same data (with the exception that the longest wavelength MIR and FIR data is included for the AG-NFITTER fits), the fitting methods and templates, specifically the inclusion of the AGN components in AGNFITTER, used are inde- pendent. AGNFITTER returns SFRs measured both from the stellar templates in the optical-UV, like FAST, as well as in the IR. In the comparison here, we use the AGNFITTER SFRs inferred from the total IR luminosities, because these are related to the total SB IR luminosity used in the following section to classify star-forming galaxies. Although not shown here, there is good agreement between the IR and optical SFRs derived by AGNFITTER, with a few extremely large optical SFRs, like those observed in the FAST SFRs in Fig. 7. Stellar-template-derived optical-UV emission SFRs are prone to significant dust extinction and can be less reliable.

Star Formation
The radio power-redshift distribution (Fig. 4) shows that we are mostly sensitive to high power sources at intermediate redshifts, while at low redshifts the opposite is true. If the radio emission is driven by star formation alone, then at radio powers P 150 MHz 10 25 W Hz −1 (corresponding to P 1.4 GHz 10 24 W Hz −1 , for a spectral index of −0.7), the required star formation rate is in excess of 25 M yr −1 (Condon 1992). This is not unreasonable for starforming galaxies at these intermediate redshifts, so we may expect some contamination of our RL AGN sample by star-forming galaxies. To explore this further, we consider the total SB IR luminosi- ties, L IR , defined as the total SB IR luminosity, L IR , integrated over the rest-frame wavelength range 1 < λ < 100 µm. Since this L IR is calculated from the fitted rest-frame component, no k-correction is needed. The values of q IR = L IR /L 150 are plotted as a function of redshift and radio power in Fig. 8, where the radio luminosities have been k-corrected assuming a spectral index of −0.7. As expected, below a redshift of ∼ 1, most of the sources lie on the FIRradio correlation and their radio emission can be attributed to star formation alone. The opposite is true at higher redshifts, but there remain some sources near the FIR-radio correlation particularly in the lower radio power range 10 25 P 150 MHz 10 26 W Hz −1 . We therefore consider galaxies with values within 2σ of the FIR-radio correlation as star-forming and exclude them from the subsequent RL AGN analysis. We use the FIR-radio correlation of Calistro Rivera et al. (2017), q IR = 1.72 (1 + z) −0.22 (σ = 0.529), which is based on LOFAR and Herschel measurements.

HERG/LERG separation
We aim to differentiate between HERGs ('cold mode' or 'radiative mode' sources), and LERGs ('hot mode' or 'jet mode' sources) based on their broadband SED information. In the remainder of the paper we use the nomenclature of HERGs and LERGs for succinctness. The AGNFITTER parameters of interest for our purposes are the disentangled host galaxy and AGN luminosities that together contribute to the emission at MIR wavelengths. The total MIR emission is the sum of the AGN torus luminosity L ν,T O , as well as the stellar emission L ν,GA and the reprocessed cold/warm dust emission L ν,SB . To allow for comparison, we define the integrated luminosities, L T O , L GA , L SB , and L BB , over the respective templates in a single rest-frame wavelength range 1 < λ < 30 µm. From these integrated luminosities we calculate the value, which quantifies the fraction of MIR emission that arises from an AGN compared with that from the stellar component of the host galaxy, independent of the MIR star-forming component. The MIR AGN emission we consider here mainly arises from the torus component, but we also include the BB component, which despite peaking at optical/uv wavelengths can be comparable to the GA component at these MIR wavelengths for strong, unobscured Type I AGN.
We specifically do not consider in the denominator of this ratio the SB component, as this arises from the reprocessed cold/warm dust in star-forming galaxies and can dominate the MIR part of the SED in moderately strong AGN in star-forming galaxies. This allows us to identify AGN across a range in star-formation rates. Essentially the ratio of eq. 4 is a measure of the strength of the MIR AGN emission compared to the stellar mass of the galaxy, traced by the GA component. The error of f AGN is calculated by propagating the errors on L T O and L GA given by AGNFITTER. We expect that HERGs will have significant contribution to the IR emission from the torus (or obscuring structure) and thus have large f AGN values, while LERGs will have little or no such contribution and low f AGN values. Similarly, we define the quantity to quantify the fraction of IR-emission due to star formation relative to that from the galaxy, independent of the AGN emission. The error of f SF is calculated by propagating the errors on L SB and L GA given by AGNFITTER.
To investigate how these values correspond to classifications based on spectroscopy, we cross-matched sources in the full optical photometry catalogue with sources from the SDSS Data Release 12 spectroscopic sample (DR12; Alam et al. 2015), using a simple nearest neighbour match within 1 arcsec. This yielded a sample of 2 315 sources for which we have an SDSS spectral classification and the same set of broadband photometry used for the radio sources in this paper. As this is an SDSS-selected sample, it is restricted in redshift to z < 0.3. Using the spectroscopic redshifts 2 from SDSS we fitted the broadband SEDs in AGNFITTER. The resulting distribution of the derived f AGN and f SF values is shown in Fig. 9, separated by their SDSS spectral classification, which is based on the optical emission lines in the SDSS spectra (see Bolton et al. 2012). There is some overlap at intermediate values, but the respective populations generally occupy different regions. It is clear that most of the 387 sources classed as 'AGN' in SDSS DR12 (with spectral class either 'QSO' or 'AGN') have values of f AGN close to one, indicating the presence of an excess of MIR emission from a torus as expected from these optical AGN. Similarly, the 217 sources identified as 'star-forming' (SDSS spectral class either 'STARFORMING', or 'STARBURST') have values of f SF close to one. Finally, the 1 711 remaining galaxies (SDSS spectral class 'GALAXY', i.e. those without any significant spectral lines) generally have small f AGN and f SF values, consistent with them being quiescent. The population of galaxies with large f SF values could be explained by a lack of signal-to-noise in the spectral lines necessary to meet the requirements to be classified as star-forming.
Since LERGs are expected to have no contribution from torus of accretion disk emission, we expect LERGs to have small values of f AGN . We consider the radio sources from the z < 0.3 Best & Heckman (2012, hereafter BH12) sample, which are separated into HERGs and LERGs based on emission line diagnostics. Given 2 Using our derived photometric redshifts for these sources instead of their SDSS spectroscopic redshifts in the AGNFITTER fits yields very similar results, and only increases the scatter slightly for this low redshift optically bright sample. the space density of sources in the BH12 catalogue, it is not unexpected that we find only LERGs within the ∼ 9 deg 2 Boötes field. All these LERGs have f AGN 0.1. While there is no strict boundary separating these sources, four per cent of normal galaxies have f AGN > 0.25, and one per cent of optical AGN have f AGN < 0.25. Based on this we define a separation of f AGN = 0.25 and in what follows we classify radio sources above this value as HERGs and below as LERGs. Finally, it should be noted that not all misclassifications are necessarily due to the faults of the SED fitting. It is possible that the spectral classifications here miss Type II obscured AGN.
We have further considered how this classification works for the 875 radio sources in the Boötes field for which we have optical spectra from AGES (Kochanek et al. 2012). The advantage of this test is that it probes sources fainter than the SDSS sample above. Similarly to the SDSS spectral classification, we use a BPT (Baldwin et al. 1981) classification for SF galaxies. While this limits the sample to z < 0.35 (366 sources) where the relevant emission lines lie within the AGES spectral coverage, it provides a clean separation between star-forming galaxies and AGN. We measured the strength and width of the from Kauffmann et al. (2003) for sources with SNR> 3 in all these lines to classify SF galaxies. The remaining sources are assumed to be RL AGN and are separated into HERGs and LERGs based on the strength and equivalent width of the [O III]λ 5007Å line (Pracy et al. 2016;Best et al. 2014;Best & Heckman 2012 Fig. 10 shows the distribution of the derived f AGN and f SF values for thse sources. As for the SDSS sample, a similar trend is seen, in that the SF galaxies generally have f SF close to one while HERGs have high f AGN compared to the lower f AGN values of LERGs. One per cent of LERGs have f AGN > 0.25 while 14 per cent of HERGs have f AGN < 0.25. It should be noted that the SF/LERG/HERG classification shown in this plot pertains only to the classification based on the optical spectra of these sources and as such the radio emission may be a result of either star-formation or AGN activity, thus most of the sources with high f SF values lie on the FIR-radio correlation (cf. Section 6.1).
Finally, we investigated a single cut in 22 µm flux separating LERGs and HERGs (Gürkan et al. 2014). Similar to (Janssen et al., in preparation), we find a large spread in 22 µm flux for both HERGs and LERGs classified either by their f AGN values or their optical spectra. However HERGs do generally have higher 22 µm fluxes than LERGs.

Radio AGN at intermediate redshifts
The aim of this paper is to study the population of RL AGN at intermediate redshifts. The radio power-redshift plot (Fig. 4) shows that at low redshifts, z 0.3, the radio-optical sample is dominated by low luminosity radio sources and contains very few high power sources, while at higher redshifts we can only probe high power sources. For this reason we cannot use this sample to directly compare high luminosity sources at both low and high redshift. The wide LOFAR surveys will provide the areal coverage needed for such a comparison. The rms in the radio map varies across the field of view (see Williams et al. 2016) between 100-250 µJy beam −1 , meaning that at a given redshift the lowest-power sources can only be detected over a smaller area. We do not make a cut on radio flux density, but account for incompleteness resulting from the varying detection area later (Sections 7.1.4 and 7.2). From the P − z plane it is clear that there is increasing imcompleteness above z = 2 and that at this redshift we can observe sources only with radio powers above P 150 MHz ≥ 10 25 W Hz −1 . Therefore, to compare the same sources across redshifts in this Section, we study only the high power (P 150 MHz ≥ 10 25 W Hz −1 ) sources at intermediate redshifts 0.5 ≤ z < 2. The final sample consists of 624 sources, which we divide into three redshift intervals: (i) 0.5 ≤ z < 1.0 (134 sources), (ii) 1.0 ≤ z < 1.5 (262 sources), (iii) 1.5 ≤ z < 2.0 (228 sources).
These numbers reflect the number of sources in each bin based on their best photometric redshifts. We note that the final redshift bin may be incomplete below P 150 MHz 10 25.5 W Hz −1 .

Local Reference Sample
As a local comparison sample we use the catalogue compiled by BH12. This matched radio-optical catalogue was constructed from the seventh data release (DR7; Abazajian et al. 2009) of the Sloan Digital Sky Survey (SDSS) spectroscopic sample and the NRAO Very Large Array (VLA) Sky Survey (NVSS; Condon et al. 1998) and the Faint Images of the Radio Sky at Twenty centimetres (FIRST; Becker et al. 1995). The optical data includes parameters from the value-added spectroscopic catalogues (VASC) created by the Max Plank Institute for Astrophysics and Johns Hopkins University (MPA-JHU) group 3 (Brinchmann et al. 2004). This includes information from the imaging data such as magnitudes and sizes (York et al. 2000), as well as derived properties including the stellar mass based on fits to the photometry (Kauffmann et al. 2003). The spectroscopy also provides D n 4000 (Balogh et al. 1999) which, like galaxy colour, provides a guide to the stellar population age. BH12 separated the sources into star-forming galaxies and RL AGN (7302 sources), the latter of which are further sub-divided into HERGs and LERGs, based on their optical photometric and spectroscopic parameters. Noting the different observed radio frequency, we select sources with P 1.4 GHz > 10 24 W Hz −1 , broadly comparable to P 150 MHz > 10 25 W Hz −1 , assuming a spectral index of α = −0.7. This local radio-optical sample consists of 3736 radio sources between 0.01 < z ≤ 0.3.

HERG/LERG composition
The distribution of f AGN values for our intermediate redshift and high power, P 150 MHz ≥ 10 25 W Hz −1 , samples is plotted in Fig. 11, where we show the distribution for all radio sources, and within each of the three redshift intervals. There is a maximum in the overall distribution for sources with 0.9 < f AGN < 1, a second maximum for sources with 0.1 < f AGN < 0.2, and a minumum around f AGN ≈ 0.5. The two highest redshift intervals both show the growing peak at f AGN ≈ 1, which suggests that there are more 'strongly AGN-dominated' sources and fewer 'AGN-free' sources in these intervals than at 0.5 < z ≤ 1.0.
The number of HERGs and LERGs in each redshift interval are given in Table 1, including those for the local reference sample. The percentage of HERGs and LERGs within each redshift interval is given relative to the total number of radio sources in that interval. Here we have considered the sources that lie on the FIR-radio correlation (see Section 6.1) as SF galaxies. It can be seen from these numbers that the fraction of HERGs rises between z ≈ 0 and 0.5 ≤ z < 1.0, and then again between 0.5 ≤ z < 1.0 Figure 11. Distribution of the fraction of IR AGN emission, f AGN , defined by equation 4 for the full sample (black), and the three redshift intervals: 0.5 ≤ z < 1.0 (orange), 1.0 ≤ z < 1.5 (red), and 1.5 ≤ z < 2.0 (dark red). The dotted vertical line shows the value of f AGN = 0.25 that we use to separate HERGs and LERGs. and 1.0 ≤ z < 1.5. The fraction of LERGs on the other hand is strongly declining between all redshift intervals. The spectroscopic completeness for both AGN types in the first redshift interval is ∼ 30 per cent, which for LERGs drops to below ∼ 2 per cent in the two higher redshift intervals, while for HERGs it only drops to ∼ 15 per cent.

Colour-mass distribution
We now consider the distribution of the RL AGN, both HERGs and LERGs, in colour-mass space. This is plotted in Fig. 12 4 for both optical and radio sources where we plot the 0.1 (u − r) colour (defined in Section 4.2) against stellar mass for both optical and radio sources in each of the four redshift intervals. The f AGN values for the radio sources are colour-coded in the 2-d distribution and the 1-d distributions of both stellar mass and colour are shown for the optical and radio sources as well as separately for the HERGs and LERGs. Here we use the value f AGN = 0.25 to separate the HERGs and LERGs. In comparing the local and higher redshift samples, we note that the parameters used for the HERG/LERG separation are different. However, they provide a qualitative comparison for the distribution of the radio and optical source populations in colourmass space. Given the use of photometric parameters we expect 4 A preliminary version of this figure was presented in Williams et al.
Table 2. Two-sided Kolmogorov-Smirnov statistics in comparing the HERG and LERG distributions in colour and mass within each redshift interval.
0.5-1.0 0.62 1.1 · 10 −10 0.57 4.3 · 10 −9 1.0-1.5 0.46 2.5 · 10 −6 0.48 9.6 · 10 −7 1.5-2.0 0.44 4.8 · 10 −3 0.61 2.3 · 10 −5 there to be some fraction of catastrophic outliers. In particular, a few points at very high stellar masses, notably in the highest redshift bin, could be a result of poorly determined photometric redshifts and fitted masses. The colour-mass distributions of optical and radio sources are clearly different, which is not unexpected. At all redshifts radio sources tend to be more massive and redder compared to the parent galaxy population. The properties of HERGs and LERGs are also different, in that the HERGs span a wider range of stellar masses 10 9 < M * /M < 10 11.5 and colours. Going from the lower redshift bin, 0.5 ≤ z < 1.0, to the higher redshift bins, we see that the colour distribution of LERGs changes from showing a clear peak at 0.1 (u − r) > 2.5, to becoming flatter where red and blue galaxies approximately contribute the same. Similarly for HERGs, we see an increasing contribution of very blue objects, 0.1 (u − r) < 0.1, with redshift. In general, though, LERGs are always more likely to be hosted by massive red galaxies.
We compute the two sample Kolmorogorov-Smirnov two sided test statistics, and in all cases can reject the null hypothesis that the two samples are randomly drawn from the same distribution. The Kolmogorov-Smirnov statistics and p-values are given in Table 2. In the highest redshift bin in these plots the radio source population is slightly incomplete due to the varying rms in the LO-FAR map (see Section 7.1).

Radio-loud fraction
The mass-dependence of the RL fraction can be an indicator of the accretion mode of the RL AGN, largely because of the different dependence of the fuelling source (hot vs. cold gas) on stellar mass (Best et al. 2006) and the relationship between black hole and stellar mass for elliptical galaxies. As this is not a volume-limited sample, following Janssen et al. (2012) and Williams & Röttgering (2015), we use the RL fraction defined by: where the sets A and R are, respectively, all galaxies and all radio sources in a given bin, defined by the parameters of mass (x) and accretion mode (y) using the classification in Section 7.1.2 for radio sources. The maximum accessible volume over which each source can be observed, V i , is determined by both the minimum and maximum distance at which a given source would be included in the sample as a result of the selection criteria: where V max and V min are the volumes enclosed within the observed sky area out to the maximum and minimum distances respectively. The minimum accessible volume is a result of the lower redshift limit in a given bin. The maximum accessible volume is determined by the flux limits of the optical (< 24 mag) and radio rms Figure 12. Colour, 0.1 (u − r), versus stellar mass (main panels). The density of optical sources with m I < 24 mag is plotted in black, in log units, and the radio sources are plotted in colour. The subpanels show the normalised stellar mass (top) and colour distributions (right) for optical sources (grey) and radio sources (black). The HERG ( f AGN > 0.25) and LERG ( f AGN < 0.25) distributions are shown in cyan and magenta respectively, normalised to the total number of radio sources. In the sub-panels, the HERG distributions are multiplied by a factor of 10 for visibility. Four redshift intervals are plotted separately: (a) the local BH12, 0.01 ≤ z < 0.3, spectroscopic sample where HERGs and LERGs are classified spectroscopically; (b-d) the three higher redshift samples from the Boötes field where the fill-colour of the radio points indicates their f AGN values.
map as well as the upper limit of the redshift bin. Following Hardcastle et al. (2016), the radio V max is calculated as d max dA. The completeness function is determined from an rms map created by masking the LOFAR rms map by the optical coverage area, which excludes regions around bright optical sources. The total sky area of the masked map is 9.27 deg 2 . For the optical, we use the rest frame I-band magnitude determined from INTERREST to compute the V max . In order to take into account the uncertainties on the photometric redshifts, we consider the full photometric redshift PDFs from Duncan et al. (in prep) in the calculation of the RL fractions. We do this by generating 100 realisations where the redshifts for each source are randomly drawn from their respective PDFs 5 . We compute for each realisation the RL fraction described above. Uncertainties are calculated as the statistical Poissonian errors. We Table 3. Slope of the RL fraction as a function of mass for HERGs and LERGs in the three redshift intervals.
z HERGs LERGs 0.5-1.0 0.55 ± 0.23 2.31 ± 0.29 1.0-1.5 0.67 ± 0.18 2.18 ± 0.13 1.5-2.0 1.17 ± 0.20 1.79 ± 0.33 then take the median of the RL fraction realisations and the errors given are based on the 16th and 84th percentiles. This also takes into account which sources are included in each redshift interval. The RL fraction for all radio sources, and separately for HERGs and LERGs is shown in Fig. 13. It can be seen that the RL fraction for HERGs at all masses increases with redshift. It is possible the RL fraction in the highest redshift bin is slightly underestimated due to incompleteness between 10 25 < P 1.4 GHz < 10 25.5 W Hz −1 (cf. Section 7.1). The rising fraction of HERGs in massive galaxies becomes similar to the fraction of LERGs in the most massive galaxies at the highest redshifts. These results are consistent with those from many of the earlier radio surveys, which suggest that the most powerful radio galaxies (P 1.4 GHz 10 26 W Hz −1 ) at z 1 are predominantly HERGs hosted by the most massive (M * 10 12 M ) galaxies (e.g. Eales et al. 1997;Jarvis et al. 2001;Seymour et al. 2007;Fernandes et al. 2015). The dependence on stellar mass is flatter for HERGs than it is for LERGs at all redshifts. Table 3 gives the value of the slope fitted to log f RLlog M * /M over the range 10 9 < M * /M < 10 12 . The shape of the RL fraction as a function of stellar mass for LERGs remains almost the same steep function, f RL ∼ M 2.0−2.5 * . This is very similar to that observed in the local Universe (Tasse et al. 2008b;Simpson et al. 2013). Importantly, this suggests that the fueling mechanism for LERGs remains the same out to redshifts of about 2, a prediction made by B14.

Luminosity functions for HERGs and LERGs
The luminosity functions (LFs) of the different radio AGN populations were first constructed by Best & Heckman (2012). Since then, B14 have looked at the evolution of these LFs out to z ∼ 1. Recently, Pracy et al. (2016) studied their evolution in finer detail using a much larger sample of over 5000 sources with spectroscopic redshifts and classifications, again covering the range 0.005 < z < 0.75. Within this redshift range LERGs show little evolution, explained by the very slow evolution of available massive elliptical host galaxies. However, these are predicted to drop off more strongly after z > 1 (see B14). In this paper we probe this evolution out to z ∼ 2, making use of photometric redshifts and the photometric AGN classification described in the previous section.
In order to take into account the uncertainties on the photometric redshifts, we consider the full photometric redshift PDFs from Duncan et al. (in prep) in the calculation of the LFs, in a similar way to the RL fraction. We do this by generating 100 realisations of the LFs, where in each realisation the redshifts for each source are randomly drawn from their respective PDFs 6 . We compute for each realisation the individual radio LFs for the SF sources and the two populations of AGN in the standard way, using the ρ = Σ i 1/V i formalism (Schmidt 1968;Condon 1989) where V max is calcuated the same way as in Section 7.1.4. Uncertainties are calculated as the statistical Poissonian errors. We then take the median of the LF realisations and the errors given based on the 16th and 84th percentiles. Given the large area covered, we assume that the effects of cosmic variance are negligible even for very massive galaxies.
The local LFs, split into AGN and SF categories, have been well studied at higher frequencies (Best et al. 2005a;Mauch & Sadler 2007;Prescott et al. 2016), and recently also with LOFAR (Hardcastle et al. 2016). Although not plotted here our SF and total AGN LFs show very good agreement with those of Hardcastle et al. (2016). The AGN LFs have been studied for the two accretion modes to a lesser extent only more recently (Best & Heckman 2012;Pracy et al. 2016), with the reasons for the observed differences largely attributed to different selections and classification schemes (Pracy et al. 2016). Due to the relatively small areal coverage of the sample presented here, it is not very well suited to probing the AGN population in the local universe; however, we have included the 0.01 < z < 0.3 LFs for the various populations here (Fig. 14) as a comparison and as test of the ability of the photometric methods to separate the populations. For comparison we show the LFs from Pracy et al. (2016) and BH12 and find our LFs to be broadly consistent with these.
At intermediate redshifts we can probe the AGN LFs across a broad range in radio power, up to P 150 MHz ≈ 10 28 W Hz −1 . The LFs for the RL AGN classified as HERGs and LERGs are shown in Fig. 15 and tabulated in Table 4, for three redshift intervals 0.5 < z ≤ 1.0, 1.0 < z ≤ 1.5, and 1.5 < z ≤ 2.0. The SF LFs are also included. From these plots it is clear that the HERG LFs show no statistically significant evolution across this redshift range -the LFs differ by less than 2σ . This is in contrast to the strong evolution seen between z < 0.3 and 0.5 < z < 1.0 . Despite the lack of evolution there is a suggestion that at the high power end (P 150 MHz 10 26 W Hz −1 ) the space density peaks in the second redshift interval 1.0 < z ≤ 1.5. This is only marginally significant as the space density above P 150 MHz > 10 26 W Hz −1 in the second interval is only ∼ 2σ above that in the first interval and ∼ 1.5σ above that in the third interval. Nevertheless, it is consistent with what is known about the space density evolution of the RL AGN population as a whole. The space density of the highest power (P 1.4 GHz 10 25 W Hz −1 ) radio sources is well known to peak at z ≈ 1 (Dunlop & Peacock 1990;Rigby et al. 2011Rigby et al. , 2015. In contrast, the space density of LERGs is strongly declining significantly. The evolution of the LERG population is of particular interest. In Fig. 16 we plot only the LERG LFs along with the models of B14 at the relevant redshifts. These are models fitted to the LERG LFs out to z < 1. For more details of these models we refer the reader to B14 and mention here the key points. The 'model 1' models are for pure density evolution, where the evolution of LERGs is driven purely by the strongly declining population of quiescent galaxies, with a possible time lag between the formation of a quiescent galaxy and the onset of radio activity (variant b). The 'model 2' models are similar but include a luminosity evolution of the sources, i.e. the luminosities of the radio sources systematically increase with redshift. The third type of model, 'model 3' includes a contribution of HERG-type sources. The LERG LFs do suggest the 'model 2a' or 'model 2b' are prefered over some of the other models. Figure 13. The fraction of galaxies hosting a radio source (RL fraction) as a function of stellar mass for a radio-power cut-off of P 150 MHz > 10 25 W Hz −1 in the three redshift bins, for HERGS (left) and LERGs (right). The errors are determined from Poisson statistics.

CONCLUSION
We have presented a catalogue of radio sources from the 150-MHz LOFAR observations of the Boötes field (Williams et al. 2016) matched to optical sources from the multi-band photometry catalogue (Brown et al. 2007) available for the Boötes field, which have photometric redshifts from Duncan et al. (in prep). We have performed SED fitting using AGNFITTER (Calistro Rivera et al. 2016), including the FIR flux measurements, to ascertain the relative contribution of AGN and galaxy emission in the MIR (1 < λ < 30 µm) and used this to classify radio galaxies as HERGs or LERGs. We have used a well-defined sub-sample of 624 highpower ( P 150 MHz > 10 25 W Hz −1 ) radio sources with redshifts between 0.5 ≤ z < 2 to study the RL fraction as a function of galaxy mass. We have derived the radio LFs separately for HERGs and LERGs. Our key findings can be summarised as the following: (i) the fraction of HERGs is significantly higher at z = 2 compared to z = 0; (ii) LERGs are generally found in passive galaxies, while the host galaxies of HERGs are both star forming and passive; (iii) LERGs tend to be hosted by more massive galaxies than HERGs and are predominantly found in redder galaxies, but at higher redshifts can also be found in bluer galaxies; (iv) the fraction of galaxies hosting LERGs is a strong function of stellar mass like it is in the local Universe, suggesting these sources are still fuelled by accretion of hot gas; (v) at moderate redshifts, 0.5 < z < 2.0, the LF of LERGs undergoes strong negative evolution consistent with the decline of quiescent galaxies in the Universe.
While we have used the full photometric redshift PDFs to take into account the uncertainties on their estimates, this kind of panchromatic SED decomposition would benefit greatly from spectroscopic redshifts, in particular for the sources with broad photometric redshift distributions. We have used a local comparison sample at low redshift, there is no directly comparable sample to track the evolution from z ≈ 0. In the near future, the full LO-FAR surveys (Shimwell et al. 2017), both wide and deep, combined with spectroscopic information from WEAVE-LOFAR (Smith et al. 2016), providing both redshifts and spectral classification, will al- Figure 15. The radio LFs for our three intermediate redshift samples, separated into SF galaxies, HERGs and LERGs. For comparison and as a guide for the eye, the z < 0.3 (dotted) and 0.5 < z < 1.0 (dashed) HERG and LERG LFs fitted by B14 are shown in all panels. We also show the simple model SF LFs from Novak et al. (2017) at each redshift interval. low us to address both these issues and will allow this to be done with greater accuracy and in more detail for larger samples.    Figure B1. Examples of good quality SED fits (with likelihood ≈ −1) with significant star formation contribution (high f SF values) for HERGs (left) and LERGs (right) in the three redshift bins 0.5 < z ≤ 1.0 (top), 1.0 < z ≤ 1.5 (middle) and 1.5 < z ≤ 2.0 (bottom). In all cases ten realisations from the parameters' posterior probability distributions are plotted giving an indication of the uncertainties in the fitted components. These show the total SED (red) and the individual components: the AGN torus (TO; purple), the starburst (SB; green), the galaxy (GA; yellow) and the blue bump (BB; blue). The red points show the total SEDs integrated across the filter bandpasses and the black points with errorbars show the observed luminosities. Figure B2. Examples of good quality SED fits (with likelihood ≈ −1) with minimal star formation contribution (low f SF values) for HERGs (left) and LERGs (right) in the three redshift bins 0.5 < z ≤ 1.0 (top), 1.0 < z ≤ 1.5 (middle) and 1.5 < z ≤ 2.0 (bottom). In all cases ten realisations from the parameters' posterior probability distributions are plotted giving an indication of the uncertainties in the fitted components. These show the total SED (red) and the individual components: the AGN torus (purple), the starburst (green), the galaxy (yellow) and the blue bump (blue). The red points show the total SEDs integrated across the filter bandpasses and the black points with errorbars show the observed luminosities. Figure B3. Examples of poorer SED fits (with likelihood ≈ −20) for HERGs (left) and LERGs (right) in the three redshift bins 0.5 < z ≤ 1.0 (top), 1.0 < z ≤ 1.5 (middle) and 1.5 < z ≤ 2.0 (bottom). In all cases ten realisations from the parameters' posterior probability distributions are plotted giving an indication of the uncertainties in the fitted components. These show the total SED (red) and the individual components: the AGN torus (purple), the starburst (green), the galaxy (yellow) and the blue bump (blue). The red points show the total SEDs integrated across the filter bandpasses and the black points with errorbars show the observed luminosities.