© ESO 2021

1 Introduction

Combining extremely high sensitivity with a wide field-of-view, the Low Frequency Array (LOFAR; van Haarlem et al. 2013) offers an unprecedented capability for performing large statistical surveys of the radio sky. The LOFAR Surveys Key Science Project (Rottgering 2010) is undertaking a set of tiered surveys over the northern sky at 120–168 MHz. The first release of data from the all-sky LOFAR Two-metre Sky Survey (LoTSS; Shimwell et al. 2017) presented 424 deg² to an average sensitivity of 71 μJy beam⁻¹ at 150 MHz¹. Complementary to the wide area LoTSS data, the LoTSS Deep Fields First Data Release (LoTSS DR1) reaches radio continuum sensitivities comparable to or deeper than the deepest surveys currently available, and over orders of magnitude wider areas (> 50 deg² ; see Tasse et al. 2021 and Sabater et al. 2021 - Paper I and Paper II, respectively). Located in some of the best-studied northern extragalactic survey fields – Boötes, European Large Area Infrared Space Observatory Survey-North 1 (ELAIS-N1, or EN1 hereafter), and the Lockman Hole (LH) – the LoTSS Deep Field data reach a current rms sensitivity of ~32, 20, and 22 μJy beam⁻¹ at 150 MHz,respectively, sufficient to detect radio-quiet AGN and extremely star-forming galaxies out to the highest redshifts (z > 5).

Extracting the maximum scientific value from the LOFAR radio continuum observations requires robust identification of the host galaxies of radio sources, alongside knowledge of the source redshifts, to extract intrinsic physical properties for both the radio sources (e.g. physical size, luminosity) and their host galaxies. By design, the LOFAR deep fields are located in regions of the sky thatcontain extensive ancillary imaging data, from ultraviolet (UV) all the way to far-infrared (FIR). Kondapally et al. (2021; hereafter Paper III) present new multi-wavelength optical to mid-IR photometry catalogues for all three fields alongside careful cross-identification with the LOFAR radio source population (with host identifications for ≳ 97% of radio sources). However, the wide range of intrinsic host properties in radio continuum selected samples, including both active galactic nuclei (AGN) and extreme star-forming galaxies, means that even with extensive ancillary data, deriving reliable photometric redshifts (photo-zs) can be non-trivial (Norris et al. 2019).

In Duncan et al. (2018a, hereafter D18a), we demonstrated that photo-z estimates for radio continuum sources obtained by combining multiple template-fitting estimates can be both more precise and more reliable than those using just one single template library. However, there remained key subsets of the AGN population for which template-based photo-zs could not produce satisfactory results (e.g. IR or X-ray selected AGN at 1 < z < 3).

In Duncan et al. (2018b, hereafter D18), we built upon the previous template-fitting method of D18a by incorporating additional machine-learning-based photo-z estimates trained for different AGN subsets. When combined together within a hierarchical Bayesian (HB) combination framework, the resulting consensus photo-z estimates improve on the performance of either individual method (specifically template-fitting or machine learning, see also Cavuoti et al. 2017; Fotopoulou & Paltani 2018, for other successful approaches). This combined “hybrid” approach was successfully applied to shallow optical and mid-IR photometry over very wide fields as part of LoTSS DR1 (Shimwell et al. 2019; Williams et al. 2019). The resulting photo-zs – presented in Duncan et al. (2019a, hereafter D19) – provided precise and reliable photo-zs (≈ 3% scatter and < 2% outlier fraction; OLF) out toz ~ 0.8 for LOFAR radio sources with host-dominated optical-infrared spectral energy distributions (SEDs; specifically, sources for which there is no strong evidence for a dominant AGN contribution to the SED at X-ray, optical, or IR wavelengths). Beyond z ~ 0.8, the shallow optical data leads to incomplete samples and highly uncertain redshift estimates for this population in LoTSS DR1. Thanks largely to the extensive training samples available for machine learning estimates, reliable photo-z estimates for quasar-like sources in D19 extended out to z ≲ 3.

In the case of the Boötes field, a test field for both D18a and D18b, we have already demonstrated that our hybrid approach can produce high quality photo-zs across a broad range of source types and redshifts. However, Boötes is currently unique among the LoTSS Deep Fields with respect to its larger sample of high quality spectroscopic redshifts both for calibration of template fits and for training of machine learning photo-z estimates. Neitherthe EN1 nor LH fields contain the representative samples of AGN spectroscopic redshifts necessary to apply the hybrid method of D18b with the same effectiveness as achieved in Boötes. Until large samples of spectroscopic redshifts provided by the forthcoming WEAVE-LOFAR survey (Smith et al. 2016), a modified approach is therefore required if we wish to maximise the quality of photo-zs across all three LOFAR Deep Fields and fully exploit the extraordinary sample of faint radio continuum sources they provide.

In this paper we present photometric analysis of the matched aperture multi-wavelength optical to mid-IR data in each of the three LoTSS Deep Fields. The primary aim of the paper is to provide optimal photo-z estimates for all sources in the fields to enable the exploitation of the deep radio continuum observations. Additionally, given that all three fields have deep near-infrared (NIR) or mid-infrared (mid-IR) photometry required for robust SED fitting, we provide stellar mass estimates for the subset of the optically detected sources for which we can make reliable measurements – enabling a wide range of science and providing a valuable reference sample for studies of the faint LOFAR radio population.

The remaining sections of this paper are set out as follows. In Sect. 2 we summarise the properties of the multi-wavelength photometry and spectroscopic samples used for photo-z analysis in this study. Next, in Sect. 3 we outline the key changes to the hybrid photo-z method used in this analysis. In Sect. 4 we analyse the precision and accuracy of the resulting photo-zs across all three fields as a function of redshift, magnitude, source type and radio luminosity. Section 5 then presents the method used for deriving stellar mass estimates for the optical sample, along with derivation of mass completeness limits and tests demonstrating their overall quality and limitations. In Sect. 6 we then briefly explore the physical properties of the LoTSS Deep Fields radio source population based on our derived redshifts and stellar masses. Finally, Sect. 7 presents a summary of our work. Throughout this paper, all magnitudes are quoted in the AB system (Oke & Gunn 1983) unless otherwise stated. We also assume a Λ Cold Dark Matter cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3 and Ω_Λ = 0.7.

2 Data

2.1 Photometry

Optical to mid-IR imaging, photometry and radio source identification for each of the LOFAR deep fields are presented inPaper III, to which we refer the reader for full details. In summary, available deep photometry in the EN1 and Lockman Hole fields has been mosaicked onto a common pixel and flux scale. For each field, forced aperture photometry in all bands has then been performed based on detections in two stacked χ² S/N images (optical to NIR and mid-IR) and the resulting sets of catalogues have been merged to produce a consistent multi-wavelength catalog. In the Boötes field, existing forced aperture photometry catalogues based on detections in I and Spitzer Infrared Array Camera (IRAC; Fazio et al. 2004) 4.5 μm images have been merged following the same procedure as used for the other two fields. Aperture corrections based on either curve of growth analysis for sources in the field (EN1 and LH) or the analytic point spread function (PSF; Boötes) have then been calculated to provide total flux estimates in each filter.

For all three fields we make use of 3′′ apertures for optical to NIR bands and 4′′ apertures for IRAC (due to the lower resolution in the IRAC imaging). As described in Paper III, the galactic extinction (E(B− V)) for each source has been calculated based on the Milky Way extinction map of Schlegel et al. (1998), queried using the dustmaps Python package (Green 2018). Filter-dependent extinction factors are then calculated by convolving the respective filter response curves with the Milky Way dust extinction law of Fitzpatrick (1999). In Fig. 1 we visually summarise the wavelength coverage and sensitivity of the resulting photometry catalogues used for the photo-z analysis in this paper.

In addition to the processing steps above that are presented in Paper III, for the photo-z estimation and SED fitting work in this paper we include further processing as follows. Due to the combination of large areas, the number of unique bands and depth of the photometry, small numbers of spurious datapoints are inevitable. Therefore, as a final additional step before photo-z estimation, we automatically filter the photometric datasets for spurious datapoints. Specifically, we filter for unphysical colours or extreme datapoints that are likely to be caused by artefacts such as stellar diffraction spikes, cosmic rays, cross-talk etc. Measurements that lie >2.5 mag above or <1 mag below those at both adjacent filters at shorter and longer wavelengths, indicative of extreme excess or deficit, are excluded. The large magnitude cuts are chosen to not exclude strong emission or absorption line sources. Similarly, measurements that result in a > 5 magnitude jump or drop in colour between consecutive filters are also excluded. These cuts are designed to be conservative, excluding clearly unphysical colours whilst not affecting large but physical colours such as Lyman break or strong emission line features. We note that the full optical catalogues released in Paper III have not been processed inthis way as users may wish to apply different cuts. However, the exact photometric catalogues used for photo-z estimation (including outlier identification) are available through the LOFAR Surveys Data Release site².

There likely remain a number of spurious measurements for which more detailed analysis would be required for identification (e.g. iterative fits excluding individual photometry points, cf. Chung et al. 2014). However, given the large number of available bands and the overall high quality of imaging available in these fields, we do not expect that our photo-z estimates (or science results derived from them) are significantly affected by any remaining anomalous datapoints.

Fig. 1 Illustration of the wavelength coverage and depths for the photometric datasets used for photo-z estimation in the three LoTSS Deep Fields. The violin plot for each filter shows the distribution of magnitudes for 5σ sources. The widths of the violin plots are indicative of the area available at a given depth within each field and are scaled by the fraction of sources with 5σ detections compared to the total size of the catalog. For example, the deeper Spitzer IRAC photometry of the smaller SERVS regions (at 3.6 and 4.5 μm) are evident in the EN1 and LH plots.

2.2 Multi-wavelength classifications

We broadly classify all sources in the photometric samples using the following additional criteria:

Infrared AGN are identified using the Spitzer IR colour and monotonically increasing mid-IR SED criteria presented by Donley et al. (2012). As we are primarily concerned with identifying only the robustly selected AGN sources, we also require > 5σ detections in all four IRAC bands. For the subset of the optical catalogues which satisfy our primary quality cuts, we find 0.7, 0.3 and 0.3% of the catalogues are found to satisfy the IR AGN criteria in Boötes, EN1, and LH, respectively.

