The LOFAR Two-metre Sky Survey: Deep Fields data release 1. V. Survey description, source classiﬁcations, and host galaxy properties

Source classiﬁcations, stellar masses, and star-formation rates are presented for ≈ 80 000 radio sources from the ﬁrst data release of the Low Frequency Array Two-metre Sky Survey (LoTSS) Deep Fields, which represents the widest deep radio surv e y ev er undertaken. Using deep multi-wavelength data spanning from the ultraviolet to the far-infrared, spectral energy distribution (SED) ﬁtting is carried out for all of the LoTSS Deep host galaxies using four different SED codes, two of which include modelling of the contributions from an active galactic nucleus (AGN). Comparing the results of the four codes, galaxies that host a radiative AGN are identiﬁed, and an optimized consensus estimate of the stellar mass and star-formation rate for each galaxy is derived. Those galaxies with an excess of radio emission o v er that e xpected from star formation are then identiﬁed, and the LoTSS Deep sources are divided into four classes: star-forming galaxies, radio-quiet AGN, and radio-loud high-excitation and low-excitation AGN. Ninety-ﬁve per cent of the sources can be reliably classiﬁed, of which more than two-thirds are star-forming galaxies, ranging from normal galaxies in the nearby Universe to highly-starbursting systems at z > 4. Star-forming galaxies become the dominant population below 150-MHz ﬂux densities of ≈ 1 mJy, accounting for 90 per cent of sources at S 150MHz ∼ 100 μ Jy. Radio-quiet AGN comprise ≈ 10 per cent of the o v erall population. Results are compared against the predictions of the SKADS and T-RECS radio sky simulations, and improvements to the simulations are suggested.


