The LOFAR Two-metre Sky Survey: The nature of the faint source population and SFR-radio luminosity relation using Prospector

Spectral energy distribution (SED) fitting has been extensively used to determine the nature of the faint radio source population. Recent efforts have combined fits from multiple SED-fitting codes to account for the host galaxy and any active nucleus that may be present. We show that it is possible to produce similar-quality classifications using a single energy-balance SED fitting code, Prospector, to model up to 26 bands of UV$-$far-infrared aperture-matched photometry for $\sim$31,000 sources in the ELAIS-N1 field from the LOFAR Two-Metre Sky Survey (LoTSS) Deep fields first data release. One of a new generation of SED-fitting codes, Prospector accounts for potential contributions from radiative active galactic nuclei (AGN) when estimating galaxy properties, including star formation rates (SFRs) derived using non-parametric star formation histories. Combining this information with radio luminosities, we classify 92 per cent of the radio sources as a star-forming galaxy, high/low-excitation radio galaxy, or radio-quiet AGN and study the population demographics as a function of 150 MHz flux density, luminosity, SFR, stellar mass, redshift and apparent $r$-band magnitude. Finally, we use Prospector SED fits to investigate the SFR$-$150 MHz luminosity relation for a sample of $\sim$$133,000~3.6~\mu$m-selected $z<1$ sources, finding that the stellar mass dependence is significantly weaker than previously reported, and may disappear altogether at $\log_{10} (\mathrm{SFR}/M_\odot~\mathrm{yr}^{-1})>0.5$. This approach makes it significantly easier to classify radio sources from LoTSS and elsewhere, and may have important implications for future studies of star-forming galaxies at radio wavelengths.


INTRODUCTION
Understanding the evolution of galaxies from the early stages of the Universe to the present is one of the key goals in modern astrophysics.Two critical aspects of the evolutionary pathway of galaxies are the build-up of stellar mass over cosmic time (i.e., the star formation history or SFH; see review by Madau & Dickinson 2014), and feedback from active galactic nuclei (AGN; see reviews by McNamara & Nulsen 2007;Fabian 2012;Heckman & Best 2014;Hardcastle & Croston 2020).Radio observations are particularly important for conducting simultaneous studies of both star formation and nuclear activity in galaxies, being virtually unaffected by the dust obscuration that plagues shorter wavelengths (e.g.Condon & Ransom 2016;Smith et al. 2016).Radio emission in galaxies is a combination ★ E-mail: s.das2@herts.ac.uk of two principal components (e.g.Condon 1992): thermal free-free emission (which contributes significantly only at frequencies above ∼ 1 GHz), and non-thermal synchrotron radiation.The synchrotron component dominates at low frequencies, and is powered by (i) the acceleration of cosmic rays by supernova remnants of short-lived massive stars (Blumenthal & Gould 1970;Condon 1992;Condon & Ransom 2016) and (ii) jets and lobes produced by the central supermassive black hole (SMBH; see reviews by Bridle & Perley 1984;Padovani et al. 2017;Hardcastle & Croston 2020) if the galaxy hosts an AGN.
AGN activity can be separated into two fundamental modes based on the efficiency of matter accretion onto the SMBH (see Hardcastle et al. 2007;Best & Heckman 2012;Janssen et al. 2012;Heckman & Best 2014;Mingo et al. 2014;Tadhunter 2016;Hardcastle & Croston 2020, and references therein).Radiatively efficient AGN activity (also known as "radiative-mode" or "quasar-mode") is known to form over at least 25 deg2 , a factor of 2 − 4× deeper compared to the LoTSS wide survey, with extensive wide-field UV, optical, and infrared coverage.Of particular relevance to this work, Best et al. (2023, hereafter B23) performed SED fitting of about 80,000 radiodetected LoTSS deep fields sources using four codes; Magphys (da Cunha et al. 2008), Bagpipes (Carnall et al. 2018), Cigale (Noll et al. 2009;Boquien et al. 2019), and AGNfitter (Calistro Rivera et al. 2016).There are two principal differences in these codes.Firstly, the Cigale and AGNfitter models account for the contribution of AGN to the SED while Magphys and Bagpipes do not.Secondly, Magphys, Bagpipes, and Cigale assume energy balance between dust attenuation and emission, while the version of AGNfitter used in B23 does not.By comparing the results from these different SED fitting codes, B23 identified galaxies hosting radiative-mode AGN and derived optimised consensus estimates of the stellar mass and SFR.Furthermore, they identified galaxies with radio luminosity exceeding that expected on the basis of the star formation processes alone.Using this information, B23 classified the radio sources as either SFGs, RQ AGN, LERGs, or HERGs.
Energy balance SED fitting represents the state-of-the-art for estimating the physical properties of galaxies from photometry at high-redshift (and subsequently classifying them).However, due to the computational resources required, and because obtaining scientifically-usable results from even a single code is not trivial, using multiple SED fitting codes to perform galaxy classifications may not be optimal.
In this work, we explore the viability of using a single SED fitting code to classify galaxies on the basis of their multi-wavelength photometry (following the approach adopted by B23) and use our new SED fitting results to revisit the stellar mass-dependence previously reported in the SFR-150 MHz 2 radio luminosity relation (e.g., Gürkan et al. 2018, hereafter G18, andSmith et al. 2021, hereafter S21).To do this, we use the SED fitting code Prospector (Leja et al. 2017;Johnson et al. 2021) as (i) it can model varying contributions of nebular emission and AGN to galaxy SEDs, and (ii) it is capable of modelling multi-wavelength photometry and spectroscopy on an equal footing.Even though Prospector can not model the radio SEDs of galaxies (cf.e.g.Cigale, Dey et al. 2022, Magphys, da Cunha et al. 2015, and Prospect, Robotham et al. 2020;Thorne et al. 2023), these aforementioned attributes are particularly useful in the context of this work, since radio-selected sources are known to be rich in both AGN and emission lines (e.g.Netzer 1990;Simpson et al. 2006).While the current work will focus on the photometric classification (since spectra exist for only a minority of the sources in the LoTSS deep fields; e.g.Sabater et al. 2019;Duncan et al. 2021, Drake et al. in preparation), the ability to model spectra simultaneously with Prospector will become particularly useful in the coming years as we build up complete spectroscopy for sources in the LoTSS deep fields with the WEAVE-LOFAR survey (Smith et al. 2016).
This paper is structured as follows: in Section 2, we describe the datasets used in this work.In Section 3, we discuss the process of SED fitting using Prospector.In Section 4, we outline the criteria that we use to identify sources hosting radiative-mode and/or radioloud AGN, summarise our classification scheme, and look into the demographics of the different galaxy classes as a function of a number of physical parameters.We also use our parameter estimates to revisit the relation between star formation rate and 150 MHz radio luminosity.In Section 5, we look into the stellar mass-dependent calibration of the SFR-150 MHz radio luminosity relation.Finally, in Section 6 we summarise the key findings of this work.

DATA
Of the four LoTSS deep fields, the ELAIS-N1 region has the most sensitive 150 MHz data in LoTSS DR1, making it ideal for studying faint star-forming galaxies.ELAIS-N1 also benefits from some of the deepest wide-field optical, near-IR and mid-IR surveys.In this work, we focus on the 6.7 deg 2 of this field where the best ancillary data are available.Kondapally et al. (2021) generated 26 bands of aperture-matched photometry by combining aperture-and galactic extinction-corrected fluxes from UV to the mid-IR wavelengths with far-IR photometry measured using the XID+ tool (Hurley et al. 2017;McCheyne et al. 2022).UV-mid IR fluxes come from the following: Spitzer Adaptation of the Red-sequence Cluster Survey (SpARCS) -band (Wilson et al. 2009;Muzzin et al. 2009); Panoramic Survey Telescope and Rapid Response System (PanSTARRS) , , , ,  bands (Kaiser et al. 2010;Chambers et al. 2016); , , , , , and the narrow-band NB921 data from Hyper-Suprime-Cam (HSC; Miyazaki et al. 2012Miyazaki et al. , 2018) ) Subaru Strategic Program (HSC-SSP) public DR1 (Aihara et al. 2018); -and -band data from the UK Infrared Deep Sky Survey (UKIDSS) Deep Extragalactic Survey (DXS) DR10 (Lawrence et al. 2007); Infrared Array Camera (IRAC) 3.6-8m data from the Spitzer Wide-Area Infrared Extragalactic Survey (SWIRE; Lonsdale et al. 2003), and 3.6-4.5mdata from the Spitzer Extragalactic Representative Volume Survey (SERVS; Mauduit et al. 2012).Duncan et al. (2021) combined state-of-the-art template fitting and machine learning techniques to obtain photometric redshift (photo-, or  photo ) information for all the optical/near-IR selected sources in the LoTSS deep fields.These photo-s complemented the spectroscopicallycompiled redshifts (spec-, or  spec ) which were available for a small minority (approximately 0.2 per cent) of the deep field sources.Far-IR photometry used in this work included 24 m data from the Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 2004) instrument on Spitzer, 250, 350 and 500 m data from the Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010) instrument, and 100 and 160 m data from Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010).The latter two were obtained as part of the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012).Duncan et al. (2021) flagged sources with prior identification as AGN based on optical spectra, X-ray properties, or mid-IR colour cuts.We refer the reader to Kondapally et al. (2021), McCheyne et al. (2022) and references therein for further details.This multiwavelength source catalogue contains ∼ 2.1 million objects with redshifts up to  ≈ 7.