Optical AGN: were identified through cross-matching the optical catalogue with the Million Quasar catalogue compilation of optical AGN, which is primarily based on Sloan Digital Sky Survey (SDSS; Alam et al. 2015) and other literature catalogues (Flesch 2015). Additionally, sources flagged as AGN based on spectroscopic observations are included in the optical AGN sample. For the subset of the optical catalogues which satisfy our primary quality cuts, we find consistently that ≈ 0.1% of the optical sources in all three fields satisfy the optical AGN selection criteria.

X-ray AGN: in the Boötes field where deep X-ray observations are present over a large area, X-ray AGN were identified by cross-matching the positions of sources in our catalogue with the X-Böotes Chandra survey of NDWFS (Kenter et al. 2005). We calculate the X-ray-to-optical flux ratio, X∕O = log₁₀(f_X∕f_opt), based on the I band magnitude following Brand et al. (2006) For a source to be selected as an X-ray AGN, we require that an X-ray source has X∕O > −1 or an X-ray hardness ratio >0.8 (Bauer et al. 2004). In total, we find 2811 X-ray AGN in Boötes based on these criteria. In EN1 and LH, bright X-ray sources were identified based on the Second ROSAT All-Sky Survey (2RXS; Boller et al. 2016) and the XMM-Newton Slew Survey (XMMSL2)³, as in LoTSSDR1. X-ray sources were matched to their IR counterparts using the published AllWISE (Wide-field Infrared Survey Explorer, WISE; Wright et al. 2010) cross-matches of Salvato et al. (2017), with sources then matched to the deep fields photometricdataset using the corresponding AllWISE source positions.

We note here that as in previous works, these broad sample selections are designed to identify clear AGN dominated SEDs for the purposes of optimising the photo-z analysis. The classifications are not intended to be complete samples of the AGN population within the fields. Subsequent studies will combine the multi-wavelength photometry with radio and FIR information to provide robust source classifications for the LOFAR detected population (see Best et al., in prep.; Paper V).

2.3 Spectroscopic redshift samples

Spectroscopic redshifts for sources in Boötes are taken from a compilation of observations within the field comprising primarily the results of the AGN and Galaxy Evolution Survey (AGES; Kochanek et al. 2012) spectroscopic sample, with additional redshifts provided by a large number of smaller surveys in the field including Lee et al. (2012, 2013, 2014), Stanford et al. (2012), Zeimann et al. (2012, 2013) and Dey et al. (2016). Included in this sample are a number of proprietary redshifts used for spectroscopic redshift training only (M. Brown, priv. comm.). In the release catalogues provided in this paper, we incorporate publicly available spectroscopic redshifts compiled as part of the Herschel Extragalactic Legacy Project (HELP, PI: S. Oliver)⁴.

Our spectroscopic redshift samples for the EN1 and LH fields are also based on the HELP compilations. The majority of spectroscopic redshifts in these compilations originate from the SDSS spectroscopic sample (Alam et al. 2015), with additional data from a number of smaller spectroscopic follow-up campaigns (namely Berta et al. 2007; Swinbank et al. 2007a,b; Lacy et al. 2013).

For all three fields, the spectroscopic samples are matched to the optical photometry catalogues using a simple nearest neighbour match with a maximum radius of 1′′. Table 1 summarises the spectroscopic redshift samples available for each field and the corresponding redshift and magnitude distributions are presented in Figs. 2 and 3, respectively. The disparity in currently available spectroscopic samples between Boötes and the other two deep fields is evident, with much larger samples of low-redshift galaxies available and significant AGN samples that extend to higher redshift and fainter magnitudes than available in EN1 or LH. However, we can see from Fig. 3 that the available spectroscopic sources are biased towards brighter magnitudes in all fields. Nevertheless, a key goal of this work is to provide consistent photo-z estimates across the three fields in all parameter space – ideally of comparable quality.

3 Photometric redshift estimates

We estimatephoto-zs for the full matched aperture optical catalogues presented in Paper III following a modified version of the hybrid approach presented inD18a and D18b and applied to LoTSS DR1 in D19. For full details of the motivation behind the hybrid approach, as well as details regarding the HB combination method we refer the reader to those papers. In this paper we focus solely on the modifications or changes to the method presented in D19 that are specific to this work.

Fig. 2 Redshift distributions for the spectroscopic redshift training and test samples available in each of the fields (see Table 1 for total numbers). While Boötes has a much greater number of sources with zspec available, these are largely limited to z < 1 for galaxies but extend over a greater redshift range for the known AGN.

Fig. 3 Normalised optical magnitude distributions for the galaxy (blue histogram) and AGN (red histrogram) spectroscopic redshift samples in comparison to the full parent photometric sample. For EN1 and LH we plot the r-band magnitude distribution, while for Boötes we show the nearest equivalent wavelength, R. The offset in the magnitude distributions highlights the importance of high quality template estimates in the regime where spectroscopic training samples are not available.

3.1 Template fitting

Due to the limited availability of spectroscopic training samples in two of the deep fields (EN1 and Lockman Hole) and the scientific focus on higher redshift where spectroscopic redshifts will always be limited, the primary modifications to our hybrid photo-z method are designed to maximise the accuracy and reliability of our template-fitting based estimates for the AGN population. Our aim is to reduce the potential need for machine learning estimates and provide consensus photo-zs for AGN in EN1 and LH that are of comparable quality to those produced in Boötes (D18b). For the template fitting photo-z estimates, we use updated versions of all three template libraries that either improve the wavelength coverage or broaden the representation of different AGN SED types. The three updated libraries now used for template fitting are:

Updated EAZY models

In previous work, the “default” EAZY template set (version 1.3) was found to produce photo-z estimates with the smallest scatter and OLF for star-forming or quiescent galaxy populations D18a. Here we make use of the new template set derived from the Flexible Stellar Population Synthesis code (FSPS; Conroy et al. 2009; Conroy & Gunn 2010). While designed to reproduce the same representative combination of stellar emission dominated galaxy SEDs, the revised FSPS templates now include dust reprocessed emission in the mid-IR. The emission from dust extends the range of rest-frame wavelengths that can reliably be used to constrain the photo-z estimate, which is particularly valuable in the deep field datasets that includes photometry out to 8 μm.

Extended Atlas Library

The second template set, the “Atlas of Empirical SEDs” (Brown et al. 2014) was recently extended to incorporate a wide range of the AGN and quasi-stellar object (QSO) populations not previously present in the original library (Brown et al. 2019). Additionally, the new Atlas of AGN SEDs incorporates a range of combinations of Seyfert type AGN spectra and different underlying host stellar populations (with varying AGN to host contributions). As illustrated in Brown et al. (2019), the inclusion of a fully representative range of AGN SEDs (and combinations of AGN + host galaxy) leads to significant improvement in the photo-z statistics for AGN dominated galaxies – the area for which machine learning estimates were found to provide the greatest improvement in the hybrid method (D18b).

Revised “XMM-COSMOS” Team templates

Finally, we also make use of a new iteration of the “XMM-COSMOS” galaxy and AGN SEDs (Polletta et al. 2007; Salvato et al. 2008, 2011) as presented in Ananna et al. (2017). Optimised for significantly larger survey fields than previous iterations (comparable to the total LOFAR Deep Fields coverage), the implementation in Ananna et al. (2017) includes more luminous quasar SED types than previous versions – as would be expected in the larger survey volumes than the deep pencil-beam surveys for which the Salvato et al. (2011) library was optimised.

The implementation of all three libraries within our template fitting method is the same as their corresponding versions in D18a and D18b. One key exception is the incorporation of new template library specific rest-frame model uncertainties within the fitting process that we outline below.

3.1.1 Photometric zero-point offsets

The inclusion of small magnitude offsets, or zero-point offsets, to the observed photometry of some datasets has been shown to improve photo-z estimates from template fitting (e.g. see Dahlen et al. 2013). While typically small (≲ 10%), these additional offsets can often substantially reduce the overall scatter or OLF for photo-z estimates. In Hildebrandt et al. (2012), detailed comparison between the resulting photo-z performancefor different levels of photometry precision indicates that zero-point flux offsets serve largely to correct for PSF effects. Given the nature of the photometry used in the work, we would therefore expect that the inclusion of zero-point offsets during template fitting will be beneficial.

Zero-point offsets for all template sets are derived for each photometric dataset following the method outlined in D18a. In summary, for 50% of the spectroscopic redshift subset (with 50% retained for validation/testing), the template set is fit to the observed photometry with the redshift fixed to the true redshift and the corresponding zeropoint offset is then calculated from the median offset between the observed and fitted flux values for sources with S∕N > 3 in that band.

As previously seen for the Boötes photometry in D18a, the inclusion of the zero-point offsets during template fitting leads to substantial improvement in photo-z quality (as tested using the 50% of spectroscopic sources not included in the derivation of zero-point offsets). Similar improvements are observed for both EN1 and LH. For reference, we present the derived zero-point offsets for all template library and field combinations in Appendix B. Within a given field, we find good agreement between the largest zero-point corrections derived for the different template sets (e.g. the UKIRT Infrared Deep Sky Survey, UKIDSS, K-band data in EN1 and LH are consistently found to require a correction of ~10%). We conclude that for the datasets employed in this analysis, the zero-point corrections are likely correcting for small offsets in the photometry itself, rather than large systematic errors in the templates themselves (we note that the three different template sets encompass both model and empirically derived templates).

3.1.2 Template model uncertainties

Incorporating an estimate of the rest-frame model uncertainties during photo-z template fitting has been shown to significantly improve the accuracy of the resulting estimates (Brammer et al. 2008; Dahlen et al. 2013). As part of our revised photo-z template fitting procedure we therefore estimate the rest-frame uncertainties for each template library from the spectroscopic training sample following the method outlined by Brammer et al. (2008). For the Boötes field (chosen due to its large available spec-z sample), for each template set we fit the templates to the spectroscopic sample (including the zero-point flux correction) while fixing the fitto the known redshift and then measuring the distribution of flux residuals.

