Cosmic evolution of low-excitation radio galaxies in the LOFAR two-metre sky survey deep fields Open Access

ABSTRACT

1 INTRODUCTION

It has been over two decades since the discovery of the tight correlation between the mass of the central supermassive black hole (SMBH) M_•, and the velocity dispersion (σ), of the bulge component of a galaxy (the M_• − σ relation; Ferrarese & Merritt 2000; Merritt & Ferrarese 2001), suggesting a co-evolution between the SMBH and its host galaxy. Accretion of matter on to an SMBH can power an active galactic nucleus (AGN), producing vast amounts of energy, which, when efficiently coupled to the surrounding gas, can disrupt star formation in the host galaxy or offset cluster cooling; this process is known as AGN feedback. It is widely accepted that almost all massive galaxies harbour an SMBH at their centres (e.g. Kormendy & Ho 2013), with observations showing that feedback from AGN plays a key role in shaping the observed galaxy populations in the local Universe (see reviews by McNamara & Nulsen 2007; Cattaneo et al. 2009; Fabian 2012; Heckman & Best 2014). This wealth of observational work is complemented by galaxy formation models (see review by Somerville & Davé 2015) and cosmological simulations (e.g. Bower et al. 2006; Croton et al. 2006; Kaviraj et al. 2017; Davé et al. 2019), which often require feedback from AGN to successfully reproduce the observed local galaxy properties, such as the high-mass end of the galaxy luminosity function (e.g. Bower et al. 2006) and the observed bi-modality in galaxy colours (e.g. Cattaneo et al. 2006), and to prevent runaway cooling flows in galaxy clusters (e.g. Peterson & Fabian 2006; Cattaneo et al. 2009).

Detailed characterization of the AGN in the local Universe has revealed that AGN can be split into two classes, which show differences in their host galaxy properties: ‘radiative-mode’ and ‘jet-mode’ AGN (e.g. Allen et al. 2006; Hardcastle, Evans & Croston 2007; Best & Heckman 2012; Janssen et al. 2012; Heckman & Best 2014; Mingo et al. 2014; Tadhunter 2016; Ching et al. 2017b; Hardcastle & Croston 2020). Radiative-mode AGN accrete matter on to the black hole in a radiatively efficient manner; such an accretion flow results in the formation of an optically thick, geometrically thin, accretion disc, and a dusty obscuring structure (e.g. Shakura & Sunyaev 1973); these AGN can sometimes launch jets. The photoionization of gas by photons from the accretion disc results in high-excitation emission lines in their optical spectra (e.g. Best & Heckman 2012). In contrast, jet-mode AGN are fuelled by radiatively inefficient accretion on to the black hole (Narayan & Yi 1994, 1995) and do not display signs of an optically thick accretion disc or torus structure or AGN activity at other wavelengths but are efficient at producing bipolar relativistic jets.

Radio observations provide a powerful technique for identifying AGN based on the synchrotron emission from the jets that result from the acceleration of relativistic charged particles by the central engine; these sources are appropriately called ‘radio-loud’ AGN. The radio regime offers the only way to identify and study the jet-mode AGN class, and also allows the selection of radiative-mode AGN in a manner that is unaffected by dust extinction. Due to the nature of the emission line properties of these sources, historically, the radio-detected radiative-mode AGN are known as high-excitation radio galaxies (HERGs), while the jet-mode AGN are referred to as low-excitation radio galaxies (LERGs).

Detailed characterization of the host galaxy properties of these two AGN classes has revealed that LERGs are typically hosted in massive red galaxies with old stellar populations and very massive black holes (e.g. Tasse et al. 2008; Smolčić 2009; Best & Heckman 2012; Janssen et al. 2012; Mingo et al. 2014; Ching et al. 2017a; Whittam et al. 2018; Williams et al. 2018), whereas the HERGs tend to be hosted in less massive, bluer galaxies with recent star-formation activity (Heckman & Best 2014). Compared to HERGs, LERGs are also found to have typically lower radio luminosities (e.g. Best & Heckman 2012; Best et al. 2014; Pracy et al. 2016), and are more often found in rich group or cluster environments (e.g. Best 2004; Tasse et al. 2008; Gendre et al. 2013; Sabater, Best & Argudo-Fernández 2013; Croston et al. 2019). These results are consistent with the idea that the radio-loud AGN activity in these two classes of AGN is triggered by different mechanisms related to the available fuel source; in this picture, LERGs are thought to be fuelled by the cooling hot gas from haloes present in their massive host galaxies, whereas the HERGs tend to accrete efficiently from cold gas, from either external (e.g. mergers or interactions) or internal (e.g. bars or instabilities) processes, with the more plentiful gas supply leading to the formation of an accretion disc (e.g. Allen et al. 2006; Hardcastle et al. 2007; Heckman & Best 2014).

However, in reality, the picture is not quite so simple; in recent years, improved models for the nature of the accretion flow from the hot phase suggest that under certain conditions, thermal instabilities in the hot gas medium can result in the formation of cold filaments of gas that then ‘rain’ down on to the black hole (e.g. Sharma et al. 2012; Gaspari, Ruszkowski & Oh 2013; Gaspari, Brighenti & Temi 2015; Gaspari, Temi & Brighenti 2017). This means that accretion on to the SMBHs in LERGs may also be from gas in a ‘cold’ state. Indeed there is now growing observational support for this scenario, with observations of a handful of LERGs being fuelled by cold gas (e.g. Tremblay et al. 2016; Ruffa et al. 2019b).

It has been suggested that the Eddington-scaled accretion rate on to the black hole, rather than the gas origin alone, may determine the formation of a radiatively efficient/inefficient AGN (e.g. Best & Heckman 2012; Best et al. 2014; Hardcastle 2018; Hardcastle & Croston 2020). An analogy of this model can be drawn to the different spectral states of X-ray black hole binaries, where a switch from the ‘soft’ (analogous to radiative-mode AGN) to the ‘hard’ (analogous to the jet-mode AGN) state is governed by a change in the nature of the accretion flow occurring at a few per cent of the Eddington rate (see reviews by Remillard & McClintock 2006; Yuan & Narayan 2014). In such a model, fuelling by cooling of hot gas generally occurs at low accretion rates in massive galaxies which host massive black holes, leading to a LERG, whereas low-mass galaxies with abundant cold gas supply are likely to result in higher (Eddington-scaled) accretion rates, and hence in the formation of an HERG. However, due to the stochastic nature of accretion processes, accretion from cold gas can occur at low Eddington rates leading to a LERG. Moreover, Whittam et al. (2018) found a difference in the Eddington-scaled accretion rate distribution of LERGs and HERGs out to z ∼ 0.7, but with considerable overlap in these distributions in their deeper radio data compared to that found by Best & Heckman (2012) in the local Universe. Therefore, the origin of the differences in these two populations and the precise physical mechanisms that trigger the different modes of AGN in different galaxies are not well understood.

In the local Universe, the incidence of LERGs shows a strong dependence on the stellar mass (⁠|${\propto}M_{\star }^{2.5}$|⁠) and black hole mass (⁠|${\propto}M_{\bullet }^{1.6}$|⁠; Best et al. 2005a; Janssen et al. 2012). In these massive systems, the time-averaged energetic output from the jets, in the form of mechanical energy, is found to balance the radiative cooling losses from the hot gas (e.g. Best et al. 2006; McNamara & Nulsen 2007; Hardcastle et al. 2019). Therefore, feedback from LERGs plays an important role in the evolution of massive galaxies as the medium into which the jet energy is deposited also forms the eventual fuel source for the AGN, thus providing the conditions for a self-regulating AGN feedback cycle. Understanding the cosmic evolution of the LERG population and the host galaxies in which they reside is crucial in marking their role in early galaxy evolution and testing this picture established from detailed observations in the local Universe.

Best et al. (2014) were the first to study the evolution of the HERG and LERG luminosity functions, separately, out to higher redshifts (z ∼ 1) using a combination of various radio-AGN data sets and spectroscopic information. They found that for the LERGs, the space density decreases mildly with redshift at low-luminosities (⁠|$L_{\rm {1.4\, GHz}} \lesssim 10^{24}\, \rm {W\, Hz^{-1}}$|⁠) but increases with redshift at higher radio luminosities. They developed various models to explain the observed evolution but concluded that deeper, higher redshift data were required to distinguish between the models. Similarly, Pracy et al. (2016) constructed the evolving HERG and LERG luminosity functions (at z < 0.75) using a sample of 5000 radio-AGN derived from cross-matching the Faint Images of the Radio Sky at Twenty-cm (FIRST; Becker, White & Helfand 1995) survey at 1.4 GHz and the Sloan Digital Sky Survey (SDSS; York et al. 2000). They also found that the LERG population shows little-to-no evolution over this redshift. More recently, Butler et al. (2019) investigated the evolution of a sample of 1729 LERGs using 2.1 GHz observations of the XMM extragalactic survey field (XXL-S) using the Australia Telescope Compact Array (ATCA). They found that the space densities of the LERGs showed weak evolution out to z ∼ 1.3. At low radio frequencies, Williams et al. (2018) studied the evolution of HERGs and LERGs, classified using photometry alone, between 0.5 < z ≤ 2 using 150 MHz LOw Frequency ARray (LOFAR; van Haarlem et al. 2013) observations of the Boötes field. Over this redshift range, they found that the LERG population shows a strong decline in space densities that is consistent with the decline in the space densities of the expected hosts, i.e. the massive quiescent galaxies.

In this paper, we present the first robust measurement of the LERG luminosity functions out to z ∼ 2.5 and study the cosmic evolution of their host galaxy properties using a sample of over 10 000 LERGs constructed by combining deep radio observations from LOFAR with deep, wide-area multi-wavelength photometry. The radio observations come from the first data release of the LOFAR Two-metre Sky Survey Deep Fields (LoTSS-Deep; Duncan et al. 2021; Kondapally et al. 2021; Sabater et al. 2021; Tasse et al. 2021), covering three extragalactic fields and forming the deepest radio continuum survey to date at low frequencies. Our combination of deep and wide radio and multiwavelength data sets, covering ∼25 deg², is ideal for probing much fainter luminosities out to higher redshifts than previous studies (e.g. Best et al. 2014; Pracy et al. 2016; Williams et al. 2018; Butler et al. 2019) and also allows a better sampling of the bright end of the luminosity function compared to previous deep observations over small areas (e.g. Smolčić et al. 2017a,b), while limiting the effects of cosmic variance. We focus on the evolution of the LERG population as our deep radio data set is particularly well-suited to sample the low-luminosity AGN population, which is dominated by LERGs; this allows us to characterize the evolution of this population in unprecedented detail.

The paper is structured as follows. In Section 2 we describe the radio and multiwavelength data sets used, and the selection of the sample of LERGs. Section 3 describes the construction of the luminosity functions, completeness simulations, and the evolution of the radio-AGN luminosity functions. Section 4 presents the evolution of the LERG luminosity functions and comparisons with literature. Section 5 investigates the incidence of LERGs as a function of stellar mass, split into those hosted by quiescent galaxies and star-forming galaxies. Section 6 presents the LERG luminosity functions split into those hosted by star-forming and quiescent galaxies and models the evolution of the quiescent LERG population. We draw our conclusions in Section 7.

Throughout this work, we use a flat ΛCDM cosmology with Ω_M = 0.3,  Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹, and a radio spectral index α = −0.7 (where S_ν ∝ ν^α). Where quoted, magnitudes are in the AB system (Oke & Gunn 1983), unless otherwise stated.

2 DATA

In this paper, we utilize a large sample of radio-detected sources combined with other multiwavelength data sets from LoTSS-Deep Data Release 1 in the European Large-Area ISO Northern-1 (ELAIS-N1), Lockman Hole, and Boötes fields.