Radio Observations
The 150 MHz radio data in the ELAIS-N1 deep field were taken with the LOFAR High Band Antenna (HBA) over a period of many years (totalling 164 hours of on-source integration time in DR1) and includes observations from the Dutch LOFAR stations, reaching spatial resolution (equivalent to the full width at half maximum of the restoring beam) of 6 arcsec with an RMS below 30 Jy beam −1 .Away from the bright sources, the central region of ELAIS-N1 reaches RMS as deep as 20 Jy beam −1 .Sabater et al. (2021) performed calibration and imaging, and catalogued the extracted sources.Owing to the excellent multiwavelength data coverage in this field, Kondapally et al. (2021) identified optical/near-IR counterparts for over 97 per cent of the radio-detected sources.We refer the reader to Sabater et al. (2021) and Kondapally et al. (2021) for details on source extraction and radio-optical cross-matching.

Datasets for SED fitting
In this work we use two different photometric samples.
(i) Radio-selected sample -The final cross-matched catalogue in the ELAIS-N1 field contains 31,610 radio sources (Kondapally et al. 2021).We refer to this subset of sources as the "radio-selected sample".Redshift information (spectroscopic and/or photometric) is available for 30,538 of these sources (Duncan et al. 2021).Out of these 30,538 sources, 30,470 satisfied the flag cuts recommended in Kondapally et al. (2021).The flag cuts ensured that these sources lay in the overlapping coverage area of PanSTARRS -band, UKIDSS-DXS -band, and Spitzer-SWIRE 4.5m channel, while excluding regions around bright stars as masked in the Spitzer data.
(ii) SWIRE-selected sample -The majority of sources in the multiwavelength catalogue were not individually catalogued at 150 MHz by Sabater et al. (2021).Excluding these sources risks biasing the results towards the minority of radio-bright sources.Therefore we also follow the S21 sample definition and select sources with SWIRE 3.6m flux density exceeding 10 Jy.This roughly corresponds to the 5 detection threshold in the SWIRE 3.6m data and is in line with SWIRE webpage3 recommendations.S21 adopt the SWIRE 3.6m selection as these data are extremely sensitive and are less susceptible to the influence of dust obscuration and sample biases than samples identified at shorter wavelengths.We restrict this sample to  < 1, where Duncan et al. (2021) photo- estimates are the most reliable, bringing the sample size to 145,162 sources.Finally, we exclude sources at  < 0.05 to avoid potential IRAC sources with extended radio emission larger than the LoTSS beam size of 6 arcsec.More than 99.5 per cent of sources at  < 1 and over 98 per cent of sources at  < 0.2 have an FWHM smaller than 6 arcsec (see S21), therefore this redshift cut does not significantly impact our work.The resulting source sample contains 133,544 sources and is referred to as the "SWIRE-selected sample".
We use the Sabater et al. (2021) 1 indicates the sample size as a function of 150 MHz luminosity, with the 150 MHz non-detected sources assigned to the (numerically dominant) lowest luminosity bin if their maximum-likelihood luminosity estimates were lower than this.The individually undetected sources (with S/N < 5, represented by the blue data points below the lower limit of the purple distribution in Figure 1) are numerically dominant and together can provide significant analytical insights, making it critical to include these sources in our analysis.

SED FITTING
Prospector (Leja et al. 2017;Johnson et al. 2021) is an SED fit-ting code that uses the Flexible Stellar Populations Synthesis (FSPS;Conroy et al. 2009;Johnson et al. 2023) code to model photometric and/or spectroscopic data and infer stellar population properties.It uses the Padova isochrone table (Bertelli et al. 1994;Girardi et al. 2000;Marigo et al. 2008) and the MILES spectral library (Sánchez-Blázquez et al. 2006).The dynamic nested sampling tool Dynesty (Speagle 2020) is used to estimate Bayesian posteriors and evidences.We adopt the Initial Mass Function (IMF) from Kroupa (2001).Since Prospector does not support evolving metallicities (in contrast to, e.g., Prospect; Thorne et al. 2022b), and even non-evolving metallicities remain poorly constrained based on photometry alone (e.g.Smith & Hayward 2018), we set the stellar metallicity at solar.While Thorne et al. (2022b) have shown that the stellar masses and SFRs are not systematically offset by this choice of metallicity, we intend to revisit this choice in future works once WEAVE-LOFAR spectra are available.
Source redshifts have been fixed at spec- for sources where they are available (5.1 per cent of the radio-detected sample and 1.5 per cent of the SWIRE-selected sample).For the remaining sources, they were fixed at the photo- estimates from Duncan et al. (2021).Following B23, neither redshift uncertainties, nor those sources lacking redshift information were considered further.Applying the flags recommended in Kondapally et al. (2021), we ensure that only those sources with the most reliable photo- are included in our work. 4e use the two-component Charlot & Fall (2000) dust attenuation model, which consists of separate birth cloud and diffuse dust screens, to account for the differential attenuation of starlight from younger and older stellar populations.We use the Kriek & Conroy (2013) attenuation curve which is parametrised by the diffuse dust -band optical depth, the diffuse dust attenuation index, and the ratio of the birth cloud optical depth to the diffuse dust screen optical depth, while energy balance is assumed between dust attenuation and emission (see Leja et al. 2017, for further discussions).The three parameter Draine & Li (2007) dust emission templates are used to describe the shape of the IR SED.Mid-IR emission from the AGN torus is described by the CLUMPY models (Nenkova et al. 2008a,b).These models are parametrised by the 4 − 20m AGN luminosity (which is expressed as a fraction of the total stellar bolometric luminosity) and the optical depth of the AGN torus.Nebular lines and continuum emission are calculated through the implementation of the CLOUDY photoionization code (Ferland et al. 2013), wherein the ionizing sources are represented by the stellar populations generated using FSPS (Byler et al. 2017).This is parametrised by the intensity of the ionizing spectrum, which is denoted by a dimensionless ionization parameter (and kept free in the model), and the gas-phase metallicity.While Leja et al. (2017) demonstrated that gas-phase metallicity is not adequately constrained for individual galaxies using photometry, works such as Bellstedt et al. (2020), Bellstedt et al. (2021) and Thorne et al. (2022a) have highlighted the benefits of modelling galaxies using an evolving metallicity (e.g. for recovering the best estimates of the cosmic star formation history, the CSFH).Since that functionality is not available in Prospector, we have made the pragmatic decision to set the gas-phase metallicity equal to the solar metallicity.Gas-phase metallicity Equal to log( ⊙ ) Table 1.Description of the Prospector parameters used in the final model.† The Prospector models use the total stellar mass formed in a source as a free parameter, which is larger than the currently surviving stellar mass.We use the latter for analysis purposes throughout this work.
Star formation histories have traditionally been modelled with a number of functional forms (such as constant or exponentiallydeclining SFHs, commonly referred to as 'parametric' SFHs; e.g.Buat et al. 2008;Maraston et al. 2010;Gladders et al. 2013;Simha et al. 2014;Carnall et al. 2018, 2019, andreferences therein).Parametric SFHs often struggle to model complex SFH features, such as episodes of star-burst activity, sudden quenching, and rejuvenation (Simha et al. 2014;Smith & Hayward 2015;Carnall et al. 2019;Leja et al. 2019).Non-parametric SFHs promise a way out by not assuming a specific functional form for the SFH.In its simplest form, the lifetime of the galaxy is divided into discrete time bins and the stellar mass formed in each bin is fitted (e.g., Cid Fernandes et al. 2005;Ocvirk et al. 2006;Tojeiro et al. 2007;Leja et al. 2017).Iyer & Gawiser (2017), Leja et al. (2019), Lower et al. (2020), andCiesla et al. (2023) have shown that non-parametric SFH models provide the flexibility to describe the full extent of complex SFHs, potentially leading to more robust estimates of galaxy properties and outperforming parametric SFH models, usually at the expense of a few extra free parameters.Given the inability of photometric data to fully constrain SFHs of galaxies (e.g.Smith & Hayward 2015) the choice of prior distribution plays an important role (e.g.Leja et al. 2019).
In this work, we adopt the Prospector non-parametric "continuity" prior, which directly fits for Δ log 10 (SFR) between neighbouring time bins using a Student's -distribution with scale factor  = 0.3 and two degrees of freedom ( = 2).This prior discourages abrupt changes in SFR between adjacent bins, but remains flexible enough to fit starbursts, star-forming and quenched galaxies (Leja et al. 2019;Johnson et al. 2021;Haskell et al. 2023b, Das et al. in preparation).We divide the SFH up such that the first two and the last time bins are kept the same for all sources, covering 0 <   < 10 Myrs, 10 <   < 100 Myrs, and 0.85  <   <   , respectively.Here,   represents the age of the Universe at the object's redshift, and   the lookback time.The remaining period between 100 Myrs and 0.85  is evenly spaced in logarithmic time.The number of bins is selected such that log 10 (Δ  /GYr) > 0.02 for each bin.The number of bins ranges from nine (for the lowest redshift sources) to six (for the highest).Following Smith et al. (2012), an SED fit is deemed acceptable if the best-fit  2 is below the 99 per cent confidence threshold for the given number of photometric bands5 .
Tests show that with this setup, Prospector produces good agreement with the spectroscopic SFRs (obtained from dust extinction and aperture corrected H emission line measurements) for the small minority of sources in the MPA-JHU catalogue (Brinchmann et al. 2004, see Appendix A1), and a plausible distribution of derived galaxy properties in the SFR-stellar mass plane relative to the socalled galaxy "main sequence" (e.g. the presence of a main-sequence and distinct clumps of starbursts and quiescent galaxies, as expected from works such as Noeske et al. 2007 andSchreiber et al. 2015), which is shown in Appendix A2.We intend to further investigate the agreement between Prospector and spectroscopic SFRs in a future work.Table 1 summarises the key Prospector model parameters and the prior distributions that we adopt for this work.
In Figure 2, we present the SED output and the recovered SFH for two radio-selected galaxies: one exhibiting strong MIR emission, due to the presence of a radiative-mode AGN, and another without such emission.Both of these sources were fit twice -once with an AGN component included in the fitting process, and once with any possible AGN contribution neglected.It is evident that the source with a strong MIR emission is better fitted when an AGN template is included in the model, however the impact of the inclusion of AGN on the recovered SFHs is more complex (visible in the top right panel of Figure 2).On the other hand, for the source without an obvious MIR excess, the SED fits as well as the recovered SFH from the two sets of fits are virtually indistinguishable whether an AGN component is included or not.