The measured rest-frame residuals for each template library and the corresponding 1σ ranges are illustrated in Fig. 4. The resulting template error functions for each library are then shown in Fig. 5. Consistent with previous measurements of the photo-z template errors (Brammer et al. 2008), we find that rest-frame UV and mid-IR wavelengths have the highest model uncertainties. We note however that model uncertainties estimated by this method are representative of the average uncertainty across the full template setgiven the spectroscopic sample available. As demonstrated by Brown et al. (2019), the model uncertainties for known AGNsamples can be substantially larger at certain rest-frame wavelengths. However, a full Bayesian template fitting framework that can account for variance associated with each individual template (see e.g. Leistedt et al. 2019) is not practical for this work.

In our implementation of the rest-frame model uncertainties for this study, we also identify and correct for a subtle but important statistical error in the version of the Brammer et al. (2008) photo-z code used. Specifically, by default EAZY simply stores the best-fit χ² at each step in the chosen redshift grid (). The redshift likelihood is then typically taken as (e.g. Dahlen et al. 2013; Finkelstein 2014), which may then be convolved with a redshift prior (e.g. a magnitude prior). However, when redshift-dependent model errors are included in the fitting, the simplifying assumption of dropping the normalisation in the Gaussian likelihood term is no longer valid since it is no longer a constant. This has the effect of the P(z) becoming biased towards redshift ranges where the template error itself is maximised. For bright sources with strong colour features, this has minimal effect, but for fainter sources this can significantly bias the P(z). Correctly accounting for the redshift-dependent model errors is trivial, with the maximum likelihood at a given redshift, , being: (1)

where for a given filter, i, σ_i (z) is the total error (i.e. ), f_i,obs the observed flux and the model flux of the best-fitting template at that redshift. Taking the logarithm of both sides, the log-likelihood can be separated into two terms: (2)

where the right-hand term can be rewritten as (3)

with the best-fit χ² at each redshift step of the template fitting (as calculated and stored by EAZY). When deriving the full photo-z posterior, we therefore calculate the first term of the right-hand side of Eq. (3) for every source after the fact and incorporate alongside magnitude priors.

For all three fields, we use an IRAC 4.5μm prior derived from the Boötes spectroscopic sample following the method outlined in Sect. 5.1.1 of D18a. The use of the magnitude prior for the photo-z estimates was found to improve overall statistical performance at z < 2 (reduced scatter), the key area of interest for deep fields.

Fig. 4 Rest-frame residuals for the template fits using each of the three SED libraries employed in this work. Background grey points correspond to individual data points (i.e. one for each fitted filter per source) while the solid coloured lines correspond to the median residual within a given wavelength bin and the shaded region corresponds to the 1σ distribution (16 and 84th percentiles).

Fig. 5 Rest-frame model uncertainties derived for each of the template libraries used during our analysis (solid coloured lines). The model uncertainties are calculated by subtracting in quadrature the scaled average fractional error (black dotted line) from the median absolute rest-frame residuals (coloured symbols). However, the contribution of the average error is negligible.

3.2 Gaussian process redshift estimates

Following previous iterations of our hybrid redshift methodology, we incorporate machine learning estimates derived using the Gaussian process redshift code GPZ (Almosallam et al. 2016a,b). As in D18b, when training the GPZ classifiers, we employ the magnitude and colour-based weighting scheme (based on the method presented in Lima et al. 2008) to benefit from GPZ’s cost-sensitive learning features. We train GPZ using 25 basis functions and allowing variable covariances for each basis function (i.e. the “GPVC” of Almosallam et al. 2016a). Finally, we also follow the practices outlined in Sect. 6.2 of Almosallam et al. (2016a) and allow pre-processing of the input data to normalise or de-correlate the features (also known as “sphering” or “whitening”).

A key change to our method from previous implementations is the choice of input magnitudes. While the use of standard logarithmic magnitudes has produced good photo-z results in our previous efforts (D19), they are however not suited to datasets containing large numbers of non-detections (i.e. flux measurements with S∕N < 2) or low S/N measurements. In multi-wavelength forced photometry catalogues such as those employed here where depth varies significantly between filters (Fig. 1), non-detections in individual bands are inevitable and therefore a number of sources may have zero or negative flux measurements and hence are undefined in standard logarithmic magnitudes. In many cases, most notably high-redshift galaxies, those flux limits provide valuable colour (and hence redshift) information.

We therefore make use of asinh magnitudes (Lupton et al. 1999, sometimes referred to as “luptitudes”) for GPZ analysis, since they are able to incorporate zero or negative flux measurements and therefore allow us to train using all available measurements⁵. For a given flux, f (with a flux zeropoint f₀ = 3631 Jy), we define the asinh magnitudes as (4)

where the softening parameter, b, for each band is derived from the median flux uncertainty on ≈ 5σ sources across the field (specifically, 4.95−5.05σ)⁶. For high S/N measurements, the choice of flux zeropoint ensures that the asinh magnitudes are equal to traditional AB magnitudes. As demonstrated in Buchs et al. (2019), photo-z estimates from asinh magnitudes are typically not sensitive to the details of the softening parameter used if the data are of similar depth across the field. For individual photometric bands within the deep fields optical datasets, this assumption is valid due to the relatively homogeneous depth within our datasets (the heterogeneity within the fields is largely a result of different spatial coverage of different bands).

Due to the differing range and sources (i.e. telescopes or surveys) of optical to mid-IR photometric imaging across the three fields, GPZ must be trained separately for each of the corresponding spectroscopic datasets. The primary GPZ classifier for each field is trained on optical sources that do not satisfy any of the AGN selection criteria – corresponding to the significant majority of both the training sample and photometric catalog. In Boötes, we also use additional GPZ classifiers trained on the IR, optical and X-ray selected AGN subsets. However, due to the lack of comparable quality X-ray imaging, in the EN1 and LH fields we use additional IR and optical AGN GPZ estimates only. As illustrated in D18a, sources can satisfy multiple AGN selection criteria. The resulting spectroscopic training samples for the different AGN selections are as follows: the IR, Optical and X-ray AGN sample sizes are 1936, 1750 and 1307 (2714 in total) for Boötes, while for EN1 and LH the IR and Optical AGN sample sizes are 346 and 563 (593), and 579 and 1152 (1182), respectively. For all training, each input sample was split at random into training, validation and test samples consisting of 80% (training), 10% (validation), and 10% (test) of the full sample, respectively.

The photometric bands used for GPZ training for each subset and field are listed in Table 2. The exact choice of bands used was based on balancing the maximum number of bands (and hence information) with the fraction of each field covered by the respective bands. We note that the relative merit of including specific bands in the GPZ estimates varies between fields due to the variation in depth and coverage of different bands (for example, less improvement was gained by including IRAC 5.8 and 8.0μm photometry for the galaxies sample in EN1 – see Table 2).

3.3 Calibration of photo-z uncertainty

As outlined in Sect. 3.1.2, a key goal of the modifications to our template fitting method is to incorporate model uncertainties. Nevertheless, additional calibration of the uncertainties on the resulting photo-z is still necessary. To quantify the over- or under-confidence of our photo-z estimates, we follow the method outlined in Sect. 3.3.1 of D18b (and originally proposed by Wittman et al. 2016) and calculate the distribution of threshold credible intervals, c, where the spectroscopic redshift intersects the redshift posterior. For perfectly accurate estimates of the uncertainties, the cumulative distribution of credible intervals, , should follow a straight 1:1 relation, i.e. a quantile-quantile (or Q − Q) plot. Curves that fall below this 1:1 relation indicate that there is overconfidence in the photo-z errors (i.e. the P(z)s are too sharp) while curves that fall above indicate under-confidence.

Following D18a, we scale the uncertainties on the template fitting estimate for a source i, such that (5)

where α(m) is a magnitude-dependent function following: (6)

with α(m) being constant (α_η) below some characteristic apparent magnitude, m_η, and following a simple linear relation above this magnitude (Ilbert et al. 2009). For Boötes, EN1, and LH we use the I, r_PS and r-band optical asinh magnitudes, respectively for calculating the magnitude dependence of the error scaling and assume a characteristic magnitude of m_η = 16. The parameters α_η and κ are then fit using the EMCEE Markov chain Monte Carlo (MCMC) fitting tool (Foreman-Mackey et al. 2013) to minimise the Euclidean distance between the measured and ideal distributions.

When calibrating the uncertainties produced by GPZ, we calculate the threshold credible interval following (D19) and scale the uncertainties as a function of magnitude in the same manner as for the template estimates: (7)

where α(m_i) follows the same functional form as Eq. (6). As above, the parameters for α(m_i) are optimised through MCMC minimisation of the difference between the measured and ideal distributions.

As highlighted in Fig. 3, the available spectroscopic sample is biased towards bright optical magnitudes. Even though our photo-z uncertainties are calibrated as a function of magnitude, the biased sample could result in improved photo-z uncertainty accuracy for the brightest sources at the expense of the faint population. In order to prevent the optimisation of the photo-z uncertainties being dominated by the most populous magnitude ranges, we calibrate the uncertainties using a magnitude balanced subset of the total spec-z population following the approach presented in D19. For both the AGN and galaxy samples separately, a subset of each training sample is created by randomly selecting up to 750 sources in magnitude bins (of width = 1 mag) over the range covered by the spectroscopic redshift subset. Calibration of the uncertainties is then done on two-thirds of this subsample, with the other one-third retained for testing. We note that this balancing of the spec-z sample is most important for the Boötes field where the larger samples of bright magnitude selected spec-zs could significantly bias the optimisation.

3.4 Hierarchical bayesian combination

To produce the final consensus redshift prediction for a given source, we use the HB combination method outlined by D18a (based on the method presented in Dahlen et al. 2013) and subsequently extended to hybrid GPZ + template estimates in D18b. In summary, the HB combination produces a consensus redshift prediction, P(z) from a set of n individual predictions while marginalising over the probability that any individual P(z) is incorrect. Hyper-parameters for the fraction of measurements that are bad, f_bad, and the relative covariance between the different estimates β, are optimised using training data to ensure that the posterior redshift distributions more accurately represent the redshift uncertainties.