2.1 Radio data

The radio data in the three fields were obtained from observations taken with the LOFAR High Band Antenna (HBA) with high spatial resolution (6 arcsec) and a central observing frequency of 146 MHz for the ELAIS-N1 field and 144 MHz for the Lockman Hole and Boötes fields.¹ Multi-epoch observations were used to build deep radio images, which cover 68 deg² in each field down to the 30 per cent power point of the primary beam. The radio images, with total integration times of 168, 112, and 80 h, reach RMS noise levels of 20, 22, and 32 |$\mu$|Jy beam⁻¹ in the central regions of the ELAIS-N1, Lockman Hole, and Boötes fields, respectively. The full details of the radio calibration and imaging are described by Tasse et al. (2021) for Lockman Hole and Boötes, and by Sabater et al. (2021) for ELAIS-N1. Source extraction was performed on the Stokes I images in each field using Python Blob Detector and Source Finder (pybdsf; Mohan & Rafferty 2015), with source catalogues extracted out to the 30 per cent power point of the primary beam.

2.2 Multiwavelength data

The three LoTSS Deep Fields contain deep, wide-area panchromatic photometry; this rich multiwavelength data set makes these fields ideal for determining robust photometric redshifts and performing spectral energy distribution (SED) fitting.

The full details of the available multiwavelength data, including the depth and areas covered by each survey, are given by Kondapally et al. (2021). In summary, ultraviolet (UV) data in all fields come from the Galaxy and Evolution Explorer (GALEX) space telescope (Morrissey et al. 2007). The optical data in the three fields comes from a combination of the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS-1; Kaiser et al. 2010) Medium Deep Survey, the Hyper-Suprime-Cam Subaru Strategic Program (HSC-SSP) survey data release 1 (Aihara et al. 2018), the Spitzer Adaptation of the Red-sequence Cluster Survey (SpARCS; Muzzin et al. 2009; Wilson et al. 2009), Red Cluster Sequence Lensing Survey (RCSLenS; Hildebrandt et al. 2016), and the NOAO Deep Wide Field Survey (NDWFS; Jannuzi & Dey 1999). In both ELAIS-N1 and Lockman Hole, near-infrared (NIR) data comes from the UK Infrared Deep Sky Survey (UKIDSS) Deep Extragalactic Survey (DXS; Lawrence et al. 2007). The mid-infrared (MIR) data in all fields comes from the Spitzer space telescope (Werner et al. 2004). In ELAIS-N1 and Lockman Hole this consists of observations from the Spitzer Wide-area Infra-Red Extragalactic (SWIRE; Lonsdale et al. 2003) survey and the Spitzer Extragalactic Representative Volume Survey (SERVS; Mauduit et al. 2012). In Boötes, the mid-infrared data come from Spitzer observations of the NDWFS field (Eisenhardt et al. 2004; Ashby et al. 2009; Decadal IRAC Boötes Survey).

In ELAIS-N1 and Lockman Hole, we make use of new forced, matched-aperture multiwavelength catalogues with UV to mid-infrared photometry generated by Kondapally et al. (2021), providing more robust and, at some wavelengths, deeper catalogues than those previously available in the literature. The full details of the catalogue creation process are provided by Kondapally et al. (2021); in summary, sources were detected using deep χ² images that incorporated information across optical-MIR wavelengths, with fluxes extracted in circular apertures and corrected to total fluxes by performing aperture corrections based on sources typical of distant galaxies. In the Boötes field, we made use of existing point spread function (PSF) matched I-band and 4.5  |$\mu$|m band selected catalogues (Brown et al. 2007, 2008) to generate a combined multiwavelength catalogue, as described by Kondapally et al. (2021).

Additional far-infrared (FIR) data come from the FIR-deblended catalogues generated by the Herschel Extragalactic Legacy Project (HELP; Shirley et al. 2021) and McCheyne et al. (2022). This data set compiles observations from Spitzer-MIPS (Rieke et al. 2004) at 24 |$\mu$|m, and imaging from the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010; with photometry at 100 and 160 |$\mu$|m), and the Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010; with bands at 250, 360, and 520 |$\mu$|m) instruments on-board the Herschel Space Observatory (Pilbratt et al. 2010). The observations from Herschel come from the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012).

2.3 Multiwavelength counterparts and value-added catalogues

The new multiwavelength catalogues generated in these fields were used to identify the counterparts of the LOFAR-detected sources in the LoTSS Deep Fields. The host-galaxy identification process was limited to the regions of each field with the best available multiwavelength coverage,² totalling ∼26 deg².

The host-galaxy identification method is described in detail by Kondapally et al. (2021). In summary, counterparts were identified using a combination of the colour-based statistical likelihood ratio (LR) method (de Ruiter, Willis & Arp 1977; Sutherland & Saunders 1992; Williams et al. 2019) and a visual classification scheme. The LR method is suitable for radio sources with a secure radio position; a decision tree was developed to select such sources. The LR method cross-matching achieved a reliability and completeness of >97.5 per cent across the three fields. Sources not selected for LR cross-matching were identified using two main visual inspection workflows with any issues with source associations, for example, blending of distinct physical sources or association of core and lobe components of an AGN, also being corrected at this step. The first visual workflow used was LOFAR Galaxy Zoo (LGZ; Williams et al. 2019), a Zooniverse-based framework, where the consensus decision from five volunteers is used to form the source associations and identifications. Sources that required additional inspection after LGZ or were otherwise more complex, for example, radio blends, were finally classified by a single expert in a workflow designed with more functionality than LGZ. The cross-matching process resulted in the identification of the host-galaxies for >97 per cent of the radio-detected sources across the three fields. Full details of the statistics and magnitude distributions of the counterparts are given in Kondapally et al. (2021).

Building on the robust multiwavelength photometric catalogues in each field, photometric redshifts for the full multiwavelength catalogues were presented by Duncan et al. (2021). These were generated using a hybrid approach, combining template fitting and machine learning methods, developed for the next generation of radio surveys (see Duncan et al. 2018a,b). In this paper, we use spectroscopic redshifts where available and reliable, which is the case for a small fraction of the sources (∼5, 5, and 22 per cent of all radio-sources in ELAIS-N1, Lockman Hole, and Boötes, respectively); otherwise the photometric redshifts are used.

2.4 Spectral energy distribution fitting

Spectral energy distribution (SED) fitting was performed using four different SED fitting codes for all LoTSS Deep Fields sources with counterparts that satisfy the quality cuts defined by Kondapally et al. (2021) and Duncan et al. (2021); this process is described in detail by Best et al. (in preparation).

The SED fitting was performed using this same input catalogue using each of agnfitter (Calistro Rivera et al. 2016), Bayesian Analysis of Galaxies for Physical Inference and Parameter Estimation (bagpipes; Carnall et al. 2018), Code Investigating Galaxy Evolution (cigale; Burgarella, Buat & Iglesias-Páramo 2005; Noll et al. 2009; Boquien et al. 2019), and Multiwavelength Analysis of Galaxy Physical Properties (magphys; da Cunha, Charlot & Elbaz 2008). One of the key differences between the SED fitting codes employed by Best et al (in preparation) is that agnfitter and cigale, unlike the other two routines, are able to model emission from AGN which imprint features in the mid-infrared regime in particular. This is done by incorporating models for the accretion disc and dusty torus surrounding the AGN in both SED fitting codes; for cigale, this includes a run based on the Fritz, Franceschini & Hatziminaoglou (2006) AGN models, and a separate run based on the skirtor (Stalevski et al. 2012, 2016) prescription which assumes a clumpy two-medium torus model (rather than a smooth torus structure; see Boquien et al. 2019; Yang et al. 2020); for agnfitter, this includes two-component models for the AGN, torus models (Silva, Maiolino & Granato 2004), and an accretion disc emission model (Richards et al. 2006). The other key difference between the different codes is that bagpipes, magphys, and cigale all enforce an energy balance between the UV-optical emission from starlight absorbed by dust and the re-radiated emission in the infrared; this ensures that the spectral distributions of galaxies are physically consistent. The version of agnfitter used by Best et al. (in preparation) does not enforce this, for example, to account for cases where the UV emission is spatially offset from the dust emission (e.g. Calistro Rivera et al. 2018; Cochrane et al. 2021); the most recent version allows users to enforce a strict energy balance. The inclusion of energy balance in the magphys, bagpipes, and cigale results has implications when comparing our results with previous agnfitter-based work from Williams et al. (2018); these are discussed in Appendix  A.

The outputs from the four different SED fitting routines were combined to generate consensus estimates of physical galaxy properties; the key parameters relevant for analysis in this study are the stellar masses and star-formation rates (SFRs). This process is described in detail by Best et al. (in preparation) but in summary, for sources that showed no signatures of a radiative-mode AGN, the stellar masses and SFRs were estimated by taking the average of the magphys and bagpipes results (accounting for their goodness of fits); both stellar masses and SFRs agree very well with each other, with a scatter |${\lt}0.15\, \rm {dex}$| (see Best et al., in preparation). For sources that showed signs of a radiative-mode AGN (see below; also referred to as ‘SED AGN’ colloquially hereafter), the stellar masses and SFRs were taken by averaging the two cigale runs with the Fritz et al. (2006) and skirtor AGN models, provided that a good fit was found; agnfitter was found to result in less reliable fits and was hence excluded from this step. We use these consensus values determined by Best et al. (in preparation) throughout this paper unless otherwise stated.

The ‘SED AGN’ (i.e. AGN showing emission from the accretion disc or torus), were identified using the outputs of the four SED fitting routines. First, Best et al. (in preparation) defined a diagnostic based on the f_AGN parameter, which corresponds to the ratio of MIR luminosity that arises from AGN components compared to that from the stellar galaxy component, as fitted by both agnfitter and cigale. In particular, Best et al. used the 16th percentile of the f_AGN parameter, to avoid bias from objects with large uncertainties on f_AGN. Secondly, Best et al. also considered the goodness-of-fit estimates from bagpipes and magphys (which do not include AGN components) compared to those from cigale and agnfitter; the latter two codes should find a better goodness of fit for sources with significant AGN contribution. Based on the combination of the |$f_{\rm {AGN, 16\text{th}}}$| parameters and the relative goodness of fit values, Best et al. identified the likely ‘SED AGN’. Finally, a small subset of the LOFAR-detected sources have either optical spectroscopy or X-ray detections, indicating the presence of a typical AGN (see Duncan et al. 2021); most of these were already identified as ‘SED AGN’ by the method above, but any additional sources were also added to this sample.

2.5 Selection of LERGs and HERGs

Radio continuum surveys detect both synchrotron emission associated with jets from AGN, and also emission from star-formation activity. The faint radio source population, especially at |$S_{\rm {150\, MHz}} \lesssim 1\, \rm {mJy}$|⁠, is expected to be dominated by star-forming galaxies (Wilman et al. 2008; Padovani 2016). Source classification of the LOFAR-detected sources was therefore performed to separate star-forming galaxies from different classes of AGN using the outputs from the SED fitting process. This process is described in detail by Best et al. (in preparation) and summarized below.

As low-frequency radio observations trace cosmic-ray electrons from supernovae from massive (recently formed) stars, there is a well-known relation between the radio luminosity and SFR for star-forming galaxies (e.g. Calistro Rivera et al. 2016; Gürkan et al. 2018; Smith et al. 2021, and references therein), with a correlation also seen between the far-infrared and radio luminosities (far-infrared radio correlation; FIRC). We selected radio-AGN (also known as ‘radio-excess AGN’) as sources that show excess radio emission (⁠|${\gt}0.7\, \rm {dex}$|⁠, ≈3σ) compared with that expected from star-formation processes alone using a ridgeline analysis by Best et al. (in preparation). To this sample, we added the small fraction of sources that are resolved and show extended radio emission that is associated with jets from the AGN even if they show a radio-excess |${\lt}0.7\, \rm {dex}$|⁠. The other sources without such radio-excess are largely star-forming galaxies (SFGs) but can also include ‘radio-quiet AGN’ which are known to broadly follow the same radio luminosity – star formation rate relation (e.g. Gürkan et al. 2018). Both of these groups are excluded from this analysis as we are focused on the radio-loud AGN population in this study. The luminosity functions of these populations are studied by Bonato et al. (2021) and Cochrane et al. (in preparation).