Comparison of physical parameters
The principal physical quantities of interest for this work are the current stellar mass, the SFR and specific SFR (each averaged over the last 10 and 100 Myrs), and the dust luminosity computed between 8 − 1000m.The current stellar mass refers to the surviving stellar mass after accounting for losses during evolution (e.g.due to AGB winds, supernovae) and is different from the total stellar mass formed which is a free parameter in Prospector models.Throughout this work, the term "stellar mass" will refer to the current stellar mass, unless specified otherwise.For each of these parameters, we construct marginalised PDFs by weighting the output Dynesty samples (analogous to Markov chains) by their corresponding importance weights.The median likelihood estimates and the 16 th and 84 th percentile bounds are computed from these PDFs.We compare these parameters to those obtained using the popular SED fitting code Magphys, for sources in the radio selected sample for which both Magphys and Prospector produced fits that are considered acceptable following the  2 based criterion outlined in Smith et al. (2012).As noted above, the Magphys SED fitting of LoTSS Deep Fields galaxies was done by S21 and B23, and we refer the reader to these for a detailed account of the process, highlighting here some key differences.For the stellar SEDs Magphys uses a library of models based on 50,000 exponentially declining SFHs with random bursts superposed, assuming the Chabrier (2003) IMF.It uses the two-component Charlot & Fall (2000) dust model and the principle of energy balance to model the dust emission, neglecting any possible contribution of AGN.
In Figure 3, we compare the median likelihood estimates of stellar mass, SFR (averaged over the last 100 Myrs), sSFR, and dust luminosity estimated using Magphys and Prospector.SFRs have been converted to the Kroupa (2001) IMF where necessary, using the corrections provided in Madau & Dickinson (2014).These plots are colour-coded by the 16 th percentile value of the AGN fraction estimated using Prospector, to emphasize the impact of including AGN templates in the models.The ability of the AGN fraction parameter in SED fitting codes to distinguish between AGN dust emission and dust heated by star formation is strongly dependent on the MIR photometric coverage as well as on the S/N ratio of the MIR observations (Leja et al. 2018).Therefore, it is not uncommon for estimated AGN fractions to have large associated uncertainties (e.g.Ciesla et al. 2015).Following B23, we found the 16 th percentile of the posterior of the AGN fraction to be a more reliable and robust indicator of AGN activity.Therefore, throughout this work, the term AGN fraction will be used to denote the 16 th percentile values (  AGN,16 ). Figure 3 shows that estimates from Prospector and Magphys are in good general agreement given the uncertainties (sources with high  AGN,16 aside).For the majority of sources, the stellar masses, SFRs and dust luminosities estimated by the two codes are within 0.3 dex of each other.The Prospector results show that a large number of sources with high median likelihood stellar masses ( star > 10 11 M ⊙ ), SFRs (> 10 2  ⊙ /yr), and dust luminosities ( dust > 10 12  ⊙ ), also have high AGN fractions, leading to lower Prospector estimates of stellar mass, SFR, and dust luminosity than Magphys (which treats stellar emission as the sole origin of emission).Similar trends were observed by Thorne et al. (2022a), who used a method based on Prospect AGN fractions6 to identify AGN hosts.They noted a reduction of up to 2 dex in the estimated SFRs for sources with high Prospect AGN fractions when modelled with AGN templates, as opposed to when these sources were fitted without AGN templates.Additionally, Prospector returns lower median likelihood SFRs with larger error bars than Magphys for galaxies with low SFRs (less than 1  ⊙ yr −1 ).While it may seem tempting to consider these larger error bars as a weakness of using Prospector with non-parametric SFHs, they are in fact one its major strengths, since previous works (e.g.Pacifici et al. 2023;Haskell et al. 2023a) have shown that Magphys uncertainties are likely to be underestimated for such galaxies.
To determine the degree of consistency between our work and B23, in Figure 4 we also compare the median likelihood stellar mass and SFR estimates with the "consensus" measurements obtained by B23 (derived for each galaxy by combining the results from four SED fitting codes, as discussed in Section 1).We find remarkable agreement, similar to the comparison with Magphys estimates in Figure 3, perhaps unsurprisingly since the Magphys estimates (or a median of the Magphys and Bagpipes estimates) are adopted for the large majority of sources not classified as containing AGN by B23.Nevertheless, the good agreement for the high-stellar mass and high-SFR sources where AGN are more prevalent is very encouraging.The principal apparent disagreement between the B23 consensus values and the Prospector estimates appears in the sources with log 10 ( B23 / ⊙ yr −1 ) ≲ 0, however this 'difference' is not significant once the error bars are considered.Additionally, Prospector SED fits indicate that a lower number of radio-detected sources have extremely high SFRs (149 sources with SFR > 10 3  ⊙ yr −1 ) when compared to both Magphys SED fits (684) and B23 consensus estimates (418).This potentially relieves the tension between the SFR function as derived in works such as Smit et al. (2012), Duncan et al. (2014) and Katsianis et al. (2017), and that calculated in recent studies using SFRs estimated from SED fits (which predict an excess of extremely star forming galaxies, e.g.Gao et al. 2021).

RADIO SOURCE CLASSIFICATION
In this section, we look into the classification of radio detected sources into the underlying SFG, RQ-AGN, LERG, and HERG population classes.This involves two distinct parts: (i) identifying those sources hosting a radiative-mode AGN (Section 4.1) and (ii) identifying radio-excess sources (Section 4.2).In Section 4.3, we combine the two criteria to give the overall classification scheme and discuss the characteristics of the various classes.