As in D19, GPZ estimates are evaluated on the same redshift grid as used during the template fitting procedure. If a source does not have a photo-z estimate for a given GPZ estimator (either through not satisfying the selection criteria for a given subset or lack of observations in a required band) it is assumed to have a flat redshift posterior for that specific estimator. GPZ therefore contributes no additional information to the consensus HB estimates for these sources.

For the application in this work, we assume 0 ≤ f_bad ≤ 0.05 and 0 ≤ f_bad ≤ 0.2 for the galaxy and AGN subsets, respectively, and a flat prior on the redshift distribution for “bad” estimates. However, one further improvement to our method in this work is that we optimise the hyper-parameter β (i.e. the degree of covariance between different photo-z estimates) as a function of magnitude. This change is necessary to ensure that the resulting consensus photo-z posteriors provided accurate uncertainties across all magnitude ranges. Assuming a constant β for all magnitudes results in underestimates of the uncertainties for the brightest and faintest sources.

In the following section, we present analysis of the consensus photo-z estimates for each field.

4 Photometric redshift properties

In this work we have estimated photo-zs for all sources within the optical catalogues for which estimates could be made, including parts of the fields with very limited photometriccoverage (including regions outside the area where radio to optical cross-identifications have been performed in Paper III). However, due to the heterogeneous photometry coverage within the deep fields, the resulting photo-z quality will also vary across the field. In the following analysis we examine the photo-z quality within the core regions of each field where photometric data is relatively homogeneous. These areas broadly conform to those where GPZ estimates are available and where radio source cross-matching has been performed for the LoTSS Deep Field data (Paper III). Specifically, we apply the following cuts in each fieldbased on the flags in the input photometric catalogues of Paper III: In Boötes, we require that FLAG_DEEP≠0 to remove duplicate or masked sources fromthe I-band detected catalogue sources. In EN 1, we require FLAG_OVERLAP ≥6 to restrict the analysis to the region that includes both NIR and Spitzer mid-IR photometric coverage. Due to the smaller coverage of the NIR data relative to the optical and mid-IR coverage, in LH we do not require NIR data when analysing the photo-z in this field - restricting the analysis to regions with SpARCS optical coverage (FLAG_OVERLAP ≥ 6) which overlaps with the Spitzer imaging.

Additionally, in all fields we require FLAG_CLEAN = 1 to exclude sources within the optical bright star mask. Compared to the initial combined total of ≈7.2 × 10⁶ catalogue sources available across the three fields, these cuts reduce the number of sources with the most robust photo-z estimates to a combined ≈5 × 10⁶ sources.

Figure 6 presents a qualitative illustration of the final consensus redshifts after the error calibration for all input estimates and the tuning of the Bayesian combination hyper-parameters. We show the stacked redshift posteriors as a function of spectroscopic redshift for each field, with the spectroscopic samples separated into the galaxy or host-dominated (top row) and AGN subsets (bottom row). From Fig. 6 we can see that the photo-z quality for galaxies at z < 1 is excellent for all three fields. At higher redshifts, the spectroscopic samples become much more limited (especially outside of Boötes). However, from the limited number of spec-zs available in this regime, it is clear that the performance deteriorates between 1 ≲ z ≲ 1.5 for the host-dominated population, and gets substantially worse beyond z ~ 1.5. Examining the AGN population posteriors, we find that photo-z performance looks generally good out to z > 3 in all fields – albeit with significantly larger scatter than the host-dominated population. Few spectroscopically confirmed sources at z > 4 exist in the EN1 and LH but we find that our photo-z estimates in Boötes provide accurate estimates out to z > 6 (including the confirmed radio-loud quasar at z = 6.1, McGreer et al. 2006). Given that hybrid photo-z estimates in this regime are found to be dominated by the template estimates (see D18b), we would expect the performance in the EN1 and LH datasets to therefore be comparable in quality (or potentially better in the case of EN1). Nevertheless, we caution against the use of the photo-zs above z ~ 4 without careful analysis of the individual estimates.

4.1 Overall photo-z statistics

The z_spec versus z_phot distribution shown in Fig. 6 allows us to qualitatively assess the photo-z performance. However, a more quantitative analysis is required to enable comparison between fields and for defining appropriate selection criteria for future science exploitation. When calculating photo-z statistics, we use the median of the primary redshift peak (z_1,median) following D19: see Sect. 4.1 of that paper for details on how this value is defined. For our measure of robust scatter, we then use the normalised median absolute deviation, σ_NMAD, defined as: (8)

where Δz = z_1,median − z_spec. Similarly, we define the OLF as sources where (9)

as is common in the literature (e.g. Dahlen et al. 2013).

In Table 3 we present the overall σ_NMAD and OLF statistics derived for the galaxy and identified AGN subsets. To ensure a fair comparison across the three fields, for the galaxy subset we limit the statistics to those sources with z_1,median < 1.5 while for the AGN subsample we limit to z_1,median < 4.0. We note that these cuts are defined using z_1,median (as opposed to z_spec) to better represent the performance expected for science samples that will be defined based on their photo-z alone.

We find that all three fields exhibit very similar σ_NMAD and OLF metrics for both subsets, with less than 2% scatter and OLF for galaxies or host-dominated sources. As already seen in Fig. 6 the photo-z quality for the AGN subset is significantly worse, with OLF reaching ≈ 20% or higher. Since the available spectroscopic redshifts are biased to typically brighter galaxies and lower redshifts, these statistics are not fully representative of the true photo-z accuracy and reliability for the average galaxy population.

In Figs. 7 and 8 we present the resulting robust scatter and OLF as a function of redshift and magnitude, respectively. The resulting statistics provide a quantitative confirmation of the photo-z quality visible in Fig. 6, where we find that our photo-z estimates have excellent scatter <0.04 × (1 + z) and OLF ~5% for the galaxy/host-dominated population at z < 1 across all three fields. For all three fields, there is a deterioration in the photo-z scatter and OLF with increasing redshift. However, in LH and EN1, not enough spectroscopic redshifts are available above z ~ 1 redshift to assess the quality beyond this point. In Boötes, where spec-zs are available in this regime, we find that the trend in the statistics verifies the visual interpretation of Fig. 6 – whereby the measured scatter deteriorates above z ~ 1 before becoming substantially worse at z > 1.5.

In contrast to the strong correlation with redshift observed for the galaxy subset, the measured σ_NMAD for the identified optical, IR and X-ray AGN samples exhibits negligible evolution with redshift. We also find that the OLF appears to decline as a function of redshift. This trend is likely driven by biases within the spectroscopic sample (with the z > 2 training and test samples are typically dominated by optically bright QSOs) and the presence of the strong Lyman break feature redshifting into the u-band at z ~ 3, leading to improved template fitting estimates (see also D18b, for further discussion).

We find that the measured OLF for known AGN at a given redshift or 4.5 μm magnitude in EN1 and LH is higher than for Boötes. We note that this difference is likely driven by the demographics of the respective spectroscopic redshift samples available. The EN1 and LH spec-z samples are dominated by QSOs while the Boötes sample includes a large number of sources explicitly selected on X-ray or IR AGN criteria. In D19, photo-z performance for AGN selected as optical QSOs was measured to be significantly worse than for those selected based on the other AGN criteria.

Fig. 6 Stacked redshift posterior estimates for the galaxy (or host-dominated) population (top) and for the combined AGN selected population (bottom; IR, X-ray or optically selected) for each of the three fields. The spectroscopic redshift samples available become particularly sparse at z > 1 for the EN1 and LH fields. Grey dotted and dashed line correspond to ± 0.1 × (1 + zspec) and ± 0.2 × (1 + zspec), respectively. The colour scale is illustrative and represents the summed P(z) (normalised such that ∫P(z)dz = 1) for each zspec bin, more precise photo-z therefore result in darker peaks.

Fig. 7 Robust scatter (σNMAD; upper panels) and OLF (lower panels) for the consensus photo-z estimates (z1,median) as a function of spectroscopic redshift. Shaded regions illustrate the statistical uncertainties on the respective metrics derived from bootstrap resampling.

Fig. 8 Robust scatter (σNMAD; upper panels) and OLF (lower panels) for the consensus photo-z estimates (z1,median) as a functionof 4.5 μm magnitude. Shaded regions illustrate the statistical uncertainties on the respective metrics derived from bootstrap resampling. The magnitude range plotted for each field is representative of the magnitude distribution of the spectroscopic sample and is not indicative of the relative photometric depths (see Fig. 1).

4.1.1 Comparison with existing literature photo-z estimates

In D18b we compared the quality of photo-z produced by our hybrid methodology with existing estimates in the literature, finding that our results produced better statistics than redshifts presented by Brodwin et al. (2006), with the substantial benefits of the accurate redshift posteriors. Thanks to the improved template fitting methodology in this study, the results presented in this data release represent a further improvement in photo-z quality available for this field. Specifically, when calculating photo-z quality statistics for the same subset of spectroscopic sources as used in Table 3, the photo-z presented in D18b have σ_NMAD = 0.032 and OLF = 0.028 for galaxies andσ_NMAD = 0.12 and OLF = 0.31 for the known AGN. The results presented in this paper therefore represent up to a factor of ~ 2 improvement in redshift quality averaged over the available spectroscopic population (we note that the degree of improvement in parameter space not probed by the existing spec-z sample could differ substantially from this, both positively and negatively).

In the EN1 and LH fields, there are a number of recently published photo-z catalogues to which we can directly compare our results. As part of the Hyper Suprime-Cam Subaru Strategic Program Second Public Data-release (HSC PDR2; Aihara et al. 2019), photo-z estimates derived from the HSC optical photometry (g,r,i,z,y) are available for the EN1 field for much of the area covered by the LoTSS Deep Field optical catalogues. The HSC PDR2 photo-z estimates presented by Nishizawa et al. (2020, see also Tanaka et al. 2018) are based on two different approaches: empirical estimates using the “Direct Empirical Photometric method” (DEmP, Hsieh & Yee 2014), and template based estimates following the method presented in Tanaka (2015, MIZUKI).