As discussed in the introduction, the total radio-AGN population consists of two classes of AGN: HERGs and LERGs. HERGs display strong optical emission lines indicating the presence of an accretion disc and dusty obscuring structure; these sources are consistent with being the radio-loud subset of the radiative-mode AGN population. The LERGs show powerful radio emission from the AGN but little to no evidence of having an accretion disc or obscuring structure in the IR through to the X-ray regime, and as such are not identified as AGN at other wavelengths. Using the above definitions of the two classes of AGN, in this study, we define LERGs as sources that host a ‘radio-excess AGN’ but not an ‘SED AGN’ (based on the SED fitting; see Section 2.4). Likewise, HERGs are defined as sources that are classed as both ‘radio-excess AGN’ and ‘SED AGN’.

The final sample of LERGs and HERGs used in this paper and for the construction of the luminosity functions is limited to sources that have a radio detection with peak flux density at >5σ level based on the local RMS, and excludes sources that are masked in optical bright-star masks (i.e. flag_clean =1; see Kondapally et al. 2021 for details) to ensure a clean and robust sample. In addition, we limit our analysis of the luminosity functions to 0.03 < z ≤ 2.5; at lower redshifts the volume sampled by LoTSS-Deep is small and there may be incompleteness due to the larger angular size of the nearby sources, and at higher redshifts, we would reach the limits of the multiwavelength data sets available in these fields, beyond which the source classifications of Best et al. from the SED fitting process become less secure. The above criteria result in a sample of 10 481 LERGs (of which 4563 are hosted by quiescent galaxies and 5918 are hosted by star-forming galaxies; see Section 6), and 1302 HERGs, across the three deep fields.

Fig. 1 (top) shows the 150 MHz radio luminosity as a function of redshift for the sample of LERGs from the LoTSS Deep Fields. The red line corresponds to the point source 5σ detection limit calculated based on the noise in the central region of ELAIS-N1, |$\sigma =20\, \rm {\mu }$|Jy, and a radio spectral index α = −0.7. Fig. 1 (bottom) shows the redshift distribution for the sample of LERGs (split into those hosted by quiescent and star-forming galaxies; see Section 5.1) and HERGs. We note that the ‘peaks’ seen in the redshift distribution, most prominent for the star-forming LERGs, are caused by the use of the median of the posterior photometric redshift distribution, which for very faint sources can suffer from aliasing due to gaps in the filter coverage. The full posterior distribution for such sources is smoother, and hence these ‘peaks’ seen are not physical. Neither the SED fitting nor the source classifications account for uncertainties in the photometric redshifts, however, the impact of such features on the luminosity functions can be reduced by applying an optical magnitude completeness selection (see Section 3.1). Moreover, because of the wide redshift bins used for the analysis in this paper, this small aliasing has little effect on the derived luminosity functions.

Figure 1.

Top: Radio luminosity as a function of redshift for the LERGs detected in the LoTSS Deep Fields. The colourbar with a logarithmic scale shows the number of sources across the parameter space. The red line shows the point source 5σ detection limit using |$\sigma =20\, \rm {\mu Jy}$| (achieved in the central region of ELAIS-N1) and assuming a spectral index α = −0.7. Bottom: Histogram of redshifts for the sample of LERGs split into those hosted by quiescent (red dashed line) and star-forming (blue dotted line) galaxies (see Section 5.1), along with the HERGs (orange line), selected in the LoTSS Deep fields.

Open in new tabDownload slide

3 TOTAL RADIO-AGN LUMINOSITY FUNCTIONS

3.1 Building the luminosity functions

Figure 2.

Cumulative area covered down to a given noise level in the three LoTSS Deep Fields. The area is computed after applying the optical-based bright-star masks and only including the regions with coverage from the best available multiwavelength surveys in each field (see Kondapally et al. 2021). The total area across the three fields covers ≈24.6 deg².

Open in new tabDownload slide

Figure 3.

Radio flux density completeness corrections computed using the method outlined in Section 3.2. Top: Completeness curves in the ELAIS-N1 field for different assumed source-size distributions. Curves show completeness based on applying the size-distribution of moderately bright sources (S_tot = 1−10 mJy; i.e. where completeness is high) classified as either SFGs or AGN in Best et al. (in preparation). In addition, the completeness curves resulting from the Windhorst et al. (1990) and 2 × Windhorst et al. (1990) integral angular size distribution, often used in the literature, are also shown. Bottom: Completeness curves for each of the three LoTSS Deep Fields assuming an ‘AGN’ size-distribution; this is used in the construction of the AGN LFs for each field. A table listing these adopted corrections for each field is shown in Table 1. The ‘combined’ curve is the area-weighted average of the completeness curves in the three fields, used for constructing the combined LFs across the three fields.

Open in new tabDownload slide

In constructing the LFs, we consider only the radio sources that have a host galaxy counterpart; this corresponds to 97 per cent of the radio sources in the LoTSS Deep Fields (see Kondapally et al. 2021). For the remaining 3 per cent of sources, the lack of a multiwavelength counterpart identification means that no photometric redshift was available, and hence SED fitting and subsequent source classification were not possible. Kondapally et al. (2021) performed an optical to mid-IR stacking analysis of the subset of these sources with secure radio positions (which forms the majority) and found that these sources are likely to be predominantly z > 3 radio-AGN. Therefore, we expect that this will have a negligible effect on the derived LFs in this paper as we consider only sources with z ≤ 2.5.

In the calculation above, we have assumed that the radio flux density sets the limit on the V_max for a given source. It is, however, possible that the optical/IR data set may set the limits on the maximum observable volume for a source. The proper application of optical/IR magnitude limits on the LFs is complicated by the use of χ² detection images, analysis of these would involve detailed modelling of how they would affect both the photometric redshift accuracy and the source classifications; this is beyond the scope of this paper. We have, however, considered the effect of applying an optical/IR magnitude limit in the calculation of V_max and find that the radio data set predominantly sets the limits on the V_max. The LFs constructed when also considering an optical/IR magnitude limit agree very well with those based on the radio selection effects alone. Therefore, in this paper, we consider only the radio limits when constructing the LFs.

To compute the uncertainties on the luminosity functions, we performed bootstrap sampling (random sampling by replacement) of the catalogue to generate a distribution of 1000 realizations of the luminosity function. The lower and upper 1σ uncertainties on our luminosity functions are then determined from the 16th and 84th percentiles of the bootstrap samples. For the faint luminosity bins, where the samples are large and the uncertainties computed from bootstrapping correspondingly small, the true uncertainties are likely to be dominated by other factors such as the photometric redshift errors and source classification uncertainties. Therefore, we set a minimum uncertainty of 0.03 dex in the luminosity functions reported, based on the ∼7 per cent photometric redshift outlier fraction in these fields.

3.2 Radio completeness corrections

The radio flux density completeness corrections are generated by performing simulations of inserting mock sources of various intrinsic source-size and flux density distributions into the radio image, and then recovering them using the same pybdsf parameters as used for the real sources (see Sabater et al. 2021 for the pybdsf parameters used).

We simulated mock sources at fixed total flux density intervals separated by |$0.2\, \rm {dex}$| in the range 0.4 < S < 40 mJy and, using a finer sampling interval of |$0.08\, \rm {dex}$| in the range |$90 \lt S \lt 400\, \rm {\mu }$|Jy to better probe the expected sharp decline in completeness at faint flux densities. For the Lockman Hole and Boötes fields, we used the same sampling intervals but only simulating sources down to |${\sim}110\, $| and |${\sim}190\, \rm {\mu }$|Jy, respectively, due to the slightly shallower depth of the radio data. The flux intervals used are listed in Table 1.

For each field, we simulated 120 000–150 000 mock sources to sample the full range in (total) flux density and source-size parameter space, while also obtaining robust statistics for the bright and extended rare sources. Practically, this was done by inserting 1000 mock sources with convolved sizes between 6 and 30 arcsec, where 6 arcsec is the size of the LOFAR beam, for a given (total) flux density value into the radio image and extracting the sources using pybdsf, with this step repeated many times (see Section 3.2.1 for details). The injected sources were modelled as Gaussians, although the structure of real sources may be more complex. We ensure that a mock source is placed at least twice its full width at half-maximum (FWHM; along the major axis) away from other mock sources and real radio sources to avoid source overlapping, which can complicate the process of determining if a simulated source has been detected. To define a mock source as being ‘detected’, we require a detection in the mock catalogue within a given separation between the input position and the extracted pybdsf position. For fainter simulated flux bins, the cross-match radius is set to be three times the expected rms positional uncertainty for a given SNR following Condon et al. (1998) and assuming a FWHM of 7 arcsec, the typical for a high-SNR compact source (e.g. Shimwell et al. 2022); this results in a cross-match radius of ∼3.5 arcsec for the faintest flux bins. At brighter simulated flux densities, this positional uncertainty becomes very small; we therefore set a minimum cross-match radius of 2 arcsec. These angular separation criteria were validated by examining the change in the number of cross-matches as a function of separation.

3.2.1 Source-size distributions

Completeness depends not only on the flux density but also on the size of the source as source detection is performed based on the peak flux density of a source; therefore, for a given total flux density, the peak flux density for a larger source is more likely to fall below the detection threshold than for a smaller source. However, an accurate determination of the source-size distribution of the sub-mJy radio source population at low frequencies is lacking (see work by Mandal et al. 2021 at characterizing the faint low-frequency source counts) and we must therefore make some assumptions in deriving the corrections.

We start by assuming that the observed size distribution of sources, within a flux density range that is unaffected by completeness, is an accurate description of sizes at fainter flux densities. As our work is focused on generating completeness corrections suited for the AGN subset of the radio population, we generate an ‘AGN’ source-size distribution by selecting all sources classified as radio-excess AGN or SED AGN in the total flux density range |$1\, \rm {mJy} \lt S_{\rm {150,tot}} \lt 10\, \rm {mJy}$|⁠, where we expect the sample to be largely complete. We consider only sources with sizes in the range 6–30 arcsec (along both major and minor axes); larger sizes are not used in our simulations as such sources are poorly represented by a Gaussian surface brightness profile; the small number of sources with larger sizes are all placed at 30 arcsec in our simulations. Within each simulated flux density bin, we weight the simulation output by this size distribution. To determine the completeness correction for each flux density interval, we then considered the subset of mock sources with total flux density >5σ based on the local rms, and determined the fraction of these that were detected by pybdsf with a peak flux density above 5σ (thus matching the selection criteria applied to the observed source sample, and the limits adopted for the V_max calculation). The ‘AGN’ completeness corrections obtained in this way in ELAIS-N1 are shown by the grey line in Fig. 3 (top panel). For comparison, we repeat this process using the size distribution of the star-forming galaxy (SFG) subset of the radio population with |$1\, \rm {mJy} \lt S_{\rm {150,tot}} \lt 10\, \rm {mJy}$|⁠, deriving the pink completeness curve in Fig. 3. The relatively lower completeness for the ‘AGN’ line compared to the ‘SFG’ curve is driven by the sizes of (resolved) radio-AGN being significantly larger than the SFG population; this results in a higher fraction of extended sources (with consequently lower peak flux densities for a given total flux density), which leads to a lower completeness, as expected.