Identification of radiative-mode AGN
Radiative-mode AGN are typically identified using optical emissionline based diagnostics (e.g.Baldwin et al. 1981;Veilleux & Osterbrock 1987;Kauffmann et al. 2003;Best & Heckman 2012;Tanaka 2012;Comerford et al. 2022) or mid-IR colour-colour cuts (e.g., Lacy et al. 2004;Stern et al. 2005;Donley et al. 2012;Messias et al. 2012).However, these methods are fraught with several caveats.Observations of the Baldwin-Phillips-Terlevich (BPT; Baldwin et al. 1981) emission lines are strongly affected by dust obscuration as well as by redshift effects.Thus the AGN identification methods that rely on optical emission lines are heavily biased against dust-obscured and low luminosity AGN (Padovani et al. 2017;Assef et al. 2018).Furthermore, spectroscopic data are not widely available; B23 report that only 5.1 per cent of the radio-detected ELAIS-N1 deep field sources have associated spectroscopic information, though the WEAVE-LOFAR survey (Smith et al. 2016) will soon offer complete spectroscopic coverage of these sources7 .Irrespective, selection criteria based on mid-IR colour-cuts often suffer from low signal-tonoise measurements and sample incompleteness, especially when dealing with fainter galaxies.Source samples classified using these broad colour-colour cuts are also likely to be contaminated by a large number of higher-redshift inactive galaxies (Barmby et al. 2006;Georgantopoulos et al. 2008;Gürkan et al. 2014;Donley et al. 2012;Messias et al. 2012).
SED fitting can account for the possible contribution of an AGN whilst simultaneously modelling the host galaxy stars and dust.The AGN fraction parameter in Prospector directly estimates the fractional contribution of AGN to the total mid-IR (4−20m) luminosity, and SED modelling can be used to identify AGN at lower luminosities than mid-IR colour cuts.For example, Leja et al. (2018)  sample by employing a similar method based on Prospect AGN fractions.The top panel of Figure 5 shows that the distribution of  AGN,16 for the radio-detected sample has a prominent peak below 10 −3 and a long tail above this value.Roughly 92 per cent of the sources have  AGN,16 less than 10 −3 .The tail of this distribution is primarily composed of sources that were marked as AGN by Duncan et al. (2021) on the basis of previous optical, spectroscopic, or X-ray studies.For this latter group, the distribution of  AGN,16 is bimodal, with a minimum falling between 10 −3 and 10 −2 .We identify those sources with  AGN,16 > 10 −3 as radiative-mode AGN hosts, noting that the total number of SFGs and AGN we classify in this way depends only weakly on the chosen value for  AGN within the limits of 10 −3 and 10 −2 .Leja et al. (2018) and Thorne et al. (2022a) demonstrated that incorporating AGN emission in the galaxy model produces better fits (i.e., lower  2 ) for sources with a significant AGN fraction compared to fits without AGN.Therefore to further improve our ability to identify sources hosting radiative AGN, we produce Prospector fits with and without the inclusion of an AGN component in the model.Figure 2 presents an example where a source with excess MIR emission is better fit when AGN templates are included in the model.Inclusion of AGN templates in the model led to improved fits for 66 per cent of the radio-detected sources (middle panel of Figure 5), and 96 per cent of the radio-detected sources which were previously classified as AGN by optical, spectroscopic, or X-ray studies were better fitted when AGN templates were incorporated into the fitting process.
To systematically identify the sources hosting radiative AGN, we consider the distribution of  2 No AGN / 2 AGN (where the two terms denote the goodness of the best fit SED omitting and including AGN, respectively), which is shown in the lower panel of Figure 5.For values below unity, this distribution can be reasonably well-fit by a Gaussian with standard deviation  = 0.057.The sources that show a significant improvement in fits after inclusion of AGN templates are identified as AGN.A fit is considered to have significantly im- No AGN exceeds  2 AGN by a factor of 1.17, corresponding approximately to 3.
Combining the two criteria, we classify a source as a radiativemode AGN if either of the following two conditions is met: (i) If the fit including AGN in the model is considered acceptable according to the  2 criterion described in Section 3, and  AGN,16 > 10 −3 , or, (ii) If  2 No AGN > 1.17 2 AGN and  AGN,16 > 10 −3 .We flag sources that were not fitted by Prospector (owing to either a lack of redshift information or not satisfying the recommended cuts in the LoTSS Deep fields catalogue, following Duncan et al. 2021) as unclassified.Applying the aforementioned criteria, we identify 3,925 (±58) out of the 30,470 radio-selected sources as radiative-mode AGN hosts (where the quoted uncertainties have been obtained by using bootstrap resampling from among our catalogue to generate 10,000 realisations, and calculating the median, 16 th , and 84 th percentile bounds of the resulting PDF).In this way, we have identified 82 per cent of the 503 optically-/spectroscopically-/X-ray -identified AGN and 73 per cent of the 1540 AGN selected through the Donley et al. (2012) colour cuts8 .To assess the performance of our classification scheme, we follow the radiative-AGN selection criteria laid out in B23 and use the data from their Magphys, Bagpipes, Cigale, and AGNfitter catalogues to reproduce their classifications.B23 identified 3,129 out of the 30,470 sources as radiative-mode AGN, identifying 83 per cent and 75 per cent of the optically/spectroscopically/X-ray-identified and Donley et al. (2012) AGN, respectively.Although we identify 3,925 sources as hosting radiative AGN as opposed to the 3,129 sources identified by B23 (i.e., ∼ 25 per cent more sources), 84 per cent of the radiative-AGN classified by B23 are also identified as such by our method 9 .It is clear that the results obtained by our method to identify radiativemode AGN, using a single SED fitting code, are comparable to those obtained using four SED codes in B23. 9 If we adopt  AGN, 16 > 10 −2 the corresponding value is 80 per cent.

Identification of radio-excess sources
One of the key tools underpinning our ability to produce robust estimates of SFR and classify radio sources using SED fitting is the relationship between SFR and radio emission in SFGs, thought to arise due to the acceleration of cosmic rays by supernovae.This is further underscored by the well-studied far-infrared radio correlation (FIRC) in SFGs, which has been observed to be linear over several orders of magnitude (e.g.van der Kruit 1971;de Jong et al. 1985;Helou et al. 1985;Yun et al. 2001;Jarvis et al. 2010;Bourne et al. 2011;Smith et al. 2014;Magnelli et al. 2015;Delhaize et al. 2017;Read et al. 2018;Molnár et al. 2021;Delvecchio et al. 2021;McCheyne et al. 2022).The precise nature of the FIRC, however, depends on a number of factors, such as the balance between the cosmic ray electron (CRE) escape fraction and the optical depth to UV photons (e.g.Lisenfeld et al. 1996;Bell 2003;Lacki et al. 2010).
Nevertheless, radio observations provide a potentially highly attractive means of estimating SFRs for galaxies up to high redshifts, and have enabled direct studies of the correlation between SFR and radio luminosity (e.g.Condon 1992;Cram et al. 1998;Bell 2003;Brinchmann et al. 2004;Garn et al. 2009;Kennicutt et al. 2009;Davies et al. 2017;Tabatabaei et al. 2017).Free-free emission is a significant contributor to the radio continuum at frequencies in the GHz range and above (Condon 1992;Murphy et al. 2011).In contrast, low frequency radio observations are largely unaffected by free-free emission.Therefore radio surveys at low frequencies, such as LoTSS, are ideal for studying the relationship between radio luminosity and SFR.
Works  2022) have used the sensitive LoTSS observations to study the so-called "SFR- 150 MHz relation".At 150 MHz, the primary source of radio emission in galaxies is the nonthermal synchrotron radiation arising either due to stellar processes or accretion activities within the central engine, including jets (Padovani et al. 2017).The latter are commonly exhibited by radio-loud active galactic nuclei (RL AGN; e.g.Urry & Padovani 1995;Wilson & Colbert 1995;Best et al. 2005;Harrison et al. 2014;Heckman & Best 2014;Hardcastle et al. 2019;Das et al. 2021).RL AGN can be identified as those sources with a radio luminosity that significantly exceeds what is expected from stellar origin alone (Hardcastle et al. 2016;Williams et al. 2018;Hardcastle et al. 2019;Smolčić et al. 2017b;Calistro Rivera et al. 2017).The SFR- 150 MHz relation can be used to predict the radio emission expected from star formation processes alone, if the SFRs are known.In this section, we use the method of S21 to determine the SFR- 150 MHz relation for our SWIRE-selected galaxy sample10 .The main steps of this approach are summarised below, though we refer the reader to S21 for a full discussion of the method.
It is crucial to account for the 150 MHz properties of the sources irrespective of whether or not they are detected in the LoTSS catalogue.Failing to do so risks biasing the results towards the radiobright sources.We therefore use the Prospector SED fits for the SWIRE-selected sources.This sample, as demonstrated in Section 2.3, includes sources that were not individually catalogued at 150 MHz by Sabater et al. (2021).It is therefore key to getting a true measure of the low frequency radio source population, given that the source counts of radio sources at low frequencies are dominated by low luminosity SFGs (e.g.Wilman et al. 2008;de Zotti et al. 2010;Williams et al. 2016;Calistro Rivera et al. 2017;Williams et al. 2019;Mandal et al. 2021;Hale et al. 2023).Excluding objects flagged as AGN by Duncan et al. (2021) and avoiding those for which Prospector fails to produce acceptable fits, we are left with 120,307 SWIRE-selected galaxies for our analysis.We first obtain the median likelihood SFRs (averaged over 100 Myr) along with asymmetric uncertainties derived using the 16 th , 50 th and 84 th percentiles of the derived SFR PDFs, which are equivalent to median ±1 in the limit of normally distributed data.Next, we use 150 MHz flux densities and uncertainties from the Sabater et al. ( 2021) catalogue where they exist, and supplement these for the remaining sources with flux densities (and associated RMS uncertainties) measured from the pixel corresponding to the source position in the ELAIS-N1 maps from LoTSS Deep fields DR1.For each source, we then construct a two-dimensional PDF in the SFR- 150 MHz parameter space with logarithmic axes: 70 equally spaced logarithmic SFR bins between −3 < log 10 (/ ⊙ yr −1 ) < 3 and 180 equally spaced log-bins between 17 < log 10 ( 150 MHz /W Hz −1 ) < 26.This PDF is populated for each source by creating 100 random samples from a normal distribution based on the median-likelihood Prospector SFR estimates, accounting for the asymmetric uncertainties.
For the 150 MHz luminosities we use 100 independent samples from a symmetric normal distribution in linear space, centred on the luminosity obtained from the flux density estimate (either from the catalogue, or from the pixel measurements) assuming the best redshift and a spectral index of  = 0.7.The standard deviation on the luminosity distribution is obtained by scaling the uncertainties on the flux density measurements such that the signal-to-noise on the derived luminosity is equal to that on the measured flux densities.The 2D PDFs for each source are then summed, to generate a PDF that accounts for the uncertainties in both the SFR and  150 MHz .To ensure that sources with low S/N and negative radio luminosities were included (since not doing so would censor the distribution and lead to skewed median-likelihood luminosity measurements), we follow S21 and arbitrarily assign samples with log 10 ( 150 MHz ) < 17 to the lowest  150 MHz such that each source is equally represented in the 2D stacked PDF.Similarly, we arbitrarily assign samples with log 10 (/ ⊙ yr −1 ) < −3 to the lowest log SFR bin.We then calculate the 16 th , 50 th , and 84 th percentiles of the  150 MHz distribution in each log SFR bin, with the 50 th percentile values denoting the median likelihood SFR- 150 MHz relation.
As discussed at length in S21, the SFR- 150 MHz relation has significant dispersion, which we must account for in order to identify those sources with a radio excess.To do this, we consider the difference between the 50 th and 84 th percentiles 11 of the  150 MHz 11 We have elected to use the 50 th and 84 th percentiles rather than the 16 th and 50 th percentiles since the distribution is asymmetric, with increased scatter distribution for log 10 (/ ⊙ yr −1 ) > 1.5, finding  = 0.289 dex.This is in good agreement with the results of Cochrane et al. (2023), who found a scatter of 0.3 dex by comparing the shapes of SFR functions for radio-and SED-estimated SFRs.We specifically avoid including sources with SFRs less than 1.5 ⊙ yr −1 to mitigate the impact of large uncertainties associated with the best-fit relation at low SFRs (e.g. Figure 6).Using this information, we identify sources that exceed the best-fit SFR- 150 MHz relation by 0.867 dex (equivalent to a 3 excess) as having significant radio-excess (this criterion is indicated by the blue dot-dashed line in Figure 7).Additionally, as seen in Section 3.1, SED fitting cannot accurately measure the SFRs for sources with very low SFRs, and thus the radio-excess classifications of these sources are likely to be unreliable.To avoid this issue, B23 classified 0.4 per cent of their sources which had SFRs < 0.01 ⊙ yr −1 as radio-excess only if their radio luminosities exceeded that expected for a source with SFR = 0.01 ⊙ yr −1 by 0.7 dex.If the radio luminosity of such a source was below that value but still surpassed the radio-excess threshold determined by the SED-fitted SFR estimate, then it was designated as Unclassifiable.We follow B23 and similarly mark 0.1 per cent of our sample as Unclassified.Among the 30,470 sources in the radio-selected sample, we identify 5,250 ± 66 sources as having a radio excess (B23 identified 4,786 sources).Uncertainty in the count of radio excess sources is calculated using bootstrapping, in a similar manner that was adopted for radiative-mode AGN identification in Section 4.1.The median likelihood SFR estimates from Prospector are lower than the consensus SFRs estimated by B23 for 92 per cent of the sources where the radio-excess classifications differ.In these cases, the radio-excess thresholds calculated for the Prospector results are lower than the thresholds based on B23's results.We therefore classify about 10 per cent more sources as radio-excess than B23.
We mark sources with unacceptable fits according to the  2 threshon the bottom side due to the lower signal to noise ratio of the fainter sources in the 150 MHz maps.old criterion as "Unclassified", since the Prospector model does not apply and we cannot use the SFR estimates to quantify the expected 150 MHz luminosity.This results in a higher number of sources being flagged as unclassified compared to B23; we have 1,532 compared to 240 (these numbers differ from those shown in Table 2 since there the numbers include sources irrespective of whether they have redshifts or are flagged in the cuts recommended by Duncan et al. 2021).Our classifications are in agreement with B23 for 77 per cent of the sources they classified as radio-excess and 90 per cent of the sources they identified as not having a radio excess.While there are non-negligible areas of implied disagreement among the radio-loud classification (e.g. the ∼ 23 per cent of the B23 radio-excess sources not identified using our method which are not apparent in Figure 7 since the B23 consensus estimates are not plotted), this type of classification is inherently uncertain: we do not know what the "correct" answer should be for these sources (e.g.Drake et al. in preparation).
Irrespective, the good general agreement in the radio-AGN classification between the two works offers significant encouragement.