To directly compare our results with the HSC PDR2 estimates, we cross-match the catalogues with a maximum separation of 1′′ and calculate the photo-z quality statistics for all three estimates (this work, HSC DeMP and HSC MIZUKI) for the subset of spectroscopic sources with measurements in both our results and at least one of the HSC catalogues. The resulting statistics, presented in Table 4, demonstrate that the estimates presented in this work are again comparable to or better than the existing literature. While the overall statistics for the galaxy/host-dominated AGN population offer only a small improvement on the excellent scatter and OLF of the HSC estimates, our estimates for the known AGN population offer significantly better precision and reliability. This improvement is even starker when restricting the analysis to the higher redshift AGN population (z > 1, for which we provide statistics in Table 4 in parentheses), where our results have lower σ_NMAD and OLF by up to a factor of ~3 compared with those available from HSC PDR2.

Additionally, Pforr et al. (2019) present template based photo-z estimates for 18 deg² of multi-wavelength optical photometry in the Spitzer Extragalactic Representative Volume Survey (SERVS; Mauduit et al. 2012) over five fields, including subsets of the EN1 and LH fields. Pforr et al. (2019) quote robust scatter and outlier statistics calculated using the same definitions as in this paper for their full spectroscopic test sample (and the subset with the best available photometry), finding σ_NMAD = 0.042 (0.037) and OLF = 0.105 (0.028) in the EN1 field and σ_NMAD = 0.067 (0.03) and OLF = 0.205 (0.048) in LH. As the photo-z estimates inPforr et al. (2019) are primarily designed for galaxies (with no AGN templates included in the fitting), we can compare these values with the statistics for the galaxy population presented in Table 3. It is clear that the results presented in this work represent a substantial improvement in overall photo-z quality.

Our improved photo-z precision and reliability for the EN1 and LH fields highlights the benefits of not just our photo-z methodology but also the advantages of the full aperture matched photometry across the UV to mid-IR regime provided by Paper III. In particular, the availability of NIR and mid-IR photometry is crucial for providing reliable estimates for the AGN population at higher redshifts (a key area of scientific interest for the LoTSS Deep Fields). Furthermore, the additional colour information provided by forced photometry across all available photometric bands is vital, offering huge improvements in photo-z estimates orSED fitting when compared to the cross-matched catalogues employed by Pforr et al. (2019, see also Nyland et al. 2017).

4.1.2 Photo-z properties for the LOFAR population

Finally, we explore the quality of the consensus photo-z estimates as a function of their radio properties. Due to the larger spectroscopic sample required to bin in multiple properties, we explore the statistics only in the Boötes field. However, given the similar performance and redshift or magnitude trends observed across the three fields we would expect any observed trends to hold for all fields. Figure 9 presents the σ_NMAD and OLF as a function of both spectroscopic redshift and 150 MHz radio luminosity. When converting from observed flux density to rest-frame radio luminosity, we assume an average spectral slope of α = −0.7 for all sources – consistent with the typical slope observed in previous studies (Calistro Rivera et al. 2017)

In D19 we observed a clear evolution in the scatter and OLF of the radio source population with z_spec, whereby the photo-z properties of the highest redshift sources are significantly worse than for sources with similar radio luminosity at low redshift. However, within a given spectroscopic redshift bin, we found no evidence for any significant trend with radio luminosity in either the scatter or OLF.

Although exhibiting a noisier evolution than observed in D19, our photo-z estimates follow the same overall trend as a function of redshift and radio power. Specifically, for a fixed redshift we do not observe a strong correlation in σ_NMAD or OLF as a function of radio luminosity. Averaged over all redshifts, the measured scatter and OLF increase with increasing radio power. However, this trend is clearly driven by the redshift distribution of spectroscopic sources available at a given radio power.

4.2 Accuracy of the photo-z uncertainties

In Fig. 10 we illustrate the accuracy of the final calibrated redshift posteriors for the AGN and galaxy subsets. Shown in both plots are the cumulative distribution () of threshold credible intervals, c, for bins of apparent optical magnitude (coloured lines). As outlined in Sect. 3.3, uncertainty calibration was performed as a function of I band for Boötes and the deepest available r-band for EN1 and LH. For all three LoTSS Deep fields, the uncertainties for the galaxy spectroscopic sample are very well calibrated, with the measured uncertainties lying close to the desired 1:1 relation across the full magnitude range.

For theAGN population, the uncertainties are also generally very well calibrated across the available magnitude ranges. However, in Boötes we find that the wings of our photo-z posteriors overestimate the true uncertainties for the very brightest (I < 19) and faintest (I > 23) AGN population. In these regimes the available spectroscopic training sample is much more limited and the optimisation of the uncertainties is therefore still dominated by the more numerous AGN types. When compared to the results obtained in D19, we find that the photo-z uncertainties presented in this work are more accurate (with closer to the desired 1:1 trend). However, we caution that this does not necessarily mean the photo-z posteriors at a given magnitude are more precise.

As with the training of the machine learning estimates, we note that the large spectroscopic training samples available in the forthcoming WEAVE-LOFAR spectroscopic survey will allow for significant improvements in the calibration of the photo-z uncertainties in future studies. Nevertheless, the analysis presented in this section verifies that the photo-zs provided forthe LoTSS Deep Field release are of a high quality and suitable for scientific exploitation.

Fig. 9 Robust scatter (σNMAD; top) and OLF (bottom) for the consensus photo-z estimate as a function of spectroscopic redshift and 150 MHz radio continuum luminosity in the Boötes field. Top and side panels: trends averaged over all redshifts and luminosities, respectively. For a cell to be plotted we require a minimum offive galaxies - for some redshifts (or luminosities) we only have the number statistics available to plot the statistics averaged over all luminosities (or redshift).

5 Stellar mass estimates

Thanks to the extensive panchromatic photometry available in each field (including deblended photometry for mid-IR to FIR), detailed physical properties can be derived for the LOFAR source population through SED fitting. Future papers will exploit this information to derive robust source classifications, AGN accretion modes and star-formation rates for the faint radio source population, however running such codes for the full optical catalogues (> 10⁶ sources) is impractical. Nevertheless, estimating the stellar masses of the full optically selected population is critical for understanding how radio AGN affects galaxy evolution. Here, we present stellar mass estimates (and rest-frame optical colours) using a simpler grid-based SED fitting approach that scales to the massive samples available across the three deep fields.

Fig. 10 Cumulative distribution of threshold credible intervals, c, () plot for the final consensus photo-z estimates for the galaxy (or host-dominated AGN) population (top) and the optical, infrared and X-ray selected AGN population (bottom) within each of the deep fields. Coloured lines represent the distributions in bins of apparent optical magnitude (I-band for Bootes, r-band for EN1and LH). Lines that fall above the 1:1 relation illustrate under-confidence in the photo-z uncertainties (uncertainties overestimated) while lines under illustrate over-confidence (uncertainties underestimated).

5.1 Spectral energy distribution fitting

Stellar masses are estimated using the PYTHON-based SED fitting code previously used by Duncan et al. (2014, 2019b). Composite stellar populations are generated using the stellar population synthesis models of Bruzual & Charlot (2003) for a Chabrier (2003) initial mass function (IMF), with the following assumptions:

Star-formation histories

Recent studies have shown that the widely used assumption of exponentially declining star-formation histories can lead to biases in the ages and star-formation histories derived from photometric SED fits, resulting in biased stellar mass estimates (see e.g. Pacifici et al. 2016; Carnall et al. 2019a; Leja et al. 2019, and references therein). Carnall et al. (2018) show that the more complex double power-law parametrisation provides sufficient flexibility to accurately describe the star-formation histories of a wide range of possible formation and quenching mechanisms. As a compromise between the tractability of fitting large samples and the optimal prior assumption on star-formation histories, we therefore define a grid of SFH based on the double power-law model with the priors on the range of power-law slopes and turnover ages taken from Carnall et al. (2019b). Specifically, we use a logarithmically spaced grid of 7 values from 0.1 ≤ α_i ≤ 1000 for both power-law slopes () and 7 turnover ages from 0.1 to 1 times the age at observation (totalling 343 star-formation history models for every age step). For a given source the time since the onset of star-formation is defined as the time since a fixed formation redshift of z_f = 20 at the closest redshift step in a grid from 0 < z < 1.5.

Stellar metallicity

We assume constant stellar metallicities at fixed values of .

Nebular emission

Due to the redshift range being probed and the available constraints on the rest-frame NIR SEDs of our target sample, nebular emission lines are not expected to have a significant effect on the inferred stellar masses (unlike at z > 3 where they can have a significant effect, e.g. Stark et al. 2013; Schenker et al. 2013). Nevertheless, a simple prescription for nebular emission is included in the model SEDs allowing for escape fractions of . Details of the assumed emission line ratios for Balmer and metal lines, as well as the nebular continuum prescription can be found in Duncan et al. (2014). As in previous studies we make the simplifying assumption that the gas-phase stellar metallicity equals the stellar metallicity.

Dust attenuation

We incorporate dust attenuation following the two-component dust model of Charlot & Fall (2000),which includes contributions from birth clouds (for stellar populations younger than 10⁷ yr) and an additional screen component for all ages. The wavelength dependence on the extinction curve and fraction of attenuation from the ISM are fixed at A(λ) ∝ λ^−0.7 and 0.3, respectively – as originally presented by Charlot & Fall (2000, to which we refer the reader for formal definitions).

After convolution of the model spectral energy distributions with the photometric filters of each field at each redshift in the grid, the model grid is fit to the observational datapoints for all optical catalogue sources with a photo-z z_1,median < 1.5 – the regime for which photo-z can be considered reliable for the galaxy population. Fitting to the observed photometry is done using a simple least-squares fit to all available photometric bands for a given source (see Sect. 2.5 of Duncan et al. 2019b, for details). However, rather than using simply the normalisation (and hence stellar mass) of only the best-fit individual model, we marginalise over the full set of stellar population parameters (star-formation history, dust attenuation, metallicity) to derive a likelihood weighted distribution for the inferred stellar mass. This marginalisation implicitly assumes a flat prior on the respective stellar population parameters. For the resulting catalog, we then take the median and 1σ (16 and 84th percentiles) of the stellar mass distribution as the estimate and corresponding uncertainty for a given source. We note however that this uncertainty does not reflect the full statistical uncertainty on the stellar masses as we do not account for the uncertainty in redshift.