INTRODUCTION
Understanding the formation and evolution of galaxies requires a detailed knowledge of the baryonic processes that both drive and quench the process of star formation within galaxies across cosmic time.In this regard, the faint radio sky provides one of the most important windows on the Universe, as it offers a direct view onto three critical (and overlapping) populations of objects: star-forming galaxies, 'radio-quiet' active galactic nuclei (AGN), and low luminosity radio galaxies (e.g.Padovani 2016).
Arguably the most important observational test for any model of galaxy formation is measurements of the evolution of the cosmic star-formation rate density across cosmic time, and the distribution of that star formation amongst the galaxy population at each redshift, as a function of stellar mass, galaxy morphology, environment, and other properties.These crucial measurements require large, unbiased samples of star-forming galaxies over a wide range of redshifts.Much progress has been made in understanding the star-forming galaxy population, at least out to cosmic noon at  ∼ 2, using a variety of star-formation indicators (e.g.Madau & Dickinson 2014).The primary uncertainty is the effect of dust: by cosmic noon, around 85 per cent of the total star-formation rate (SFR) density of the Universe is dust-enshrouded (e.g.Dunlop et al. 2017), and a sub-millimetre (sub-mm) or far-infrared (far-IR) view of the Universe paints a very different picture of galaxy properties to that of a population selected at optical (rest-frame ultraviolet) wavelengths (e.g.Cochrane et al. 2021).Current far-IR surveys are limited by sensitivity to the more extreme systems, where contamination of the far-IR light by AGN emission is also a concern (e.g.Symeonidis & Page 2021).
Radio emission provides a tool to observe the activity of galaxies in a manner that is independent of dust.For sources without AGN, the low-frequency radio emission arises primarily from recent supernova explosions of massive (young) stars (see reviews by Condon 1992;Kennicutt 1998), and thus directly traces the current star-formation rate (unless sufficiently low radio frequencies are reached such that free-free absorption becomes important; e.g.Schober et al. 2017).New generation radio interferometers offer sufficient sensitivity and field-of-view to survey large samples of star-forming galaxies out to high redshifts.Crucially, they can also provide sufficient angular resolution that deep surveys are not generally affected by the source confusion that limits the capabilities of surveys with sub-mm and far-IR telescopes such as the Herschel Space Observatory, for which the vast majority of sources in deep surveys are blends (e.g.Oliver et al. 2012;Scudder et al. 2016).
Star formation within massive galaxies is widely believed to be regulated in some manner by AGN, due to the large outflows of energy associated with the growth of supermassive black holes.AGN activity occurs in two fundamental modes (e.g.see reviews by Heckman & Best 2014;Hardcastle & Croston 2020).At high accretion rates, accretion of material on to a black hole is understood to occur through a 'standard' geometrically-thin, optically-thick accretion disk (Shakura & Sunyaev 1973), in which around 10 per cent of the rest-mass energy of the accreting material is emitted in the form of radiation ('radiative' or 'quasar-like' AGN).These AGN can drive outflowing winds through thermal or radiation pressure (e.g.Fabian 2012, and references therein), which may have a substantial effect on the evolution of the host galaxy.Radiatively-efficient AGN sometimes possess powerful twin radio jets ('radio-loud' quasars or their edge-on counterparts, the 'high-excitation radio galaxies'; HERGs), and many recent works also suggest that even those that do not ('radio-quiet' AGN) frequently (or maybe even always) possess weak radio jets (Jarvis et al. 2019;Gürkan et al. 2019;Macfarlane et al. 2021;Morabito et al. 2022, and references therein).These AGN are detectable in deep radio surveys, either due to the weak radio jets or due to the star formation that can accompany the AGN activity.
At lower accretion rates, typically below about 1 per cent of the Eddington accretion rate, the nature of the accretion flow on to a supermassive black hole is believed to change: the accretion flow is thought to become geometrically thick and radiatively inefficient (Narayan & Yi 1994, 1995).A characteristic feature of these advection-dominated or radiatively-inefficient accretion flows is that most of the energy that they release is in the form of two-sided radio jets ('jet-mode' AGN; also referred to as 'low-excitation radio galaxies').These jet-mode AGN dominate the radio sky at intermediate flux densities (above a few mJy), and the radio waveband is by far the most efficient means of identifying these sources.Jet-mode AGN have been very wellstudied in the nearby Universe (e.g.Best & Heckman 2012), where it is now widely accepted that they play a critical role in the evolution of massive galaxies and clusters, providing an energy input that counter-balances the radiative cooling losses of the surrounding hot gas and thus preventing that gas from cooling and forming stars (see reviews by McNamara & Nulsen 2007;Fabian 2012;Kormendy & Ho 2013;Heckman & Best 2014;Hardcastle & Croston 2020, and references therein).Deeper radio surveys, probing the faint radio sky, enable these low-luminosity AGN to be detected and studied to higher redshifts (Best et al. 2014;Pracy et al. 2016;Williams et al. 2018;Whittam et al. 2022), and hence their role in the evolution of massive galaxies to be determined across cosmic time.
Deep radio surveys can therefore offer a unique insight into many aspects of the galaxy and AGN population.However, to extract the maximum science from deep radio surveys, it is essential that they are carried out in regions of the sky which are extremely well-studied at other wavelengths across the electromagnetic spectrum.The ancillary data are required to identify the radio source host galaxies, to estimate their redshifts, to classify the nature of the radio emission (star formation vs radiatively-efficient AGN vs jet-mode AGN) and to determine the physical properties of the host galaxies (stellar mass, star-formation rate, environment, etc).
Until recently, the state-of-the-art in wide-area deep radio surveys was the VLA-COSMOS 3 GHz survey (Smolčić et al. 2017a), which used the Very Large Array (VLA) to cover 2 deg 2 of the Cosmic Evolution Survey (COSMOS) field, arguably the best-studied degreescale extragalactic field in the sky.Smolčić et al. (2017b) investigated the multi-wavelength counterparts of the ≈10,000 radio sources detected, and provided classifications, which then allowed several further investigations of the radio-AGN and star-forming populations (e.g.Smolčić et al. 2017c;Novak et al. 2017;Delvecchio et al. 2017;Delhaize et al. 2017).Nevertheless, even the VLA-COSMOS 3 GHz survey does not have sufficient sky area to cover all cosmic environments, and may therefore suffer from cosmic variance effects, as well as having limited source statistics at the highest redshifts.The on-going MeerKAT International GigaHertz Tiered Extragalactic Exploration (MIGHTEE) 1.4 GHz survey aims to extend sky coverage at this depth to 20 deg 2 ; Heywood et al. (2022) provide an early release, with Whittam et al. (2022) deriving source classifications for 88 per cent of the ≈ 5, 000 sources with host galaxy identifications over 0.8 deg 2 in the COSMOS field.
The Low Frequency Array (LOFAR; van Haarlem et al. 2013) Two-metre Sky Survey (LoTSS) Deep Fields have a similar goal at lower frequency.The first data release (hereafter LoTSS-Deep DR1) was made public in April 2021: the radio data reach rms sensitivity levels ≈ 4 times deeper than the wider all-northern-sky LoTSS survey (Shimwell et al. 2017(Shimwell et al. , 2019(Shimwell et al. , 2022)), corresponding to approximately the same effective depth as the VLA-COSMOS 3 GHz survey (for a source with typical radio spectral index,  ≈ 0.7, where   ∝  − ) but over an order of magnitude larger sky area (Tasse et al. 2021;Sabater et al. 2021, hereafter Papers I and II respectively).An extensive optical and near-infrared cross-matching process has identified and provided detailed photometry for over 97 per cent of the ≈80,000 radio sources detected over the central regions of the target fields where the best ancillary data are available (a combined area of 25 deg 2 ; Kondapally et al. 2021, Paper III).These data have been used to provide high-quality photometric redshifts (Duncan et al. 2021, Paper IV).In this paper, the 5th of the series, these data are combined with far-IR data to carry out detailed spectral energy distribution (SED) fits to the multi-wavelength photometry from ultraviolet (UV) to far-IR wavelengths, using several different SED fitting codes.Using the results of this analysis, the radio sources are classified into their different types, and key physical parameters of the host galaxies, such as their stellar masses and star-formation rates, are determined.
The layout of the paper is as follows.In Sec. 2 the LoTSS Deep Fields survey is described: this section outlines the choice of target fields, and places the first data release in to the context of the eventual full scope of the survey.Sec. 3 then describes the data that will be used in the paper and outlines the application of the SED fitting algorithms.Sec. 4 describes how the results are used to identify the (radiative-mode) AGN within the sample.The results of the different SED fitting algorithms are compared in Sec. 5, and used to define consensus measurements for the stellar mass and star-formation rate of each host galaxy.Combining this information with the radio data, Sec.6 then describes the identification of radio-excess AGN.Sec. 7 summarises the final classifications of the objects in the sample, and investigates the dependence of these on radio flux density, luminosity, stellar mass and redshift.In Sec. 8 the results are compared against the predictions of the most widely-used radio sky simulations, and suggestions made for improvements to those simulations.Finally, conclusions are drawn in Sec. 9.The classifications derived are released in electronic form and are used for detailed science analysis in several further papers (Smith et al. 2021;Bonato et al. 2021;Kondapally et al. 2022;McCheyne et al. 2022;Mingo et al. 2022;Cochrane et al. 2023, and others).

LOFAR observations of the LoTSS Deep Fields
The International LOFAR Telescope (van Haarlem et al. 2013) is a remarkably powerful instrument for carrying out deep and wide radio surveys of the extragalactic sky, owing to its high sensitivity, high angular resolution (6 arcsec at 150 MHz when using only Dutch baselines, improving to 0.3 arcsec with the international stations included), and in particular its wide field-of-view.The primary beam full-width at half-maximum (FWHM) of the Dutch LOFAR stations is 3.8 degrees at 150 MHz, giving a field-of-view of more than 10 deg 2 in a single pointing.International stations have a larger collecting area and a correspondingly smaller beam: 2.5 deg FWHM; 4.8 deg 2 field-of-view.The LoTSS survey (Shimwell et al. 2017(Shimwell et al. , 2019(Shimwell et al. , 2022) ) is exploiting LOFAR's capabilities by observing the entire northern sky, with a target rms depth of below 100Jy beam −1 at favourable declinations (the non-steerable nature of the LOFAR antennas means that sensitivity decreases at lower elevations).Nevertheless, LoTSS only scratches the surface of the depth that radio surveys with LOFAR are capable of reaching.LoTSS provides an excellent census of the radio-loud AGN population which dominates the bright and intermediate radio sky, but samples only the brighter end of the radio-quiet AGN and star-forming galaxy populations which become dominant as the LoTSS flux density limit is approached.
The LoTSS Deep Fields provide a complementary deeper survey, aiming to reach a noise level of 10-15 Jy beam −1 over a sky area of at least 30 deg 2 .LoTSS-Deep is designed to have the sensitivity to detect Milky-Way-like galaxies out to  > 1, and galaxies with star-formation rates of 100 yr −1 to beyond  = 5 (e.g.Smith et al. 2016), as well as being able to detect typical radio-quiet quasars right out to redshift 6 (Gloudemans et al. 2021).The sky area makes it possible to: (i) sample the full range of environments at high redshifts -for example, it is expected to include 10 rich proto-clusters at  > 2; (ii) include statistically meaningful samples of rarer objects (such as  > 5 starbursts); (iii) build large enough samples of AGN and starforming galaxies (over 100,000 of each expected to be detected) to allow simultaneous division by multiple key properties, such as luminosity, redshift, stellar mass and environment.
LoTSS-Deep is being achieved through repeated 8-hr LOFAR observations of the regions of the northern sky with the highest quality degree-scale multi-wavelength data.The four target fields are the European Large Area ISO Survey Northern Field 1 (ELAIS-N1; Oliver et al. 2000), the Boötes field (Jannuzi & Dey 1999), the Lockman Hole (Lockman et al. 1986) and the North Ecliptic Pole (NEP); these are described in more detail in Section 2.3.
Table 1 outlines the anticipated final depths of each field based on awarded observing time.Scaling by depth and area from radio source counts in shallower LoTSS-Deep observations, the final LoTSS Deep Fields are expected to detect more than 250,000 radio sources within the central 35 deg 2 , overlapping the best multi-wavelength data.Figure 1 compares the sensitivity, field-of-view, and angular resolution of the LoTSS Deep Fields to other completed and on-going radio surveys.The final LoTSS Deep Fields dataset will be unrivalled in its combination of depth and area.The inclusion of the international stations will also provide an angular resolution which is unmatched by any competitor survey: indeed, at low frequencies, the LoTSS Deep Fields with international baselines will remain unique even in the era of the Square Kilometre Array (SKA).
In order to account for the smaller primary beam of the international stations, from LOFAR Observing Cycle 14 onwards the pointing positions for the LoTSS-Deep observations of the Lockman Hole, Boötes and NEP fields have been dithered around a small mosaic.The mosaics have been designed to ensure good coverage of the sky area with the best-quality multi-wavelength data, within the primary beam of the international stations, while keeping offsets small enough so that there is negligible loss of sensitivity over this region when imaging with only Dutch stations.

LoTSS-Deep DR1
This paper considers the radio source catalogues from the first LoTSS Deep Fields data release.LoTSS-Deep DR1 released the reduced LOFAR images and catalogues constructed from data taken before October 2018 (Paper I; Paper II), along with the optical/IR catalogues and host galaxy identifications (Paper III) and photometric redshifts (Paper IV).These LoTSS-Deep DR1 LOFAR observations focused on the ELAIS-N1, Boötes and Lockman Hole fields, due to the earlier availability of the multi-wavelength data in those fields.The LoTSS-Deep DR1 LOFAR images included only the data from the Dutch LOFAR stations, not the international stations, due to the additional complications associated with calibrating the long baselines and the The area of best ancillary data is defined in Paper III. Quoted rms noise levels are those at the centre of the field.The marginally lower sensitivity in Boötes compared to the other fields is due to its lower declination, and hence lower average elevation during the observations.The 'number of sources in DR1 full area' quoted is over the full catalogues presented in Paper I and Paper II, out to the 30 per cent power point of the primary beam (ie.over ∼ 25 deg 2 in each field).symbol shows what should be achievable after including the international stations (improved angular resolution, additional depth due to the extra collecting area, but smaller areal coverage due to the smaller primary beam of the international stations).Descriptions of the surveys included on the plot (listed from high to low effective rms depth) can be found in the following references: GLEAM (Wayth et al. 2015); WENSS (Rengelink et al. 1997); TGSS (Intema et al. 2017); SUMSS (Mauch et al. 2003); NVSS (Condon et al. 1998); GLEAM-X (Hurley-Walker et al. 2022); RACS (Hale et al. 2021); FIRST (Becker et al. 1995); XXL-GMRT (Smolčić et al. 2018); VLASS (Lacy et al. 2020) associated computing requirements (see e.g.Morabito et al. 2022;Sweĳen et al. 2022, for a description of recent advances towards a pipeline for international stations).The data allow an angular resolution of 6 arcsec to be achieved: higher angular resolution images will be produced in later data releases.
As shown in Table 1, the images in LoTSS-Deep DR1 already reach an rms noise level below 20Jy beam −1 at 150 MHz at the centre of the deepest field (ELAIS-N1), away from bright sources.Sensitivity decreases with primary beam attenuation towards the outer regions of the field; dynamic range effects are also present around bright sources but only a few percent of the image suffers from significantly increased noise levels due to these calibration issues (Paper I; Paper II).Over 170,000 sources are catalogued, with peak flux densities above 5 times the local rms noise, across the full radio area of the three fields; as with all radio catalogues, imcompleteness effects come in as the flux limit is approached (see Kondapally et al. 2022;Cochrane et al. 2023, for an analysis of the completeness for AGN and SFGs, respectively).More than 80,000 sources are catalogued in the central regions with the best multi-wavelength data (Paper III).As can be seen in Figure 1, LoTSS-Deep DR1 broadly matches the depth of the VLA-COSMOS 3GHz survey but over an order of magnitude larger sky area; similarly it matches the recent MeerKAT MIGHTEE Early Release (Heywood et al. 2022) in rms depth (the latter being limited by source confusion owing to its lower angular resolution), but again over larger area.

Multi-wavelength data in the LoTSS Deep Fields
ELAIS-N1, Boötes, Lockman Hole and NEP are the premier largearea northern extragalactic fields, with vast amounts of telescope time across the electromagnetic spectrum invested in observing these fields over the last two decades.Imaging at optical and near-IR wavelengths reaches 3-4 magnitudes deeper than typical all-sky surveys, allowing host galaxy identifications for over 97 per cent of the hosts of the radio sources in LoTSS-Deep DR1 (Paper III) compared to just 73 per cent using all-sky surveys in the LoTSS DR1 release (Williams et al. 2019).Other datasets, such as deep Herschel and Spitzer data in these fields, are irreplaceable, and add greatly to the scientific potential: Herschel data are a key tool to constrain obscured starformation rates, while the mid-IR wavelengths covered by Spitzer contain the diagnostic emission from the AGN torus.This range of complementary data makes these excellent fields to study not only the high-redshift AGN and luminous star-forming galaxies detected by LOFAR, but also to understand how this activity sits within the wider cosmological context of the underlying galaxy population.
As well as their combined benefit of sky area and sample size, each of the four LoTSS Deep Fields possesses unique characteristics or datasets which further enhance its specific scientific potential, whilst complementing each other.The specific data available in each field are summarised here; a more complete description of the available data in the ELAIS-N1, Lockman Hole and Boötes fields (but not NEP, as it was not included in the LoTSS-Deep DR1) can be found in Paper III, which also provides the coverage maps of each survey and the resulting catalogues.

ELAIS-N1
ELAIS-N1 has an ideal declination (+55 deg) for LOFAR observations, and is also a target field for LOFAR's Epoch of Reionisation studies (Jelić et al. 2014), providing a combined motivation for the observations.ELAIS-N1 benefits from some of the deepest widefield optical, near-IR and mid-IR imaging.It is one of the Medium Deep Fields from the Panoramic Survey Telescope and Rapid Response Sysytem (Pan-STARRS-1) survey (Chambers et al. 2016) Mauduit et al. 2012), which is around a magnitude deeper at 3.6 and 4.5m in the central 2.4 deg 2 .Longer wavelength data in the field have been taken using both Spitzer (24m data with the Multi-band Imaging Photometer for Spitzer; MIPS) and the Herschel Space Observatory, the latter as part of the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012), one of the deepest large-area Herschel surveys.HerMES observed ELAIS-N1 at 100m, 160m, 250m, 350m and 500m.

Boötes
The Boötes field is the target of some of the deepest wide-field optical imaging, in the   ,  and  filters from the NOAO Deep Wide Field Survey (Jannuzi & Dey 1999), in the -band from the zBoötes survey (Cool 2007), and in the  and  bands from the Large Binocular Telescope (Bian et al. 2013), all covering around 10 deg 2 .The same sky region has been observed in the near-IR ,  and  bands (Gonzalez et al. 2010) and using Spitzer from 3.6 to 8.0m as part of the Spitzer Deep Wide Field Survey (SDWFS; Ashby et al. 2009).Catalogues of galaxies in the Boötes field were generated by Brown et al. (2007Brown et al. ( , 2008)).Boötes has also been observed by Herschel as part of HerMES, and by Spitzer-MIPS, adding far-infrared measurements to the dataset.
In addition to this, Boötes benefits from excellent wide-field Xray coverage, including a deep Msec Chandra survey over the full 9.3 deg 2 field (Masini et al. 2020).The comparison between deep radio and deep X-ray observations opens many new scientific avenues, such as investigating the relationship between jet power and accretion rate in AGN, and determining the black hole accretion rates of starforming galaxies to investigate the co-evolution of galaxies and black holes.Boötes also possesses a vastly higher number of spectroscopic redshifts than the other northern deep fields, largely due to the AGN and Galaxy Evolution Survey (AGES; Kochanek et al. 2012): these are also very valuable for training photometric redshifts for the radio source population (e.g.Paper IV).

Lockman Hole
Located (like ELAIS-N1) at an ideal declination for LOFAR (+58 deg), the Lockman Hole is one of the regions of sky with the lowest Galactic HI column density (Lockman et al. 1986), making it ideal for extragalactic studies, especially at IR wavelengths due to its low IR background.For this reason, the Lockman Hole has been the target of some of the widest deep coverage in the optical to mid-IR bands.Optical data in the Lockman Hole has been taken by SpARCS in ,,, over 13.3 deg 2 , and by the Red Cluster Sequence Lensing Survey (RCSLenS; Hildebrandt et al. 2016)  The Lockman Hole is arguably the best-studied of the deep fields at other radio frequencies (e.g.Mahony et al. 2016;Prandoni et al. 2018;Morganti et al. 2021).The multi-frequency radio data allow detailed investigations of radio spectral shapes, identifying peaked, remnant and re-started sources, and giving a unique insight into the physics and lifecycles of radio-loud AGN (e.g.Brienza et al. 2017;Jurlin et al. 2020).

North Ecliptic Pole
The North Ecliptic Pole is an interesting field due to its location in the continuous viewing zone (CVZ) of many space telescopes, including the JWST, the eROSITA X-ray mission and Euclid.Until very recently, the multi-wavelength data quality in the NEP was inferior to the other three LoTSS Deep Fields, but this is rapidly changing.The NEP is the location of the Euclid Deep Field North which will provide deep subarcsecond near-IR imaging to depths of  = 26 over 10 deg 2 (and slightly shallower over a wider 20 deg 2 region).Such deep data will enable mass-complete samples to be defined down to ∼ 10 10  at  = 3 and normal star-forming galaxies to be detected out to  > 6.The combination of matched sub-arcsecond near-IR and radio continuum imaging (with LOFAR's international baselines) offers a unique opportunity to study the structural evolution of galaxies, for example comparing the spatial distribution of star formation (probed by LOFAR) versus stellar mass (probed by Euclid) within galaxies, to cleanly distinguish between different growth scenarios (e.g.'insideout' or 'outside-in' growth) over large samples of massive galaxies with  < 1.
Given these forthcoming datasets, a number of photometric surveys have been recently undertaken to provide matching observations at other wavelengths, including the Hawaii Two-0 survey (McPartland et al. 2023).Additionally, the Euclid/WFIRST Spitzer Legacy Survey has obtained mid-infrared imaging over the central 10 deg 2 of the field using Spitzer that is ∼0.8mag deeper than the SERVS data available in ELAIS-N1 and Lockman Hole.
As shown in Table 1, the NEP is not included in LoTSS-Deep DR1, and hence not included in the analysis of this paper, as the radio data were not available at the time of the optical cross-identification.An image from 72-hrs of data is now available and will be published by Bondi et al. (2023).Furthermore, as LOFAR observes two HBA pointings simultaneously, observations of the NEP field have included a parallel beam centred on the Abell 2255 cluster, which has also produced an ultra-deep low-frequency image of that field (Botteon et al. 2022).

Optical to mid-IR data
For the three fields presented in LoTSS-Deep DR1 (ELAIS-N1, Boötes, Lockman Hole), Paper III presented photometric catalogues from ultraviolet to far-infrared wavelengths.The reader is referred to that paper for a full description of the catalogues; here, a brief overview is provided.
For the ELAIS-N1 and Lockman Hole fields, data from UV through to mid-IR wavelengths were assembled and mosaicked on to a common pixel scale.Two combined  2 signal-to-noise images were then constructed, one by combining the optical to near-IR bands, and the other from the Spitzer 3.6 and 4.5m bands; these were treated separately due to the mis-match in angular resolution between the ground-based optical-to-near-IR and the Spitzer images.Forced aperture photometry was then performed across all bands using sources detected in each of these stacked images, and the two catalogues were merged to produce a single consistent photometric catalogue in each field.Aperture corrections were applied band-by-band based on curve-of-growth analysis for typical faint galaxies in order to provide total flux and total magnitude measurements.The photometry was corrected for galactic extinction based on the Milky Way E(B-V) extinction map of Schlegel et al. (1998) and the Milky Way dust extinction law of Fitzpatrick (1999).Uncertainties on the photometry were determined using the variations between a large number of apertures randomly placed around the fields.
For the Boötes field, forced aperture photometry catalogues already existed (Brown et al. 2007(Brown et al. , 2008) ) using magnitude-limited samples selected in the I-band and the 4.5m Spitzer band.In this case, these catalogues were used as the starting point, and were merged and corrected in a similar manner to ELAIS-N1 and Lockman Hole.In all three fields, the catalogues were then cleaned of low-significance detections (sources detected in the combined  2 image but below 3 significance in each individual band) and crosstalk artefacts, and those sources in regions around bright stars where either the cataloguing or the photometry might be unreliable were flagged, as indicated by the _ parameter.More details on all of these processes can be found in Paper III.
These photometric catalogues were then used as the basis for crossmatching with the LOFAR catalogues.Paper III outlines the selection of the studied area for which the highest-quality multi-wavelength data are available; sources within this region can be identified using the _ parameter.The cross-matching process also involved source association, such that the catalogued LOFAR sources were combined or deblended into true physical sources, where necessary.Within these defined areas, 81,951 physically distinct radio sources were catalogued over 25.65 deg 2 of sky across the three fields; optical or near-IR host galaxies were identified for over 97 per cent of these (Paper III), very much higher than the 73 per cent found for the wider LoTSS DR1 (Williams et al. 2019).
Photometric redshifts for all of the objects in the field have been presented in Paper IV.These were derived from the UV to mid-IR data by combining machine learning and template fitting approaches using a hierarchical Bayesian framework.This method is shown to provide photometric redshifts which are accurate for both galaxy populations (out to  ≈ 1.5) and sources dominated by AGN emission (out to  ≈ 4), which is important for the LOFAR sample.As part of the calibration of the photometric redshifts, small (typically <5 per cent) offsets in the zero-point magnitudes were found to improve the accuracy of the template-fit photometric redshifts.These offsets are discussed further in Section 3.3.

Far-infrared data
The addition of far-IR photometry is described by McCheyne et al. (2022), and the reader is referred to that paper for details.In summary, the far-IR fluxes were measured using XID+ (Hurley et al. 2017) which is a Bayesian tool to deblend the flux from the low resolution Herschel data into different potential host galaxies selected from optical/near-IR images.Fluxes were initially measured as part of the Herschel Extragalactic Legacy Project (HELP; Shirley et al. 2021).In HELP, an XID+ prior list of potential emitters at 24m was derived by applying a number of cuts to the optical-IR galaxy catalogue in order to select the sources most likely to be bright at 24m (those detected both at optical wavelengths and in the Spitzer 3.6-8.0mbands), and this input list was used to deblend the 24m data.Then, a second prior list was constructed from those sources with significant 24m emission (above 20Jy) and this was used to deblend the Herschel data.The posterior distributions for the fluxes derived from XID+ allow the uncertainties to be estimated.
For the LoTSS-Deep catalogue, a cross-match was first made between each LoTSS-Deep host galaxy position (or its LOFAR position if there was no host galaxy identification at optical-IR wavelengths) and the HELP catalogue.If a match was found then the HELP far-IR fluxes were assigned to the LOFAR source.If no match was found, then XID+ was re-run following the process above, but with the radio host galaxy position (or radio position in the case of no host galaxy identification) added to the prior list: this ensures that the assignment of zero flux is not simply due to the radio source having been incorrectly excluded from the prior list.

Final catalogues for spectral energy distribution fitting
In order to ensure consistency and reliability across the different spectral energy distribution (SED) fitting codes used in this paper, it was important to ensure that the input dataset was as robust as possible, and that all photometric errors were uniformly treated.
For each field, a catalogue was produced combining the (aperturecorrected and Galactic extinction corrected) fluxes from UV to mid-IR wavelengths with the far-IR fluxes determined by XID+.Next, the small zero-point magnitude corrections determined during the photometric redshift fitting were applied: these are tabulated in Appendix B of Paper IV. Specifically, the corrections derived using the extended Atlas library (referred to as 'Brown' in that paper) were applied; this template set was chosen because it extended out to the longest IRAC wavelength and also incorporated the full range of SED types expected within the LoTSS Deep Fields sample.
The photometry catalogue was then filtered to remove photometric measurements deemed to be seriously unreliable.These unreliable measurements were identified as those which were either 2.5 magnitudes lower, or 1 magnitude higher, than the value predicted by interpolating the two adjacent filter measurements.These limits were chosen, following Duncan et al. (2019), to avoid flagging any reasonable spectral emission or absorption features, or genuine breaks, while successfully identifying those measurements that are so discrepant that they could significantly influence the SED fitting.Around 1 per cent of the photometric measurements were identified in this way; these were flagged and not used in the subsequent fitting.
Finally, in order to consistently deal with any residual photometric errors due to zero-points, aperture corrections or extinction corrections, 10 per cent of the measured flux was added in quadrature to all flux uncertainties.The resultant SED input catalogues for each field are made available in electronic form through the LOFAR Surveys website (lofar-surveys.org).

Spectral Energy Distribution fitting
Many different codes exist for fitting SEDs to an array of photometric data points for galaxies and AGN.Each of these has their own advantages and disadvantages.Pacifici et al. (2023) recently carried out a detailed comparison of different codes, finding that they provide broad agreement in stellar masses, but with more discrepancies in the star formation rates and dust attenuations derived.In this paper, four different SED-fitting codes are adopted, and a comparison of the results between these is used both to derive consensus measurements for stellar masses and star-formation rates, and to assist with the classification of the radio source host galaxies.
The 'Multi-wavelength Analysis of Galaxy Physical Properties' ( ; da Cunha et al. 2008) and 'Bayesian Analysis of Galaxies for Physical Inference and Parameter EStimation' ( ; Carnall et al. 2018Carnall et al. , 2019) ) codes each use energy balance approaches to fit photometric points from the UV through to far-IR and sub-mm wavebands.Energy balance implies that the amount of energy absorbed by dust at optical and UV wavelengths is forced to match that emitted (thermally) by the dust through the sub-mm and far-IR.The and codes are built on the same fundamental templates for single stellar populations (Bruzual & Charlot 2003) but differ in their implementation, in particular with regard to the parameterisation of the star-formation histories of the galaxies, the assumed dust models, and the approach to model optimisation.For high signal-to-noise galaxies the two codes generally give broadly consistent results (see Sec. 5), which previous studies have generally shown to be accurate (e.g.Hayward & Smith 2015).However, neither nor includes AGN emission in its model SEDs, nor do they account for AGN heating effects when determining energy balance, and therefore both can give poor fits and unreliable host galaxy parameters for galaxies with significant AGN emission.
'Code Investigating GALaxy Emission' ( ; Burgarella et al. 2005;Noll et al. 2009;Boquien et al. 2019) is another broad-band SED-fitting code which uses energy conservation between the attenuated UV/optical emission and the re-emitted IR/sub-mm emission; differs from and in that it incorporates AGN models which can account for the direct AGN light contributions and the infrared emission arising from AGN heating of the dust (more recent developments also allow for predictions of X-ray emission, cf.Yang et al. 2020).The inclusion of AGN models can give a significant advantage over and when fitting the SEDs of galaxies that have a signficant AGN contribution, allowing both more robust estimation of host galaxy parameters, and a mechanism to identify and classify AGN within the sample.However, in order to allow the additional complications of AGN fitting, for equivalent (practical) run times is not able to cover the parameter space of host galaxy properties as finely as and , leading to potentially less accurate characterisation of galaxies that do not host AGN.
All of the three codes discussed above adopt the principles of energy balance.However, if the distribution of ultraviolet light is spatially disconnected from the dust emission, as is often the case for very infrared luminous galaxies, then energy balance may not be valid; indeed, Buat et al. (2019) find for a sample of 17 well-studied dust-rich galaxies that SED-based UV-optical attenuation estimates account for less than half of the detected dust emission.This issue may be particularly pronounced in the presence of AGN, if the AGN models are not comprehensive enough to properly cover the parameter space of possible AGN SEDs.To mitigate these issues, the code (Calistro Rivera et al. 2016) models the SED by independently fitting four emission components, with each independently normalised (albeit with a prior that the energy radiated in the infrared must be at least equal to the starlight energy absorbed by dust at optical/UV wavelengths): a big blue bump, a stellar population, hot dust emission from an AGN torus, and colder dust emission.
can provide superior fits for objects where energy balance breaks down, and also for objects with strong AGN components due to its superior modelling of the big blue bump.However, the lack of energy balance and the ability of the four components to vary independently can lead to aphysical solutions, or poorer constraints on the parameters of the stellar populations (although Gao et al. 2021 find broadly good agreement in measured stellar masses and SFRs between codes with and without energy balance, at least for hyperluminous infrared galaxies).
To maximise the advantages of the different techniques, the LoTSS Deep Field host galaxies were all modelled using each of , , and .Furthermore, for , two different sets of AGN models were considered: those of Fritz et al. (2006) and those of Stalevski et al. (2012Stalevski et al. ( , 2016)), the latter of which were recently incorporated into by Yang et al. (2020).The following subsections provide details of the fitting methodology in each case.
For all SED fitting, the redshift of the source is fixed at the spectroscopic redshift,  spec , for the minority of sources for which this exists (1602, 4039 and 1466 sources in ELAIS-N1, Boötes and Lockman, respectively).For the other sources, the redshift is fixed at the median of the first photometric redshift solution,  1,median .Photometric redshift errors may introduce errors on the inferred parameters, but for most sources these are anticipated to be small since the photometric redshifts are very accurate, with a median scatter of Δ/(1 + ) ∼ < 0.015 for host-galaxy dominated sources at  < 1.5 (Duncan et al. 2021).

3.4.1
The application of to the LoTSS Deep Fields sources is described by Smith et al. (2021), and so it is only briefly summarised here.The stellar population modelling adopts single stellar population (SSP) templates from Bruzual & Charlot (2003) and the two-component (birth cloud plus interstellar medium) dust absorption model of Charlot & Fall (2000), combined to produce an optical to near-IR template library of 50,000 SEDs with a range of exponentially-declining star-formation histories with stochastic bursts superposed.The dust emission is modelled using a library of 50,000 dust SEDs constructed from dust grains with a realistic range of sizes and temperatures, including polycyclic aromatic hydrocarbons.The energy balance criterion is used to combine the two sets of templates in a physically-viable manner, to produce a model for the input photometry that stretches from near-UV to sub-mm wavelengths.
determines the best-fitting SED for every source, returning the corresponding best-fit physical parameters and their marginalised probability distribution functions (PDFs).The bestfitting stellar mass and best-fitting value of the SFR over the last 100 Myr were adopted as the stellar mass and SFR respectively; the 100 Myr timescale corresponds well to that of the expected radio emission (e.g.Condon et al. 2002).For most galaxies, very similar results are obtained if a shorter period or the current instantaneous SFR are adopted instead (although results for some individual galaxies can vary significantly).The 16th and 84th percentile of the PDFs were adopted as the 1 lower and upper limits respectively.In order to determine whether the calculated parameters are reliable, the  2 value of the fit was examined: following Smith et al. (2012), fits for which the determined  2 value was above the 99 per cent confidence limit for the relevant number of photometric bands included in the fit were flagged as unreliable.As noted by Smith et al. (2021), many of the objects that fail this test are objects with strong AGN contributions.On average, 17 per cent of sources across the three fields were flagged in this way, with ELAIS-N1 giving a significantly lower fraction (10 per cent), in line with expectations that the deeper radio data in that field should result in a higher fraction of star-forming galaxies.

3.4.2
was run on the LoTSS Deep Field sources, making use of the 2016 version of the Bruzual & Charlot (2003) SSP templates for its stellar population emission.Nebular emission is computed using the photoionization code (Ferland et al. 2017), following Byler et al. (2017).
is run using each SSP template as the input spectrum.Dust grains are included using 's 'ISM' prescription, which implements a grain-size distribution and abundance pattern that reproduces the observed extinction properties for the Interstellar Medium (ISM) of the Milky Way.A Calzetti et al. (2000) dust attenuation curve is adopted.Dust emission includes both a hot dust component from HII regions and a grey body component from the cold, diffuse dust.
A wide dust attenuation prior is adopted,   = [0, 6], which gives the code the option to fit a high degree of attenuation.The absorbed energy is re-emitted at infrared wavelengths; the dust SED is controlled by three key parameters, as described by Draine & Li (2007):  min , the lower limit of the starlight intensity; , the fraction of stars at  min ; and  PAH , the mass fraction of polycyclic aromatic hydrocarbons.The priors adopted on these parameters are broad, to allow the model to fit all types of galaxies, including those that are hot and dusty (Leja et al. 2018):  min = [0, 25],  = [0, 1], and  PAH = [0, 10]., the multiplicative factor on   for stars in birth clouds, is also fitted using the prior  = [1, 5].Metallicity is allowed to vary in the range  = [0, 2.5] ,old , where  ,old denotes solar models prior to Asplund et al. (2009).
The star-formation history (SFH) is parameterised using a double power law: SFR(t) ∝ [(t/)  + (t/) − ] −1 where  is the slope in the region of falling SFR, and  is the slope in the region of rising SFR. relates to the time at which the SFR peaks.
The code outputs posterior distributions for the fitted parameters  V ,  min , , and  PAH , , the metallicity , and the SFH parameters ,  and .Posterior distributions are also derived for the physical properties of stellar mass, star-formation rate, and specific star-formation rate, with the median and the 16th and 84th percentiles being adopted as the best-fit value and the lower and upper 1 errors.The reduced  2 of the best-fitting model was also returned.Objects with a reduced  2 above 5 were flagged as unreliable; this averaged about 9 per cent of sources across the three fields, again being lowest in ELAIS-N1 and highest in Boötes.

3.4.3
was run on the LoTSS Deep Fields sources in the manner outlined in Wang et al. (2021) and Małek et al. (2023).The choices for the input components for the modelling of the stellar population largely follow those of Pearson et al. (2018) and Małek et al. (2018).Specifically, the star-formation history was adopted to be a twocomponent model, with a delayed exponentially-decaying main starforming component (SFR delayed ∝  −/ ) plus the addition of a recent starburst.The Bruzual & Charlot (2003) SSP templates were adopted for the stellar emission.The Charlot & Fall (2000) dust attenuation model is applied to the derived SEDs, and energy-balance criteria are used to determine the quantity of emission to be reemitted in the infrared.The dust emission is calculated using the dust emission model of Draine et al. (2014), which is an updated version of the Draine & Li (2007) model and describes the dust as a mixture of carbonaceous and amorphous silicate grains.
A critical difference between and / is the inclusion of an AGN component in the models.For the LoTSS Deep Fields, was run twice, using two different AGN models: the Fritz et al. (2006) model and the model of Stalevski et al. (2012Stalevski et al. ( , 2016)).Both sets of AGN models assume point-like isotropic emission from a central source, which then intercepts a toroidal dusty structure close to the AGN.Radiative transfer models are used to trace the absorption and scattering of the AGN light by the dust in the torus, and model its re-radiation by the hot dust.The main differences between the two models are that the Fritz models adopt a smooth density distribution for the dust grains and use a 1-D approach, whereas the models treat the dusty torus as a two-phase medium with higher density clumps sitting within a lower density medium and use 3-D radiative transfer.A clumpy dust distribution was suggested by Krolik & Begelman (1988) to be necessary to stop the dust grains being destroyed by the hot surrounding gas.
returns Bayesian estimates of the stellar mass and various estimates of the recent star-formation rate of the galaxy, along with estimates of the uncertainties on these parameters.In this work, the star-formation rate averaged over the last 100 Myr is adopted, as for .
also returns a determination of the AGN fraction for the galaxy (hereafter  AGN,CG−F or  AGN,CG−S for the Fritz and models), defined as the fraction of the total infrared luminosity that is contributed by the AGN dust torus component.An uncertainty on the AGN fraction is also returned; where this is larger than the measured fraction, the 1-sigma lower limit on the AGN fraction is set to zero.Finally, the reduced 2 of the best-fitting model was used to identify unreliable fits, with objects with a reduced  2 above 5 being flagged (3 per cent and 2 per cent of sources in the Fritz and models respectively).

3.4.4
provides independent parameterisations for each of the accretion disk emission (big blue bump), the hot dust torus, the stellar component and the cooler dust heated by star formation; details of the parameterisation of these four components are provided by Calistro Rivera et al. (2016).
accounts for the effects of reddening on these emission components but without energy balance constraints.
was run on the LoTSS Deep Fields sources broadly following the implementation of Williams et al. (2018) but using an expanded set of input models ( 2; Calistro Rivera et al., in prep.).The code determines the relative importance of the four components in a few key wavelength regions, as well as broader physical parameters including estimates of the star-formation rate and the stellar mass.In this work, the IR-based estimate of the SFR was the one adopted.
Following Williams et al. (2018), an AGN fraction is defined by considering the contribution of the emission components in the 1-30m wavelength range.Note that this is different to the definition used for which considers the AGN contribution to the total IR luminosity: as the AGN peaks in the mid-IR, the AGN fractions derived by will typically be larger than those of .The AGN fraction was defined as: where  Torus,1−30 ,  SB,1−30 and  Gal,1−30 are the luminosities of the hot dust torus, the cooler dust heated by recent star formation, and the stellar component of the galaxy, respectively, all between 1 and 30m.Note that this differs slightly from the definition of Williams et al. (2018) through the inclusion of the stellar component in the denominator; this avoids a high AGN fraction being determined when the mid-infrared emission is simply dominated by the light of older stars.The uncertainties on these luminosities are used to determine the 1 upper and lower limits to the AGN fraction.Finally, returns a log likelihood for the best-fit model; the ≈ 3 per cent of objects whose fits had a log likelihood below −30 were flagged as unreliable (cf.Williams et al. 2018).

IDENTIFICATION OF RADIATIVE-MODE AGN
A characteristic feature of radiative-mode AGN is a hot accretion disk, which is being obscured in certain directions by a dusty structure (the torus).These two structures give rise to a variety of physical features that can be used to identify the radiative-mode AGN.The most widely-used of these, where spectroscopic data is available, is emission line ratios (e.g.Baldwin et al. 1981, the BPT diagram): the ionising radiation from the hot accretion disk is significantly harder than that of a young stellar population, leading to stronger highexcitation forbidden lines.Spectroscopic information is available for only a small subset of the LoTSS-Deep sources (5.1, 21.1 and 4.7 per cent in ELAIS-N1, Boötes and Lockman Hole respectively, with the AGES data in Boötes producing the large difference between the fields), so this method cannot be used for the vast majority of the sources.This will change in the coming years due to the WEAVE-LOFAR survey (Smith et al. 2016, see also Sec. 9) but alternative methods are needed for AGN identification in the meantime.
The hot dusty torus emits characteristic emission that has been widely used to identify radiative-mode AGN using mid-IR colours (e.g.Lacy et al. 2004;Stern et al. 2005).Commonly-used selections consider the four Spitzer channels centred at 3.6m, 4.5m, 5.8m and 8.0m (Channels 1 to 4 respectively); the selection is based on the premise that the emission from stellar populations generally declines with increasing wavelength through the mid-IR (since the mid-IR probes redward of the rest-frame 1.6m thermal peak of the dominant sub-solar stellar population) whereas hot AGN dust shows a rising spectrum.An equivalent approach uses the WISE mid-infrared colours (e.g.Wright et al. 2010).The exact colourspace cuts are generally defined using template tracks for galaxies and AGN to select regions of colour-space dominated by AGN.Lacy et al. (2004) and Stern et al. (2005) derived the first colourcuts based on shallow Spitzer data (hereafter referred to as the Lacy and Stern regions, respectively), and these were effective in separating out AGN from the population of relatively nearby inactive galaxies.However, the broad colour regions selected in these papers are heavily contaminated by higher redshift ( > 0.5) inactive galaxies, that deeper Spitzer surveys (such as those available in the LoTSS Deep Fields) are able to detect.Donley et al. (2012) therefore defined a much tighter region of mid-IR colour space (hereafter, the Donley region) within which AGN samples display much lower contamination, but consequently are also less complete.Even in these deep datasets, however, fainter galaxies often lack measurements in one or more channels, preventing any classification by the Stern, Lacy or Donley criteria.To help overcome this, Messias et al. (2012) derived a series of redshift-dependent colour cuts based on K-band to Channel 2, Channel 2 to Channel 4, or Channel 4 to 24m flux ratios (hereafter, the Messias regions).These allow classification of a larger fraction of galaxies, but with the same issues regarding completeness and contamination.Furthermore, simple application of colour cuts takes no account of low signal-to-noise measurements which can scatter data across the colour criteria, and can also miss some types of AGN (e.g.Gürkan et al. 2014).
The wide array of data available in the LoTSS Deep Fields allows a classification scheme to be developed which uses much more than just the mid-IR colour bands.The SED fitting described in the previous section encodes all of the mid-IR spectral expectations used in the Stern, Lacy, Donley and Messias colour criteria, but combines this with additional near-IR and optical data which allow simultaneous characterisation of the host galaxy properties; the latter allows the contribution of the host galaxy to the mid-IR to be directly predicted, and thus any additional AGN contribution to be more clearly distinguished.As an indication of this, Figure 2    AGN fraction values (e.g.0.05, 0.7) appears to be a feature of the code, perhaps due to the fairly limited sampling of the grid of AGN model parameters.The plot shows that as the AGN fraction rises above ∼0.1,objects are more likely to be identified as AGN through spectroscopic or X-ray selection or the Donley mid-IR cuts, and also that the SED fitting begins to deteriorate (higher relative  2 ) for SED codes that don't include AGN components.
through optical spectra or X-ray properties are indicated in red.It can be seen that the X-ray and spectroscopically selected AGN and the objects with high AGN fractions concentrate primarily in the selected colour-space regions, especially the Donley region, but that a significant fraction of these probable AGN are also found outside of these regions.Furthermore, there are objects within the to demonstrate the results, as this is the deepest field with the best multiwavelength data.In all cases, all three deep fields show consistent results.colour-cuts (especially the broader Lacy and Stern regions) for which predicts very low AGN contributions to the mid-IR.
The use of the four SED fitting routines provides two routes to identifying the probable AGN.First, each of and provides an estimate of  AGN , the fractional AGN contribution to the mid-IR.Second, objects which have a significant AGN contribution to their SED should be poorly fitted using or (and typically better fitted using or ). Figure 3 demonstrates these effects, by showing the AGN fraction plotted against the ratio of the  2 values determined from the SED fits without AGN components compared to those with AGN components, with points colour-coded by evidence for AGN from either spectroscopic or X-ray data, or from mid-IR colour cuts.The spectroscopic and X-ray selected AGN generally show both moderate-tohigh AGN fractions and a higher  2 using / than using / .The majority of objects which lie securely within the Donley mid-IR colour-cuts show the same characteristics.Objects that lie only within the broader Stern, Lacy or Messias colour regions typically show much lower AGN fractions and the  2 value from the / fits is lower than or comparable to that from / ; they largely overlap with the 'non-AGN' that either lie outside of these colour cuts or do not have sufficiently high signal-to-noise in their mid-IR measurements for this to be determined.Nevertheless, the SED fits are able to pick out promising AGN candidates within these categories.
An examination of the AGN fractions derived by and especially by shows that many of these have quite large uncertainties, especially for fainter galaxies with fewer securely-measured photometric points.Investigations indicated that the 16 th percentile of the posterior of the AGN fraction (i.e. the 1-sigma lower limit on the AGN fraction; hereafter P16) provided a more robust indication of the presence of an AGN.The selection of radiative-mode AGN was therefore made by considering three selection criteria (see below for a discussion of how the threshold values were set): (i) whether the P16 AGN fraction from , using the AGN models, exceeded a threshold value of 0.06 (ELAIS-N1 and Lockman Hole fields) or 0.10 (Boötes field).
(ii) whether the P16 value for the AGN fraction from , as defined in Eq. 1, exceeded a threshold value of 0.15 (ELAIS-N1 and Lockman Hole fields) or 0.25 (Boötes field).
(iii) if the lower of the reduced  2 values arising from the and SED fits was both greater than unity and at least a factor  greater than the lowest of the reduced  2 values arising from the two and the SED fits.The factor  was determined to be twice the median value of the  2 ratio between the better fit from and and the best fit from and (cf. Figure 4).This evaluated to  = 1.36 for ELAIS-N1,  = 1.59 for Lockman Hole and  = 2.22 for Boötes.
An object was classified as a radiative-mode AGN if it satisfied at least two of these three criteria.In practice, this means either that it has a determined high AGN fraction from both and or it has a high AGN fraction from at least one of the two codes combined with a superior SED fit using methods which include AGN components.The selection cuts for each criterion were set by comparing the derived classifications with the spectroscopic and Xray samples and considering the locations of the classified AGN and non-AGN on mid-IR colour-colour diagrams.The threshold values selected were different for Boötes than for the other two fields.This is because the AGN fractions calculated in that field were systematically higher than those in ELAIS-N1 or Lockman Hole (e.g. a median AGN fraction of 0.037 in Boötes using the model, compared to 0.029 in each of ELAIS-N1 and Lockman), which is likely to be due to the different manner in which the photometric catalogues were constructed in Boötes (see Paper III).Setting higher thresholds in Boötes ensured a consistency of classification across the three fields (cf.Sec. 7).Finally, a small proportion of objects did not meet these criteria but had previously been identified to be an AGN based on either optical spectra or X-ray properties; these were added to the radiative-mode AGN sample (and correspond to about 3 per cent of all radiative-mode AGN).
Fig. 4 shows the LoTSS-Deep sources on different combinations of these selection criteria, with the sources that satisfy at least two criteria, and therefore are selected as radiative-mode AGN, shown in red.It can be seen that there is a broad consistency between the different criteria: most of the selected radiative-mode AGN satisfy all three criteria and therefore are secure classifications.The main addition to this is a population of sources selected as having high AGN fractions by both and but with comparable, low  2 values from the different fitting methods; these are probably sources where and are able to pick out a weak AGN through the mid-IR emission, but there is little-to-no direct AGN light through the optical to near-IR spectrum and so and are still able to provide a good fit to the majority of the spectrum.
Fig. 5 shows the selected radiative-mode AGN and non-AGN on a series of mid-IR colour-colour diagrams, compared against the evolving colours of various galaxy template models.The panels are split by redshift ranges, in order to allow a clearer comparison against the template expectations.At each redshift, the panels show the Lacy and Donley colour plots (left), the Stern colour plot (middle), and the appropriate Messias plot (right).Template SED models were drawn from the 'Galaxy SED Atlas' of Brown et al. (2014) combined with the 'AGN SED Atlas' of Brown et al. (2019).SEDs were selected from these libraries for: (i) elliptical galaxies (as expected to be seen for jet-mode AGN); (ii) star-forming galaxies; (iii) AGN (including both quasars and edge-on 'type-II' AGN); and (iv) composite spectra, produced by combining a set of Seyfert AGN spectra with host galaxy spectra, with a range of weights.
The template tracks for the different galaxy classes confirm both the motivation for, and the shortcomings of, the colour-colour selection criteria: the Donley region relatively cleanly selects AGN at  < 2.5 but is incomplete for composite systems; the Stern and Lacy regions are more complete for composite systems but contaminated, especially at the higher redshifts; the Messias cuts perform relatively well, especially at the highest redshift where the use of the 24m colour gives a clear advantage, but still have some incompleteness and contamination.The red points show the objects selected as radiative-mode AGN by the techniques outlined above.At all redshifts these broadly overlap the regions of the AGN and composite templates, extending where appropriate beyond the colour-selection limits.It is clear, however, that in the  > 2.5 redshift range there remains a significant population of objects that are not classified as AGN, and yet which lie in similar regions of colour-space to the AGN.At these redshifts, as is evident from Fig. 5, it is only the Channel 4 and 24m filters that are able to probe rest-frame wavelengths where an AGN template becomes clearly distinct from the galaxy templates, and the composites are even more difficult to distinguish.Especially with the typically low signal-to-noise of the galaxies in this highest redshift bin, the SED fitting techniques may be less reliable: although the classifications are provided for all sources, readers should treat these with caution at  > 2.5, where there may well be a degree of incompleteness in the AGN sample.

COMPARISON OF DERIVED PROPERTIES AND CONSENSUS MEASUREMENTS
Two of the most important galaxy properties to determine are the stellar mass and the star-formation rate.Each of the SED fitting codes provides an estimate of these parameters.This section discusses how these values are combined to produce consensus measurements for each source.
In brief summary, for sources which do not host an AGN, the and codes ought to provide the best measurements of mass and SFR, because these models offer a significantly broader selection of galaxy templates.Indeed, for these sources, the results from these two codes show excellent agreement in their estimates of both stellar mass (median absolute difference of just 0.09 dex) and SFR (0.14 dex).The consensus values of the stellar mass and SFR for non-AGN were therefore generally derived from the logarithmic mean of the and results.
For radiative-mode AGN, the and results are potentially unreliable as they do not include any AGN component in their SED modelling.The two runs (with the Fritz and AGN models) should be more reliable, and indeed these two agree with each other well: the median absolute difference is only 0.09 dex in stellar mass and 0.13 dex in SFR.
is found to provide less consistent results, but is valuable for the small fraction (≈ 2 per cent) of sources which are highly AGN-dominated, and for which 's superior modelling of the AGN UV emission is required.
The consensus values of the stellar mass and SFR for radiative-mode AGN were therefore typically derived from the logarithmic mean of the two results, except where failed to provide an acceptable fit, in which case the values were adopted.Sections 5.1 and 5.2 now provide (for stellar mass and SFR respectively) a much more detailed comparison of the outputs of the different SED fitting codes, along with a full description of how the generalised approach discussed above was adapted in cases where one or more of the SED codes failed to provide an acceptable fit.Readers not interested in these finer details may wish to skip to Section 6.

Consensus stellar masses
For sources which are not identified to be a radiative-mode AGN, the results from the and codes show excellent agreement in their estimates of stellar mass: where both  2019), separated into ellipticals (pink), star-forming galaxies (yellow), AGN (purple), and composites (green).As can be seen, the broad colour cuts suffer to various extents from both incompletensss and contamination.The selected AGN broadly align with the regions of colour space covered by the AGN and composite template spectra.and pass the threshold for an acceptable fit (see Section 3.4) the median absolute difference in stellar mass is just 0.09 dex, with over 90 per cent of sources agreeing within 0.25 dex; the outliers are generally the faintest sources, at low masses or high redshifts.also gives very similar values, with a median difference in stellar mass of only 0.11 dex, and over 85 per cent agreeing within 0.25 dex.
shows much lower agreement, however, with a median difference in stellar mass of 0.27 dex compared to the estimates from the other codes.This inconsistency for is likely to be associated with the lack of an energy balance in the fitting process.
For these non-AGN the consensus stellar mass was derived from the mean of the logarithm of the stellar masses derived using and , as long as both codes provided an acceptable fit to the data (≈ 86 per cent of the non-AGN, though rising to nearly 95 per cent in ELAIS-N1).If one of the two codes provided a bad fit and the other a good fit (11 per cent of cases), then the stellar mass estimate from the well-fitting code was adopted as the consensus measurement.If both codes produced fits below the acceptability threshold then the values of the two stellar mass estimates were examined: if they agreed with each other within 0.3 dex (≈ 2 per cent of cases) then it was likely that the unreliability of the SED fits was driven by some outlier points that did not invalidate the stellar mass estimates, and so the logarithmic mean of the two values was adopted as the consensus stellar mass.If the two values disagreed by more than 0.3 dex, then the stellar mass estimates of the two fits were examined as well: if the full range of all 4 stellar masses was less than 0.6 dex (≈ 0.3 per cent of cases) then the logarithmic mean of the four measurements was adopted as the consensus measurement; if the range was larger than 0.6 dex (≈ 0.6 per cent of sources) then it was deemed that no reliable stellar mass could be provided.A comparison of the consensus masses derived against the estimates from each code individually is shown by the black points in Fig. 6, confirming visually the good agreement of the and codes, broad agreement of , and larger scatter of for these sources.
For radiative-mode AGN, the two runs provide stellar mass estimates that agree well with each other: the median absolute difference is only 0.09 dex, with 90 per cent of sources within 0.3 dex.Compared to these values, as expected, the results from and show greater scatter (each 0.16 dex median difference) and also a larger fraction of outliers where the codes significantly over-estimate the mass due to AGN light being incorrectly modelled as stellar emission (cf.Fig. 6).Again, shows a larger dispersion in stellar mass measurements relative to the other codes, with a median absolute difference of 0.49 dex; this may be due to the stellar component being fitted independently without an energy balance constraint, with some stellar light perhaps being incorrectly modelled as AGN emission or vice versa, although it could also be related to the different approach to modelling the AGN emission.For these reasons, for the radiative-mode AGN, if both runs provided acceptable fits then the logarithmic mean of the stellar masses from these two runs was accepted as the consensus mass (with excluded due to its higher proportion of outliers); this was the case for just over 94 per cent of the radiative-mode AGN.Otherwise, if just one of the runs provided an acceptable fit (≈ 3 per cent of cases) then the stellar mass from that run was adopted.If neither run provided a good fit, but did, then there was a likelihood that this was a case where either energy balance was breaking down or the superior modelling of the AGN UV emission by was helping the fit; in these 2 per cent of cases, the stellar mass estimate was used.Otherwise, it was decided that no reliable stellar mass estimate was possible.
Fig. 6 shows a comparison on each mass estimate against the consensus mass derived, and illustrates the trends discussed above.The lower-right panel also compares the consensus masses against those derived in Paper IV using a grid-based SED fitting mechanism (see also Duncan et al. 2019).This comparison is interesting because the stellar masses in Paper IV are derived for all galaxies in the field, not only the radio sources, and therefore allow a comparison between the radio sources and the underlying population.In Paper IV it is argued that the stellar mass estimates are only reliable out to  ∼ 1.5, and so this is set as an upper limit for the plotted points.As can be seen, the agreement between the Paper IV stellar masses and the consensus masses derived here is very good for the non-AGN, with no significant systematic offset (< 0.1 dex) and a median scatter of 0.11 dex.The performance for AGN is slightly worse, but still good, with a median scatter of 0.23 dex.These results confirm that the Paper IV masses provide reliable measurements for the broader population that can be used in comparison against the consensus masses for the radio source population.
In this paper, no attempt is made to derive uncertainties on the consensus stellar masses for individual sources.Uncertainties arise both due to statistical errors in the individual fits and systematic effects between different SED codes.Each SED code offers an estimate of its statistical uncertainty for each source, and the difference between the stellar masses from different SED codes can be used to gauge the size of the systematic errors.Another source of error is that during the SED fitting the redshift of the source is fixed at the best photometric redshift (unless a spectroscopic redshift is available): uncertainties in the photometric redshift are likely to be a significant contributor to the mass uncertainty for any given source.Instead of calculating uncertainties for individual sources, therefore, the approach taken here is to derive characteristic uncertainties on stellar mass as a function of the galaxy's mass and redshift.The characteristic uncertainties are evaluated in Appendix A, and are found to be typically around 0.1 dex for higher mass sources at  < 2, increasing towards higher redshifts and lower masses.

Consensus SFRs
Estimation of consensus SFRs follows broadly the same principles as those of the stellar masses, in the preferred use of the and results for the non-AGN and with the results generally used for the AGN.As would be expected (cf.Pacifici et al. 2023), the agreement in SFR estimates between the different codes is not quite as good as that of stellar masses, but still strong.For non-AGN, the SFR estimates of and show systematic differences of less than 0.1 dex, with a median scatter of only 0.14 dex and over 75 per cent of cases agreeing within 0.3 dex.The measurements agree comparably well at large SFRs, but frequently provide higher SFR estimates than either or at lower SFRs.
suffers from a significant systematic offset of, on average, more than 0.3 dex higher SFRs than the other estimators.For the radiative-mode AGN, the two SFR estimations show good agreement with each other (median difference 0.13 dex).Both and systematically over-estimate the SFRs of these radiative-mode AGN, by around 0.15 dex on average.Fig. 7 provides a visual illustration of these effects.
To determine the consensus SFRs, like for stellar masses, the outputs from and are primarily considered for the non-AGN.The only significant difference in approach arises because of a small proportion of sources (around 9 per cent of all the non-AGN sources, mostly at lower SFRs) for which returns an acceptable fit, but the SFR is dramatically below that of and with an uncertainty that can be several orders of magnitude larger than the estimated value.These very low SFRs arise because of the parametric (exponentially-declining) form of the SFR history, which can lead to unrealistically-low best-fit SFRs at large ages where the e-folding time is short, but with considerable uncertainty.For these sources, the SFR estimates are found to broadly agree with the values, with both often within the 1 confidence interval of the fit.Therefore, sources for which the fit is deemed to be good, but the uncertainty on the SFR estimate is more than 5 times the estimate itself, are treated differently.In these cases, if provides an , and all give broadly consistent results for non-AGN, but differ for the AGN subset, for which the results should be more reliable.
masses show a small systematic offset compared to the other codes, and more outliers at high mass.The lower right plot examines the masses produced in Paper IV (only out to  < 1.5); these are seen to give consistent results with only slightly larger scatter.This is of interest because these stellar masses were produced for the entire galaxy population in these deep fields, not just the LoTSS-Deep host galaxies.and give broadly consistent results for non-AGN; their performance on objects identified as (radiative-mode) AGN is more mixed, but generally reasonable where the fit is not flagged as a bad fit.
's SFR estimations for non-AGN generally perform well at higher SFRs (especially with the AGN model), but over-predict the SFR in some lower-SFR galaxies.The estimated SFRs of objects selected as AGN show a high degree of consistency between the two different runs.
SFRs show more scatter and a small systematic offset compared to the other codes.
acceptable fit then the estimate is adopted as the consensus value; if it does not, but the and estimates agree within 0.5 dex then the logarithmic mean of the and values is taken as the consensus value; otherwise, the results are deemed inconsistent and no consensus SFR is derived.Other than these cases, the approach to derive consensus SFRs for the non-AGN exactly matches that for deriving stellar masses.Similarly, for the radiative-mode AGN, the approach for stellar masses using (or occasionally ) estimates is replicated for the SFRs.Fig. 7 compares the consensus SFRs against the estimates from each individual code.The spread in derived values between different codes is comparable to that in the analysis of Pacifici et al. (2023).As with stellar masses, no attempt is made to provide a source-bysource uncertainty on the consensus SFR, but Appendix A discusses the typical errors; except for the few per cent of lowest-SFR objects at each redshift (where the uncertainties increase greatly), these can be broadly approximated as Δ(SFR) ≈ 0.1(1 + ) 0.5 dex.

IDENTIFICATION OF RADIO AGN
As discussed in the introduction, star-forming galaxies show a tight correlation between their radio luminosity and their SFR2 .This relation allows the identification of sources which possess significant radio emission associated with AGN activity, as they will appear offset to larger radio luminosities than would be predicted from their SFR (cf.Delvecchio et al. 2017;Williams et al. 2018;Whittam et al. 2022).Relationships between SFR and low frequency radio luminosity have been previously derived at relatively low redshifts by Calistro Rivera et al. (2017), Brown et al. (2017), Gürkan et al. (2018) and Wang et al. (2019), and most recently by Smith et al. (2021) using the LoTSS-Deep data in ELAIS-N1.As discussed by Smith et al., in order to determine an accurate relation it is essential to properly account for non-detections, otherwise there is a risk that the derived relation will be dependent on the depth of the radio imaging, with the bias decreasing as the depth of the radio imaging increases.Smith et al. derive their relationship out to  ≈ 1 using a near-IR magnitude selected sample, finding log 10  150MHz = 22.22 + 1.06 log 10 (SFR) for the sample as a whole (where  150MHz is in units of W Hz −1 and SFR in units of  yr −1 ), based on SFRs derived using .In this paper, the use of the consensus SFRs, and the extension to higher redshifts, may be expected to lead to small changes in the best-fit relation.A suitable relation is therefore derived using a 'ridgeline' approach.In this approach, the sources are binned into different (narrow) bins in SFR, and within each bin the distribution of radio luminosities of the detected sources is examined.The peak of the distribution is identified as the ridgeline point.Provided the radio survey is sufficiently deep then, especially in the presence of a distorted distribution (the star-forming population plus a distribution of radio-excess AGN), this method should provide a more reliable value than the mean or median of the distribution of detected sources.The radio luminosities and SFRs of the LoTSS-Deep sources are shown in the upper panel of Fig. 8, along with the calculated ridgeline points, which can be well-fitted by the relation log 10 ( 150MHz /W Hz −1 ) = 22.24 + 1.08 log 10 (SFR/M yr −1 ) (2) The uncertainty on the ridgeline gradient is ±0.06, and the uncertainty on the intercept at log 10 (SFR) = 1.5 (the median value, where the errors on the gradient and intercept are uncorrelated) is ±0.07.To within 1, there is no difference in this relation between those sources classified as radiative-mode AGN or not.The relation derived from the ridgeline is fully consistent with that of Smith et al. (2021), agreeing within 0.1 dex over the full range of star-formation rates probed.
The distribution of radio luminosities below the ridgeline can be reasonably well-fitted by a Gaussian distribution of width 0.22 dex; this also holds in different bins of star-formation rate, with the Gaussian width remaining constant (to ±0.02 dex) from low to high SFR.The distribution above the ridgeline shows a much more extended tail, as expected.In ELAIS-N1 and Lockman Hole, radio-excess sources are here defined as those sources with radio luminosities exceeding the ridgeline value by 0.7 dex, corresponding to approximately 3.It should be noted that this limit corresponds to approximately 0.8 dex above the relation of Smith et al. (2021) at high SFR; these authors derived a scatter in their relation of around 0.3 dex at SFR > 10 yr −1 (at lower SFRs they measured lower scatter, but noted that this might be due to the limiting depth of the radio imaging); Cochrane et al. (2023) also derive a similar value for the scatter.Therefore, the radio-excess selection adopted here also broadly corresponds to a 3 excess relative to the Smith et al. relation.In Boötes (where the input photometry was different), it is found that the scatter in the SFR-radio relation increases towards higher redshifts, and adoption of a fixed 0.7 dex cut-off leads to an excess of radio-AGN at higher redshifts compared to the other two fields.To remedy this, in Boötes the radio excess threshold is modified slightly to (0.7 + 0.1) dex, which brings the classifications in this field in line with those in ELAIS-N1 and Lockman (cf. Fig 9).
There is a small population of radio sources with consensus SFRs well below 0.01 yr −1 .SFRs at this level cannot be accurately estimated by the SED fitting codes, and thus have large associated uncertainties.This makes a radio-excess classification based on the consensus SFR potentially unreliable for these sources.To avoid this issue, these sources were only classified as radio-excess if their radio luminosity exceeded (by 0.7 dex) that expected for a SFR of 0.01 yr −1 .If their radio luminosity was below that level, but above the radio-excess limit for their estimated consensus SFR, they were deemed to be unclassifiable in terms of radio excess (0.4 per cent of sources).
Finally, a small proportion of sources do not reach the radio-excess selection threshold, but are clearly extended or multi-component radio sources, inconsistent with simply being star-forming galaxies.Those sources which are either multi-component sources associated through the LOFAR Galaxy Zoo effort (Paper III) with a physical size in excess of 80 kpc, or single component sources with a major axis size in excess of 80 kpc and which also exceed the resolved source threshold defined in Shimwell et al. (2019) by at least a factor of 1.5, were deemed to be clearly extended.These sources were added to the radio-excess sample if they were not already included (just under 0.5 per cent of the sample).
The lower panels of Fig. 8 show the ratio of measured radio luminosity over that expected from the consensus SFR as a function of redshift (left) and stellar mass (right); the horizontal dashed lines show the expected relation for star-forming galaxies and the radioexcess threshold, and the blue circles indicate again the peak of the distribution at each redshift.It can be seen that there is a weak variation of the population distribution with redshift, but no consistent trend, and the distribution peak never moves more than 0.2 dex (< 1) from the ridgeline value.Radio excess sources are found across all redshifts.The apparent gradual decline in the ratio with increasing redshift at  > 2.5 may be due to an increasing incompleteness in the classification of radiative-AGN at these redshifts (see Sec. 4), leading to an over-estimate in the SFR of some sources.
Regarding stellar mass, it is immediately clear from the lower-right panel of Fig. 8 that the proportion of radio-excess sources increases very strongly with mass, in particular for those objects not selected to be radiative-mode AGN.This is the well-known trend that, in the local Universe, the radio-loud AGN fraction shows a very strong mass dependence (e.g.Best et al. 2005;Sabater et al. 2019).Kondapally et al. (2022) use this LoTSS-Deep sample to investigate the cosmic evolution of this trend.Fig. 8 also shows a weak variation of the peak of the distribution of observed-to-predicted radio luminosity with mass, with a consistent trend of higher mass galaxies having on average a slightly higher radio luminosity for a given SFR.This has been previously seen in the radio luminosity to SFR relation (e.g.Gürkan et al. 2018;Smith et al. 2021), but it remains unclear to what extent this is due to an intrinsic mass-dependence of the amount of radio emission arising from star formation, as opposed to the effect of a contribution from a population of radio-weak AGN, more prevalent at higher stellar masses, that fall below the selection limit for radio-excess sources.
Regardless, the variations in Fig 8 are sufficiently small (in both redshift and stellar mass) that the use of a single SFR-radio relation does not significantly affect the selection of radio-excess sources.

FINAL RADIO SOURCE CLASSIFICATIONS, AND DEPENDENCIES
In the previous sections, LoTSS-Deep sources have been identified as either radiative-mode AGN or not, and either radio-excess sources or not, with a small number of sources being unclassifiable in each case.
Here, these are combined to derive a final set of source classifications.
• Sources which are neither radiative-mode AGN nor radio-excess sources are classified simply as star-forming galaxies (SFGs).Note that this may include some quiescent galaxies (with SFRs below the stellar mass vs SFR main sequence) whose low redshift nevertheless allows the star formation to be detected by LOFAR.
• Sources which are radiative-mode AGN but which do not display a radio excess are radio-quiet AGN (RQAGN; including the radio-quiet quasars) • Sources which are not radiative-mode AGN but do display a radio excess are the population of jet-mode AGN.Traditionally these sources are referred to as low-excitation radio galaxies (LERGs) • Sources which are both radiative-mode AGN and radio-excess sources are sources such as radio-loud quasars (Type I or Type II).These are traditionally referred to as high-excitation radio galaxies (HERGs).
• Finally, any source which could not be reliably classified in either of the criteria was left as unclassified.
Table 2 shows the number of sources of each class in each field.As can be seen, the majority population in LoTSS-Deep DR1 is the star-forming galaxies: these comprise just over two-thirds of the total population, rising to over 70 per cent in the deepest field, ELAIS-N1.Radio-quiet AGN contribute nearly 10 per cent of the total, with the two radio-loud classes contributing around 18 per cent between them, mostly as LERGs.Five per cent of the sources are unclassified.Of these, around 3 per cent are the sources without host galaxy identifications or redshifts for which no SED fitting could be carried out, and the remaining 2 per cent are mostly fainter galaxies for which the SED fitting algorithms either did not provide acceptable fits or provided highly inconsistent results.
Table 3 provides the first five lines of the classification data for each source in ELAIS-N1, along with the consensus mass and SFR measurements; the full catalogues for each field are provided electronically.More extensive catalogues, including the key outputs of each SED fitting code that were used to derive these, are made available on the LOFAR Surveys website (lofar-surveys.org).
Figure 9 shows the distribution of the different classes of source as a function of various properties of the host galaxy.The top panels show the distribution with respect to the 150-MHz flux density: the left panel shows the fraction at a given flux density, and the right panel shows the cumulative fraction above a given flux density.The population is dominated by radio-loud AGN above flux densities of about a mJy.The bulk of these are the LERGs, but with the fraction of HERGs beginning to rise at the highest flux densities, where the coverage of the sample begins to run out due to lack of sky area for these rarer bright sources.This rise of the HERG population is seen even more starkly in the middle left panel, which shows the distribution as a function of radio luminosity, and is in line with expectations from the relative luminosity functions of these two populations (e.g.Best & Heckman 2012;Best et al. 2014).At lower flux densities (and below 150 MHz luminosities of around 10 25 W Hz −1 ), star-forming galaxies take over the sample and quickly become the dominant population, accounting for over 90 per cent of sources at the limiting flux density reached in ELAIS-N1 (and more than 75 per cent of the cumulative population above  150MHz ≈ 100Jy).The switch between a starformation dominated population and a radio-loud AGN dominated population occurs at around  150MHz ≈ 1.5 mJy, which is fully consistent with the switch point at higher frequency of  1.4GHz ≈ 200Jy (found by Smolčić et al. 2017b) or  1.4GHz ≈ 250Jy (found by Padovani 2016), considering the typical radio spectral index of these sources.
At all flux densities below a few mJy there is a significant population of radio-quiet AGN, accounting for just under 10 per cent of all sources over the 100 Jy to 1 mJy flux density range.This is slightly lower than the fraction found in observations at higher frequencies: early work by Simpson et al. (2006) suggested that 20 per cent of sources with 100 Jy ∼ <  1.4GHz ∼ < 300 Jy are radio-quiet AGN, while the COSMOS 3GHz work of Smolčić et al. (2017b) indicated between 15 and 20 per cent (as determined from the 70 per cent subset of their 'High Luminosity AGN' sample that shows no radio excess).The origin of this difference is not completely clear.It may be related to different implementations of the radio-loud to radio-quiet separation, but more likely is associated with the radio-quiet AGN having a flatter spectral index than star-forming galaxies (e.g.due to a greater proportional contribution of flatter-spectrum core emission) and therefore lesser prominence at the lower frequencies probed by LOFAR.Given the steepness of the radio source counts, a difference   The fraction of sources of each different class (star-forming galaxies in grey; radio-quiet AGN in purple; low-excitation radio galaxies in blue; high excitation radio galaxies in orange; unclassifiable sources in yellow) as a function of radio flux density (upper panels; left gives fraction at a given flux density, and right gives cumulative fraction above a flux density), radio luminosity (middle left), stellar mass (middle right; for sources with  < 1.8 only -see text), optical r-band magnitude (lower left) and redshift (lower right; out to a final bin of 4 <  < 6).On each plot, the solid line for each class represents the derived fraction, and the shaded region indicates the calculated uncertainty.The open symbols show the values derived from each individual field (square = ELAIS-N1; asterisk = Lockman Hole; diamond = Boötes), where there are at least 5 sources from that field in the given bin, and demonstrate the broad agreement between fields.Note that the rise of the radio-quiet AGN population at the highest stellar masses is probably an artefact of larger mass uncertainties for these sources; see text for details.
Table 3. Classification results and consensus measurements for each source.The table shows the first five sources in ELAIS-N1: full catalogues are available electronically.Columns give the full source identifier, the radio ID number, the total 150 MHz flux density (in Jy), the redshift, the final radiative-mode AGN classification (1=AGN, 0=non-AGN, −1=unclassifiable), the logarithm of the consensus stellar mass (in solar masses), the logarithm of the consensus SFR (in solar masses per year), the radio excess (in dex), a flag to indicate extended radio sources (as defined in Sec.6; 1=extended, 0=compact), the final radio-AGN classfication (1=radio-AGN, 0=no radio excess, −1=unclassifiable), and the overall classification (SFR=star-forming galaxies; RQAGN=radio-quiet AGN; LERG=low-excitation (jet-mode) radio galaxy; HERG=high-excitation (quasar-mode)  The additional panels of Fig. 9 show the distribution of source classes as a function of redshift, stellar mass and optical magnitude.Note the strong rise of unclassified sources at  < 0.1; low SFRs for these galaxies can also lead to ambiguous radio excesses, while in addition the aperture photometry and aperture corrections used for the LoTSS Deep Field photometry (Paper III) are not optimised for these low redshifts, and resulting errors will affect the SED fitting.At these redshifts, it is in any case better to use the shallower, wider-area LoTSS surveys.All populations are seen over the full range of optical magnitudes.As expected, the LERG population shows increasing importance at higher stellar masses (note that this panel only includes redshifts  < 1.8 as mass estimates become increasingly less reliable at higher redshifts).The radio-quiet AGN show a dramatically increasing importance at stellar masses above 10 11.5  , but this is likely to be an artefact, driven by larger mass uncertainties for these sources due to the potential AGN contributions to their spectra: the number of sources at these very highest masses is relatively low, and so a few sources scattered up to high masses due to wider uncertainties on their masses, or due to errors in the photometric redshifts pushing them to higher redshift (and hence higher luminosity and mass), can artificially dominate the population.Interestingly, starforming galaxies are seen across the full range of redshifts studied; this indicates that the LoTSS-Deep sample is not only able to study normal star-forming galaxies in the low and moderate redshift Universe, but also to select starbursting galaxies in the early Universe.
All of these results are broadly consistent across the three fields (indicated by the open symbols in Fig. 9).In Sec. 4, the threshold levels for selection of radiative-mode AGN were set slightly differently in Boötes than the other two fields, based on the typically higher  AGN values found for the known spectroscopic and X-ray AGN and colour-selected probable AGN.The consistency of the classifications between fields in Fig. 9 gives confidence that this variation in thresholds is indeed appropriate.The remaining variations are consistent with what might be expected from cosmic variance, and indicate the importance of combining the multiple fields in order to overcome these effects, as well as to build a large statistical sample of sources.

COMPARISONS WITH SIMULATED SKY MODELS
Radio sky simulations provide a valuable tool for predicting the populations of radio sources that will be observed in a given survey.
In addition to the planning of future radio surveys (e.g.Norris et al. 2013) or predictions of parameter constraints achievable with those (e.g.Raccanelli et al. 2012;Harrison et al. 2016), these simulations are a valuable tool in assessing the completeness of different radio surveys (e.g.Hale et al. 2023), or in generating random samples for clustering analyses (e.g.Siewert et al. 2020).The two most widely used radio sky simulations in the literature are the SKA Design Study (SKADS) Simulated Skies (Wilman et al. 2008) and the more recent Tiered Radio Extragalactic Continuum Simulation (T-RECS; Bonaldi et al. 2019).
The starting point for these simulations is the measured luminosity functions of different source populations, and their cosmic evolution, which has typically been measured out to intermediate redshifts.The luminosity functions are then extrapolated to lower luminosities (lower flux densities), evolved out to higher redshifts, and potentially converted to a different observed frequency.Comparison of the predictions of these models against new deep observations such as the LoTSS Deep Fields provides a critical test of the assumptions that go into the radio sky simulations, and an opportunity to revise and improve these.
SKADS provides simulated predictions for four different radio source populations: star-forming galaxies, radio-quiet AGN, and two populations of radio-loud AGN.The two radio-loud AGN populations represent a low-luminosity and a high-luminosity component that Wilman et al. (2008) associated with the FRI and FRII morphological sub-populations (Fanaroff & Riley 1974), but which also map reasonably well onto the LERG and HERG classifications, respectively, used in this paper.Thus, all four radio source populations can be directly compared between the SKADS simulations and the LoTSS-Deep data.The radio-loud AGN population in T-RECS is constructed from luminosity functions for steep-and flat-spectrum radio sources together with BL Lac objects: these do not map onto the radio-AGN subclasses considered here, so comparisons with T-RECS can only be made with the radio-loud AGN population as a whole.T-RECS also includes predictions for SFGs, but does not include a separate radio-quiet AGN population: instead, T-RECS assumes that the radio emission of radio-quiet AGN is dominated by the on-going star-formation and thus that the radio-quiet AGN are encompassed within the star-forming population.
For both the SKADS and T-RECS simulations, a predicted source population was extracted over a randomly-located sky area corresponding to each of the three LoTSS Deep Fields.The simulations include sources to well below the flux limits of the observation and so, to replicate the observations, the LoTSS-Deep completeness simulations of Kondapally et al. (2022) and Cochrane et al. (2023) were used to determine the probability that each simulated source would be detected, and the source was randomly included in, or excluded from, the simulated catalogue in accordance with that probability.Figure 10 shows how the resultant simulated samples compare against the LoTSS-Deep data in both flux density (left panels) and redshift (right panels).Note that the small dip in the redshift distribution of all LoTSS-Deep populations over 1.0 <  < 1.5 is due to an aliasing effect in the photometric redshifts, particularly in the ELAIS-N1 and Lockman Hole fields, probably due to the lack of H-band data; this is discussed in more depth in Cochrane et al. (2023), but is not a significant issue for the analysis in the current paper.
The upper panels of Figure 10 show the simulation vs data comparison for a simple split into the two T-RECS source populations: star-forming galaxies plus radio-quiet AGN, against radio-loud AGN (HERGs + LERGs).Note that as well as allowing a comparison against both T-RECS and SKADS, this population split is arguably the most robustly determined in the LoTSS-Deep dataset, as it depends only on the presence or absence of a radio-excess rather than the (more difficult to establish) evidence for a radiative AGN.These upper panels show that both T-RECS and SKADS describe fairly well the transition between these two populations with decreasing radio flux density.T-RECS also provides an accurate match to the redshift distribution out to redshift  ∼ 4, beyond which the simulated source counts fall below those measured in the data; it is not clear whether this is a shortcoming of the simulation, or whether the photometric redshifts of the highest redshift sources become less reliable.The SKADS simulations also match the data reasonably well out to redshift  ∼ 2, but thereafter they over-predict the number of radio-loud AGN and under-predict the star-forming galaxy population.
The lower panels of Figure 10 provide further analysis of the SKADS simulations, split into the four sub-populations.Here, signficant differences are observed between the simulated and observed datasets.First, SKADS underpredicts the number of SFGs by a factor ≈ 2 at all redshifts  ∼ > 0.2.This is a result which has previous been established (e.g Bonaldi et al. 2016;Smolčić et al. 2017a); Hale et al. (2023) use a 'modified SKADS' model where they double the number of star-forming galaxies.Second, SKADS substantially overpredicts the number of radio-quiet AGN at lower redshifts and lower flux densities compared to the observations.Although it cannot be excluded that this is due to misclassification of faint radio-quiet AGN as star-forming galaxies in the observational data, a more likely explanation is that, as discussed earlier, this is due to an assumed radio spectral index of 0.7 for the radio-quiet AGN; a flatter spectral index (or curved spectral shape due to low-frequency absorption) would lead to a lower prevalence of these sources at the low frequencies of the LoTSS-Deep data.The combination of fewer SFGs and more RQAGN gives rise to the good agreement at low redshifts in the upper panel.For the radio-loud AGN, the difference in the high redshift number counts comes primarily from an over-prediction of the LERG population; the high redshift evolution of these sources was unknown at the time of the SKADS simulations, and so was assumed to be flat beyond  ∼ 0.7; recent works (e.g Kondapally et al. 2022) show this to be a reasonable assumption out to  ∼ 2, but with indications of a decline between 2.0 <  < 2.5, suggesting a breakdown of the SKADS assumptions.
In conclusion, while the SKADS simulations have been very suc-cessful in producing simulated radio skies, datasets such as LoTSS-Deep which probe new parameter space are revealing the shortcomings in our understanding 15 years ago when those simulations were first produced.The more modern T-RECS simulations provide a better match to the current dataset, but would be enhanced by the explicit inclusion of a radio-quiet AGN dataset, since the assumption that the radio emission of these sources is entirely produced through star-formation is known not to be true (see e.g.Macfarlane et al. 2021).Furthermore, explicit separation of the radio-loud population into HERG and LERG components in T-RECS would be a valuable addition and allow more detailed comparison of the simulation performance.

SUMMARY
The LoTSS Deep Fields are the widest deep radio survey ever undertaken.The LoTSS-Deep first data release, comprising ≈80,000 radio sources, is already an order of magnitude larger than previous radio source samples at this depth.The final LoTSS-Deep sample will detect > 250, 000 radio-selected sources over a 35 deg 2 region of sky, split into four different fields to largely overcome cosmic variance.
Extensive multi-wavelength photometry from the UV to the far-IR in each field facilitates a huge range of scientific exploitation.
In this paper, a combination of four different SED fitting codes has been applied to the multi-wavelength photometry of each of the LoTSS-Deep DR1 sources.Two of the four codes ( and -) include an AGN component in their SED modelling, and these offer an estimate of the AGN contribution to the overall galaxy SED.The other two codes ( and ) do not include AGN components, but offer more comprehensive coverage of the parameter space of the stellar component, and therefore are able to provide more accurate results for galaxies without AGN contributions.By combining the AGN fractional contributions estimated by and with the relative fitting ability of these two codes compared against and , those galaxies with an AGN contribution to their SED are identified.
Consensus stellar masses and star-formation rates are determined for each galaxy.For the galaxies without AGN contributions, these are generally based on the and results, which show excellent overall agreement with each other.For those which do show an AGN contribution to their spectra, the results are primarily adopted, as is shown to provide more reliable estimates than .The consensus star-formation rates are used to determine a relationship between 150 MHz radio luminosity and star-formation rate, using a 'ridgeline' approach to minimise bias from both radio selection effects and weak radio-AGN contributions.The determined relation is log 10  150MHz = 22.24 + 1.08 log 10 (SFR), where  150MHz is in units of W Hz −1 and SFR in units of  yr −1 .This is in very good agreement with previous literature studies.Radio-excess sources are then identified as those sources which show at least 0.7 dex (corresponding to ≈ 3) more radio emission than would be expected based on the star formation rate.
Using these results, the LoTSS Deep Field sources are then classified into four classes: (i) star-forming galaxies, which show neither any evidence for an AGN in their optical/IR SED nor a radio-excess; (ii) radio-quiet AGN, which do have an AGN contribution to their optical/IR SED, but show no radio excess; (iii) low-excitation radio galaxies (jet-mode radio-AGN), which show a radio excess but no optical/IR AGN signatures; (iv) high-excitation radio galaxies which show both AGN emission in their optical/IR SED and a radio ex-  (Bonaldi et al. 2019, dot-dash lines).The upper panels show the populations split just into star-forming galaxies plus radio quiet AGN (blue) versus radio-loud AGN (green), which can be compared against both SKADS and T-RECS simulations.The lower panels compare the four sub-populations against the SKADS simulation predictions; note that the separation of the two SKADS radio-loud classes does not map precisely onto the HERG/LERG classification used in this paper, although it is reasonably similar (see text).cess.Less than 5 per cent of the sources are unable to be classified.Overall, over two-thirds of the sources in the LoTSS Deep Fields are star-forming galaxies, around 16 per cent are LERGs, just under 10 per cent are radio-quiet AGN, and 2 per cent are HERGs.The three LoTSS Deep Fields show strong agreement in their source populations, despite significant differences in the input multi-wavelength photometric data.
The star-forming galaxies dominate the population below flux densities of  150MHz ≈ 1 mJy, accounting for ≈90 per cent of the sources close to the flux limit of the deepest field,  150MHz ∼ < 100Jy.In terms of luminosity, the star-forming galaxies become the largest population below  150MHz ≈ 10 25 W Hz −1 .At higher flux densities, and higher luminosities, the LERGs are the dominant population.The proportion of HERGs begins to rise significantly at the very highest flux densities and luminosities, but the LoTSS Deep Fields do not cover enough sky area to probe the regime where these become the dominant population.
Star-forming galaxies are observed across all redshifts, ranging from normal star-forming galaxies in the nearby Universe to extreme starbursting systems at  > 4. They are also observed across a wide range of optical magnitudes and stellar masses, peaking at around 10 10.5 solar masses, typical of galaxies towards the upper end of the star-forming main sequence.The proportion of radio-quiet AGN rises noticeably towards higher redshifts; it also rises sharply towards the highest stellar masses, but this is likely to be an artefact of the steep stellar mass function coupled with larger uncertainties on the stellar masses of this population.The LERG population reaches its peak importance at redshifts 1 to 3; however, the proportion of LERGs is smaller than that of the star-forming galaxies at all redshifts, stellar masses and optical magnitudes.
The observed populations are compared against the prediction of the SKADS and T-RECS radio sky simulations.SKADS is shown to underpredict the star-forming galaxy population by a factor ≈ 2 across all redshifts.It over-predicts the proportion of radio-quiet AGN in the sample.This is likely to be due to the assumption of a radio spectral index of  = 0.7 for these sources: a flatter spectral index, as indicated by recent LOFAR observations of radio-quiet quasars, would reduce the prevalence of these sources in these lowfrequency observations.Finally, SKADS over-predicts the numbers of LERGs at redshifts  > 2, as it does not account for the negative cosmic evolution of this population at high redshift beginning to be observed in the latest datasets.T-RECS provides a good match to the star-forming and radio-loud AGN populations, but its lack of a separate radio-quiet AGN population is a significant shortcoming.
The classifications, stellar masses and SFRs derived in this paper form a vital input to many other studies using the LoTSS Deep Fields first data release (Smith et al. 2021;Bonato et al. 2021;Kondapally et al. 2022;McCheyne et al. 2022;Mingo et al. 2022;Cochrane et al. 2023, and others), and the techniques developed to derive these can be applied to future data releases of the LoTSS Deep Fields.Many advances continue to be made in the LoTSS Deep Fields that, in addition to new deeper radio data, will improve classifications still further.Over the next 5 years, the WEAVE-LOFAR survey (Smith et al. 2016) will obtain around a million optical spectra of LOFAR sources, including all sources detected in the LoTSS Deep Fields, using the new William Herschel Telescope (WHT) Enhanced Area Velocity Explorer (WEAVE) multi-object spectrograph (Jin et al. 2023).WEAVE-LOFAR will provide spectroscopic redshifts for the vast majority of the star-forming galaxies, radio-quiet AGN and HERGs (especially at lower redshifts) due to their strong emission lines, removing one of the largest uncertainties in the SED fitting.It may be possible to obtain spectroscopic redshifts for LERGs from weaker lines or continuum features, and even where this is not the case, the confirmed absence of strong emission lines and AGN features will add confidence to the reliability of the photometric redshifts.For many sources, WEAVE-LOFAR will also improve source classifications through either emission line diagnostics, or emission line to radio flux ratios (cf.Best & Heckman 2012, at lower redshifts).Future imaging of these fields at 0.3-arcsec resolution, by including the international LOFAR baselines (cf.Morabito et al. 2022;Sweĳen et al. 2022), will further improve source classification by allowing compact radio cores (AGN), kpc-scale star-forming regions, and small-scale core-jet radio sources to be distinguished by their radio morphology in these fields (Morabito et al. 2022).A comparison between the SED-determined classifications and those from high resolution radio morphology will be very interesting.
The final LoTSS-Deep sample, imaged with sub-arcsec radio resolution and coupled with high-resolution optical spectroscopy for each source, will represent an extremely powerful resource for studies of the evolution of galaxies and AGN.

APPENDIX A: UNCERTAINTIES ON STELLAR MASSES AND SFRS
As discussed in the main text, no attempt is made to derive stellar mass or SFR uncertainties on a source-by-source basis: any reader interested in individual sources can examine the results of all of the different SED codes, provided in the extended tables on the lofarsurveys.orgwebsite, and make their own assessment of the relevant systematic and statistical errors.Instead, this appendix examines typical uncertainties that can broadly be considered.
For sources not identified as radiative-mode AGN, the consensus stellar masses are derived using the and results.For each source, the difference between the outputs of the two codes provides an indication of systematic uncertainties, while the confidence intervals provided by each code give an estimate of statistical uncertainty.As a broad guide to the dominant uncertainty, for each source the higher of these two values is considered.Figure A1 then shows the median of this value for all galaxies within a given bin in redshift-mass space.As can be seen, the calculated median uncertainties are typically ∼ < 0.1 dex at lower redshifts and higher masses.As expected, they increase towards lower masses and towards higher redshifts, in both cases due to the galaxies being fainter and therefore having lower signal-to-noise photometric measurements in the SED fitting.
The uncertainties in Fig. A1 can broadly be categorised in four different ranges of parameter space, with empirical estimates of the uncertainty possible for each: • Higher mass, lower redshift: specifically log 10  * ≥ (9.7 + 2.5 log 10 (1 + )) and  ≤ 2.Here the uncertainty on stellar mass is fairly constant at Δ * ≈ 0.1 dex.
Uncertainties on the stellar masses of the radiative AGN are harder to estimate in this manner, as mass estimates are derived from the two fits, and these are likely to be subject to related systematic errors.Comparing the confidence intervals of the fits with those of the non-AGN in the same redshift-mass bin, the statistical uncertainties of the radiative-mode AGN are on average 20 per cent larger than those of the non-AGN; this sets a lower limit to the mass uncertainty estimate, although it is likely that the systematic errors will also be larger in cases where the AGN contributes significantly to the optical to near-IR spectrum.
A similar approach can be followed to estimate the typical uncertainties on the consensus SFR estimates.Fig. A2 shows the result, split into bins of SFR and redshift.In this case it is apparent that the SFR estimates are generally robust until the very lowest SFRs at any redshift are reached (at most a few per cent of objects), where the uncertainties increase dramatically.For the vast majority of the population at higher SFRs, there is no strong dependence of the SFR uncertainty (in dex) on the measured SFR, but a clear trend for the uncertainty to increase with redshift, from ≈ 0.1 dex at  ∼ 0 up to 0.15 dex by  ∼ 1 and 0.2 dex by  = 3.This can be empirically ap- proximated as Δ(SFR) ≈ 0.1× (1+) 0.5 dex.The contributions to this uncertainty from differences between codes and from the statistical uncertainties within individual codes are comparable in size.
It should be emphasized again that these empirical relations are only intended to provide a guide to the approximate stellar mass and SFR uncertainties, and do not represent reliable values on a source-by-source basis.
This paper has been typeset from a T E X/L A T E X file prepared by the author.

FIRSTFigure 1 .
Figure1.The survey depth, area and angular resolution of the LoTSS Deep Fields compared to other existing and on-going radio surveys.All survey depths are converted to a 1.4 GHz equivalent rms depth using a spectral index of  = 0.7.The black points show published surveys, and the blue points show on-going surveys.The LoTSS Deep Fields are highlighted in red.The size of each symbol indicates the angular resolution of the survey, with the symbol area proportional to the beam FWHM.For the LoTSS Deep Fields final release, the larger symbol indicates the result of including just the Dutch baselines, while the smaller symbol shows what should be achievable after including the international stations (improved angular resolution, additional depth due to the extra collecting area, but smaller areal coverage due to the smaller primary beam of the international stations).Descriptions of the surveys included on the plot (listed from high to low effective rms depth) can be found in the following references: GLEAM(Wayth et al. 2015); WENSS(Rengelink et al. 1997); TGSS(Intema et al. 2017); SUMSS(Mauch et al. 2003); NVSS(Condon et al. 1998); GLEAM-X(Hurley-Walker et al. 2022); RACS(Hale et al. 2021); FIRST(Becker et al. 1995); XXL-GMRT(Smolčić et al. 2018); VLASS(Lacy et al. 2020); Stripe82(Hodge et al. 2011); LoTSS-Wide(Shimwell et al. 2019); VLA-COSMOS 1.4 GHz(Schinnerer et al. 2007); EMU(Norris et al. 2011); MIGHTEE (including Early Science -ES; Heywood et al. 2022) SSA-13 (Fomalont et al. 2006); VLA-COSMOS 3 GHz (Smolčić et al. 2017a); VLA-SWIRE (Owen & Morrison 2008); GOODS-N (Owen 2018); VLA Frontier (Heywood et al. 2021).
shows the Stern, Lacy and Donley mid-IR colour-colour plots with the LoTSS-Deep sources in Boötes 1 colour-coded by their AGN fraction as derived by using the model.Sources classified as an AGN

Figure 2 .
Figure 2. The location of the LoTSS-Deep sources on the Lacy et al. (2004) and Donley et al. (2012, left) and on the Stern et al. (2005, right) mid-IR colour-colour classification plots (for sources with S/N>2 in all four bands), in the Boötes field.The blue dashed lines on the left-hand panel show the Lacy et al. selection criteria, and the blue solid lines show those of Donley et al.On the right-hand plot, the Stern wedge is shown by the blue dashed lines.In both plots, the greyscale colour-coding indicates the AGN fraction from the SED fitting using the AGN model.Objects confirmed to be AGN through optical spectroscopy or X-ray observations are indicated by the red circles.

Figure 3 .
Figure 3.The AGN fraction plotted against the ratio of the  2 values between SED codes that do not include AGN components (the lower value for the (MP) and (BP) fits) and those that do (the lowest of the (CG) and (AF) fits), for the LoTSS-Deep sources in Boötes.Points are colour-coded according to whether they are spectroscopic or X-ray AGN (red filled circles), or satisfy the Donley criteria (with S/N>3 in each band; blue open circles), or satisfy the broader Stern, Lacy or Messias cuts (with S/N>3; black crosses), or 'non-AGN' that either do not satisfy any cuts or have too low signal-to-noise in the mid-IR for this to be determined (green triangles).The clustering at certainAGN fraction values (e.g.0.05, 0.7) appears to be a feature of the code, perhaps due to the fairly limited sampling of the grid of AGN model parameters.The plot shows that as the AGN fraction rises above ∼0.1,objects are more likely to be identified as AGN through spectroscopic or X-ray selection or the Donley mid-IR cuts, and also that the SED fitting begins to deteriorate (higher relative  2 ) for SED codes that don't include AGN components.

Figure 4 .
Figure 4.The selection criteria used to identify radiative-mode AGN, and the relative distributions of the AGN and non-AGN thus identified.The upper-left panel compares the reduced  2 value resulting from SED models including an AGN component ( , ) against those which do not ( , ), with the blue dashed line showing selection criterion (iii).The upper right plot shows the 1-sigma lower limits (16 th percentile; P16) to the AGN fraction from and (with the AGN model), with the blue dashed lines showing selection criteria (i) and (ii).The lower plots show selection criteria (i) vs (ii) and (i) vs (iii) in the left and right panels respectively.Data shown are for ELAIS-N1.Sources are selected as radiative-mode AGN if they satisfy at least two of the three criteria (or are confirmed AGN from spectroscopic or X-ray observations): these sources are shown in red.

Figure 5 .
Figure 5. Infrared colour-colour plots for the LoTSS-Deep sources in ELAIS-N1, compared with template spectra.The sources classified as (radiative-mode) AGN are plotted in red and the non-AGN in black symbols; sources are only plotted if they have a signal-to-noise of at least 3 in each of the relevant filters.For clarity, sources (and templates) are divided into three redshift ranges: the top row is for  < 1, the middle row for 1 <  < 2.5 and the bottom row for  > 2.5.For each redshift, the left-hand plot shows the mid-IR IRAC flux ratios used for the Lacy et al. (2004, blue dashed lines) and Donley et al. (2012, blue solid lines) selections.The middle column shows the Stern et al. (2005) colour criteria, with the Stern region indicated by the blue dashed lines.The right-hand column shows the selection criteria proposed by Messias et al. (2012), combining IRAC colours with the K-band flux at the lower redshifts, and with the 24m flux at the highest redshifts.In each plot the coloured lines indicate the evolution over the specified redshift range of a selection of galaxy and AGN template spectra, from Brown et al. (2014) and Brown et al. (2019), separated into ellipticals (pink), star-forming galaxies (yellow), AGN (purple), and composites (green).As can be seen, the broad colour cuts suffer to various extents from both incompletensss and contamination.The selected AGN broadly align with the regions of colour space covered by the AGN and composite template spectra.

Figure 6 .
Figure 6.A comparison of the masses derived by the different SED fitting codes against the final consensus masses for the LoTSS Deep Field sources in ELAIS-N1., and all give broadly consistent results for non-AGN, but differ for the AGN subset, for which the results should be more reliable.massesshow a small systematic offset compared to the other codes, and more outliers at high mass.The lower right plot examines the masses produced in Paper IV (only out to  < 1.5); these are seen to give consistent results with only slightly larger scatter.This is of interest because these stellar masses were produced for the entire galaxy population in these deep fields, not just the LoTSS-Deep host galaxies.

Figure 7 .
Figure 7.A comparison of the star-formation rates derived by the different SED fitting codes against the final consensus value, for the LoTSS Deep Field sources in ELAIS-N1.andgive broadly consistent results for non-AGN; their performance on objects identified as (radiative-mode) AGN is more mixed, but generally reasonable where the fit is not flagged as a bad fit.'s SFR estimations for non-AGN generally perform well at higher SFRs (especially with the AGN model), but over-predict the SFR in some lower-SFR galaxies.The estimated SFRs of objects selected as AGN show a high degree of consistency between the two different runs.SFRs show more scatter and a small systematic offset compared to the other codes.

Figure 8 .
Figure 8. Top: the distribution of radio luminosity versus SFR for LoTSS Deep Field sources in ELAIS-N1, split into those identified as radiative-mode AGN from their SED (red points) and the sources which are not radiative-mode AGN ('SED non-AGN'; black points).Within narrow bins in SFR, the 'ridgeline' points (larger blue circles) indicate the peak of the distribution of radio luminosities.These can be well-fitted by a power-law distribution shown by the solid blue line, which is in broad agreement with literature relations (green lines).Bottom: the ratio of observed radio luminosity to that predicted from the consensus SFR based on the ridgeline fit, versus redshift (left) and stellar mass (right).The horizontal dashed lines represent the expected relation and the radio-excess threshold.Solid blue points in each plot show the peak of the distribution in narrow bins.These always lie within 0.2 dex of the expected relation.Radio-excess sources are found over the full range of redshifts, but predominantly concentrate at high stellar masses.
Figure9.The fraction of sources of each different class (star-forming galaxies in grey; radio-quiet AGN in purple; low-excitation radio galaxies in blue; high excitation radio galaxies in orange; unclassifiable sources in yellow) as a function of radio flux density (upper panels; left gives fraction at a given flux density, and right gives cumulative fraction above a flux density), radio luminosity (middle left), stellar mass (middle right; for sources with  < 1.8 only -see text), optical r-band magnitude (lower left) and redshift (lower right; out to a final bin of 4 <  < 6).On each plot, the solid line for each class represents the derived fraction, and the shaded region indicates the calculated uncertainty.The open symbols show the values derived from each individual field (square = ELAIS-N1; asterisk = Lockman Hole; diamond = Boötes), where there are at least 5 sources from that field in the given bin, and demonstrate the broad agreement between fields.Note that the rise of the radio-quiet AGN population at the highest stellar masses is probably an artefact of larger mass uncertainties for these sources; see text for details.

Figure 10 .
Figure 10.A comparison of the radio source population fractions as a function of 150 MHz flux density (left panels) and the redshift distribution of radio sources (right panels) between the LoTSS-Deep data (solid lines and shaded regions) and the simulated sky predictions from SKADS (Wilman et al. 2008, dashed lines) and T-RECS(Bonaldi et al. 2019, dot-dash lines).The upper panels show the populations split just into star-forming galaxies plus radio quiet AGN (blue) versus radio-loud AGN (green), which can be compared against both SKADS and T-RECS simulations.The lower panels compare the four sub-populations against the SKADS simulation predictions; note that the separation of the two SKADS radio-loud classes does not map precisely onto the HERG/LERG classification used in this paper, although it is reasonably similar (see text).

Figure A1 .
Figure A1.Typical uncertainties on estimates of stellar mass, as a function of both stellar mass and redshift, for non-AGN in the ELAIS-N1 field.At lower redshifts and higher stellar masses, the uncertainty is generally ∼ 0.1dex, but this increases towards lower stellar masses and for redshifts  > 2.

Figure A2 .
Figure A2.Typical uncertainties on estimates of SFR, as a function of both SFR and redshift, for non-AGN in the ELAIS-N1 field.SFR uncertainties increase dramatically for the few per cent of sources at the lowest SFRs at each redshift; above that they have little dependence on SFR, but increase gradually with increasing redshift.

Table 1 .
Status of observations and imaging in LOFAR Deep Fields, including the data released in the LoTSS Deep Fields 1st data release (LoTSS-Deep DR1).

Table 2 .
The number of sources of each class in the LoTSS-Deep DR1 dataset.
radio galaxy; Unc=unclassified.Values of −99 indicate where no measurement is available.