Final classifications
Here, following B23, we combine the radiative-mode and radioexcess AGN classifications discussed in the previous sections to assign each source to one of the five classes: SFG, RQ AGN, LERG, HERG and Unclassified.Specifically: • sources that have been identified as neither radiative-mode nor having a radio excess are classified as SFG; • those that are identified as radiative-mode AGN hosts, without a radio excess, are classified as RQ AGN; • objects not identified as radiative-mode AGN hosts but showing a significant radio excess are classified as LERG; • hosts of radiative-mode AGN with significant radio excess are classified as HERG; and finally • sources for which we are unable to obtain an acceptable Prospector SED fit remain Unclassified.
The number of sources in each category is shown in Table 2.The number of star forming galaxies, radio-quiet AGN, lowexcitation radio galaxies, and high-excitation galaxies, identified in the ELAIS N1 deep field using our classification scheme, along with uncertainties calculated using bootstrap resampling.B23 classification results are also included for comparison.the degree of overlap between the sources identified following our classification scheme and those presented by B23 is illustrated in Figure 8.
Overall, our classification method identifies similar numbers of SFGs (6 per cent lower), RQ AGN (12 per cent lower) and LERGs (2 per cent higher) relative to the results of B23.However, there are larger differences apparent in the number of HERGs and unclassified sources, where our numbers are 68 and 203 per cent larger.The former change -in the number of HERGs -results primarily from the details of the precise definition of a radio excess that we discussed in Section 4.2 relative to B23 (and shown in Figure 7).This change represents an increase from 1.6 per cent to 2.7 per cent of the radio-selected catalogue, and though this increase is significant, HERGs remain a small minority of the overall source population, perhaps highlighting the difficulty of this task.In the case of B23, a source is flagged as Unclassified if it fails to produce an acceptable fit in multiple codes.Consequently, they report significantly fewer Unclassified sources compared to our results.While we do not formally classify the sources with unacceptable Prospector fits, roughly 70 per cent of the sources would be consistent with SFGs, 10 per cent with RQ AGN, 13 per cent with LERGs, and 3 per cent with HERGs, if we were to disregard the checks for fit acceptability.Interestingly, this is in good agreement with the final classification results, indicating that the ability of Prospector to produce an acceptable fit for a source is approximately independent of the underlying population class.This could suggest that the  2 based acceptability criterion we have used on the Prospector SED fits may be conservative, and we intend to investigate this possibility in a future work.
Similarly to B23, we further examine the demographics of the different AGN classes based on stellar mass, SFR, 150 MHz radio flux density, radio luminosity, redshift, and optical -band magnitude; these results are shown in Figure 9.In each of these figures, our results are shown as the solid lines (coloured according to the classifications as indicated in the legend) and with shaded areas according to the uncertainties implied assuming Binomial statistics (e.g.Cameron 2011).The corresponding values from the B23 catalogue are shown as squares with dotted lines of the relevant colour, enabling a direct comparison. 12he upper panels of Figure 9 show the distribution of the different classes based on the median likelihood estimates of stellar mass and SFR from Prospector.SFGs are the dominant population class for most of the low-to-intermediate stellar mass ranges.At higher stellar masses, the number of AGN sees a steep rise, led by LERGs.Additionally, LERGs dominate the source fractions at SFRs below ∼ 1  ⊙ yr −1 , supporting the view that they are hosted by massive, quiescent galaxies (e.g., Tasse et al. 2008;Smolčić 2009;Best & Heckman 2012;Heckman & Best 2014;Williams et al. 2018).Overall, the stellar mass and SFR demographics in Figure 9 are similar to those observed by B23, with the principal differences attributable to the larger number of Unclassified sources in this work (which become the dominant component at the highest stellar masses and SFRs).
The remaining panels in Figure 9 depict the fractions in each class as a function of 150 MHz flux density, luminosity, redshift and -band magnitude.The demographics of each source type as a function of the 150 MHz radio flux density are generally in agreement with the findings of B23.The LERG and HERG source fractions as functions of 150 MHz radio luminosity follow similar trends as seen in B23.We observe an abundance of LERGs at the highest radio luminosities.It must be noted here that the ELAIS-N1 field may not probe enough volume to reach the very highest radio powers where HERGs dominate.Nevertheless, at the very highest radio luminosities, the HERG fraction sees a steep rise, and the fractions of LERGs and HERGs appear to be roughly equal.This steep rise in HERG fraction strongly hints at HERGs becoming the dominant class at higher luminosities, which is consistent with expectations (e.g., Best & Heckman 2012;Heckman & Best 2014;Pracy et al. 2016).The population of RQ AGN remains relatively stable at around 10 per cent below 1 mJy and sees a steady decline above that.Our findings are consistent with B23's results and are lower than the observed fractions of approximately 15-20 per cent reported in studies such as Simpson et al. (2006) and Smolčić et al. (2017b).This difference is likely because of the RQ AGN typically having flatter spectral indices than SFGs (e.g.Gloudemans et al. 2021), and therefore as noted by B23, they are likely to be less prominent at the low frequencies probed by LOFAR than at GHz frequencies and higher.The transition from a sample dominated by star formation to one dominated by radio-loud AGN occurs around 1.6 to 2 mJy, which closely aligns with B23's finding of 1.5 mJy.As noted by those authors, this is consistent with the transition occurring around 200 − 250Jy at a higher frequency of 1.4 GHz, assuming typical radio spectral indices (Padovani 2016;Smolčić et al. 2017b).
We detect radio sources of each class at all redshifts and optical -band magnitudes.At  > 1, we observe lower SFG and RQ-AGN fractions compared to B23, but an increased fraction of LERGs.Our results suggest that SFGs dominate up to redshifts  ∼ 4, and are subsequently overtaken by RQ AGN.We observe a steeper decline in the SFG fraction with redshift compared to B23, whose classifications suggest a transition from a star-formation dominated sample to an AGN dominated sample at  ∼ 4.5.We note that the reduced SFG fractions relative to B23 could be caused by (i) a larger number of sources being classified as radio-loud AGN at high redshifts, owing to Prospector estimating lower SFRs for sources with extreme SFRs (often found at very high redshifts) than B23, therefore having lower radio-excess thresholds, and (ii) a larger number of potential SFGs being flagged as 'Unclassified' in this work on account of the unacceptable Prospector SED fit  2 .The LERG fraction peaks around  ∼ 2 and subsequently decreases, while RQ-AGN, HERGs and Unclassified sources all increase at  > 3.This hints at the existence of an evolutionary relationship between a more energetic, radiative-AGN dominated source population in the early stages of the universe and the jet-mode population at cosmic noon.We draw the attention of the reader to the increase in the fraction of unclassified sources with redshift, with approximately 20 per cent of sources labelled as unclassified at  = 5 in our results, and advise caution when using our classifications at high redshifts.
In Figure 10, we show the SFR vs stellar mass plane, populated by a shaded 2D histogram showing the number of SFGs at each point colour-coded according to the legend on the right-hand side, which is shared between both the top and bottom panels.In the lower panel, the SFRs have been calculated relative to the so-called "main sequence" relation from Schreiber et al. (2015), converted to our adopted Kroupa (2001) IMF.The typical scatter of ±0.3 dex is indicated by the shaded orange region.Overlaid on the SFGs are contours delineating the locations of the RQ AGN (purple) and LERGs (blue), while the HERGs are represented by cyan circles.The contour levels have been chosen such that they include 5, 50, and 90 per cent of the sample.The SFGs lie broadly around the main sequence as expected, although the observed scatter is slightly greater than the canonical ±0.3 dex reported in literature (e.g.Tacchella et al. 2015).As seen in the bottom panel, a considerable portion of the SFG distribution lies above the main sequence, indicating the presence of a population of starbursts in the radio-detected sample.The RQ AGN occupy a similar region of parameter space as the SFGs, albeit biased towards the high stellar-mass end of the distribution, consistent with our expectations that AGN predominantly reside in higher-mass galaxies (e.g., G18 and Kondapally et al. 2022).Notably, we find a significant population of LERGs (∼ 33 per cent) to lie on or above the main sequence, consistent with the distinct populations of starforming and quiescent LERGs discussed by Kondapally et al. (2022).At higher stellar masses, LERGs are hosted across all SFRs, whereas at lower stellar masses, only star-forming galaxies host LERGs.This is in agreement with the findings of Kondapally et al. (in preparation), and supports their hypothesis that different fuelling mechanisms may be prevalent in different classes of LERGs.Lastly, we observe that 58 per cent of the HERGs lie within the typical 0.3 dex range of the main sequence, with the significant majority of HERGs falling within a 0.6 dex scatter.Our classifications thus align with works such as Best & Heckman (2012), Best et al. (2014), Pracy et al. (2016), and Kondapally et al. (2022), who noted that HERGs were typically hosted by star-forming galaxies.