In addition to the stellar mass estimate for each source, we also derive the rest-frame magnitude for all photometric filters in a given field (regardless of whether that filter was used during fitting) based on the best-fitting template. Future papers that incorporate the full panchromatic SED (including deblended 24–500 μm observations) will provide robust estimates of the star-formation rates and other key physical properties within the LoTSS Deep Field source population.

5.2 Flux zero-point and total flux corrections

During testing and calibration of the stellar mass estimates, we found that additional flux zeropoint corrections and model uncertainties similar to those used during photo-z estimation are required for accurate and reliable mass estimates (see also Tanaka 2015). Offsets to the flux zero-points for stellar mass estimates are calculated following a similar approach to that applied during photo-z template fitting. The subset of optical sources with robust spectroscopic redshifts are first fit with no flux offsets applied and with no additional errors included to account for model uncertainty. Corrections to each observed filter are then calculated based on the median ratio between the observed best-fit model fluxes in that band. The sample is then re-fit incorporating the initial zero-point corrections to derive iterative changes to the corrections, with the process repeated a final time. After the third iteration of zero-point corrections, we then estimate the residual model uncertainties for the SED fitting model grid following the same method as in Sect. 3.1.2. The resulting corrections are included in Appendix B for reference.

When the final derived zero-point corrections are applied during fitting, we observe a significant improvement in the distribution of best-fit χ² values for the full photometric samples. For example in EN1 we find a median reduction per source of Δχ² = −0.68, and a reduction in the average χ² for the whole sample by ≈2. Furthermore, we also find an increased agreement between the observed stellar mass functions (SMFs) between fields and with those available in the literature. We note that when fitting stellar masses using the zeropoint corrections derived from the photo-z estimates (specifically the combined template sets of Brown et al. 2014, 2019), we find no systematic offset in the resulting stellar masses – with a median offset of ≲0.03 dex, significantly smaller than the typical uncertainty on individual estimates.

However, due to the use of fixed apertures (and aperture corrections) within the optical photometry catalogues, we find that we underestimate the mass of the most massive galaxy population at lower redshifts (which have sizes significantly larger than those of the faint galaxies or stars used for the aperture corrections). To estimate corrections from aperture to total fluxes for brighter resolved sources we cross-match the optical photometry in the deep fields to the model-fitting photometry of the Legacy Surveys (Dey et al. 2019). For sources detected in Legacy Surveys imaging, we calculate the ratio between the LOFAR optical catalogue (Paper III) aperture corrected flux and the model-fitting photometry based total flux in the available bands (g, r, z for EN1 andLH and z for Boötes). After accounting for any small global offset between the two flux measurements based on the median flux ratio for the source types used for aperture corrections, the resulting flux ratio can be used to estimate an approximate aperture to total flux conversion for each source. In the faint limit where a source has a S/N less than 5σ or is undetected in Legacy Surveys imaging, we assume that sources in the LOFAR optical catalogue are unresolved and the aperture corrected fluxes provide a reliable estimate of the total flux. In the stellar mass catalogue presented in this data release we provide the stellar mass and rest-frame magnitude values as measured, alongside total flux corrections where available.

5.3 Stellar mass completeness

When performing studies of different galaxy or AGN populations across redshift and environments in flux limited surveys, understanding the associated completeness limits is essential for minimising biases and defining appropriate volume limited samples. We are therefore interested in the mass completeness limits of our surveys in order to enable reliable studies in mass-selected samples.

We estimate the detection completeness limit in the K (or K_s) and IRAC 4.5 μm bands empirically by fitting a power law distribution to the observed number counts for > 5σ sources in a regime where the catalogues are known to be complete. The 90 and 50% completeness limits are then defined as the magnitude at which the observed number counts are 0.9 and 0.5× the expected number counts predicted by the power law distributions. Although the true magnitude completeness limits will vary as a function of intrinsic source size and morphology, completeness estimates derived from simple power-law fits have been shown to be in good agreement with completeness estimates derived from detailed simulations (through injection and detection of fake sources within the imaging, see e.g. Guo et al. 2013). We estimate the 90% magnitude completeness limits to be K or K_s = 20.44, 21.78, and 21.87 for Boötes, EN1, and LH, respectively. While for the IRAC 4.5 μm band we estimate 90% magnitude completeness limits of 21.33, 21.19, and 21.18 respectively, based on the SDWFS observation in Boötes and the SWIRE observations in the other two fields.

To estimate the mass completeness limit associated with the measured detection completeness limit, we derive the limiting stellar mass-to-light ratio as a function of redshift in two ways. Firstly, following the empirical method presented by Pozzetti et al. (2010), in small redshift bins we take the 20% faintest sources above the 90% magnitude completeness limit. The measured stellar masses for this sample are scaled to the magnitude limit () and the mass completeness limit derived from 95th percentile of the scaled mass distribution. This measurement is repeated for both the observed stellar masses and K (K_s) and IRAC 4.5 μm magnitude limits in each field.

Our second derivation of the mass completeness limit is based on a maximally old stellar population representative of the highest mass-to-light ratio expected within the sample. Using the same stellar populations models as used for SED fitting, we model a dust-free exponentially declining star-formation history with an e-folding time of 50 Myr and a formation redshift of z = 20. For each redshift bin, the stellar population is convolved with the K (K_s) and IRAC 4.5 μm filter response curves with the stellar population age set by the time since the formation redshift. The mass completeness limit corresponding to the maximally old population is then derived by scaling the observed template flux at each redshift bin to the corresponding magnitude completeness limit.

In Fig. 11 we find that the two approaches for estimating the limiting mass-to-light ratio at each redshift yield almost identical results. Selecting based on either near- or mid-IR, all three fields are complete to at 1.4 < z ≤ 1.5 - sufficient to probe below the knee of the galaxy SMF at this redshift. At z < 1, the deeper NIR observations provided by UKIDSS Deep Extragalactic Survey (DXS) mean that the EN1 and LH fields are complete to significantly lower masses if using K band to select samples. We note that by restricting observations to the Spitzer SERVS region of the EN1 and LH fields (≈ 1 magnitude deeper than SWIRE), mass complete samples at z ~ 1.5 could likely reach ~0.4 dex lower in mass at the expense of total volume probed.

5.4 Stellar mass validation

To validate our stellar mass estimates and verify the consistency both between the three LOFAR deep fields and with the literature, we estimate the galaxy SMF in each field. After masking regions affected by bright stars and restricting to regions with both optical and NIR or mid-IR coverage, the resulting survey areas in each field are 8.63, 6.5 and 7.24 deg² for Boötes, EN1, and LH, respectively (we note that these areas correspond to the region where our selection criteria have been satisfied, rather than the full areas quoted in Paper III). For the purposes of this comparison, we construct our estimate of the galaxy SMF through simple calculation of the volume density in bins of stellar mass. We do not include additional corrections to weight for the reduced volume for sources near the detection limit (i.e. the classical 1∕V_max method) or for the effects of Eddington bias on the high mass bins.

For each redshift range, we limit to the 90% mass completeness limit based on the deepest available band over a large area, specifically the UKIDSS DXS K for the EN1 and LH fields, and the IRAC 4.5 μm band for Boötes. We exclude known IR, X-ray and optical AGN that are likely to result in poor SED fits or biased stellar mass estimates. Additionally, we apply a cut on the best-fit χ² to exclude sources with χ²∕N (where N equals the number of filters used during fitting) to exclude the 5% worst fits indicative of unidentified stars or AGN within the sample.

Our derived SMF estimates for four redshift bins (0.15 < z ≤ 0.3, 0.3 < z ≤ 0.5, 0.5 < z ≤ 1.0, and 1.0 < z ≤ 1.5) are presented in Fig. 12 alongside a selection of published SMFs from the literature that probe comparable redshifts (Pozzetti et al. 2007; Muzzin et al. 2013; Leauthaud et al. 2016; Wright et al. 2018; Leja et al. 2020, McLeod et al., in prep.). We plot error bars for our binned SMF estimates that represent the uncertainty only from Poisson noise and the approximate cosmic variance based on the volume probed for each field. Our cosmic variance uncertainties are calculated following the prescription and code presented in Moster et al. (2011), with the standard deviation on the number counts at a given mass interpolated from the mass ranges produced by the code.

As noted above, we find that total flux corrections are necessary to produce a broad agreement between the three fields and with observed SMFs and those of the literature. However, due to the different assumptions used for aperture corrections between Boötes (analytic corrections based on the PSF), EN1 and LH (empirical corrections), we include total flux corrections at z < 0.5 for the EN1 and LH fields, while for Boötes we include corrections out to z < 1. We find that corrections are not required at higher redshift, with the three fields providing results that are consistent with the population becoming largely unresolved.

At z > 0.3, where the fields probe a large representative volume, our simple galaxy SMFs are in excellent agreement both with published literature values and with self-consistent results across all three fields. In the lowest redshift bin, 0.15< z ≤ 0.3, we observe an increased scatter between the SMF estimates of the three different fields, with EN1 having a consistently higher normalisation than the other two fields. Given the significantly smaller volume probed in this bin (~ 3.4× less than at 0.3 < z ≤ 0.5), larger cosmic variance between fields is to be expected. While the estimated cosmic variance uncertainties do not fully account for the offset between fields in all redshift bins, the overall normalisation and shape of our observed SMFs are still broadly consistent with the range published in the literature.

Finally, as with the photo-zs, we can compare our stellar mass estimates on a source by source basis with those available in the literature. Specifically, we compare our EN1 masses with those presented in the HSC PDR2 (Nishizawa et al. 2020, see Sect. 4.1.1), which assume the same Chabrier (2003) IMF as used in this work. For sources with in our catalog, we find a median offset of 0.09 dex between ourstellar masses and those of the HSC-DeMP (Hsieh & Yee 2014) methodology, with a corresponding robust scatter of 0.28 dex. For the HSC-MIZUKI template based estimates (Tanaka 2015), we find a median offset of just 0.01 dex and robust scatter of 0.32 dex. We note that the typical difference between our estimates and either of the HSC estimates is smaller than the difference between the two HSC estimates themselves when calculated for the same galaxy sample (a median offset of 0.14 dex).