In subsequent analysis, we use the ‘AGN’ completeness corrections in the construction of the luminosity functions. The ‘AGN’ completeness corrections for all three LoTSS Deep fields are shown in the bottom panel of Fig. 3 and listed in Table 1. Using ‘AGN’ based corrections, the observations in LoTSS Deep Fields reach a completeness of 50 per cent (and 90 per cent) at |$150\, \rm {\mu Jy}$| (⁠|$1.3\, \rm {mJy}$|⁠), |$209\, \rm {\mu Jy}$| (⁠|$1.4\, \rm {mJy}$|⁠), and |$289\, \rm {\mu Jy}$| (⁠|$1.4\, \rm {mJy}$|⁠), in ELAIS-N1, Lockman Hole, and Boötes, respectively. We note that in all fields, the completeness does not reach 100 per cent; this is largely due to the source finding algorithm struggling to detect simulated sources placed in the higher noise, and lower dynamic range regions near bright genuine sources. We generate a ‘combined’ completeness curve for use in constructing the LFs by performing an area-weighted average of the completeness curves in the three fields (black line in Fig. 3, bottom).

3.2.2 Application of the completeness corrections

The completeness corrections were applied, following equation (4), by linearly interpolating the corrections at the flux density of each source. We applied a maximum completeness correction of a factor of 10 as any larger corrections are likely not reliable. Then, to determine the point where the completeness corrections to the LFs are too large to be reliable, we recalculated the luminosity functions but this time without applying any completeness corrections (i.e. by setting C_radio[S(z)] = 1 in equation 4); in our analysis we do not plot or list the space densities for the luminosity bins where the difference between the data points with and without the corrections is larger than 0.3 dex. To account for uncertainties in the completeness corrections (e.g. the lack of knowledge of the true source-size distribution), we add 25 per cent of the completeness correction in quadrature to the error obtained from bootstrap sampling at each luminosity bin. This value is motivated by experiments of how the completeness correction varies with different bootstrap samples drawn from the simulations.

To take full advantage of the LoTSS Deep Fields, we calculate the LF in each of the three LoTSS Deep Fields separately, and build a combined LF across the three fields covering |${\sim}25\,$|deg² to obtain more robust number statistics across the full luminosity range and to limit the effects of cosmic variance.

3.3 The local radio-AGN luminosity function

Although the LoTSS Deep Fields cover a relatively small volume at low redshifts, comparison of the low redshift LF against previous measurements in the literature is useful. Using the methods outlined above, we have built the local 0.03 < z ≤ 0.3 luminosity function for the radio-excess AGN (939 sources) in the LoTSS Deep Fields, which is shown in Fig. 4 and listed in Table 2. The combined luminosity function across the three fields is shown by pink filled circles, with the shaded pink region showing the 1σ uncertainties. The luminosity functions for each individual deep field are also shown by pink crosses, triangles, and squares (with their respective error bars) for ELAIS-N1, Lockman Hole, and Boötes, respectively. We note that the ELAIS-N1 data points are consistently higher than the other fields, an effect that is also seen in the K-band number counts (see fig. 2 of Kondapally et al. 2021), likely due to large-scale structure within the field as the volume probed at these low redshifts is relatively small. The Boötes data points are offset lower than the other fields, which at high luminosities is likely due to cosmic variance effects, and at low luminosities may also be from incompleteness. Any differences in the radio flux density calibration between the three fields may also introduce an offset in the LFs.

Figure 4.

Local (0.03 < z ≤ 0.3) radio-excess AGN luminosity functions for LoTSS Deep Fields built with a bin width of |$\Delta \log L = 0.3\,$|dex. The combined luminosity function is shown in pink filled circles (for 939 sources) with shaded regions showing the 1σ uncertainties computed from bootstrapping. The luminosity functions for each field individually are also shown. We find fairly good agreement with the Best & Heckman (2012) and Smolčić et al. (2017b) results from GHz surveys, scaled to 150 MHz using a spectral index α = −0.7.

Open in new tabDownload slide

Fig. 4 also shows comparison between our local radio-AGN luminosity functions and other previous studies. We find good agreement with the local radio-AGN luminosity function from Smolčić et al. (2017b) built using deep radio imaging from the VLA-COSMOS 3 GHz Large project (Smolčić et al. 2017a). They identified radio-AGN using a 3σ redshift-dependent radio-excess selection based on the FIRC following Delvecchio et al. (2017). Their LFs were computed over 0.1 < z < 0.4 (105 sources), and for illustration in Fig. 4 are shifted to 150 MHz using a spectral index α = −0.7. Best & Heckman (2012) presented a large sample of radio-detected AGN drawn from the combination of FIRST and NVSS with SDSS spectroscopic sample data. We find good agreement with the luminosity function of Best & Heckman (2012), shifted to |$\rm {\mathit{ L}_{150\, MHz}}$| using a spectral index α = −0.7, in particular at |$\rm {\mathit{ L}_{150\, MHz}} \sim 10^{23.5} - 10^{25}\, \rm {W\, Hz^{-1}}$|⁠. We are however not able to sample enough volume to probe significantly above the break in the luminosity function. Best & Heckman (2012) also found their LF to be in good agreement with other determinations (e.g. Mauch & Sadler 2007; Pracy et al. 2016).

Also shown in Fig. 4 (green triangles) is the radio-AGN luminosity function from the shallower but wider LoTSS Data Release 1 (LoTSS-DR1) from Sabater et al. (2019). Their data, covering |${\gt} 400\,$| deg², is better suited to sample higher radio luminosities. We find good agreement at moderate luminosities, however our luminosity function is consistently offset to higher space densities by |${\sim}0.5\,$|dex, especially at |$\rm {\mathit{ L}_{150\, MHz}} \lesssim 10^{23.5}\, \rm {W\, Hz^{-1}}$|⁠. We find that the Best & Heckman (2012) LFs appear to turn over at around the same point in luminosity. There are a few possible reasons for the difference with both of those studies. First, neither Best & Heckman (2012) nor Sabater et al. (2019) apply any completeness corrections; we re-derive the Sabater et al. luminosity function by simply applying a 10σ radio flux density cut, which increases the space densities by |${\sim}0.2\,$| dex at |$L_{\rm {150\, MHz}} \lesssim 10^{22.3}\, \rm {W\, Hz^{-1}}$|⁠. We also re-derived our LF without applying completeness corrections; this reduces the difference with the Sabater et al. LF down to the level of field-to-field variations seen across the LoTSS Deep Fields which are a result of cosmic variance and radio flux density calibration differences. Secondly, the Best & Heckman (2012) and Sabater et al. (2019) radio-AGN samples were defined by combining their radio data with the SDSS main galaxy spectroscopic sample; sources like quasars or radio-quiet quasars are missing from this sample, likely biasing the luminosity function to lower space densities. Moreover, their classification scheme was tuned to select only ‘radio-loud AGN’, which excludes Seyfert-like AGN, biasing their luminosity function at low luminosities. Finally, the median redshift of the Best & Heckman (2012) sample is z_med ∼ 0.16, and that of Sabater et al. (2019) is z_med = 0.14, whereas the median redshift for the LoTSS Deep Fields sample (in this redshift bin) is z_med = 0.21; any cosmic evolution of this population can therefore contribute to the difference in space densities observed. We also note that this difference is not simply due to a misclassification of sources between radio-AGN and SFGs, as a similar offset (i.e. higher space densities of SFGs in the LoTSS Deep Fields) is also found by Cochrane et al. (in preparation) when comparing the local SFG luminosity function from the LoTSS Deep Fields with LoTSS-DR1.

3.4 Evolution of the radio-AGN LFs

Fig. 5 shows the redshift evolution of the total radio-AGN LF in the LoTSS Deep fields over the redshift range 0.4 < z ≤ 2.5. The combined LoTSS Deep LF is shown as pink circles, with the LFs from individual fields shown as open pink symbols (using the same symbols as in Fig. 4). For comparison, we show the parametric fit to the local radio-AGN luminosity function by Best et al. (2014) in each panel, shifted to 150 MHz using α = −0.7. We again compare our results with the evolution of the total radio-AGN population presented in Smolčić et al. (2017b), shifted to |$L_{\rm {150\, MHz}}$| using α = −0.7, shown by green symbols, with our redshift bins chosen to match their analysis. As is evident from Fig. 5, we find excellent agreement with their results at all radio luminosities, out to z ∼ 2.5; this gives us further confidence that our source classification method for separating radio-AGN from star-forming galaxies is appropriate out to high redshifts. We do however note the disagreement seen at all luminosities in the 1.3 < z ≤ 1.7 bin, which may be driven by the COSMOS field being overdense at these redshifts (see Duncan et al. 2021; McLeod et al. 2021), where the larger volume probed by LoTSS-Deep allows a more robust determination of the LF.

Figure 5.

Cosmic evolution (across 0.4 < z ≤ 2.5) of the total radio-AGN luminosity functions at 150 MHz for the combined LoTSS Deep sample (pink circles). The LFs in individual LoTSS Deep fields are shown in open pink symbols (same symbols as in Fig. 4). The dashed line in each panel shows the parametrized form of the local radio-AGN LF by Best et al. (2014), converted to 150 MHz using a spectral index α = −0.7 to guide the eye. The total radio-AGN luminosity functions from Smolčić et al. (2017b) and Butler et al. (2019) in green and cyan open points, respectively, converted to 150 MHz using α = −0.7, are also shown. Overall, we find good agreement with the results from both Butler et al. (2019) and Smolčić et al. (2017b) across redshift.

Open in new tabDownload slide

Also shown in Fig. 5 are the radio-AGN LFs of Butler et al. (2019) across 0.3 < z ≤ 0.6, 0.6 < z ≤ 0.9, and 0.9 < z ≤ 1.3 in the closest redshift bins to our LFs. These LFs were compiled using data from the Australian Telescope Compact Array (ATCA) 2.1 GHz observations of the XXL-S field (Butler et al. 2018) and the radio-AGN sample was selected based on radio-source luminosity, morphology, spectral indices, and radio-excess emission based on the FIRC. We also find good agreement with their LFs where available, but our results probe fainter luminosities.

4 COSMIC EVOLUTION OF THE LERG LUMINOSITY FUNCTIONS

The total radio-AGN luminosity functions presented in Section 3.4 contain a mixture of both the LERG and the HERG populations. From studies of the nearby Universe, these two populations are expected to evolve differently with the LERGs dominating the space densities at low luminosities and the HERGs dominating at high radio luminosities (e.g. Best & Heckman 2012); the LERG population is also particularly interesting for radio-AGN feedback cycle considerations. Previous studies (e.g. Best et al. 2014; Pracy et al. 2016; Williams et al. 2018; Butler et al. 2019) have attempted to model the evolution of the LERGs; however, the current LoTSS Deep Fields data set, with vastly greater numbers of LERGs resulting from the combination of deep radio and multiwavelength data sets over |${\sim}25\, \rm {deg}^{2}$|⁠, allows us to study the cosmic evolution of this population in unprecedented detail. The HERGs form a minority of the radio-AGN population in LoTSS-Deep DR1 (see Fig. 1), and moreover, have typically higher radio luminosities than the LERGs (e.g. Best et al. 2014; Pracy et al. 2016); the present LoTSS-Deep sample therefore does not allow us to robustly constrain the LF and the cosmic evolution of the HERGs. In this paper, we list the HERG LFs in Table 3, which are also plotted in Fig. B1, but focus our analysis on the LERG population only in the rest of the paper. Subsequent data releases covering wider areas that better sample the bright end of the luminosity function will enable detailed analysis of the evolution of the HERGs.