STELLAR MASS DEPENDENCE IN THE SFR-150 MHZ RADIO LUMINOSITY RELATION
As discussed in Section 4.2, LoTSS data have been extensively used to study the SFR- 150 MHz relation.Brown et al. (2017), Wang et al. (2019), andHeesen et al. (2022), found a super-linear relationship (slope greater than unity) between SFR and 150 MHz radio luminosity.G18 found an upturn towards larger radio luminosity at low SFRs, which is in disagreement with calorimetric arguments, and in the opposite direction to the downturn reported by Bell (2003).G18, Read et al. (2018), and S21 also found evidence for a strong stellar mass dependence in the SFR- 150 MHz relation.In this section, we seek to address the question of stellar mass-dependence using the SFR and stellar mass estimates from Prospector.

Determining stellar mass dependence using Prospector SED fits
For this analysis, we focus on the  < 1 SWIRE-selected galaxy sample with acceptable Prospector SED fits.We excluded sources that showed a significant AGN contribution (either radiative mode or radio excess, as defined in Section 4.3).We also removed the radio non-detected sources that were flagged as AGN using spectroscopy, X-ray detections, or mid-IR colour-colour cuts by prior studies (Duncan et al. 2021).After applying these cuts, we were left with a sample of 85,162 predominantly star-forming and passive galaxies.Following S21, we extend the stacked PDF method described in Section 4.2 to three dimensions.100 samples are generated for each source along the SFR, stellar mass, and  150 MHz axes.We adopt skewed error distributions in log space for SFR and stellar mass using the 16 th , 50 th , and 84 th percentile Prospector estimates. 150 MHz is sampled from a normal distribution in linear space centred at the maximum likelihood estimate and with standard deviation equal to the RMS.
We then create a three-dimensional histogram of the samples, summing over all the sources to create a stacked PDF.We use 50 equally spaced logarithmic stellar mass bins between 7.5 < log 10 (/ ⊙ ) < 11.8, 60 equally spaced logarithmic SFR bins between −3 < log 10 (/ ⊙ yr −1 ) < 3, and 180 equally spaced logarithmic  150 MHz bins between 17 < log 10 ( 150 MHz /W Hz −1 ) < 26 and calculate the median  150 MHz in each SFR and stellar mass bin.As before, the samples with  150 MHz < 10 17 W Hz −1 and SFR < 10 −3  ⊙ yr −1 are arbitrarily assigned to the lowest  150 MHz and log SFR bins, respectively.
Similar to G18 and S21, we find clear evidence of stellar mass dependence in the SFR- 150 MHz relation (Figure 11).In all the panels of Figure 11, only the bins populated by a minimum of 15 galaxies are displayed, taking into account that each source is sampled 100 times.The top-left panel of this figure shows the median-likelihood  150 MHz as a function of SFR in different stellar mass bins indicated by the colour bar to the right.Interestingly, this plot indicates that the relationship between SFR and  150 MHz becomes progressively sublinear as we move towards higher stellar mass bins.This observation may suggest the influence of non-calorimetric processes, such as the rapid diffusion of CREs within star-forming regions, as discussed by Berkhuĳsen et al. (2013) and Heesen et al. (2019), or perhaps the presence of undiagnosed low-luminosity contamination by AGN (as discussed in S21).Furthermore, in contrast to S21 our Prospector results indicate that the stellar mass dependence (as apparent by the difference between the most-and least-150 MHz luminous stellar mass bins), appears weaker at log 10 (/ ⊙ yr −1 ) > 0.5 than for the lower SFR galaxies.
The top-right panel of Figure 11 displays the median-likelihood  150 MHz as a function of stellar mass in different SFR bins.It is immediately clear that there is an increase in  150 MHz with stellar mass (consistent with S21 and Heesen et al. 2022).However, in keeping with expectations on the basis of the top-left panel, the gradient of the stellar mass-radio luminosity relation appears to flatten with increasing stellar mass (though the increasingly small stellar mass range probed for the highest SFRs precludes a definitive statement based on visual inspection alone).
To investigate these trends further, we fit our data using the stellar mass dependent form of the SFR- 150 MHz relation from G18: (2) We use emcee with 16 walkers and 10,000 iterations to calculate the uncertainties associated with these fits.The following best fit values are obtained: log 10  1 = 22.083 ± 0.004,  = 0.778 ± 0.004, and  = 0.334 ± 0.006.S21 used simulations to reveal small biases in these estimates, and suggested that the uncertainties are likely to be underestimated (with true uncertainties around   = 0.011,  log 10  1 = 0.016, and   = 0.037; correcting for these biases will be discussed at length in Shenoy et al. in prepara-tion).By way of comparison, S21 obtained best fit estimates of log 10  1 = 22.218±0.016, = 0.903±0.012,and  = 0.332±0.037,while G18 found log 10  1 = 22.13 ± 0.01,  = 0.77 ± 0.01, and  = 0.43 ± 0.01.While the exponent of SFR dependence () seen in our results is close to that found by G18, the exponent of the stellar mass dependence differs.On the other hand, we see a similar degree of stellar mass dependence to S21, but the SFR dependence differs between the two works.To determine whether incompleteness in our sample selection could be playing a role in our findings, we followed S21 and repeated our analyses using a 95 per cent stellar mass-complete sample of 64,838 sources derived using the method of Pozzetti et al. (2010), finding that our results are statistically unchanged.We therefore conclude that stellar mass completeness does not significantly impact our results.
In the bottom-left panel of Figure 11, we plot the stellar massadjusted radio luminosity (obtained by dividing the  150 MHz data shown in the top panels by the best fit stellar mass-dependent term in Equation 2) in an attempt to reveal the SFR- 150 MHz relationship with the stellar mass effects accounted for.At log 10 (/ ⊙ yr −1 ) > 0.5, the  150 MHz values in different stellar mass bins virtually collapse into a single line, while at lower SFRs there is still significant scatter, to a similar degree seen in S21.The differences observed between our results and those of S21, as well as the inferred stellar mass dependence in the SFR- 150 MHz relationship, may be partially attributed to the assumed SFH and the various explicit or implicit priors used in the SED fitting process.While the non-parametric continuity prior used in this work has been shown to be flexible enough for modelling a diverse range of galaxy SEDs (Leja et al. 2019;Johnson et al. 2021;Haskell et al. 2023b, Das et al. in preparation), it is still by no means perfect and may not be uniformly appropriate at all stellar mass scales (e.g.lower mass galaxies are expected to have burstier SFHs; Hopkins et al. 2023).Moreover, the apparent 'break' in the SFR- 150 MHz relation might suggest the requirement for additional physics.