Given the good agreement and consistency between our three fields at the redshifts of interest, we are therefore confident, that with appropriate sample selection and quality cuts, our stellar mass estimates are suitable for robust quantitative studies across all three fields. Furthermore, the large samples and volumes probed by the combined deep fields dataset over a wide range of cosmic history offer reliable reference samples for studying radio properties as a function of history, mass, and environment.

Fig. 11 Observed stellar mass distribution as a function of redshift for the three LoTSS Deep Fields. The background density plot shows the mass distribution of sources brighter than the 90% magnitude limit of the most sensitive reference band available in that field (4.5 μm for Boötes and K for EN1 and LH). Solid lines represent the 90% mass completeness limits, , derived empirically (Pozzetti et al. 2010) and dashed lines show the corresponding limits derived from model stellar populations. The blue completeness curves are derived based on the K or Ks magnitude limit, while red curves are derived from the widest area 4.5 μm limit (i.e. SWIRE for EN1 and LH).

Fig. 12 Estimated galaxy SMFs in each LoTSS deep field based on the stellar mass estimates presented in this data release (black data points), with a sample of results from the literature over the same redshift ranges. For the SMF estimates from this work, we only plot only the mass ranges that are estimated to be 90% complete (see Sect. 5.3). Total flux corrections to the inferred stellar masses are applied for all fields at z < 0.5 and for Boötes at z < 1.

Fig. 13 Distribution of LoTSS Deep Field sources as a function of redshift, zBest (Z_BEST in our catalogues), and radio power L150MHz assuming a constant spectral index of α = −0.7. The colour scale illustrates the combined number of sources across the three deep fields with a logarithmic scale. The 5σ luminosity limit for the wide area LoTSS DR1 (Shimwell et al. 2019) data is plotted in red for comparison.

6 Radio source properties

The combination of exquisite radio continuum observations and extensive multi-wavelength data covering > 20 deg² provided by the LoTSS Deep Fields data release offers the potential for a wide range of studies of the radio population – ranging from studies of the cosmic star-formation history through to obscuration free studies of the black hole accretion history (e.g. Smolčić et al. 2017a; Novak et al. 2017). In addition to the SED fitting presented in this paper, extensive detailed physical modelling of the radio population that incorporates additional FIR information (e.g. Calistro Rivera et al. 2016) can offer further insight into the physical properties of the faint radio source population. Building upon the multi-wavelength datasets presented in Paper III and the photo-zs presented here, Best et al. (in prep.; Paper V) combine results from multiple SED fitting tools to derive consensus source classifications for the radio population. Further papers will utilise these classifications and the associated detailed SED modelling to study various open questions regarding AGN and galaxy evolution.

We therefore do not present a detailed analysis of the LoTSS Deep Field radio source population in this study. However, to demonstrate the scientific potential offered by the LoTSS Deep Fields release, we present some basic properties of the radio source population derived from the catalogues presented in this work. In Fig. 13, we show the overall distribution of the LoTSS Deep Field sources as a function of both redshift (based on our best available estimate, z_Best, where spectroscopic redshifts are used when available and photo-zs otherwise) and of the inferred low-frequency radio luminosity, L_150MHz.

The colour scale in Fig. 13 illustrates the total number of sources per cell across all three fields, with no additional correction for radio source completeness. However, we do apply a cut on the measured photo-z, restricting the sample to z_Best < 4 and also exclude sources affected by bright stars. Clearly illustrated by this figure is the broad dynamic range offered by the LoTSS Deep Fields, which offer large samples of low power radio sources (log₁₀(L_150MHz∕W Hz⁻¹) <25) while also probing a sufficient volume to detect statistical samples of the high power radio source population.

In Fig. 14, we present a different view of the radio source population, showing the mean log stellar mass in cells of redshift and radio luminosity. We limit the analysis to z_Best < 1.5, where stellar mass estimates are available, and exclude sources with identified AGN components likely to bias the stellar mass estimates. At low redshift and low radio power, the average stellar mass is significantly below M^⋆ at these redshift ranges (M^⋆ ≈ 10¹¹ M_⊙, see Fig. 12); this is consistent with the expectation that the radio continuum population at these luminosities is dominated by star-forming galaxies (Best et al. 2005; Smolčić et al. 2017b). As radio luminosity increases, the radio population is hosted by increasingly massive galaxies. This apparent evolution is likely driven by two trends: the first being the transition from the regime where star-forming galaxies dominate the radio source population to where AGN dominate, and the second being the strong correlation between stellar mass and radio AGN activity for radio AGN (Sabater et al. 2019).

We note the important caveat that this analysis excludes the radio AGN that are hosted by sources with significant dust torus emission (aswould be detected by IRAC colour criteria) as well as radio-quiet quasars. The trends observed are therefore not representative of the full evolution in radio AGN host properties over this redshift range. Nevertheless, Fig. 14 illustrates the potential of the catalogues presented in this work, which when combined with additional source classification and complex SED fitting can offer a detailed picture of the evolution of the radio source population.

Fig. 14 Average stellar mass (taken as the mean in log10 space) of the LoTSS Deep radio population as a function of redshift, zBest (Z_BEST in our catalogues), and radio power L150MHz – assuming aconstant spectral index of α = −0.7. The 5σ luminosity limit for the wide area LoTSS DR1 (Shimwell et al. 2019) data is plotted in red for comparison (with the solid line representing the redshift range with high quality photo-z for the comparable galaxy population D19). The typical M⋆ for the galaxy SMF at these redshifts is marked on the colour scale for reference (white dashed line).

7 Summary

In this paper we present details of photo-z and stellar mass estimates produced for the LoTSS Deep Fields DR1. Photo-zs are estimated for all optical sources within the three Deep Fields (Boötes, ELAIS-N1 and LH), totalling over 5 million estimates across a combined ~ 25 deg² after appropriate optical quality cuts. Building on previous work our photo-z method combines multiple template fitting and empirical training based estimates to produce a consensus redshift prediction with well-calibrated photo-z uncertainties.

Based on the available spectroscopic training and test sample in each field, the resulting consensus photo-zs have robust scatter ranging from σ_NMAD = 0.016 to 0.02 for galaxies and/or host dominated AGN sources, and from σ_NMAD = 0.064 to 0.07 for identified AGN sources. Our estimated OLF for the corresponding subsets range from 1.5 to 1.8% and 18 to 22%, respectively. Similar to previous studies we find that the photo-z quality is a function of both optical magnitude and spectroscopic redshift.

Exploring the photo-z quality of the LoTSS Deep Field radio source populations we find that there is no strong trend in photo-z quality as a function of radio luminosity (for a fixed redshift), reproducing trends observed in previous studies (Duncan et al. 2019a). However, there is clear deterioration in photo-z quality as a function of redshift for a given radio luminosity that we attribute to selection effects in the spectroscopic sample and/or intrinsic evolution within the radio population.

In the redshift range for which we find our photo-z estimates to be reliable for host-dominated SEDs (z < 1.5), we exploit the extensive wavelength coverage to produce consistent stellar mass estimates across the three fields, as well as estimates of the corresponding stellar mass completeness. From simple estimations of the galaxy SMFs within each field and comparison with the literature, we validate our stellar masses provide reliable and self consistent estimates suitable for statistical studies across all three fields.

The catalogue presented in this work builds upon both the LOFAR radio data presented in Papers I and II, and the optical catalogues and radio-optical cross identification presented in Paper III. All data produced in this work are made available for public release to enable the full exploitation of the LoTSS Deep Field survey for a wide range of scientific goals.

Acknowledgements

K.J.D. and H.J.A.R. acknowledge support from the ERC Advanced Investigator programme NewClusters 321271. R.K. acknowledges support from the Science and Technology Facilities Council (STFC) through an STFC studentship via grant ST/R504737/1. P.N.B. and J.S. are grateful for support from the UK STFC via grant ST/R000972/1 M.B. acknowledges support from INAF under PRIN SKA/CTA FORECaST and from the Ministero degli Affari Esteri della Cooperazione Internazionale – Direzione Generale per la Promozione del Sistema Paese Progetto di Grande Rilevanza ZA18GR02. R.A.A.B. acknowledges support from the Glasstone Foundation. M.J.H. acknowledges support from the UK STFC via grant ST/R000905/1. M.J.J. acknowledges support from the UK Science and Technology Facilities Council [ST/N000919/1] and the Oxford Hintze Centre for Astrophysical Surveys which is funded through generous support from the Hintze Family Charitable Foundation. M.K.B. and A.W. acknowledge support from the National Science Centre, Poland under grant no. 2017/26/E/ST9/00216. K.M. has been supported by the National Science Centre (UMO-2018/30/E/ST9/00082). I.P. acknowledges support from INAF under the SKA/CTA PRIN “FORECaST” and the PRIN MAIN STREAM “SAuROS” projects. This paper is based (in part) on data obtained with the International LOFAR Telescope (ILT) under project codes LC0 015, LC2 024, LC2 038, LC3 008, LC4 008, LC4 034 and LT10 01. LOFAR (van Haarlem et al. 2013) is the Low Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, which are owned by various parties (each with their own funding sources), and which are collectively operated by the ILT foundation under a joint scientific policy. The ILT resources have benefitted from the following recent major funding sources: CNRS-INSU, Observatoire de Paris and Université d’Orléans, France; BMBF, MIWF-NRW, 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 Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation. This work is based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l’Univers of the Centre National de la Recherche Scientifique (CNRS) of France and the University of Hawaii. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. RCSLenS data processing was made possible thanks to significant computing support from the NSERC Research Tools and Instruments grant program. This work is based in part on observations made with the Spitzer Space Telescope, which was operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA.

Appendix A Catalogue description

The contents of the catalogue added by this work are as follows:

• Z_BEST -Best available redshift estimate

• Z_BEST_SOURCE - The source of the best available redshift, Z_BEST where 1 corresponds to spectroscopic redshift and 0 corresponds to the photo-z presented in this work