In Fig. 6, we show the evolution of the LERG LF for the LoTSS deep fields in four redshift bins (0.5 < z ≤ 1.0, 1.0 < z ≤ 1.5, 1.5 < z ≤ 2.0, and 2.0 < z ≤ 2.5), each spanning 3–4 decades in luminosity. The combined LoTSS deep fields LFs, also tabulated in Table 3, are shown as pink circles. The LFs for individual fields are also shown in pink, using the same symbols for each field as in Fig. 4. In each panel, the black dotted line shows the parametrized form of the local (z < 0.3) LERG luminosity function determined by Best et al. (2014), scaled to 150 MHz using a spectral index α = −0.7. We find that the LERG population shows modest evolution in the range 0 < z < 0.75, especially at high luminosities, although the details of this depend on the assumed spectral index. Beyond z = 0.75, we see a relatively mild evolution in the LFs with redshift out to z ∼ 2.5.

Figure 6.

Cosmic evolution (across 0.5 < z ≤ 2.5) of the LERG luminosity functions at 150 MHz for the combined LoTSS Deep sample (pink circles), showing relatively mild evolution across the redshift bins examined. Also shown are the LFs in individual LoTSS deep fields (in pink; same symbols as in Fig. 4) and the parametric fit to the local LERG LF (black dotted line) from Best et al. (2014), scaled to 150 MHz assuming a spectral index α = −0.7. The pure density evolution (PDE), pure luminosity evolution (PLE), and combined luminosity and density evolution (LDE) models from fits to individual redshift bins are shown as orange (dash-dotted), blue (dotted), and dark blue (dashed) lines, respectively (see Section 4.1 for details). In the 0.5 < z ≤ 1.0 bin, the LERG LF from Best et al. (2014, green triangles), scaled to 150 MHz using a spectral index α = −0.7, shows good agreement with our results. The LFs from Williams et al. (2018), reaching z ∼ 2 (yellow squares in the three lowest redshift bins), have systematically lower space densities compared to our results. Further tests show that the differences are largely driven by our improvements in the source classification method used.

Open in new tabDownload slide

We compare our results to the LERG LF derived by Best et al. (2014) at 0.5 < z ≤ 1 in Fig. 6, shifted to |$L_{\rm {150\, MHz}}$| using a spectral index α = −0.7. Best et al. (2014) compiled catalogues of radio-detected AGN from eight surveys, covering different radio-depths and areas, to obtain a sample of 211 radio-loud AGN. They used spectroscopic information to classify their AGN sample into LERGs and HERGs, representing the first study of the evolution of the two modes of AGN, separately, out to z ∼ 1. The LoTSS-Deep LF shows good agreement with their resulting luminosity function (green triangles), with our data probing fainter luminosities.

Also in Fig. 6, we show the LERG LFs computed by Williams et al. (2018), who studied the cosmic evolution (0.5 < z ≤ 2) of a sample of 1224 LOFAR-detected sources within the Boötes field (albeit with shallower radio data). They used the combination of the FIR-radio correlation of Calistro Rivera et al. (2016) and SED fitting via agnfitter to perform source classifications. Radio sources that lie 2σ away from the FIR-radio correlation (Calistro Rivera et al. 2017) were identified as radio-loud AGN. Then, the various AGN and galaxy models fitted by agnfitter were used to quantify the fraction of MIR emission arising from the AGN compared to the galaxy component (f_AGN); sources with f_AGN > 0.25, indicative of significant emission from the torus, were classified as HERGs, whereas the remainder are classified as LERGs, which are expected to show little to no torus emission. Their final sample contains 243 LERGs (and 398 HERGs) within 0.5 < z ≤ 2, with their resulting LF shown as yellow squares in Fig. 6.

We note that the Williams et al. LFs appear systematically offset to lower space densities than our data set in all redshift bins, and also appear offset by |${\sim}0.4\, \rm {dex}$| at moderate to faint luminosities (⁠|$L_{\rm {150\, MHz}} \lt 10^{26}\, \rm {W\, Hz^{-1}}$|⁠) compared to Best et al. (2014). We have performed tests to investigate the sources of this discrepancy, which are detailed in Appendix  A. First, we require a 0.7 dex (⁠|${\sim}\, 3\sigma$|⁠) radio-excess over the radio luminosity versus SFR relation, whereas Williams et al. (2018) used a 2σ cut (with σ = 0.529) on the redshift-dependent FIRC of Calistro Rivera et al. (2017); therefore the radio-excess AGN classification scheme of Williams et al. (2018) is more conservative than our selection and is expected to result in systematically lower space densities of radio-excess objects. Secondly, we find that improvements in the input models for agnfitter since the analysis of Williams et al. (2018) result in a change in source classification for a significant number of sources, in particular at higher redshifts. Finally, in our analysis, the SFR used for identifying radio-excess AGN is estimated from the consensus measurements from four SED fitting codes by Best et al. (in preparation), whereas Williams et al. used estimates from agnfitter alone. Best et al. (in preparation) show that estimates of SFRs (and infrared luminosities) from agnfitter are systematically higher than those estimated from other SED fitting routines. As detailed in our analysis in Appendix  A, the combination of differences in the SED fitting and source classification criteria results in the apparent discrepancy with Williams et al. (2018).

4.1 Modelling the cosmic evolution of the LFs

We first considered a luminosity and density evolution (LDE) model where we allowed ρ^⋆(z) and L_⋆(z) to be free parameters: both of these parameters were then fitted for each redshift bin using the emcee Markov Chain Monte Carlo (MCMC) fitting routine (Foreman-Mackey et al. 2013), with the resulting fitted ρ^⋆(z) and L_⋆(z) shown in Table 4. This LDE model can also be expressed in terms of a (1 + z) parametrization such that the evolution of the normalization is expressed as |$\rho ^{\star }(z) = \rho ^{\star }(0) \times (1 + z)^{\alpha _{\rm {D}}}$|⁠, and, the redshift evolution of the characteristic luminosity, L^⋆(z), is given as |$L^{\star }(z) = L^{\star }(0) \times (1 + z)^{\alpha _{\rm {L}}}$|⁠. Here, α_D and α_L correspond to the density and luminosity evolution parameters, respectively. This can be used to both illustrate how strong the evolution is, and whether it follows a simple trend with redshift. For these cases, we used the Mauch & Sadler normalization |$\rho ^{\star }(0) = \frac{1}{0.4} 10^{-5.5} \log \rm {L^{-1}}\, \rm {Mpc^{-3}}$| (converted to per log luminosity) and the characteristic luminosity |$L^{\star }(0) = 10^{25.27}\, \rm {W\, Hz^{-1}}$| (scaled to 150 MHz assuming α = −0.7) along with their two slopes. The corresponding values for the evolution parameters for this LDE model, along with the reduced chi-squared goodness of fit values, |$\chi ^{2}_{\nu }$|⁠, are also listed in Table 4.

For each redshift bin, we also considered a pure luminosity evolution (PLE) model (i.e. we set α_D = 0) and a pure density evolution (PDE) model (i.e. we set α_L = 0). The resulting best-fitting α_D and α_L parameters (as appropriate for each model) for each redshift bin are listed in Table 4, along with the |$\chi ^{2}_{\nu }$| values. The resulting best-fitting LFs for all three models at each redshift bin are also shown in Fig. 6. We find that both the PLE and PDE models are not sufficient at describing the evolution of the LERGs, in particular at high radio luminosities, as also noted by the resulting |$\chi ^{2}_{\nu }$| values in Table 4. The LDE model is better able to match the high-luminosity end of the LF. The LDE fits indicate that at least out to z ∼ 2, there is a mild increase in the characteristic luminosity and a corresponding mild decrease in the characteristic space density with increasing redshift. This appears to reverse in the highest redshift bin; however, it should be noted that there is some degeneracy between the fitted L_⋆(z) and ρ_⋆(z), especially due to the lack of a prominent break in the LFs.

Finally, we considered how uncertainties in the source classification criteria employed affect our results on the LERG luminosity functions; this is discussed in detail in Appendix  B. In summary, given the good agreement between our LFs and literature results of the total radio-AGN LFs across redshift, we studied the effects of varying the selection of ‘SED AGN’ on the LERG to HERG classification. We find that the evolution of the LERG population is largely insensitive (at the |${\lesssim}0.2\,$|dex level) to the exact threshold used for separating the two modes of AGN, even under a very conservative definition of HERGs; this gives us confidence in our source classification criteria adopted and on the resulting evolution of the LERG LFs.

To understand the origin of the mild evolution seen in the luminosity functions of the LERGs, we investigate the incidence of the LERGs as a function of stellar mass and redshift in different types of host galaxies in Section 5.

5 PREVALENCE OF AGN ACTIVITY WITH STELLAR MASS AND STAR-FORMATION RATE

It is well-known that the prevalence of radio-AGN increases strongly with stellar mass in the local Universe (e.g. Best et al. 2005b; Smolčić et al. 2009; Janssen et al. 2012; Sabater et al. 2019) and early studies out to z ∼ 1 suggest that this trend holds at earlier cosmic time (e.g. Tasse et al. 2008; Simpson et al. 2013). In this paper, we study the redshift evolution of LERG activity by measuring how the fraction of galaxies that host a LERG varies with stellar mass across redshift.

For the LERG sample, we used the consensus stellar masses derived from SED fitting (see Section 2.4). For the underlying galaxy population, we used the stellar masses (50th percentile of the posterior distribution) computed by Duncan et al. (2021) using a grid-based SED fitting method (see also Duncan et al. 2014, 2019) for the full multiwavelength catalogue in each field (resulting in a total of ∼1.8 million sources after satisfying the redshift and stellar mass completeness limits; see below). Duncan et al. (2021) validated their stellar masses for the population as a whole using comparison with literature galaxy stellar mass functions, and also showed them to be in good agreement on a source-by-source basis (in ELAIS-N1, where comparison was possible); Best et al. (in preparation) also find no systematic offset between the Duncan et al. and the consensus stellar masses (although with a scatter of 0.11 dex for non-AGN, and 0.23 dex for AGN). As detailed by Duncan et al. (2021), the photometric redshifts, and hence the derived stellar masses, are found to be reliable at z < 1.5 for host-galaxy dominated sources; we therefore restrict our analysis in this section to z < 1.5. To avoid biasing our results by stellar mass incompleteness, we restrict analysis to masses above the 90 per cent stellar mass completeness limits estimated by Duncan et al. (2021) for each field separately (given the varying depths of the multiwavelength data). We estimated the stellar mass completeness limit in each of the five redshift bins as the stellar mass above which a source would be detected over the full redshift bin, and simply removed sources with stellar masses below this completeness limit from the analysis.

In Fig. 7, we present the fraction of all galaxies that host a LERG as a function of stellar mass (hereafter; LERG stellar mass fractions) in five redshift bins (0.3 < z ≤ 0.5, 0.5 < z ≤ 0.7, 0.7 < z ≤ 1.0, 1.0 < z ≤ 1.2, 1.2 < z ≤ 1.5), shown by different coloured lines, for three equally spaced radio luminosity limits (shown in different panels) from |$L_{\rm {150\, MHz}} \ge 10^{24}\, \rm {W\, Hz^{-1}}$| to |$L_{\rm {150\, MHz}} \ge 10^{25}\, \rm {W\, Hz^{-1}}$|⁠. The lowest radio luminosity limit chosen corresponds to a 5σ detection (based on noise in the deepest regions) at z ∼ 1.5, and is broadly comparable to the 1.4 GHz luminosity limit (⁠|$L_{\rm {1.4\, GHz}} = 10^{23}\, \rm {W\, Hz^{-1}}$|⁠, assuming α = −0.7) typically used in similar studies in the local Universe (e.g. Best et al. 2005b; Janssen et al. 2012). To account for the lack of a volume-limited sample, we weight each source by a 1/V_max factor based on the maximum volume that a source could be detected in at above 5σ based on the local radio RMS level. The error bars show Poisson uncertainties (following Gehrels (1986) for stellar mass bins with |${N} \rm \lt 10$|⁠). In total, this results in 2722, 1776, and 835 LERGs (within 0.3 < z ≤ 1.5) with |$L_{\rm {150\, MHz}} \ge 10^{24}\, \rm {W\, Hz^{-1}}$|⁠, |$L_{\rm {150\, MHz}} \ge 10^{24.5}\, \rm {W\, Hz^{-1}}$|⁠, and |$L_{\rm {150\, MHz}} \ge 10^{25}\, \rm {W\, Hz^{-1}}$|⁠, respectively.