The influence of passive galaxies
Our findings in Section 5.1 point towards a more pronounced dependence on stellar mass in the SFR- 150 MHz relationship at log 10 (SFR) < 0.5 ⊙ yr −1 .To investigate a possible cause, we repeat our previous analysis after excluding any quiescent galaxies, which we identify as follows.For each source, we calculate the SFR expected based on the parametrisation of the star-formation ratestellar mass relation from Schreiber et al. (2015), and identify those 7,081 galaxies with Prospector SFR lower than < 0.6 dex below the relation as passive galaxies13 .The results appear similar -and the best-fit parameters obtained using Equation 2 to fit the data after correcting for biases (log 10  1 = 22.085 ± 0.004,  = 0.783 ± 0.005, and  = 0.338 ± 0.006) are consistent with the values derived in Section 5.1.This is not unexpected since the star-forming galaxy population is vastly numerically dominant, and we are using medianlikelihood  150 MHz values in each bin.Excluding the 7,081 quiescent galaxies in this way results in an apparently reduced scatter in the 150 MHz radio luminosity-SFR relation.The apparent reduction in scatter, however, can be attributed purely to the smaller sample size (since fewer bins meet the minimum occupancy threshold of 15 galaxies, especially at low SFRs) rather than a difference in the relationship between luminosity and SFR for passive galaxies, for example.

CONCLUSIONS
In this work, we have used Prospector (Leja et al. 2017;Johnson et al. 2021) to estimate the physical properties of sources over 6 deg 2 of the ELAIS-N1 field.Prospector offers substantial advantages compared to other SED fitting codes, namely, the inclusion of AGN on an even footing with the stellar population synthesis, the implementation of non-parametric star formation histories, the modelling of nebular emission lines and the use of a dynamic nested sampling framework to efficiently sample multi-modal parameter spaces and produce realistic uncertainties alongside the source properties.The main conclusions of this work are as follows: (i) We successfully produce radio source classifications for 92 per cent of the 31,610 150 MHz sources in the ELAIS-N1 field of the LoTSS first data release.Each source is classified as either a star-forming galaxy, radio-quiet AGN, high-excitation radio galaxy (HERG) or low-excitation radio galaxy (LERG).
(ii) We have shown that our classifications are of a quality which is similar to those produced using four codes in a previous work.This was possible because the inclusion of AGN models in Prospectoron an even footing with the stellar population synthesis -enables us to build on previous work (e.g.Best et al. 2023) by identifying sources that are the likely hosts of radiative mode AGN, as well as those for which there is a radio excess above that expected on the basis of the Prospector SFR.
(iii) We used Prospector to estimate the physical properties of 133,000 3.6 m-selected sources, and revisited the relationship between a source's SFR and its radio luminosity.
(v) We found that, using Prospector to fit exactly the same input photometry as used by previous works (e.g.S21), the form of the stellar mass dependence in the SFR- 150 MHz relation differs from what had previously been observed using Magphys for the SED fitting.The Prospector results indicate a far-tighter relation, with lower scatter especially at log 10 (/ ⊙ yr −1 ) > 0.5 when the stellar mass effects are accounted for.This work has focused on SED modelling using only photometric data, and it is clear that even using all available codes and every photometric band there can still be significant doubt over what the "correct" classification is for any particular radio source.The most reliable radio source classifications require spectroscopy (e.g.Best & Heckman 2012, Drake et al. in preparation).In this context, it is timely that the WEAVE instrument (Dalton 2016;Jin et al. 2023) has started collecting data, and the initial batch of over a million spectra of LOFAR-selected sources will be available in the coming months as part of the WEAVE-LOFAR survey (Smith et al. 2016).This will ultimately include every source in the ELAIS-N1 field brighter than ∼100 Jy.Prospector is uniquely placed to deal with the WEAVE-LOFAR data since it can deal with photometry and spectra on an equal footing.We are hopeful that including the WEAVE spectra in our modelling will allow us to put stronger constraints on the physical parameters (e.g.SFRs) obtained for ELAIS-N1 sources, and significantly improve the accuracy of galaxy classifications.There are great opportunities for this kind of galaxy analysis given the range of new and forthcoming data sets in the era of massivelymultiplexed spectroscopic surveys such as 4MOST (de Jong et al. 2019), Subaru-PFS (Greene et al. 2022) and MOONS (Cirasuolo & MOONS Consortium 2016), for studying new samples from e.g.Euclid, LMT/TolTEC (Pope et al. 2019;Wilson et al. 2020) as well as the SKA (Rawlings et al. 2004).

APPENDIX A: ASSESSMENT OF PROSPECTOR DERIVED SFRS A1 Comparison with spectroscopically measured SFRs
To determine the reliability of SFRs obtained using SED fitting, we compare Prospector (Johnson et al. 2021) SFR estimates with those obtained from H emission line measurements.The latter, which accounts for the possible presence of dust, is widely regarded as the gold standard tracer of SFR up to 10 Myrs (see Calzetti 2013;Tacchella et al. 2022a).We use the SFRs measured from the dust extinction and aperture corrected Sloan Digital Sky Survey (SDSS DR8; Abazajian et al. 2009) line fluxes provided by the Max Planck Institute for Astrophysics and the Johns Hopkins University (MPA-JHU) group (see Brinchmann et al. 2004).We match the optical positions of sources in our SWIRE-selected sample with the MPA-JHU catalogue using a nearest-neighbour algorithm with a maximum 1 arcsec separation, finding 241 common sources.Our analysis is limited to 96 sources not flagged as stars by the Schlegel classification (Bolton et al. 2012), with significant (> 3) H detection (allowing reliable dust corrections) and acceptable Prospector fits.For this test, we average Prospector SFRs over 10 Myrs to enable a likefor-like comparison given the expected timescale associated with the Balmer line emission (e.g.Kennicutt & Evans 2012).Figure A1 shows excellent consistency between the Prospector and Balmerline SFRs.We obtain a best-fit relation of log 10 ( Prospector ) = log 10 ( H ) − (0.03 ± 0.03), as well as a reduced  2 ≈ 1 (where  =  Prospector −  H ;   represents the Prospector and H SFR uncertainties added in quadrature).