• Z_SPEC -Literature Spectroscopic Redshift

• Z_SOURCE - Source of the spectroscopic redshift

• Z_QUAL - Spectroscopic redshift quality, where flag Q = 3 means probable, Q ≥ 4 means reliable. Lower reliability redshifts have not been included.

• AGN_ZSPEC - Spectroscopic AGN/QSO flag where provided

• z1_median - Median of the primary redshift peak above 80% HPD CI

• z1_min -Lower bound of the primary 80% HPD CI peak

• z1_max -Upper bound of the primary 80% HPD CI peak

• z1_area - Integrated area of the primary 80% HPD CI peak

• z2_median - Median of the secondary redshift peak (if present) above 80% HPD CI

• z2_min - Lower bound of the secondary 80% HPD CI peak

• z2_max - Upper bound of the secondary 80% HPD CI peak

• z2_area - Integrated area of the secondary 80% HPD CI peak

• nfilt_eazy - Number of filters included in EAZY template fit|

• nfilt_atlas - Number of filters included in Atlas+AGN template fit

• nfilt_ananna - Number of filters included in Ananna et al. (2017) template fit

• chi_r_best - χ² / nfilt for best-fit galaxy or AGN template (any library)

• chi_r_stellar - χ² / nfilt for best-fit stellar template

• stellar_type - Stellar type of best-fit stellar template

Also included for all sources are the multi-wavelength AGN classifications used during photo-z estimation.

• AGN - Sources flagged by any one of optAGN, IRAGN or XrayAGN

• optAGN - Flag indicating whether source is included in Million Quasar catalogue compilation (Flesch 2015), where 1 means a source is included. Sources flagged as AGN based on their spectroscopic redshifts are also flagged.

• IRAGN - Source satisfies Donley et al. (2012) IR AGN selection criteria.

• XrayAGN- Source has X-ray counterpart.

For the Boötes field, where sources have been matched to the X-Böotes Chandra survey of NDWFS (Kenter et al. 2005), we provide the additional associated values:

• XrayFlux_0.5-2 [10⁻¹⁴ erg cm⁻² s⁻¹] XBoötes Soft X-ray Flux

• XrayHardness - XBoötes X-ray Hardness Ratio

For sources with either Z_SPEC or z1_median<1.5, we provide our estimate of the galaxy stellar mass as well as the rest-frame magnitudes for the individual best-fitting SED model. The catalogue columns associated with the additonal SED fitting are:

• zmodel -Model grid redshift used in stellar mass fit

• chi_best - χ² for the single best-fit model SED (for Z_BEST)

• Mass_median - 50th percentile of the marginalised stellar mass posterior (for Z_BEST), in units of log₁₀(M∕M_⊙)

• Mass_l68 - 16th percentile of the marginalised stellar mass posterior (for Z_BEST), in units of log₁₀(M∕M_⊙)

• Mass_u68 - 84th percentile of the marginalised stellar mass posterior (for Z_BEST), in units of log₁₀(M∕M_⊙)

• Nfilts -Number of photometric bands included in stellar mass fit.

• ap_to_model_[x] - Estimated to aperture to total flux correction derived from the model fitting photometry of the Legacy Survey in a band, x⁷.

• ap_to_model_err_[x] - Statistical uncertainty on ap_to_model_[x] derived from the combined flux uncertainties of the LoTSS Deep Field optical and Legacy Surveys flux measurement.

Finally, for each photometric band used during stellar mass estimation we provide the rest-frame magnitude in that filter for the best-fit SED template:

• [x]_rest Rest-frame magnitude in a given band, x, for best-fit SED (for Z_BEST)

Appendix B Zeropoint flux corrections

The inclusion of the zero-point offsets during template fitting has been demonstrated to lead to substantial improvement in photo-z quality. As outlined in Sect. 3.1, zero-point offsets for the photo-z estimates for each template library and dataset were derived following the method outlined in D18a. Similarly, zero-point corrections for stellar mass estimates are calculated following the method described in Sect. 5.2. Here we provide the multiplicative correction factors for all template sets in each of the LoTSS Deep Fields.

References

All Tables

All Figures

Fig. 1 Illustration of the wavelength coverage and depths for the photometric datasets used for photo-z estimation in the three LoTSS Deep Fields. The violin plot for each filter shows the distribution of magnitudes for 5σ sources. The widths of the violin plots are indicative of the area available at a given depth within each field and are scaled by the fraction of sources with 5σ detections compared to the total size of the catalog. For example, the deeper Spitzer IRAC photometry of the smaller SERVS regions (at 3.6 and 4.5 μm) are evident in the EN1 and LH plots.

In the text

	Fig. 2 Redshift distributions for the spectroscopic redshift training and test samples available in each of the fields (see Table 1 for total numbers). While Boötes has a much greater number of sources with zspec available, these are largely limited to z < 1 for galaxies but extend over a greater redshift range for the known AGN.
In the text

Fig. 3 Normalised optical magnitude distributions for the galaxy (blue histogram) and AGN (red histrogram) spectroscopic redshift samples in comparison to the full parent photometric sample. For EN1 and LH we plot the r-band magnitude distribution, while for Boötes we show the nearest equivalent wavelength, R. The offset in the magnitude distributions highlights the importance of high quality template estimates in the regime where spectroscopic training samples are not available.

In the text

	Fig. 4 Rest-frame residuals for the template fits using each of the three SED libraries employed in this work. Background grey points correspond to individual data points (i.e. one for each fitted filter per source) while the solid coloured lines correspond to the median residual within a given wavelength bin and the shaded region corresponds to the 1σ distribution (16 and 84th percentiles).
In the text

	Fig. 5 Rest-frame model uncertainties derived for each of the template libraries used during our analysis (solid coloured lines). The model uncertainties are calculated by subtracting in quadrature the scaled average fractional error (black dotted line) from the median absolute rest-frame residuals (coloured symbols). However, the contribution of the average error is negligible.
In the text

Fig. 6 Stacked redshift posterior estimates for the galaxy (or host-dominated) population (top) and for the combined AGN selected population (bottom; IR, X-ray or optically selected) for each of the three fields. The spectroscopic redshift samples available become particularly sparse at z > 1 for the EN1 and LH fields. Grey dotted and dashed line correspond to ± 0.1 × (1 + zspec) and ± 0.2 × (1 + zspec), respectively. The colour scale is illustrative and represents the summed P(z) (normalised such that ∫P(z)dz = 1) for each zspec bin, more precise photo-z therefore result in darker peaks.

In the text

	Fig. 7 Robust scatter (σNMAD; upper panels) and OLF (lower panels) for the consensus photo-z estimates (z1,median) as a function of spectroscopic redshift. Shaded regions illustrate the statistical uncertainties on the respective metrics derived from bootstrap resampling.
In the text

Fig. 8 Robust scatter (σNMAD; upper panels) and OLF (lower panels) for the consensus photo-z estimates (z1,median) as a functionof 4.5 μm magnitude. Shaded regions illustrate the statistical uncertainties on the respective metrics derived from bootstrap resampling. The magnitude range plotted for each field is representative of the magnitude distribution of the spectroscopic sample and is not indicative of the relative photometric depths (see Fig. 1).

In the text

Fig. 9 Robust scatter (σNMAD; top) and OLF (bottom) for the consensus photo-z estimate as a function of spectroscopic redshift and 150 MHz radio continuum luminosity in the Boötes field. Top and side panels: trends averaged over all redshifts and luminosities, respectively. For a cell to be plotted we require a minimum offive galaxies - for some redshifts (or luminosities) we only have the number statistics available to plot the statistics averaged over all luminosities (or redshift).

In the text

Fig. 10 Cumulative distribution of threshold credible intervals, c, () plot for the final consensus photo-z estimates for the galaxy (or host-dominated AGN) population (top) and the optical, infrared and X-ray selected AGN population (bottom) within each of the deep fields. Coloured lines represent the distributions in bins of apparent optical magnitude (I-band for Bootes, r-band for EN1and LH). Lines that fall above the 1:1 relation illustrate under-confidence in the photo-z uncertainties (uncertainties overestimated) while lines under illustrate over-confidence (uncertainties underestimated).

In the text

Fig. 11 Observed stellar mass distribution as a function of redshift for the three LoTSS Deep Fields. The background density plot shows the mass distribution of sources brighter than the 90% magnitude limit of the most sensitive reference band available in that field (4.5 μm for Boötes and K for EN1 and LH). Solid lines represent the 90% mass completeness limits, , derived empirically (Pozzetti et al. 2010) and dashed lines show the corresponding limits derived from model stellar populations. The blue completeness curves are derived based on the K or Ks magnitude limit, while red curves are derived from the widest area 4.5 μm limit (i.e. SWIRE for EN1 and LH).

In the text

Fig. 12 Estimated galaxy SMFs in each LoTSS deep field based on the stellar mass estimates presented in this data release (black data points), with a sample of results from the literature over the same redshift ranges. For the SMF estimates from this work, we only plot only the mass ranges that are estimated to be 90% complete (see Sect. 5.3). Total flux corrections to the inferred stellar masses are applied for all fields at z < 0.5 and for Boötes at z < 1.

In the text

Fig. 13 Distribution of LoTSS Deep Field sources as a function of redshift, zBest (Z_BEST in our catalogues), and radio power L150MHz assuming a constant spectral index of α = −0.7. The colour scale illustrates the combined number of sources across the three deep fields with a logarithmic scale. The 5σ luminosity limit for the wide area LoTSS DR1 (Shimwell et al. 2019) data is plotted in red for comparison.

In the text

Fig. 14 Average stellar mass (taken as the mean in log10 space) of the LoTSS Deep radio population as a function of redshift, zBest (Z_BEST in our catalogues), and radio power L150MHz – assuming aconstant spectral index of α = −0.7. The 5σ luminosity limit for the wide area LoTSS DR1 (Shimwell et al. 2019) data is plotted in red for comparison (with the solid line representing the redshift range with high quality photo-z for the comparable galaxy population D19). The typical M⋆ for the galaxy SMF at these redshifts is marked on the colour scale for reference (white dashed line).

In the text