Figure 7.

Fraction of galaxies hosting a LERG as a function of stellar mass in five redshift bins between 0.3 < z ≤ 1.5 in the LoTSS deep fields. Panels from left to right show this for increasing radio luminosity limits of |$L_{\rm {150\, MHz}} \ge 10^{24}\, \rm {W\, Hz^{-1}}$|⁠, |$L_{\rm {150\, MHz}} \ge 10^{24.5}\, \rm {W\, Hz^{-1}}$|⁠, and |$L_{\rm {150\, MHz}} \ge 10^{25}\, \rm {W\, Hz^{-1}}$|⁠. The error bars represent Poisson uncertainties, following Gehrels (1986) for |$N \rm\lt 10$|⁠. The fraction of galaxies that host a LERG shows a steep stellar mass dependence, in particular at high radio luminosity limits; this relation remains at least out to z ∼ 1.5. We also find evidence of this relation flattening at lower radio luminosity limits suggesting the presence of an additional fuelling mechanism which is investigated in Section 5.1.

Open in new tabDownload slide

For the high radio-power LERGs, corresponding to a radio luminosity limit of |$10^{25}\, \rm {W\, Hz^{-1}}$| (right-hand panel of Fig. 7), we find that the fraction of galaxies that host a LERG has a steep dependence on the stellar mass, approaching ∼10 per cent at the highest stellar masses; these results agree with observations in the local Universe (e.g. Best et al. 2005a; Janssen et al. 2012) and show no signs of redshift evolution. As the radio luminosity limit is lowered (i.e. going from the right-hand panels to left in Fig. 7), we observe an increase in the fraction of LERGs at lower stellar masses at all redshifts, and in particular at higher redshifts, resulting in a shallower dependence of the overall LERG fraction with stellar mass. Physically, the strong dependence on stellar mass is expected from arguments based on the fuelling mechanism of the LERGs (e.g. Best et al. 2006; Hardcastle et al. 2007). The flattening of this relation, in particular at lower radio luminosities and stellar masses where star-forming galaxies are expected to dominate the galaxy stellar mass function, suggests the presence of a star-forming LERG population with a potentially different fuelling mechanism. We investigate this in more detail in Section 5.1 by considering the dependence of AGN activity on both stellar mass and star-formation rate.

5.1 Prevalence of LERGs in star-forming and quiescent galaxies

Fig. 8 shows the SFR-M_⋆ relation for the LERG sample in ELAIS-N1 in the same five redshift bins as in Fig. 7 (0.3 < z ≤ 1.5), split into those hosted by quiescent and star-forming galaxies (red and blue points, respectively), separated using the criterion from equation (8). The black dashed line corresponds to the sSFR < 0.2/t_H(z) threshold used to separate the two populations; this results in 1202 quiescent and 1113 star-forming LERG host galaxies within 0.3 < z < 1.5. We note that the vast majority of the star-forming LERGs lie on or below the main-sequence. An investigation of the rest-frame (u - r) colours (uncorrected for dust reddening) reveals that the median colours for both the star-forming and quiescent LERGs are consistent with being close to the green valley (e.g. Schawinski et al. 2014), as is also evident from Fig. 8. Recent work by Mingo et al. (2022) found a tail of LERGs with high sSFRs based on a study of the resolved radio-loud AGN within LoTSS-Deep; this is consistent with our study that finds the existence of a significant population of LERGs hosted by star-forming galaxies.

Figure 8.

The location of LERGs in the ELAIS-N1 field on the SFR-stellar mass (M_⋆) plane in five redshift bins (0.3 < z ≤ 1.5). The red and blue points correspond to quiescent and star-forming host galaxies of LERGs, respectively. The shaded contours show the SFR-M_⋆ distribution (from magphys) for the underlying IRAC flux-selected sample from Smith et al. (2021) in each redshift bin. The dashed black line corresponds to the sSFR selection threshold adopted for identifying quiescent galaxies (i.e. f_sSFR = 1/5 in equation 8), as described in Section 5.1. Similarly, the green dotted lines, above and below this line correspond to sSFR thresholds based on f_sSFR = 1/3 and f_sSFR = 1/10, respectively. The vertical lines correspond to the 90 per cent stellar mass completeness limits for each redshift bin in ELAIS-N1; sources below this are removed from the analysis of LERG stellar mass fractions in Fig. 9.

Open in new tabDownload slide

The green dotted lines above and below this dividing line correspond to f_sSFR = 1/3 and f_sSFR = 1/10, respectively, in equation (8); we use these variations on the standard sSFR threshold to test the robustness of our results on the prevalence of LERGs (see below). Vertical dashed lines show the stellar mass completeness limits in ELAIS-N1 in each redshift bin that are applied when generating the LERG stellar mass fraction plot. The grey shaded contours in each redshift bin show the SFR-M_⋆ relation for the underlying parent galaxy population in ELAIS-N1; these are drawn such that the outermost contour encompasses 99 per cent of all the underlying sources. For this underlying population, we use the F|$_{\rm {3.6\, {\mu }m}} \gt 10\, \mu$|Jy IRAC selected sample of Smith et al. (2021) in ELAIS-N1, consisting of 183 399 sources for which SED fitting was performed using magphys. We limit the following analysis to the ELAIS-N1 field as a similar magphys SED fitting output for the underlying population is not available in the other two fields. Limiting to the sources with the best available multi-wavelength coverage, good χ² SED fits, and within the chosen redshift range (i.e. 0.3 < z ≤ 1.5) results in a final sample of 140 754 sources in the underlying population that are used in this analysis. We separate this underlying population into star-forming and quiescent galaxies using the same sSFR threshold as that used for the LERG population. As evident from the panels in Fig. 8, the f_sSFR = 1/5 threshold is found to be an appropriate division between star-forming and quiescent galaxies at all redshifts.

Fig. 9 shows the fraction of quiescent galaxies (of a given mass) that host a LERG (left-hand panel), and similarly, the fraction of star-forming galaxies hosting a LERG (right-hand panel), for the same five redshift bins as in Fig. 7, with a radio luminosity limit of |$L_{\rm {150\, MHz}} \ge 10^{24}\, \rm {W\, Hz^{-1}}$|⁠. We note that this radio luminosity limit is not used to define the LERG sample, but rather the luminosity range over which the LERG stellar mass fraction analysis is carried out. The black dashed line also plotted in each panel shows the |$f_{\rm {LERG}} \approx 0.01 (M_{\star }/10^{11}\, \rm {M_{\odot }})^{2.5}$| relationship found for the prevalence of jet-mode AGN with stellar mass in the local Universe (Best et al. 2005a). The exact value of the normalization in this relationship depends on the luminosity limit and hence on the spectral index adopted; the line drawn in Fig. 9 is for illustration purposes only. The fraction of quiescent galaxies hosting a LERG agrees well with this steep stellar mass dependence found in the local Universe, showing essentially no evolution with redshift. In contrast, the LERG stellar mass fractions for the star-forming host galaxies show a weaker dependence on stellar mass, with this relation evolving with redshift such that the fraction increases with increasing redshift at a fixed stellar mass, with the prevalence of these LERGs in star-forming hosts reaching comparable levels to the quiescent hosts at high redshifts.

Figure 9.

Fraction of quiescent galaxies that host a LERG (left-hand panel), and the fraction of star-forming galaxies that host a LERG (right-hand panel) as a function of stellar mass with |$L_{\rm {150\, MHz}} \ge 10^{24}\, \rm {W\, Hz^{-1}}$| in ELAIS-N1. The number in each panel lists the number of quiescent and star-forming LERGs, respectively, within 0.3 < z ≤ 1.5, for the given radio luminosity limit. The error bars represent Poisson uncertainties, following Gehrels (1986) for |$N \rm \lt 10$|⁠. The black dashed line shows the |$f_{\rm {LERG}} \approx 0.01 (M_{\star }/10^{11}\, \rm {M_{\odot }})^{2.5}$| relation found from studies of LERGs in the local Universe; the normalization of this relation however depends on the radio luminosity limit, and hence also on the assumed spectral index. We find that the quiescent LERG stellar mass fractions show a steep stellar mass dependence with essentially no evolution with redshift out to z ∼ 1.5. The star-forming LERGs show a much shallower dependence on stellar mass, with an increase in the fraction at higher redshifts.

Open in new tabDownload slide

Fig. 10 shows the results of this power-law fitting process for the quiescent (red) and star-forming (blue) LERG populations in the same five redshift bins used in Fig. 9. The left-hand panel of Fig. 10 shows the redshift evolution of the power-law slope, β, and right-hand panel shows the redshift evolution of the normalization at |$M_{\star } = 10^{11}\, \rm {M_{\odot }}$| for both populations. The error bars show the 1σ uncertainties in the fitted parameters. We omit data points where the uncertainties are of order or larger than the magnitude of the fitted parameters due to a poor fit; this in particular affected the lowest redshift bin of the star-forming subset due to large uncertainties in the LERG stellar mass fractions. In each panel, the solid lines correspond to f_sSFR = 1/5, our adopted threshold (see Section 5.1 for details) for separating star-forming (blue) and quiescent (red) host-galaxies used throughout this study. Also shown are the results from fitting the LERG stellar mass fractions that were derived based on variations of the f_sSFR selection thresholds equal to 1/3 and 1/10 (see also Fig. 8). The black dashed line in the left-hand panel corresponds to the power-law slope (β = 2.5) found in studies of jet-mode AGN population in the local Universe (e.g. Best et al. 2005a; Janssen et al. 2012).

Figure 10.

Results of power-law fits to the star-forming and quiescent LERG stellar mass fractions with a functional form given by equation (9). Left-hand panel shows the evolution of the power-law slope for quiescent (red lines) and star-forming (blue lines) hosts and the right-hand panel similarly shows the evolution in the normalization of the LERG stellar mass fractions at |$M_{\star } = 10^{11}\, \rm {M_{\odot }}$|⁠. The solid line (with crosses) corresponds to the results obtained from using the f_sSFR = 1/5 threshold adopted throughout this paper for separating star-forming and quiescent galaxies (see Section 5.1). The dashed and dotted lines correspond to results obtained from using variations of the f_sSFR threshold adopted. The black horizontal line in the left-hand panel corresponds to the slope obtained in the local Universe (e.g. Best et al. 2005a; Janssen et al. 2012). We find a clear difference in the slope of the LERG stellar mass fractions between the quiescent and star-forming hosts of LERGs. The steep power-law slope of the quiescent LERGs agrees well with the local Universe measurement at all redshifts, consistent with these AGN being fuelled from the hot gas haloes in massive galaxies; the flatter relation for the star-forming LERGs suggests that a different fuelling mechanism, associated with accretion at low levels from the cold gas within these galaxies, may be present. Confidence for this link with the cold gas also comes from the increase in the normalization of the LERG stellar mass fractions (right-hand panel) for this sub-group at earlier times when gas fractions were higher.

Open in new tabDownload slide