A2 Star forming main sequence
To further assess the performance of Prospector, in Figure A2 we show a comparison between the location of our SWIREselected sample in the stellar mass-SFR plane, derived using the Prospector (top) and Magphys (bottom) SED fits.We focus on the 120,307 sources from the complete SWIRE-selected sample for which Prospector produced acceptable fits.To produce the PDFs, we follow an approach similar to the one used in Section 4.2, wherein we generate 100 samples from the stellar mass and SFR distributions for each source using 50 bins evenly spaced between 7 < log 10 ( ★ /  ⊙ ) < 12, and −3 < log 10 (/ ⊙ yr −1 ) < 3, respectively.Sources with log 10 (/ ⊙ yr −1 ) < −3 were arbitrarily assigned to the lowest log SFR bins.The PDFs shown are derived by summing the individual stellar mass-SFR PDFs over the whole sample.
In both panels we have overlaid a shaded region enclosing the location of the "main sequence" (hereafter MS) expected for star forming galaxies at 0 <  < 1 using the Schreiber et al. (2015) parameterisation, and converted to our adopted IMF.In each case, a distinct MS is evident, accompanied by a region of quiescent galaxies (i.2015) parametrisation of the main sequence of star forming galaxies between 0 <  < 1, converted to our adopted Kroupa (2001) IMF. et al. 2007;Elbaz et al. 2007;Whitaker et al. 2012).Together with the results in Appendix A1, these results offer significant encouragement that our Prospector SED fitting is producing realistic results.
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.(Left panel) 150 MHz radio luminosity coverage of the radiodetected and the  < 1 SWIRE-selected source samples used in this work.(Right panel) Histogram showing the number count of sources binned by 150 MHz luminosity for the two data samples.Sources with radio luminosity less than 10 17 W Hz −1 (including radio non-detections) were arbitrarily assigned to the lowest radio luminosity bin at log 10 ( 150 /W Hz −1 ) = 17.The blue data points below the lower limit of the purple distribution represent sources with < 5 150 MHz radio flux measurements.Including these sources is key to leveraging the low S/N measurements in a statistical manner.

Figure 2 .
Figure 2. SED fits (left-hand panels) and recovered star formation histories (right-hand panels) for two example radio sources.ILTJ161528.25+544449.9, at a redshift of 1.8, with an estimated log 10 (  AGN,16 ) = 0.32 (indicating the likely presence of radiative-mode AGN), is shown in the top panels.Meanwhile, ILTJ161146.91+553623.4,at  = 0.4 with log 10 (  AGN,16 ) = −4.33 (implying a likely lack of radiative-mode AGN activity), is displayed in the bottom panels.In left-hand panels the photometry is shown with black dots and error bars, while solid lines represent SED fits randomly drawn from the Dynesty chains.In all panels, red lines indicate fits with an AGN component included, while blue lines refer to fits without AGN.
flux densities and uncertainties for the 11 per cent of the SWIRE-selected sample that have 150 MHz counterparts in the Kondapally et al. (2021) catalogue.For the remaining sources, we follow S21 and use the flux densities measured at the 150 MHz image pixel coordinates corresponding to each source.For unresolved sources, these values represent the maximum likelihood estimate of the integrated flux density.Figure 1 illustrates the different parameter spaces occupied by the two source samples, with the Radio-selected sample shown in purple, and the SWIRE-selected sample shown in light blue.The marginal histogram to the right of Figure

Figure 3 .
Figure 3. Median likelihood estimates of stellar mass (top-left), star formation rate averaged over the last 100 Myrs (top-right), specific star formation rate averaged over the last 100 Myrs (bottom left), and dust luminosity (bottom-right) with associated uncertainties from Magphys (plotted along the axis) compared with those obtained from Prospector (plotted along the -axis).The sources in each plot are colour coded by the logarithm of the 16 th percentile of AGN fraction estimated by Prospector.At the bottom of each panel, we plot the difference between the Magphys and Prospector estimates (Δ = log 10 ( Magphys,median ) − log 10 ( Prospector,median ), where  represents the physical parameter estimates).In the insets we show the distribution of this difference (Δ), with the blue histogram representing sources with log 10 (  AGN,16 ) < −2 and the orange filled histogram representing sources with log 10 (  AGN,16 ) > −2.

Figure 4 .
Figure 4. Comparison of the median likelihood estimates of stellar mass (left-hand panel) and the average star formation rate in the last 100 Myrs (right-hand panel) from Prospector with the consensus estimates from B23.The error bars along the -axis represent the uncertainties in the Prospector estimates; however, uncertainties are not available for the consensus estimates from B23.The difference between the B23 consensus and Prospector estimates (as the Δ parameter, see Figure 3) are plotted at the bottom of each panel.Histograms depicting the distribution of the Δ parameters are shown in the insets.The blue histograms represent sources with log 10 (  AGN,16 ) < −2, while the orange filled histograms represent sources with log 10 (  AGN,16 ) > −2.

Figure 5 .
Figure 5. (Top panel) Histogram showing the distribution of log 10 (  AGN,16 ).(Middle panel) Comparison of  2 arising from Prospector fits with AGN component included ( 2 AGN ) along the -axis against the  2 from fits without AGN ( 2 No AGN ).Sources classified as AGN using previous optical/spectroscopic/X-ray studies are represented as red points.(Bottom panel) Histogram showing the distribution of  2No AGN / 2 AGN .A majority of sources classified as AGN in previous works produce better fits when AGN templates are included in the physical model.

Figure 6 .
Figure 6.Stacked two dimensional heatmap showing the SFR- 150 MHz plane, constructed by summing over the PDFs of sources in the SWIREselected sample.The thick red curve represents the median likelihood SFR- 150 MHz relation.The accompanying thin red curves indicate the 16 th and 84 th percentiles of the  150 MHz distribution at each SFR, respectively.The blue and green solid lines represent the single slope and broken power law relations from G18.The black and magenta solid lines represent the best fit relations from S21 and B23 respectively.

Figure 7 .
Figure 7. Maximum likelihood 150 MHz radio luminosities of the radiodetected sample plotted as a function of the Prospector estimated median SFRs.The best-fit SFR- 150 MHz relation (obtained from the SWIREdetected sample and therefore including radio non detections) and the radioexcess cut-off used in our work (blue) and in B23 (orange) are overplotted.

Figure 8 .
Figure8.Venn diagrams showing the number of sources in each class (starforming galaxies, radio-quiet AGN, low-or high-excitation radio galaxy) that have been identified by our method (pink), by B23 (yellow) and by both (overlap region).

Figure 9 .Figure 10 .
Figure 9. Sources classified as SFG (in black), RQ-AGN (orange), LERG (blue), HERG (purple), and the sources marked unclassified (in green), represented as the fraction of total sources.We show the distribution of sources in each galaxy class as a function of (top-left panel) Prospector median likelihood stellar mass estimates, (top-right panel) Prospector median likelihood SFR estimates, (centre-left panel) 150 MHz radio flux density, (centre-right panel) 150 MHz radio luminosity, (bottom-left panel) redshift, and (bottom-right panel) optical -band magnitude.The dotted lines with open squares represent the fractions of sources in each galaxy class from B23.In the top-left and top-right panels, we plot the B23 source fractions against the B23 consensus stellar mass and SFR estimates, respectively.The distribution of the whole sample as a function of the respective quantities is shown in the background, as a normalised grey histogram with peak being equal to unity.

Figure 11 .
Figure 11.The relationship between radio luminosity at 150 MHz and the dependence on stellar mass and star formation rate (SFR), derived for the 85,162  < 1 SWIRE-selected sources excluding known/likely AGN.Stellar mass and SFRs are estimated through Prospector SED fits.In the top left-hand panel, the median likelihood 150 MHz radio luminosity -SFR plot is shown for different stellar mass bins.The different coloured lines represent distinct stellar mass bins, as indicated in the colour-bar to the right.The top right-hand panel illustrates the median likelihood stellar mass vs 150 MHz radio luminosity plots for different SFR bins, with colours again relative to the bar to the right.Similarly, in the bottom left-hand panel, the stellar mass-adjusted radio luminosity -SFR plot is displayed for varying stellar mass bins, using the best fit values from Equation 2. Finally, the bottom right-hand panel exhibits the SFR-adjusted radio luminosity vs stellar mass plot for different SFR bins.

Figure A1 .Figure A2 .
Figure A1.Comparison of Prospector estimated median likelihood SFRs (averaged over 10 Myrs) with the Balmer corrected H SFRs.The solid black line denotes the ideal 1:1 relation.The red dashed line represents the best-fit deviation of the Prospector estimated SFRs from the H SFRs, while the bottom panel shows the difference of the Prospector SFR estimates of the individual sources from H SFRs, in units of the (propagated) uncertainties.

Table 2 , while
INSU, Observatoire de Paris and Université d'Orléans, France; BMBF, MIWFNRW, MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The Science and Technology Facilities Council, UK; Ministry of Science and Higher Education, Poland; The Istituto Nazionale di Astrofisica (INAF), Italy.This research made use of the Dutch national e-infrastructure with support of the SURF Cooperative (e-infra 180169) and the LOFAR e-infra group.The Jülich LOFAR Long Term Archive and the German LOFAR network are both coordinated and operated by the Jülich Supercomputing Centre (JSC), and computing resources on the supercomputer JUWELS at JSC were provided by the Gauss Centre for Supercomputing e.V. (grant CHTB00) through the John von Neumann Institute for Computing (NIC).This research made use of the University of Hertfordshire high performance computing facility and the LOFAR-UK computing facility located at the University of Hertfordshire and supported by STFC (ST/V002414/1), and of the Italian LOFAR IT computing infrastructure supported and operated by INAF, and by the Physics Department of Turin University (under an agreement with Consorzio Interuniversitario per la Fisica Spaziale) at the C3S Supercomputing Centre, Italy.For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.