We find that the quiescent and star-forming LERG populations show distinct power-law slopes, with the quiescent LERG stellar mass fractions having a steep slope, close to the local redshift relation, that remains roughly constant with redshift. Moreover, if a stricter definition of ‘quiescence’ is used (i.e. f_sSFR = 1/10; dotted line), we find that the slope of our quiescent LERG population at higher redshifts agrees even better with the local relation. The star-forming LERG stellar mass fractions show a much shallower slope, which is interestingly similar to the relation of the radiative-mode AGN population in the local Universe (e.g. Janssen et al. 2012). As is evident from Fig. 10, the LERG stellar mass fractions resulting from the two alternative selection thresholds for quiescent/star-forming galaxies agree well with the values derived from our adopted ‘quiescence’ criteria (within 1σ based on our uncertainties); this demonstrates that our results are robust to changes in how quiescent and star-forming host-galaxies are selected.

Looking at the evolution of normalization of the power law at |$M = 10^{11}\, \rm {M_{\odot }}$| in Fig. 10 (right-hand panel), we find that for the quiescent LERGs, the normalization stays roughly constant out to z ∼ 1.5. In contrast, the star-forming subset shows a strong increase in the normalization of the LERG stellar mass fraction at higher redshifts, increasing by a factor of ∼4 by z ∼ 1, showing hints of higher prevalence at |$M = 10^{11}\, \rm {M_{\odot }}$| compared to the quiescent hosts at these redshifts.

We note that our use of a radio-excess >3σ criterion may miss low-luminosity radio-AGN, hosted in star-forming galaxies in particular, if the star-formation rate is sufficiently high that a jet with a power close to our adopted radio luminosity limit does not produce a significant radio excess. This may be particularly relevant in massive galaxies, due to the correlation between mass and star-formation rate for star-forming galaxies. Therefore this may impact the results determined in Fig. 9, leading to an artificial flattening of the relation with stellar mass for the star-forming LERGs. We have investigated applying a higher radio-luminosity limit, |$L_{\rm {150\, MHz}} \ge 10^{25}\, \rm {W\, Hz^{-1}}$| (with broader redshift bins 0.5 < z ≤ 1 and 1 < z ≤ 1.5 to account for fewer sources); this ensures a complete sample as radio jets with such high powers will dominate over any star-formation, and thus be selected as radio-excess AGN. Using this higher radio luminosity limit does result in a steeper relation with stellar mass for the star-forming LERGs. We attempted to fit a power law to these star-forming LERG stellar mass fractions, however, such a high radio luminosity limit removes a large fraction of the LERGs from the sample, leading to poorer statistics and hence we were unable to converge on a fit for the lowest redshift bin. For the 1 < z ≤ 1.5 bin, we find a power law slope β = 1.37 ± 0.57, which is ∼2σ away from the power-law slope of β = 2.5 found for the quiescent LERGs in the local Universe. This gives us further confidence that the shallower stellar mass dependence observed for the star-forming LERGs is not purely driven by a luminosity selection effect. It is also important to note that this limit only selects LERGs with the most powerful jets, rather than the ‘typical’ LERG population which we are more interested in. Furthermore, in future work, more robust measurements of the star-forming LERG population will be enabled by upcoming data sets from the WEAVE-LOFAR survey (Smith et al. 2016), which will allow a selection of radio-AGN using emission line diagnostics, overcoming the limitations of the radio-excess selection, and data from sub-arcsecond LOFAR imaging (e.g. Sweijen et al. 2022) which will enable more robust identification of jetted AGN hosted by star-forming galaxies.

5.2 The nature of star-forming and quiescent LERGs

Physically, the differences observed between the quiescent and star-forming hosts of LERGs suggests that the LERGs in these galaxies are fuelled by different mechanisms. For the quiescent LERGs, we observe a steep stellar mass dependence on the LERG fraction, similar to that found in the local Universe (e.g. Best et al. 2005a; Janssen et al. 2012), where the LERG population is dominated by red, quiescent host galaxies. This strong dependence of radio-AGN activity with stellar mass for the LERG population has been shown to be coupled to the cooling rate of hot gas from haloes in massive elliptical galaxies, supporting the idea that these LERGs are fuelled by the hot phase of the intergalactic medium (e.g. Allen et al. 2006; Best et al. 2006; Hardcastle et al. 2007). Since the observed power-law slope (and normalization) for this quiescent LERG population remains roughly constant out to z ∼ 1.5, this suggests that the higher-redshift quiescent LERG population is also fuelled by accretion from the hot medium, as in the local Universe, and that this mechanism has been in place at least since z ∼ 1.5.

In contrast, the LERG activity in star-forming galaxies shows a much flatter dependence with stellar mass, suggesting that a different physical mechanism may be responsible for fuelling the AGN in these galaxies. Given that this population is undergoing a recent episode of star-formation, likely fuelled by the availability of cold gas (e.g. via mergers or secular processes), we may expect this cold gas to provide a fuel source for the black hole as well. Cold gas accretion is often thought to lead to a radiative-mode AGN (which are widely associated with star-forming galaxies), but these star-forming LERGs are understood to be fuelled by a radiatively inefficient accretion flow; hence, the cold gas associated with the on-going star-formation must also be capable of fuelling the black hole. This is consistent with the idea that it is the Eddington-scaled accretion rate on to the black hole that determines the accretion mode (i.e. radiatively efficient/inefficient mode) and not necessarily the fuel source, such that it is possible to achieve accretion of cold material at low Eddington rates, resulting in a radiatively inefficient mode of AGN (e.g. see discussion by Best et al. 2014; Hardcastle 2018; Hardcastle & Croston 2020). In such a scenario, we would expect the HERGs to be the higher accretion rate extension of this star-forming LERG population.

Other evidence in support of this picture comes from the observation of enhanced LERG activity in interacting galaxies (e.g. Sabater et al. 2013; Gordon et al. 2019) and that in lower mass galaxies (⁠|$M_{\star } \lesssim 10^{11}\, \rm {M_{\odot }}$|⁠), the merger fraction of LERGs is higher than that of a control sample (by ∼20 per cent; Gao et al. 2020); these results could be due to the supply of cold gas brought in to the galaxy due to the interaction/merger events. Indeed there is evidence for large reservoirs of cold molecular gas (⁠|${\sim}10^{8} - 10^{10}\, \rm {M_{\odot }}$|⁠) in the host galaxies of LERGs (e.g. Okuda et al. 2005; Ocaña Flaquer et al. 2010; Ruffa et al. 2019a), along with evidence for the cold gas playing a role in the fuelling of the AGN in some systems (e.g. Tremblay et al. 2016; Maccagni et al. 2018; Ruffa et al. 2019b). However, it is also worth noting that while the presence of cold gas within the host galaxy is associated with increased AGN activity (e.g. Vito et al. 2014; Aird, Coil & Georgakakis 2019; Carraro et al. 2020), it does not necessarily imply that this gas reaches and accretes on to the central engine; the cold gas in some systems is mostly found in a relaxed rotating disc (e.g. North et al. 2019) and as such would be expected to lead to only relatively low accretion rates.

The present data are not sufficient for us to determine the physical processes underpinning the AGN fuelling and the direct link with the cold gas; instead detailed characterization of the multiphase medium and its kinematics for matched samples of the three sub-groups of the radio-excess AGN are required. None the less, the trends of shallower stellar mass dependence of the star-forming LERG stellar mass fractions (compared to the quiescent LERGs) and the increase in their normalization at earlier times is consistent with the idea that a different mechanism is responsible for fuelling the LERG activity in these more star-forming galaxies, likely associated with the cold gas component.

6 THE LUMINOSITY FUNCTIONS OF THE QUIESCENT LERGS

Given that we see this significant population of LERGs being hosted in star-forming galaxies at higher redshifts in Section 5, it is informative to study how the luminosity functions of both of these populations evolve with redshift in order to understand the evolution of the LERG population.

Fig. 11 shows the evolution of the total LERG LFs (black points, as shown in Fig. 6), but also split into the subsets that are hosted by quiescent (red triangles) and star-forming (blue squares) galaxies. The error bars represent the 1σ uncertainties with the LFs also tabulated in Table 3. Across the four redshift bins, the LFs are constructed in total for 3406 quiescent and 5024 star-forming LERGs. We find that LERGs hosted by star-forming and quiescent galaxies, are both found across all radio luminosities. At z < 1, the total LERG population is primarily hosted in quiescent galaxies at almost all radio luminosities, with the star-forming population reaching comparable space densities only at the faint end. At higher redshifts, we find a steady decline in the space density of the quiescent LERGs; simultaneously, the space density of the star-forming subset increases steadily with redshift, with a switch-over in the dominant population occurring by z ∼ 1. This star-forming LERG population dominates over the quiescent LERGs across all radio luminosities at z > 1.5. We note that this trend is not simply due to an increase in the number of star-forming galaxies at early epochs as our analysis of the LERG stellar mass fractions (see Figs 9 and 10), which accounts for this effect, shows an increase in the prevalence of this population too. The combination of the evolution of these two sub-groups can thus explain the little evolution seen in the total LERG LFs (see Section 4).

Figure 11.

The evolving LFs for the LERGs in the combined LoTSS deep fields split into the subset hosted by quiescent and star-forming galaxies for the same redshift bins as in Fig. 6. The space densities of the quiescent LERGs dominate the total LERG population at z < 1 and decrease steadily with redshift, however the star-forming LERGs begin to dominate beyond z ∼ 1; the combination of the evolution of these two sub-groups explains the evolution of the total LERG population seen in Section 4. For the quiescent LERGs, we show a broken power-law fit to the data in the 0.5 < z ≤ 1 bin, with the shaded region forming the 1σ uncertainty on the best fit. At higher redshifts, the dashed line and hatched regions represent the best fit and 1σ uncertainties from a 2-parameter LDE fit, whereas the solid line shows an LDE fit where the density evolution is fixed based on the evolution of the expected host galaxies, and can be seen to provide a good match to the data within the 1σ region at each redshift.

Open in new tabDownload slide

We also considered how uncertainties in the separation of LERGs and HERGs affect the LFs of the quiescent and star-forming LERGs; this analysis is detailed in Appendix  B. In summary, we modified and applied additional criteria for selecting ‘SED AGN’ (i.e. those that show signs of MIR emission from an AGN, typical for HERGs) and reconstructed the LFs of the HERGs and the two sub-groups of LERGs. We found that at all redshifts, the quiescent LERGs are a robust population, largely unaffected by the exact threshold used to separate HERGs from the LERGs. Similarly, the star-forming LERGs also show little change in their LFs, even for quite extreme changes in classification criteria, with changes being largely consistent within 2σ of the errors.

Here, we focus on modelling the evolution of the quiescent LERG LFs as this population is particularly important for radio-mode feedback considerations. To do this, we first fitted the 0.5 < z ≤ 1.0 quiescent LERG LFs as a broken power law (with all four free parameters), with the best-fitting model shown as a solid red line in the 0.5 < z ≤ 1 panel of Fig. 11 and the shaded region representing the 1σ uncertainties on the best fit. The fitted parameters are also tabulated in Table 5.

Fitting all parameters for a broken power law independently at higher redshifts for the quiescent LERGs is not possible due to the lack of any strong break in the LFs.³ Instead, we fixed the two power law slopes at the values derived from the fit at 0.5 < z ≤ 1, and adopted an LDE model by fitting directly for ρ_⋆ and L_⋆ at each redshift.⁴ The resulting fitted model and 1σ uncertainties are shown by the dashed lines and hatched regions in Fig. 11. Due to the lack of an obvious break in the LFs, there is a large degeneracy between the evolution of ρ_⋆ and L_⋆, as can be visualized in Fig. 12 which shows a 2D contour plot of the posterior distribution of the LDE model fit in L_⋆(z) and ρ_⋆(z) parameter space for different redshift bins. For the 0.5 < z ≤ 1 redshift bin, the contours shown are from the broken power law fit described above.

Figure 12.

2D contour plot showing the posterior distribution of ρ_⋆(z) and L_⋆(z) for the quiescent LERGs. For 0.5 < z ≤ 1, the contours are based on a broken power-law fit, while at higher redshifts, the contours show the posterior distributions from a 2-parameter LDE fit to the LFs. The large degeneracy in these parameters can be better constrained when the evolution of ρ_⋆(z) is fixed based on the evolution of the expected host galaxies of quiescent LERGs as shown by the dashed lines (see Fig. 13); based on this estimate, the shaded vertical regions correspond to the 1σ uncertainties on the evolved ρ_⋆(z).

Open in new tabDownload slide

Under a simple model, we can consider the evolution of the quiescent LERG LFs to be dictated by the availability of their expected host galaxies (i.e. the massive quiescent galaxies; e.g. Best et al. 2014), which can help reduce the degeneracy of the two-parameter LDE fit. To determine the evolution of the space densities of the hosts of the quiescent LERGs, we follow the method outlined by Best et al. (2014). Best et al. used various determinations of the quiescent stellar mass functions from the literature (from Domínguez Sánchez et al. 2011; Baldry et al. 2012; Ilbert et al. 2013; Moustakas et al. 2013; Muzzin et al. 2013) and combined this with the prevalence of LERGs as a function of stellar mass to predict the total space density of the expected LERG hosts. We use these data compiled by Best et al. (2014), but scale these by a factor of 0.018/0.01 based on the normalization of quiescent LERG stellar mass fractions found in this study (where 0.018 is the mean normalization at |$10^{11}\, \rm {M_{\odot }}$| across 0.3 < z ≤ 1.5, and 0.01 is the normalization derived by Best et al. 2005a, at their 1.4 GHz radio luminosity limit). We also remove the small scale factors that Best et al. (2014) introduced to line up the different data sets. To this analysis, we add additional points based on a more recent measurement of the quiescent stellar mass function by McLeod et al. (2021).

Fig. 13 shows the results of this analysis, where we find that the space density of the available hosts remains relatively flat out to z ∼ 0.75, beyond which there is a sharp decline with redshift. We then determined the evolution of the space density of the expected hosts at the midpoint of our LF redshift bins (i.e. at z = 0.75, 1.25, 1.75, 2.25) by interpolating and computing the mean of the six data sets (shown by black solid points in Fig. 13). This evolution of the hosts was then used constrain the LDE model by scaling the ρ_⋆(z = 0.75) from the broken power-law fit to predict ρ_⋆(z) at higher redshifts, as shown by the dashed lines in the 2D contour plot in Fig. 12 for each redshift bin, with the shaded regions representing the uncertainties in ρ_⋆(z = 0.75) from the broken power-law fit. To test this explicitly, we used the host density evolution determined above to fix the value of ρ_⋆(z) at higher redshifts, and fitted directly for L_⋆(z) only for each of the three 1.0 < z ≤ 2.5 bins. The resulting best-fitting models are shown as red solid lines for the corresponding redshift bins in Fig. 11; these are found to show a good match to the data, lying within the 1σ hatched region at each redshift.

Figure 13.

The evolution of the space density of the expected host-galaxies of the quiescent LERGs from literature data sets at z < 4. The McLeod et al. (2021) data points were determined by combining the stellar mass function of quiescent galaxies with the prevalence of quiescent LERGs with stellar mass (⁠|$f_{\rm {AGN(M_{\star })}} = 0.018 (M_{\star }/10^{11}\, \rm {M_{\odot }})$|⁠; see the text) and integrating in stellar mass (down to |$M_{\star } = 10^{9}\, \rm {M_{\odot }}$|⁠) to obtain the space density of the expected host galaxies. The other data points (from Domínguez Sánchez et al. 2011; Baldry et al. 2012; Ilbert et al. 2013; Moustakas et al. 2013; Muzzin et al. 2013) were complied by Best et al. (2014, taken from their fig. 6) and scaled to match the normalization in |$f_{\rm {AGN}}\rm {(M_{\star })}$| found in this study (see the text). The black filled points show the mean of the data sets obtained from interpolating at the midpoints of our four LF redshift bins.

Open in new tabDownload slide

These results agree well with studies of the total LERG population out to z ≲ 1 that find a mild evolution (e.g. Best et al. 2014; Pracy et al. 2016); at these redshifts, the total LERG population is dominated by the quiescent hosts which show little evolution (see Fig. 13). At higher redshifts, while the total LERG population still shows weak evolution, as also found by Butler et al. (2019, out to z < 1.3), the quiescent LERG population declines strongly with increasing redshift in line with the expected host galaxies, while the characteristic luminosity of the sources increases; this is counterbalanced by an increasing star-forming LERG population at higher redshifts.

7 CONCLUSIONS

In this paper, we have used the first data release of the LoTSS deep fields to study the evolution of the radio-AGN population out to z ∼ 2.5. The LoTSS deep fields represents the deepest low-frequency radio continuum survey to date, covering 25 deg² across the ELAIS-N1, Lockman Hole, and Boötes fields. The depth of the radio and multiwavelength data set, in combination with the wide area coverage, makes this an ideal survey to constrain the evolution of the low-luminosity radio-AGN out to high redshifts.

The total radio-AGN luminosity functions were constructed out to z ≤ 2.5 using a sample of 11 783 radio-AGN; this shows good agreement across all redshifts with the literature LFs determined from GHz surveys. We separated the radio-AGN into LERGs and HERGs via the identification of accretion disc and torus features in galaxy SEDs; the presence of such signatures is indicative of an HERG. This results in a sample of 10 481 LERGs and 1302 HERGs across the deep fields (within z ≤ 2.5), with the LERGs dominating over the HERGs across all luminosities probed by this study. Using this sample, we have then presented the first robust measurement of the cosmic evolution of the LERG luminosity functions out to z ∼ 2.5. We find relatively mild evolution within the redshift bins (0.5 < z ≤ 2.5) examined in this study; this evolution is found to be driven by a very different evolution of those LERGs hosted by quiescent galaxies and those hosted by star-forming galaxies. We therefore investigated these two populations separately.

The main results of our study for the quiescent LERGs are as follows:

• The quiescent LERG LFs dominate the space densities of LERGs across essentially all luminosities at z ≲ 0.75.

• The incidence of LERG activity in quiescent galaxies shows a strong stellar mass dependence, which can be well-described at all redshifts as |$f_{\rm {LERG,Q}} = c (M_{\star } / 10^{11}\, \rm {M_{\odot }})^{\beta }$|⁠, with a power-law slope of β ≈ 2.5 and a normalization of c ≈ 0.018 (for a luminosity limit of |$L_{\rm {150\, MHz}} \ge 10^{24}\, \rm {W\, Hz^{-1}}$|⁠). This relation for the quiescent LERGs also shows no evolution with redshift (between 0.3 < z ≤ 1.5).

• The steep dependence of LERG activity with stellar mass, and the lack of a redshift evolution suggest that the LERGs in these quiescent galaxies are fuelled by cooling of hot gas from haloes, consistent with local Universe observations, at least over the past |${\sim}9\, \rm {Gyr}$|⁠.

• The quiescent LERGs show a strong negative evolution in their LFs beyond z ∼ 0.75, which is consistent with the characteristic space density evolving in accordance with the availability of the expected host galaxies, while there is also an increase in the characteristic luminosity with increasing redshift.

The main results of our study for the star-forming LERGs are as follows:

• The incidence of LERG activity in star-forming galaxies shows a much weaker stellar mass dependence (β ≈ 1.5 consistently across redshift) than that found for the quiescent LERGs in this study, and for LERGs in the local Universe, and is instead comparable to that found for the HERGs locally.

• The incidence of LERGs in star-forming galaxies (at fixed stellar mass; i.e. the normalization c) increases by nearly a factor of four over the range 0.3 < z ≤ 1.5, reaching a higher prevalence than that of the quiescent LERGs at the highest redshifts.

• The star-forming LERGs are a minority of the overall LERG population at z < 1, but space densities of the star-forming LERGs increase with increasing redshift, such that they become the dominant population at z ≳ 1.5 over all radio luminosities.

In the local Universe, LERGs have been argued to be primarily fuelled by cooling of hot gas, typically at low accretion rates, from their massive host galaxies, whereas the incidence of HERGs show a weaker stellar mass dependence, attributed to fuelling via cold gas, typically at high accretion rates. In our study, we find a flatter stellar mass dependence of LERGs hosted by star-forming galaxies (broadly similar to the local HERGs), which increase in prevalence at earlier cosmic times, when gas fractions were higher; this suggests that a different fuelling mechanism, likely associated with the cold gas, is responsible for triggering LERG activity within these star-forming galaxies compared to the LERGs in quiescent galaxies. Therefore, it is possible that the radio-AGN activity in the star-forming LERGs are triggered by a similar mechanism to the HERGs in the local Universe, however, the present data set is not sufficient to fully understand the nature of fuelling of the star-forming LERGs. Detailed characterization of the molecular gas and other host galaxy properties of matched samples of the different subsets of the radio-AGN population are required to understand the physical processes driving the observed differences, and whether the cold gas does indeed play a role in fuelling the LERGs in star-forming galaxies.

The upcoming WEAVE spectrograph in the near future will obtain spectra for all of the radio sources detected in these LOFAR Deep Fields as part of the WEAVE-LOFAR survey (Smith et al. 2016); this will not only enable confirmation of the photometric redshifts, but also more robust host galaxy properties and source classification using emission-line diagnostics out to at least z ∼ 1. Combined with the large samples of AGN detected in these deep fields, this will allow for detailed characterization of the radio-AGN population as a function of multiple parameters simultaneously. Future data releases of the LoTSS deep fields will present deeper and wider radio and multiwavelength data, providing yet further constraints on the faint-end of the AGN luminosity functions.

ACKNOWLEDGEMENTS

We thank the anonymous referee for their helpful comments that have improved the paper. 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. For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.

RK acknowledges support from the Science and Technology Facilities Council (STFC) through an STFC studentship via grant ST/R504737/1. PNB and JS are grateful for support from the UK STFC via grants ST/R000972/1 and ST/V000594/1. We thank James Aird and Seb Oliver for useful discussions on the paper. RKC acknowledges support from the Flatiron Institute which is supported by the Simons Foundation. MJH and DJBS acknowledge support from STFC grant ST/V000624/1. MB 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. KJD acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 892117 (HIZRAD). KM has been supported by the National Science Centre (UMO-2018/30/E/ST9/00082). CLH acknowledges support from the Leverhulme Trust through an Early Career Research Fellowship. BM acknowledges support from the UK STFC under grants ST/R00109X/1,ST/R000794/1, and ST/T000295/1. IP acknowledges support from INAF under the SKA/CTA PRIN ‘FORECaST’ and the PRIN MAIN STREAM ‘SAuROS’ projects. HJAR acknowledges support from the ERC Advanced Investigator programme NewClusters 321271. WLW acknowledges support from the CAS-NWO programme for radio astronomy with project number 629.001.024, which is financed by the Netherlands Organisation for Scientific Research (NWO). This research made use of astropy, a community-developed core Python package for astronomy (Astropy Collaboration 2013, 2018) hosted at http://www.astropy.org/, and of matplotlib (Hunter 2007).

DATA AVAILABILITY

The data set used in this paper is derived from the LoTSS Deep Fields Data Release 1 (Duncan et al. 2021; Kondapally et al. 2021; Sabater et al. 2021; Tasse et al. 2021). The images and catalogues are publicly available from https://lofar-surveys.org/deepfields.html. The luminosity function data presented are available in a machine-readable format at the following repository: https://github.com/rohitk-10/AGN_LF_Kondapally22. Other results presented in the paper are available upon reasonable request to the corresponding author.

Footnotes

REFERENCES