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

R. Kondapally, P.N. Best, R.K. Cochrane, J. Sabater, K.J. Duncan, M.J.

Hardcastle, P. Haskell, B. Mingo, H.J.A. Röttgering, D.J.B. Smith, et al.

► To cite this version:

R. Kondapally, P.N. Best, R.K. Cochrane, J. Sabater, K.J. Duncan, et al. Cosmic evolution of low-excitation radio galaxies in the LOFAR two-metre sky survey deep fields. Monthly Notices of the Royal Astronomical Society, 2022, 513 (3), pp.3742-3767. 10.1093/mnras/stac1128. hal-03658876

HAL Id: hal-03658876

https://hal.science/hal-03658876

Submitted on 23 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

Rohit Kondapally<sup>1,★</sup> Philip N. Best,<sup>1</sup> Rachel K. Cochrane<sup>1,2,3</sup> José Sabater,<sup>1</sup> Kenneth J. Duncan<sup>1</sup>,

Martin J. Hardcastle<sup>1,4</sup> Paul Haskell,<sup>4</sup> Beatriz Mingo<sup>1,5</sup> Huub J. A. Röttgering,<sup>6</sup> Daniel J. B. Smith<sup>1,4</sup>,

Wendy L. Williams,<sup>6</sup> Matteo Bonato,<sup>7,8,9</sup> Gabriela Calistro Rivera<sup>1,10</sup> Fangyou Gao,<sup>11,12</sup>

Catherine L. Hale<sup>1</sup>,<sup>1</sup> Katarzyna Małek,<sup>13,14</sup> George K. Miley,<sup>6</sup> Isabella Prandoni<sup>7</sup> and Lingyu Wang<sup>11,12</sup>

Accepted 2022 April 14. Received 2022 April 13; in original form 2022 January 18

## ABSTRACT

Feedback from low-excitation radio galaxies (LERGs) plays a key role in the lifecycle of massive galaxies in the local Universe; their evolution, and the impact of these active galactic nuclei on early galaxy evolution, however, remain poorly understood. We use a sample of 10 481 LERGs from the first data release of the LOFAR two-metre Sky Survey Deep Fields, covering ~25 deg<fnref n="2" />, to present the first measurement of the evolution of the radio luminosity function (LF) of LERGs out to z ~ 2.5; this shows relatively mild evolution. We split the LERGs into those hosted by quiescent and star-forming galaxies, finding a new dominant population of LERGs hosted by star-forming galaxies at high redshifts. The incidence of LERGs in quiescent galaxies shows a steep dependence on stellar mass out to z ~ 1.5, consistent with local Universe measurements of accretion occurring from cooling of hot gas haloes. The quiescent-LERGs dominate the LFs at z < 1, showing a strong decline in space density with redshift, tracing that of the available host galaxies, while there is an increase in the characteristic luminosity. The star-forming LERG LF increases with redshift, such that this population dominates the space densities at most radio-luminosities by z ~ 1. The incidence of LERGs in star-forming galaxies shows a much weaker stellar-mass dependence, and increases with redshift, suggesting a different fuelling mechanism compared to their quiescent counterparts, potentially associated with the cold gas supply present in the star-forming galaxies.

**Key words:** accretion, accretion discs – galaxies: active – galaxies: evolution – galaxies: jets – galaxies: luminosity function, mass function – radio continuum: galaxies.

## 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 ($\sigma$), of the bulge component of a galaxy (the $M_*$ – $\sigma$ 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

<footnote n="1"><sup>1</sup>SUPA, Institute for Astronomy, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="2"><sup>2</sup>Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA</footnote>

<footnote n="3"><sup>3</sup>Harvard-Smithsonian Center for Astrophysics, 60 Garden St, Cambridge, MA 02138, USA</footnote>

<footnote n="4"><sup>4</sup>Centre for Astrophysics Research, University of Hertfordshire, College Lane, Hatfield AL10 9AB, UK</footnote>

<footnote n="5"><sup>5</sup>School of Physical Sciences, The Open University, Walton Hall, Milton Keynes MK7 6AA, UK</footnote>

<footnote n="6"><sup>6</sup>Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the Netherlands</footnote>

<footnote n="7"><sup>7</sup>INAF-Istituto di Radioastronomia, Via Gobetti 101, I-40129 Bologna, Italy</footnote>

<footnote n="8"><sup>8</sup>Italian ALMA Regional Centre, Via Gobetti 101, I-40129 Bologna, Italy</footnote>

<footnote n="9"><sup>9</sup>INAF-Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, I-35122 Padova, Italy</footnote>

<footnote n="10"><sup>10</sup>European Southern Observatory, Karl-Schwarzschild-Straße 2, D-85748 Garching bei München, Germany</footnote>

<footnote n="11"><sup>11</sup>SRON Netherlands Institute for Space Research, Landleven 12, NL-9747 AD Groningen, the Netherlands</footnote>

<footnote n="12"><sup>12</sup>Kapteyn Astronomical Institute, University of Groningen, Postbus 800, NL-9700 AV Groningen, the Netherlands</footnote>

<footnote n="13"><sup>13</sup>National Centre for Nuclear Research, ul. Pasteura 7, PL-02-093 Warszawa, Poland</footnote>

<footnote n="14"><sup>14</sup>Aix Marseille Univ. CNRS, CNES, LAM, Marseille 13013, France</footnote>

<footnote>*E-mail: rohitkondapally@gmail.com</footnote>
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 \sim 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_{\star}^{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 \sim 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_{1.4\,GHz} \lesssim 10^{24} \, \mathrm{W}\,\mathrm{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 \sim 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 \le 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 \sim 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 $\sim 25$ deg$^2$, 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 $\Lambda$CDM cosmology with $\Omega_M = 0.3$, $\Omega_\Lambda = 0.7$, and $H_0 = 70$ km s$^{-1}$ Mpc$^{-1}$, and a radio spectral index $\alpha = -0.7$ (where $S_v \propto v^\alpha$). 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.<fnref n="1" /> Multi-epoch observations were used to build deep radio images, which cover 68 deg$^2$ 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 µJy beam$^{-1}$ 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 $\chi^2$ 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 $l$-band and 4.5 µm band selected catalogues (Brown et al. 2007, 2008) to generate a combined multiwavelength catalogue, as described by Kondapally et al. (2021).

<footnote n="1">¹Hereafter, we refer to the central frequencies for all fields as 150 MHz, for simplicity.</footnote>
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 µm, and imaging from the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010; with photometry at 100 and 160 µm), and the Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010; with bands at 250, 360, and 520 µ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,<fnref n="2" /> totalling ~26 deg<fnref n="2" />.

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 <0.15 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

<footnote n="2">²See FLAG_OVERLAP criteria described in Kondapally et al. (2021, their table 5)</footnote>
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_{AGN,16th}$ 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_{150MHz} \lesssim 1$ 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 ($>0.7$ dex, $\approx 3\sigma$) 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 $<0.7$ 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\sigma$ 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

<figure><img src="image_3.png" /><figcaption>**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\sigma$ detection limit using $\sigma = 20$ µJy (achieved in the central region of ELAIS-N1) and assuming a spectral index $\alpha = -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.</figcaption></figure>

for details) to ensure a clean and robust sample. In addition, we limit our analysis of the luminosity functions to $0.03 < z \le 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\sigma$ detection limit calculated based on the noise in the central region of ELAIS-N1, $\sigma = 20$ µJy, and a radio spectral index $\alpha = -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.

# 3 TOTAL RADIO-AGN LUMINOSITY FUNCTIONS

## 3.1 Building the luminosity functions

We calculate the rest-frame radio luminosity of our sources by assuming that the radio spectrum is described by a simple power law in frequency, $\nu$, with $S_\nu \propto \nu^\alpha$, where $S_\nu$ is the radio flux density at frequency $\nu$ and $\alpha = -0.7$ is the assumed spectral index throughout this study (e.g. Calistro Rivera et al. 2017; Murphy et al. 2017). The radio luminosity $L_\nu$ can be computed using the radio flux density $S_\nu$ as

$$L_{\nu} = \frac{4\pi D_{L}^{2}(z)}{(1+z)^{1+\alpha}} S_{\nu}, \qquad (1)$$

where $D_L$ is the luminosity distance to the source and the $(1+z)^{-(1+\alpha)}$ term accounts for the radio spectrum k-correction.

The luminosity functions (LFs) were built using the standard $1/V_{max}$ technique (Schmidt 1968; Condon 1989), which weights each source in the sample by the maximum volume $V_{max}$ that the source could be observed in, given the potential redshift range, and still satisfy all selection effects to be included in the sample. This is particularly important for surveys like the LoTSS Deep Fields where the RMS in the radio images (and hence the $5\sigma$ flux density limit) varies as a function of the position in the field due to primary beam effects, increased RMS around bright sources, and facet-to-facet variations in radio data calibration. The luminosity function $\rho(L, z)$ gives the number of sources per unit co-moving volume observed per unit of log luminosity, and is given by

$$\rho(L, z) = \frac{1}{\Delta \log L} \sum_{i=1}^{N} \frac{1}{V_{\text{max},i}}, \qquad (2)$$

where $\Delta \log L = 0.3$ is the luminosity bin width in log-space, with the sum calculated over each source $i$ in a given luminosity and redshift bin.

The $V_{max,i}$ for a given source is then computed as

$$V_{\text{max},i} = \int_{z_{\text{min}}}^{z_{\text{max}}} V(z) \theta(S, z) dz, \qquad (3)$$

where $V(z)dz$ is the comoving volume across the whole sky between redshift $z$ and $z + dz$, $\theta(S, z)$ is the fractional sky coverage of the LoTSS Deep Fields that accounts for the non-uniform radio-map noise and radio flux density incompleteness, and $S$ is the radio flux density that a source $i$ with a given intrinsic luminosity would have at redshift $z$. In practice, we evaluate this integral numerically, with a step size of $\Delta z = 0.0001$ between $z_{min}$ and $z_{max}$, the minimum and

<figure><img src="image_4.png" /><figcaption>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 $\approx 24.6$ deg$^2$.</figcaption></figure>

maximum range of the redshift bin. At each redshift, $\theta(S, z)$ for each source $i$ is given as

$$\theta(S, z) = \frac{\Omega[S(z)]}{4\pi} \times C_{\text{radio}}[S(z)], \qquad (4)$$

where $\Omega[S(z)]$ is the solid angle (in units of sr) of the survey area in which a source $i$ with a flux density $S$ can be detected at $5\sigma$ given the non-uniform radio image noise, and $C_{\text{radio}}[S(z)]$ is the radio flux density completeness correction at flux density $S$. The $\Omega[S(z)]$ term is computed by performing a linear interpolation using the cumulative area versus noise plot for each field, shown in Fig. 2; these areas are obtained after limiting to the best-available multiwavelength regions in each field and applying optical-based bright star masking (see Kondapally et al. 2021). Completeness corrections are required to account for the faint undetected radio sources close to the survey detection limit, otherwise leading to an underestimate of the space density; these are discussed in detail in Section 3.2.

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 \le 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 $\chi^2$ 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\sigma$ 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 dex in the range $0.4 < S < 40$ mJy and, using a finer sampling interval of 0.08 dex in the range $90 < S < 400$ µ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 ~110 and ~190 µ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.

<table><caption>**Table 1.** Completeness corrections determined for the LOFAR Deep Fields calculated using a source-size distribution derived from the observed sources classified as AGN by Best et al. (in preparation) with $1$ mJy $< S_{150, tot} \le 10$ mJy in each field. The completeness curves for each field are shown in Fig. 3 (bottom).</caption><thead><tr><th>Flux density (mJy)</th><th>ELAIS-N1</th><th>Lockman Hole</th><th>Boötes</th></tr></thead><tbody><tr><td>0.09</td><td>0.143</td><td>—</td><td>—</td></tr><tr><td>0.11</td><td>0.221</td><td>—</td><td>—</td></tr><tr><td>0.13</td><td>0.351</td><td>0.193</td><td>—</td></tr><tr><td>0.16</td><td>0.553</td><td>0.319</td><td>—</td></tr><tr><td>0.19</td><td>0.68</td><td>0.43</td><td>0.211</td></tr><tr><td>0.23</td><td>0.742</td><td>0.57</td><td>0.322</td></tr><tr><td>0.28</td><td>0.778</td><td>0.664</td><td>0.486</td></tr><tr><td>0.33</td><td>0.802</td><td>0.751</td><td>0.624</td></tr><tr><td>0.4</td><td>0.813</td><td>0.779</td><td>0.717</td></tr><tr><td>0.63</td><td>0.849</td><td>0.841</td><td>0.816</td></tr><tr><td>1.01</td><td>0.884</td><td>0.872</td><td>0.871</td></tr><tr><td>1.59</td><td>0.921</td><td>0.907</td><td>0.912</td></tr><tr><td>2.52</td><td>0.951</td><td>0.941</td><td>0.952</td></tr><tr><td>4.0</td><td>0.961</td><td>0.959</td><td>0.971</td></tr><tr><td>6.34</td><td>0.973</td><td>0.967</td><td>0.982</td></tr><tr><td>10.05</td><td>0.975</td><td>0.974</td><td>0.981</td></tr><tr><td>15.92</td><td>0.978</td><td>0.973</td><td>0.98</td></tr><tr><td>25.24</td><td>0.984</td><td>0.978</td><td>0.989</td></tr><tr><td>40.0</td><td>0.987</td><td>0.978</td><td>0.985</td></tr></tbody></table>

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$ mJy $< S_{150, tot} < 10$ 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\sigma$ based on the local rms, and determined the fraction of these that were detected by PYBDSF with a peak flux density above $5\sigma$ (thus matching the selection criteria applied to the observed source sample, and the limits adopted for the $V_{max}$ calculation). The ‘AGN’
<figure><img src="image_5.png" /><figcaption>**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 \times$ 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.</figcaption></figure>

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$ mJy $< S_{150, tot} < 10$ 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.

To confirm the robustness of the above approach, we consider also the flux-dependent angular size distribution based on GHz surveys commonly used in the literature. Windhorst, Mathis & Neuschaefer

(1990) describe the integral angular size distribution as

$$h(> \psi) = \exp \left[ -\ln 2 \left( \frac{\psi}{\psi_{\text{med}}} \right)^{0.62} \right], \quad (5)$$

where $h(> \psi)$ is the integral angular size distribution for sources with angular sizes larger than $\psi$ (in arcsec) at a given flux density and $\psi_{med}$ is the median angular size at the given flux density. Windhorst et al. (1990, 1993) proposed a relationship between $\psi_{med}$ and the radio flux at 1.4 GHz, $S_{1.4 GHz}$.

$$\psi_{\text{med}} = 2(S_{1.4\text{GHz}} / \text{mJy})^{0.3} \text{arcsec}, \quad (6)$$

which was converted to 150 MHz using a spectral index $\alpha = -0.7$. They also considered a potential floor in this relationship at a size of 2 arcsec, however given that this needs to be convolved with the 6 arcsec LOFAR beam, this makes no significant difference to the observed size distribution. Using equations (5) and (6), $h(> \psi)$ was computed for each flux density interval and the simulation outputs were weighted by this to determine the completeness corrections. The resulting correction in the ELAIS-N1 field is shown in Fig. 3 (top panel; blue line). Other low radio frequency studies in the literature (e.g. Mahony et al. 2016; Williams et al. 2016; Retana-Montenegro et al. 2018) find good agreement between the LOFAR data and the Windhorst et al. (1990) size distributions if the median angular size relation (equation 6) is scaled by a factor of two. We therefore also compute the completeness corrections for this relation, also shown in the top panel of Fig. 3 (green line), with the larger median sizes resulting in lower completeness. We find a good agreement between the $2 \times$ Windhorst et al. (1990) and the ‘AGN’ curves, and similarly between the Windhorst et al. (1990) and the ‘SFG’ curves.

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 µJy (1.3 mJy), 209 µJy (1.4 mJy), and 289 µJy (1.4 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_{\text{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.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Local (0.03 < z ≤ 0.3) radio-excess AGN luminosity functions for LoTSS Deep Fields built with a bin width of Δ 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.</figcaption></figure>

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 ~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.

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

<table><caption>**Table 2.** The local (0.03 < z ≤ 0.3) radio-excess AGN LF for LoTSS deep fields.</caption><thead><tr><th>log L<sub>150 MHz</sub><br>log (W Hz<sup>-1</sup>)</th><th>log ρ<br>log (Mpc<sup>-3</sup> log L<sup>-1</sup>)</th><th>N</th></tr></thead><tbody><tr><td>21.15</td><td>-3.07<sup>+0.18</sup><sub>-0.21</sub></td><td>6</td></tr><tr><td>21.45</td><td>-3.00<sup>+0.11</sup><sub>-0.12</sub></td><td>20</td></tr><tr><td>21.75</td><td>-3.21<sup>+0.09</sup><sub>-0.10</sub></td><td>34</td></tr><tr><td>22.05</td><td>-3.33<sup>+0.07</sup><sub>-0.08</sub></td><td>67</td></tr><tr><td>22.35</td><td>-3.42<sup>+0.06</sup><sub>-0.07</sub></td><td>135</td></tr><tr><td>22.65</td><td>-3.61<sup>+0.05</sup><sub>-0.05</sub></td><td>187</td></tr><tr><td>22.95</td><td>-3.77<sup>+0.04</sup><sub>-0.04</sub></td><td>176</td></tr><tr><td>23.25</td><td>-3.98<sup>+0.04</sup><sub>-0.04</sub></td><td>119</td></tr><tr><td>23.55</td><td>-4.34<sup>+0.06</sup><sub>-0.07</sub></td><td>54</td></tr><tr><td>23.85</td><td>-4.42<sup>+0.06</sup><sub>-0.07</sub></td><td>47</td></tr><tr><td>24.15</td><td>-4.82<sup>+0.08</sup><sub>-0.13</sub></td><td>19</td></tr><tr><td>24.45</td><td>-4.74<sup>+0.07</sup><sub>-0.11</sub></td><td>23</td></tr><tr><td>24.75</td><td>-4.85<sup>+0.09</sup><sub>-0.11</sub></td><td>18</td></tr><tr><td>25.05</td><td>-4.99<sup>+0.12</sup><sub>-0.11</sub></td><td>13</td></tr><tr><td>25.35</td><td>-4.90<sup>+0.10</sup><sub>-0.12</sub></td><td>16</td></tr><tr><td>25.95</td><td>-5.63<sup>+0.22</sup><sub>-0.48</sub></td><td>3</td></tr></tbody></table>

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 L<sub>150 MHz</sub> using a spectral index α = −0.7, in particular at L<sub>150 MHz</sub> ~ 10<sup>23.5</sup> W Hz<sup>-1</sup>. 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 >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 ~0.5 dex, especially at L<sub>150 MHz</sub> ~ 10<sup>23.5</sup> W Hz<sup>-1</sup>. 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 ~0.2 dex at L<sub>150 MHz</sub> ~ 10<sup>22.3</sup> W Hz<sup>-1</sup>. 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
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Cosmic evolution (across $0.4 < z \le 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 $\alpha = -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 $\alpha = -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.</figcaption></figure>

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} \sim 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 \le 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 $\alpha = -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_{150\,MHz}$ using $\alpha = -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 \sim 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 \le 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.

Also shown in Fig. 5 are the radio-AGN LFs of Butler et al. (2019) across $0.3 < z \le 0.6$, $0.6 < z \le 0.9$, and $0.9 < z \le 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 ~25 deg², 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 \le 1.0$, $1.0 < z \le 1.5$, $1.5 < z \le 2.0$, and $2.0 < z \le 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 $\alpha = -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 \sim 2.5$.

We compare our results to the LERG LF derived by Best et al. (2014) at $0.5 < z \le 1$ in Fig. 6, shifted to $L_{150\text{MHz}}$ using a spectral index $\alpha = -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 \sim 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 \le 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\sigma$ 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_{\text{AGN}}$); sources with $f_{\text{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 \le 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 ~0.4 dex at moderate to faint luminosities ($L_{150\text{MHz}} < 10^{26} \text{ 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 (~$3\sigma$) radio-excess over the radio luminosity versus SFR relation, whereas Williams et al. (2018) used a $2\sigma$ cut (with $\sigma = 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

The radio luminosity function of AGN is generally modelled as a broken power law of the form

$$ \rho(L) = \frac{\rho^*}{(L^*/L)^\beta + (L^*/L)^\gamma}, \qquad (7) $$

where $\rho^*$ is the characteristic space density, $L^*$ is the characteristic luminosity, and, $\beta$ and $\gamma$ are the bright and faint end slopes, respectively (e.g. Dunlop & Peacock 1990). As evident from comparison with the $0.5 < z \le 1$ Best et al. (2014) LF in Fig. 6, the area covered by the first data release of the LoTSS deep fields is not sufficient to probe bright enough luminosities to constrain the bright-end slope of the LERG LF, at all redshifts. Therefore, we modelled the evolution of the LERG population as the luminosity evolution and density evolution of the local LF such that in equation (7), $\rho^*(z)$ accounts for the redshift evolution of the normalization and $L^*(z)$ is the redshift evolution of the characteristic luminosity. For this process, we assumed that the shape of the LF remains constant by fixing the bright-end slope $\beta = -1.27$, and the faint-end slope $\gamma = -0.49$, as found by a broken power-law fit to the local radio-AGN LF by Mauch & Sadler (2007). Although the Mauch & Sadler LF includes both the LERGs and HERGs, the faint-end slope, which is the key parameter to constrain as our data do not probe the bright end well, will be dominated by the LERGs; the Mauch & Sadler (2007) faint-end slope provides a better match to our data set than, for example, the jet-mode AGN LF of Best et al. (2014). We also found that a broken power-law fit to the radio-AGN LF of Sabater et al. (2019) gives a faint-end slope consistent with the Mauch & Sadler (2007) value, and hence gives similar results in modelling the evolution below.

We first considered a luminosity and density evolution (LDE) model where we allowed $\rho^*(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 $\rho^*(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^*(z) = \rho^*(0) \times (1+z)^{\alpha_D}$, and, the
<table><caption>**Table 3.** Luminosity functions of the LERGs and HERGs in the LoTSS deep fields in the range $0.5 < z \le 2.5$, illustrated in Fig. 6. Also shown are the LFs in the same redshift bins, for the subsets of the LERG population that are hosted by quiescent galaxies and by star-forming galaxies, respectively; these resulting LFs are discussed in more detail in Section 6, with the separation of quiescent and star-forming host galaxies of the LERGs described in Section 5.1. Space densities for bins with $N < 2$ are not shown. A machine-readable version of these data is available at https://github.com/rohitk-10/AGN_LF_Kondapally22.</caption><thead><tr><th>z</th><th>log L<sub>150 MHz</sub><br>log (W Hz<sup>-1</sup>)</th><th>log ρ<sub>LERG</sub><br>log (Mpc<sup>-3</sup> log L<sup>-1</sup>)</th><th>N<sub>LERG</sub></th><th>log ρ<sub>LERG,Q</sub><br>log (Mpc<sup>-3</sup> log L<sup>-1</sup>)</th><th>N<sub>LERG,Q</sub></th><th>log ρ<sub>LERG,SF</sub><br>log (Mpc<sup>-3</sup> log L<sup>-1</sup>)</th><th>N<sub>LERG,SF</sub></th><th>log ρ<sub>HERG</sub><br>log (Mpc<sup>-3</sup> log L<sup>-1</sup>)</th><th>N<sub>HERG</sub></th></tr></thead><tbody><tr><td rowspan="12">0.5 &lt; z &le; 1.0</td><td>23.55</td><td>-4.07<sup>+0.07</sup><sub>-0.07</sub></td><td>428</td><td>-4.26<sup>+0.07</sup><sub>-0.07</sub></td><td>279</td><td>-4.53<sup>+0.08</sup><sub>-0.08</sub></td><td>149</td><td>-5.15<sup>+0.10</sup><sub>-0.11</sub></td><td>35</td></tr><tr><td>23.85</td><td>-4.22<sup>+0.05</sup><sub>-0.05</sub></td><td>756</td><td>-4.48<sup>+0.05</sup><sub>-0.05</sub></td><td>411</td><td>-4.56<sup>+0.05</sup><sub>-0.05</sub></td><td>345</td><td>-5.29<sup>+0.07</sup><sub>-0.07</sub></td><td>66</td></tr><tr><td>24.15</td><td>-4.39<sup>+0.03</sup><sub>-0.03</sub></td><td>719</td><td>-4.69<sup>+0.03</sup><sub>-0.03</sub></td><td>361</td><td>-4.69<sup>+0.03</sup><sub>-0.04</sub></td><td>358</td><td>-5.39<sup>+0.06</sup><sub>-0.06</sub></td><td>72</td></tr><tr><td>24.45</td><td>-4.70<sup>+0.03</sup><sub>-0.03</sub></td><td>387</td><td>-5.01<sup>+0.03</sup><sub>-0.03</sub></td><td>190</td><td>-4.99<sup>+0.03</sup><sub>-0.03</sub></td><td>197</td><td>-5.41<sup>+0.05</sup><sub>-0.05</sub></td><td>76</td></tr><tr><td>24.75</td><td>-4.93<sup>+0.03</sup><sub>-0.03</sub></td><td>242</td><td>-5.15<sup>+0.03</sup><sub>-0.04</sub></td><td>144</td><td>-5.32<sup>+0.04</sup><sub>-0.04</sub></td><td>98</td><td>-5.63<sup>+0.06</sup><sub>-0.07</sub></td><td>48</td></tr><tr><td>25.05</td><td>-5.15<sup>+0.03</sup><sub>-0.04</sub></td><td>151</td><td>-5.31<sup>+0.04</sup><sub>-0.04</sub></td><td>104</td><td>-5.65<sup>+0.05</sup><sub>-0.07</sub></td><td>47</td><td>-5.91<sup>+0.08</sup><sub>-0.09</sub></td><td>26</td></tr><tr><td>25.35</td><td>-5.28<sup>+0.04</sup><sub>-0.04</sub></td><td>111</td><td>-5.42<sup>+0.04</sup><sub>-0.05</sub></td><td>82</td><td>-5.87<sup>+0.08</sup><sub>-0.10</sub></td><td>29</td><td>-6.33<sup>+0.11</sup><sub>-0.15</sub></td><td>10</td></tr><tr><td>25.65</td><td>-5.41<sup>+0.04</sup><sub>-0.05</sub></td><td>83</td><td>-5.51<sup>+0.05</sup><sub>-0.06</sub></td><td>67</td><td>-6.13<sup>+0.10</sup><sub>-0.12</sub></td><td>16</td><td>-6.43<sup>+0.10</sup><sub>-0.20</sub></td><td>8</td></tr><tr><td>25.95</td><td>-5.60<sup>+0.06</sup><sub>-0.06</sub></td><td>54</td><td>-5.72<sup>+0.06</sup><sub>-0.07</sub></td><td>41</td><td>-6.22<sup>+0.11</sup><sub>-0.15</sub></td><td>13</td><td>-6.56<sup>+0.12</sup><sub>-0.18</sub></td><td>6</td></tr><tr><td>26.25</td><td>-5.88<sup>+0.07</sup><sub>-0.10</sub></td><td>29</td><td>-5.94<sup>+0.07</sup><sub>-0.08</sub></td><td>25</td><td>-6.74<sup>+0.22</sup><sub>-0.48</sub></td><td>4</td><td>-6.86<sup>+0.22</sup><sub>-0.48</sub></td><td>3</td></tr><tr><td>26.55</td><td>-6.06<sup>+0.09</sup><sub>-0.10</sub></td><td>19</td><td>-6.20<sup>+0.11</sup><sub>-0.15</sub></td><td>14</td><td>-6.64<sup>+0.15</sup><sub>-0.22</sub></td><td>5</td><td>-6.50<sup>+0.11</sup><sub>-0.24</sub></td><td>7</td></tr><tr><td>26.85</td><td>-6.74<sup>+0.18</sup><sub>-0.30</sub></td><td>4</td><td>-6.86<sup>+0.22</sup><sub>-0.48</sub></td><td>3</td><td>-</td><td>-</td><td>-6.74<sup>+0.18</sup><sub>-0.30</sub></td><td>4</td></tr><tr><td>27.15</td><td>-6.86<sup>+0.22</sup><sub>-0.48</sub></td><td>3</td><td>-7.04<sup>+0.18</sup><sub>-0.30</sub></td><td>2</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>27.45</td><td>-6.74<sup>+0.18</sup><sub>-0.30</sub></td><td>4</td><td>-</td><td>-</td><td>-6.86<sup>+0.22</sup><sub>-0.48</sub></td><td>3</td><td>-6.56<sup>+0.18</sup><sub>-0.18</sub></td><td>6</td></tr><tr><td rowspan="12">1.0 &lt; z &le; 1.5</td><td>24.15</td><td>-4.48<sup>+0.07</sup><sub>-0.07</sub></td><td>396</td><td>-4.84<sup>+0.08</sup><sub>-0.08</sub></td><td>166</td><td>-4.73<sup>+0.07</sup><sub>-0.07</sub></td><td>230</td><td>-6.02<sup>+0.12</sup><sub>-0.15</sub></td><td>14</td></tr><tr><td>24.45</td><td>-4.66<sup>+0.04</sup><sub>-0.04</sub></td><td>545</td><td>-5.16<sup>+0.05</sup><sub>-0.05</sub></td><td>170</td><td>-4.82<sup>+0.04</sup><sub>-0.04</sub></td><td>375</td><td>-5.70<sup>+0.07</sup><sub>-0.07</sub></td><td>51</td></tr><tr><td>24.75</td><td>-4.87<sup>+0.03</sup><sub>-0.03</sub></td><td>403</td><td>-5.43<sup>+0.04</sup><sub>-0.04</sub></td><td>111</td><td>-5.01<sup>+0.03</sup><sub>-0.03</sub></td><td>292</td><td>-5.77<sup>+0.06</sup><sub>-0.07</sub></td><td>50</td></tr><tr><td>25.05</td><td>-5.15<sup>+0.03</sup><sub>-0.03</sub></td><td>224</td><td>-5.63<sup>+0.05</sup><sub>-0.05</sub></td><td>74</td><td>-5.33<sup>+0.03</sup><sub>-0.04</sub></td><td>150</td><td>-5.92<sup>+0.06</sup><sub>-0.08</sub></td><td>38</td></tr><tr><td>25.35</td><td>-5.42<sup>+0.04</sup><sub>-0.05</sub></td><td>127</td><td>-5.80<sup>+0.06</sup><sub>-0.06</sub></td><td>53</td><td>-5.65<sup>+0.05</sup><sub>-0.05</sub></td><td>74</td><td>-6.18<sup>+0.07</sup><sub>-0.11</sub></td><td>22</td></tr><tr><td>25.65</td><td>-5.56<sup>+0.04</sup><sub>-0.04</sub></td><td>94</td><td>-5.99<sup>+0.07</sup><sub>-0.08</sub></td><td>35</td><td>-5.77<sup>+0.05</sup><sub>-0.05</sub></td><td>59</td><td>-6.53<sup>+0.11</sup><sub>-0.15</sub></td><td>10</td></tr><tr><td>25.95</td><td>-5.77<sup>+0.05</sup><sub>-0.06</sub></td><td>58</td><td>-6.09<sup>+0.07</sup><sub>-0.09</sub></td><td>28</td><td>-6.06<sup>+0.07</sup><sub>-0.10</sub></td><td>30</td><td>-6.58<sup>+0.12</sup><sub>-0.18</sub></td><td>9</td></tr><tr><td>26.25</td><td>-5.98<sup>+0.07</sup><sub>-0.08</sub></td><td>36</td><td>-6.20<sup>+0.08</sup><sub>-0.10</sub></td><td>22</td><td>-6.39<sup>+0.09</sup><sub>-0.16</sub></td><td>14</td><td>-6.46<sup>+0.12</sup><sub>-0.12</sub></td><td>12</td></tr><tr><td>26.55</td><td>-6.27<sup>+0.08</sup><sub>-0.10</sub></td><td>19</td><td>-6.59<sup>+0.12</sup><sub>-0.18</sub></td><td>9</td><td>-6.54<sup>+0.12</sup><sub>-0.18</sub></td><td>10</td><td>-6.84<sup>+0.15</sup><sub>-0.22</sub></td><td>5</td></tr><tr><td>26.85</td><td>-6.25<sup>+0.10</sup><sub>-0.12</sub></td><td>20</td><td>-6.70<sup>+0.11</sup><sub>-0.24</sub></td><td>7</td><td>-6.43<sup>+0.12</sup><sub>-0.11</sub></td><td>13</td><td>-7.07<sup>+0.22</sup><sub>-0.48</sub></td><td>3</td></tr><tr><td>27.15</td><td>-6.37<sup>+0.10</sup><sub>-0.13</sub></td><td>15</td><td>-6.94<sup>+0.18</sup><sub>-0.30</sub></td><td>4</td><td>-6.50<sup>+0.11</sup><sub>-0.14</sub></td><td>11</td><td>-7.25<sup>+0.18</sup><sub>-0.30</sub></td><td>2</td></tr><tr><td>27.45</td><td>-6.94<sup>+0.18</sup><sub>-0.30</sub></td><td>4</td><td>-7.25<sup>+0.18</sup><sub>-0.30</sub></td><td>2</td><td>-7.25<sup>+0.18</sup><sub>-0.30</sub></td><td>2</td><td>-</td><td>-</td></tr><tr><td>27.75</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-6.85<sup>+0.15</sup><sub>-0.22</sub></td><td>5</td></tr><tr><td rowspan="12">1.5 &lt; z &le; 2.0</td><td>24.45</td><td>-4.61<sup>+0.08</sup><sub>-0.08</sub></td><td>306</td><td>-5.31<sup>+0.09</sup><sub>-0.09</sub></td><td>64</td><td>-</td><td>-</td><td>-6.11<sup>+0.13</sup><sub>-0.15</sub></td><td>13</td></tr><tr><td>24.75</td><td>-4.70<sup>+0.04</sup><sub>-0.04</sub></td><td>581</td><td>-5.52<sup>+0.06</sup><sub>-0.06</sub></td><td>85</td><td>-4.77<sup>+0.04</sup><sub>-0.04</sub></td><td>496</td><td>-5.83<sup>+0.07</sup><sub>-0.07</sub></td><td>45</td></tr><tr><td>25.05</td><td>-4.83<sup>+0.03</sup><sub>-0.03</sub></td><td>516</td><td>-5.79<sup>+0.05</sup><sub>-0.06</sub></td><td>57</td><td>-4.88<sup>+0.03</sup><sub>-0.03</sub></td><td>459</td><td>-5.65<sup>+0.05</sup><sub>-0.06</sub></td><td>79</td></tr><tr><td>25.35</td><td>-5.08<sup>+0.03</sup><sub>-0.03</sub></td><td>314</td><td>-5.93<sup>+0.05</sup><sub>-0.07</sub></td><td>44</td><td>-5.14<sup>+0.03</sup><sub>-0.03</sub></td><td>270</td><td>-5.83<sup>+0.05</sup><sub>-0.06</sub></td><td>55</td></tr><tr><td>25.65</td><td>-5.44<sup>+0.04</sup><sub>-0.04</sub></td><td>143</td><td>-6.37<sup>+0.09</sup><sub>-0.09</sub></td><td>17</td><td>-5.49<sup>+0.04</sup><sub>-0.04</sub></td><td>126</td><td>-6.08<sup>+0.07</sup><sub>-0.09</sub></td><td>33</td></tr><tr><td>25.95</td><td>-5.66<sup>+0.05</sup><sub>-0.05</sub></td><td>89</td><td>-6.35<sup>+0.09</sup><sub>-0.11</sub></td><td>18</td><td>-5.76<sup>+0.05</sup><sub>-0.06</sub></td><td>71</td><td>-6.38<sup>+0.09</sup><sub>-0.12</sub></td><td>17</td></tr><tr><td>26.25</td><td>-5.88<sup>+0.05</sup><sub>-0.06</sub></td><td>54</td><td>-6.53<sup>+0.11</sup><sub>-0.15</sub></td><td>12</td><td>-5.99<sup>+0.06</sup><sub>-0.07</sub></td><td>42</td><td>-6.53<sup>+0.12</sup><sub>-0.12</sub></td><td>12</td></tr><tr><td>26.55</td><td>-6.08<sup>+0.07</sup><sub>-0.08</sub></td><td>34</td><td>-6.53<sup>+0.10</sup><sub>-0.14</sub></td><td>12</td><td>-6.27<sup>+0.09</sup><sub>-0.11</sub></td><td>22</td><td>-6.66<sup>+0.12</sup><sub>-0.18</sub></td><td>9</td></tr><tr><td>26.85</td><td>-6.20<sup>+0.08</sup><sub>-0.09</sub></td><td>26</td><td>-6.66<sup>+0.14</sup><sub>-0.20</sub></td><td>9</td><td>-6.39<sup>+0.09</sup><sub>-0.11</sub></td><td>17</td><td>-6.84<sup>+0.18</sup><sub>-0.18</sub></td><td>6</td></tr><tr><td>27.15</td><td>-6.34<sup>+0.10</sup><sub>-0.10</sub></td><td>19</td><td>-6.92<sup>+0.18</sup><sub>-0.18</sub></td><td>5</td><td>-6.48<sup>+0.09</sup><sub>-0.11</sub></td><td>14</td><td>-6.67<sup>+0.12</sup><sub>-0.18</sub></td><td>9</td></tr><tr><td>27.45</td><td>-6.84<sup>+0.12</sup><sub>-0.18</sub></td><td>6</td><td>-</td><td>-</td><td>-6.84<sup>+0.12</sup><sub>-0.30</sub></td><td>6</td><td>-7.14<sup>+0.22</sup><sub>-0.48</sub></td><td>3</td></tr><tr><td>27.75</td><td>-6.67<sup>+0.12</sup><sub>-0.18</sub></td><td>9</td><td>-7.32<sup>+0.18</sup><sub>-0.30</sub></td><td>2</td><td>-6.78<sup>+0.12</sup><sub>-0.24</sub></td><td>7</td><td>-7.32<sup>+0.18</sup><sub>-0.30</sub></td><td>2</td></tr><tr><td>28.05</td><td>-7.14<sup>+0.22</sup><sub>-0.48</sub></td><td>3</td><td>-</td><td>-</td><td>-7.14<sup>+0.22</sup><sub>-0.48</sub></td><td>3</td><td>-</td><td>-</td></tr><tr><td rowspan="12">2.0 &lt; z &le; 2.5</td><td>24.75</td><td>-4.90<sup>+0.07</sup><sub>-0.07</sub></td><td>210</td><td>-5.84<sup>+0.11</sup><sub>-0.12</sub></td><td>23</td><td>-4.95<sup>+0.08</sup><sub>-0.07</sub></td><td>187</td><td>-6.19<sup>+0.13</sup><sub>-0.15</sub></td><td>13</td></tr><tr><td>25.05</td><td>-4.98<sup>+0.04</sup><sub>-0.04</sub></td><td>344</td><td>-6.23<sup>+0.08</sup><sub>-0.09</sub></td><td>19</td><td>-5.00<sup>+0.04</sup><sub>-0.04</sub></td><td>325</td><td>-5.93<sup>+0.07</sup><sub>-0.07</sub></td><td>39</td></tr><tr><td>25.35</td><td>-5.14<sup>+0.03</sup><sub>-0.03</sub></td><td>269</td><td>-6.67<sup>+0.13</sup><sub>-0.18</sub></td><td>8</td><td>-5.16<sup>+0.03</sup><sub>-0.03</sub></td><td>261</td><td>-5.85<sup>+0.06</sup><sub>-0.07</sub></td><td>53</td></tr><tr><td>25.65</td><td>-5.54<sup>+0.04</sup><sub>-0.04</sub></td><td>115</td><td>-6.90<sup>+0.13</sup><sub>-0.30</sub></td><td>5</td><td>-5.56<sup>+0.04</sup><sub>-0.04</sub></td><td>110</td><td>-6.03<sup>+0.06</sup><sub>-0.08</sub></td><td>38</td></tr><tr><td>25.95</td><td>-5.92<sup>+0.06</sup><sub>-0.07</sub></td><td>50</td><td>-7.02<sup>+0.18</sup><sub>-0.30</sub></td><td>4</td><td>-5.96<sup>+0.06</sup><sub>-0.06</sub></td><td>46</td><td>-6.39<sup>+0.09</sup><sub>-0.12</sub></td><td>17</td></tr><tr><td>26.25</td><td>-6.06<sup>+0.07</sup><sub>-0.08</sub></td><td>37</td><td>-7.15<sup>+0.15</sup><sub>-0.24</sub></td><td>3</td><td>-6.10<sup>+0.07</sup><sub>-0.09</sub></td><td>34</td><td>-6.55<sup>+0.10</sup><sub>-0.12</sub></td><td>12</td></tr><tr><td>26.55</td><td>-6.55<sup>+0.10</sup><sub>-0.12</sub></td><td>12</td><td>-7.33<sup>+0.18</sup><sub>-0.30</sub></td><td>2</td><td>-6.63<sup>+0.11</sup><sub>-0.15</sub></td><td>10</td><td>-6.73<sup>+0.14</sup><sub>-0.20</sub></td><td>8</td></tr><tr><td>26.85</td><td>-6.46<sup>+0.10</sup><sub>-0.13</sub></td><td>15</td><td>-</td><td>-</td><td>-6.46<sup>+0.11</sup><sub>-0.10</sub></td><td>15</td><td>-6.73<sup>+0.10</sup><sub>-0.20</sub></td><td>8</td></tr><tr><td>27.15</td><td>-6.74<sup>+0.14</sup><sub>-0.20</sub></td><td>8</td><td>-</td><td>-</td><td>-6.80<sup>+0.15</sup><sub>-0.24</sub></td><td>7</td><td>-7.17<sup>+0.22</sup><sub>-0.48</sub></td><td>3</td></tr><tr><td>27.45</td><td>-7.17<sup>+0.22</sup><sub>-0.48</sub></td><td>3</td><td>-</td><td>-</td><td>-7.17<sup>+0.22</sup><sub>-0.48</sub></td><td>3</td><td>-7.17<sup>+0.22</sup><sub>-0.48</sub></td><td>3</td></tr><tr><td>27.75</td><td>-7.17<sup>+0.22</sup><sub>-0.48</sub></td><td>3</td><td>-</td><td>-</td><td>-7.17<sup>+0.22</sup><sub>-0.48</sub></td><td>3</td><td>-</td><td>-</td></tr><tr><td>28.05</td><td>-7.17<sup>+0.22</sup><sub>-0.48</sub></td><td>3</td><td>-</td><td>-</td><td>-7.17<sup>+0.22</sup><sub>-0.48</sub></td><td>3</td><td>-</td><td>-</td></tr><tr><td>28.35</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-7.34<sup>+0.30</sup><sub>-0.30</sub></td><td>2</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Cosmic evolution (across $0.5 < z \le 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 $\alpha = -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 \le 1.0$ bin, the LERG LF from Best et al. (2014, green triangles), scaled to 150 MHz using a spectral index $\alpha = -0.7$, shows good agreement with our results. The LFs from Williams et al. (2018), reaching $z \sim 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.</figcaption></figure>

<table><caption>**Table 4.** Results of modelling the evolution of the LERG LFs using a pure luminosity evolution (PLE), pure density evolution (PDE), and the luminosity and density evolution (LDE) models, relative to the local relation from Mauch & Sadler (2007), as detailed in Section 4.1. The best-fitting LFs for each model are also shown in Fig. 6.</caption><thead><tr><th rowspan="2">z</th><th colspan="3">PLE</th><th colspan="3">PDE</th><th colspan="3">LDE</th></tr><tr><th>log<sub>10</sub>L<sub>*</sub>(z)</th><th>α<sub>L</sub></th><th>χ<sub>ν</sub><sup>2</sup></th><th>log<sub>10</sub>ρ<sup>*</sup>(z)</th><th>α<sub>D</sub></th><th>χ<sub>ν</sub><sup>2</sup></th><th>log<sub>10</sub>ρ<sup>*</sup>(z)</th><th>α<sub>D</sub></th><th>α<sub>L</sub></th><th>χ<sub>ν</sub><sup>2</sup></th></tr></thead><tbody><tr><td>0.5 &lt; z &le; 1.0</td><td>25.55<sup>+0.02</sup><sub>-0.02</sub></td><td>1.14<sup>+0.07</sup><sub>-0.07</sub></td><td>6.62</td><td>-4.92<sup>+0.01</sup><sub>-0.01</sub></td><td>0.76<sup>+0.05</sup><sub>-0.05</sub></td><td>9.47</td><td>-5.39<sup>+0.07</sup><sub>-0.07</sub></td><td>25.98<sup>+0.12</sup><sub>-0.11</sub></td><td>-1.19<sup>+0.28</sup><sub>-0.30</sub></td><td>2.93<sup>+0.45</sup><sub>-0.49</sub></td><td>5.59</td></tr><tr><td>1.0 &lt; z &le; 1.5</td><td>25.50<sup>+0.02</sup><sub>-0.02</sub></td><td>0.66<sup>+0.05</sup><sub>-0.05</sub></td><td>9.19</td><td>-4.93<sup>+0.01</sup><sub>-0.01</sub></td><td>0.48<sup>+0.04</sup><sub>-0.04</sub></td><td>12.44</td><td>-5.63<sup>+0.08</sup><sub>-0.09</sub></td><td>26.24<sup>+0.14</sup><sub>-0.13</sub></td><td>-1.50<sup>+0.23</sup><sub>-0.24</sub></td><td>2.76<sup>+0.36</sup><sub>-0.39</sub></td><td>5.56</td></tr><tr><td>1.5 &lt; z &le; 2.0</td><td>25.79<sup>+0.02</sup><sub>-0.02</sub></td><td>1.18<sup>+0.04</sup><sub>-0.04</sub></td><td>8.28</td><td>-4.68<sup>+0.01</sup><sub>-0.01</sub></td><td>0.96<sup>+0.03</sup><sub>-0.03</sub></td><td>15.33</td><td>-5.52<sup>+0.08</sup><sub>-0.08</sub></td><td>26.35<sup>+0.12</sup><sub>-0.11</sub></td><td>-0.95<sup>+0.17</sup><sub>-0.18</sub></td><td>2.47<sup>+0.26</sup><sub>-0.27</sub></td><td>6.28</td></tr><tr><td>2.0 &lt; z &le; 2.5</td><td>25.58<sup>+0.02</sup><sub>-0.02</sub></td><td>0.60<sup>+0.04</sup><sub>-0.04</sub></td><td>3.14</td><td>-4.83<sup>+0.01</sup><sub>-0.01</sub></td><td>0.54<sup>+0.03</sup><sub>-0.03</sub></td><td>3.70</td><td>-5.13<sup>+0.12</sup><sub>-0.11</sub></td><td>25.61<sup>+0.14</sup><sub>-0.14</sub></td><td>-0.05<sup>+0.23</sup><sub>-0.22</sub></td><td>0.67<sup>+0.27</sup><sub>-0.27</sub></td><td>3.45</td></tr></tbody></table>

redshift evolution of the characteristic luminosity, $L^*(z)$, is given as $L^*(z) = L^*(0) \times (1+z)^{\alpha_L}$. Here, $\alpha_D$ and $\alpha_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^*(0) = \frac{1}{0.4}10^{-5.5} \log L^{-1} \text{Mpc}^{-3}$ (converted to per log luminosity) and the characteristic luminosity $L^*(0) = 10^{25.27} \text{W Hz}^{-1}$ (scaled to 150 MHz assuming $\alpha = -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$, are also listed in Table 4.

For each redshift bin, we also considered a pure luminosity evolution (PLE) model (i.e. we set $\alpha_D = 0$) and a pure density evolution (PDE) model (i.e. we set $\alpha_L = 0$). The resulting best-fitting $\alpha_D$ and $\alpha_L$ parameters (as appropriate for each model) for each redshift bin are listed in Table 4, along with the $\chi^2$ 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_v^2$ 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 \sim 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 $\rho_*(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 \sim 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 $\sim 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 \le 0.5$, $0.5 < z \le 0.7$, $0.7 < z \le 1.0$, $1.0 < z \le 1.2$, $1.2 < z \le 1.5$), shown by different coloured lines, for three equally spaced radio luminosity limits (shown in different panels) from $L_{150\,MHz} \ge 10^{24} \, \mathrm{W Hz}^{-1}$ to $L_{150\,MHz} \ge 10^{25} \, \mathrm{W Hz}^{-1}$. The lowest radio luminosity limit chosen corresponds to a $5\sigma$ detection (based on noise in the deepest regions) at $z \sim 1.5$, and is broadly comparable to the $1.4 \, \mathrm{GHz}$ luminosity limit ($L_{1.4\,GHz} = 10^{23} \, \mathrm{W Hz}^{-1}$, assuming $\alpha = -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\sigma$ based on the local radio RMS level. The error bars show Poisson uncertainties (following Gehrels (1986) for stellar mass bins with $N < 10$). In total, this results in 2722, 1776, and 835 LERGs (within $0.3 < z \le 1.5$) with $L_{150\,MHz} \ge 10^{24} \, \mathrm{W Hz}^{-1}$, $L_{150\,MHz} \ge 10^{24.5} \, \mathrm{W Hz}^{-1}$, and $L_{150\,MHz} \ge 10^{25} \, \mathrm{W Hz}^{-1}$, respectively.

For the high radio-power LERGs, corresponding to a radio luminosity limit of $10^{25} \, \mathrm{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 $\sim 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

It is well known that star-forming galaxies occupy a well-defined sequence in the SFR-M$_*$ plane, known as the ‘main-sequence’ of star-formation (e.g. Whitaker et al. 2012; Speagle et al. 2014; Schreiber et al. 2015; Tasca et al. 2015), with quiescent galaxies lying below this sequence. Therefore, the ratio of the star-formation rate and the stellar mass of a galaxy, known as the specific star-formation rate (sSFR) can be used to identify quiescent galaxies (e.g. Fontana et al. 2009; Pacifici et al. 2016; Merlin et al. 2018). In this study, using the consensus SFRs and stellar masses derived from SED fitting (see Section 2.4), we select sources as quiescent galaxies if they satisfy the condition

$$ \mathrm{sSFR} < f_{\mathrm{sSFR}}/t_{\mathrm{H}}(z), \qquad (8) $$

where $t_{\mathrm{H}}(z)$ is the age of the Universe at redshift $z$ of the source, and $f_{\mathrm{sSFR}} = 1/5$ defines the threshold in sSFR. Galaxy samples selected by this criterion have been found to show good agreement with traditional UVJ rest-frame colour selected samples (e.g. Pacifici et al. 2016; Carnall et al. 2018, 2020).
<figure><img src="image_9.png" /><figcaption>**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<sub>150 MHz</sub> ≥ 10<sup>24</sup> W Hz<sup>-1</sup>, L<sub>150 MHz</sub> ≥ 10<sup>24.5</sup> W Hz<sup>-1</sup>, and L<sub>150 MHz</sub> ≥ 10<sup>25</sup> W Hz<sup>-1</sup>. The error bars represent Poisson uncertainties, following Gehrels (1986) for N < 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.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** The location of LERGs in the ELAIS-N1 field on the SFR-stellar mass (M<sub>*</sub>) 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<sub>*</sub> 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<sub>sSFR</sub> = 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<sub>sSFR</sub> = 1/3 and f<sub>sSFR</sub> = 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.</figcaption></figure>

Fig. 8 shows the SFR-M<sub>*</sub> 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<sub>H</sub>(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.

The green dotted lines above and below this dividing line correspond to $f_{\text{sSFR}} = 1/3$ and $f_{\text{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_{3.6\mu m} > 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 $\chi^2$ SED fits, and within the chosen redshift range (i.e. $0.3 < z \le 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_{\text{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_{150\,\text{MHz}} \ge 10^{24}\,\text{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_{\text{LERG}} \approx 0.01(M_*/10^{11}\,\text{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.

To quantitatively investigate the trends in the evolution of the LERG stellar mass fractions seen in Fig. 9, we parametrized the LERG stellar mass fractions as a power law of the form

$$f_{\text{LERG}}(M_*) = c \left( \frac{M_*}{10^{11} M_\odot} \right)^\beta \quad (9)$$

such that $f_{\text{LERG}}(M_*)$ is the LERG fraction at mass $M_*$, c is the normalization at $M_* = 10^{11} M_\odot$, and $\beta$ is the power-law slope. We then fitted this power-law form to the LERG stellar mass fractions shown in Fig. 9 at each redshift individually, for both the star-forming and quiescent LERG populations.

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, $\beta$, and right-hand panel shows the redshift evolution of the normalization at $M_* = 10^{11} M_\odot$ for both populations. The error bars show the $1\sigma$ 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_{\text{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_{\text{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 ($\beta = 2.5$) found in studies of jet-mode AGN population in the local Universe (e.g. Best et al. 2005a; Janssen et al. 2012).

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_{\text{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\sigma$ 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} M_\odot$ in Fig. 10 (right-hand panel), we find that for the quiescent LERGs, the normalization stays roughly constant out to $z \sim 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 $\sim 4$ by $z \sim 1$, showing hints of higher prevalence at $M = 10^{11} M_\odot$ compared to the quiescent hosts at these redshifts.

We note that our use of a radio-excess $>3\sigma$ 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_{150\,\text{MHz}} \ge 10^{25}\,\text{W Hz}^{-1}$ (with broader redshift bins $0.5 < z \le 1$ and $1 < z \le 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
<figure><img src="image_11.png" /><figcaption>**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_{150\,\text{MHz}} \ge 10^{24} \,\text{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 \le 1.5$, for the given radio luminosity limit. The error bars represent Poisson uncertainties, following Gehrels (1986) for $N < 10$. The black dashed line shows the $f_{\text{LERG}} \approx 0.01(M_*/10^{11} 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 \sim 1.5$. The star-forming LERGs show a much shallower dependence on stellar mass, with an increase in the fraction at higher redshifts.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**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_* = 10^{11} M_{\odot}$. The solid line (with crosses) corresponds to the results obtained from using the $f_{\text{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_{\text{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.</figcaption></figure>

sample, leading to poorer statistics and hence we were unable to converge on a fit for the lowest redshift bin. For the $1 < z \le 1.5$ bin, we find a power law slope $\beta = 1.37 \pm 0.57$, which is $\sim 2\sigma$ away from the power-law slope of $\beta = 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 \sim 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 \sim 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_* \lesssim 10^{11} 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 (~$10^8 - 10^{10} 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\sigma$ 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 \sim 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).

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\sigma$ 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 \le 1.5$
<figure><img src="image_13.png" /><figcaption>**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 \sim 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 \le 1$ bin, with the shaded region forming the $1\sigma$ uncertainty on the best fit. At higher redshifts, the dashed line and hatched regions represent the best fit and $1\sigma$ 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\sigma$ region at each redshift.</figcaption></figure>

<table><caption>**Table 5.** Results of modelling the evolution of the quiescent LERGs. We fitted a broken power law to the $0.5 < z \le 1$ LF. Using this broken power-law fit (i.e. constant slopes), the LF at higher redshifts is modelled with a luminosity and density evolution. To reduce the degeneracy between parameters, the evolution of $\rho_*(z)$ was fixed based on the evolution of the expected host galaxies of the quiescent LERGs (see Fig. 13), with $L_*(z)$ fitted as a free parameter. The best-fitting models are also shown in Fig. 11.</caption><thead><tr><th>z</th><th>log<sub>10</sub>ρ*(z)</th><th>log<sub>10</sub>L<sub>*</sub>(z)</th><th>β</th><th>γ</th><th>χ<sup>2</sup><sub>ν</sub></th></tr></thead><tbody><tr><td>0.5 &lt; z &le; 1.0</td><td>-6.37<sup>+0.17</sup><sub>-0.23</sub></td><td>27.03<sup>+0.39</sup><sub>-0.22</sub></td><td>-2.88<sup>+1.11</sup><sub>-1.76</sub></td><td>-0.55<sup>+0.03</sup><sub>-0.03</sub></td><td>3.57</td></tr><tr><td>1.0 &lt; z &le; 1.5</td><td>-6.92</td><td>27.33<sup>+0.04</sup><sub>-0.04</sub></td><td>...</td><td>...</td><td>3.39</td></tr><tr><td>1.5 &lt; z &le; 2.0</td><td>-7.28</td><td>27.91<sup>+0.06</sup><sub>-0.06</sub></td><td>...</td><td>...</td><td>2.28</td></tr><tr><td>2.0 &lt; z &le; 2.5</td><td>-7.61</td><td>27.93<sup>+0.12</sup><sub>-0.12</sub></td><td>...</td><td>...</td><td>5.89</td></tr></tbody></table>

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 \le 1$ panel of Fig. 11 and the shaded region representing the $1\sigma$ 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 \le 1$, and adopted an LDE model by fitting directly for $\rho_*$ and $L_*$ at each redshift.⁴ The resulting fitted model and $1\sigma$ 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 $\rho_*$ and $L_*$, as can be visualized in Fig. 12 which shows

³Interestingly, in Fig. 11, it can be seen that all redshifts, the LF could potentially be fit better with a shallower faint-end slope, with then an upturn in the space densities of the quiescent LERGs at the faintest luminosities. At the lowest redshift bin ($0.5 < z \le 1$), the upturn starts to appear at $L_{150\,MHz} \approx 10^{24.5} \, \text{W Hz}^{-1}$; this is typically below the limit probed by previous surveys. If this upturn is real, then it might suggest that a double Schechter function may be more appropriate than a broken power law for modelling the luminosity functions of the quiescent LERGs in particular. However, the present data set does not have enough data points at the faint end beyond the upturn to provide constraints on such a model; deeper radio data, for example from subsequent data releases of the LoTSS Deep Fields, will allow us to better sample the faint end of the quiescent LERGs.

⁴We also examined the PLE and PDE models, but neither is able to match the data well: the PLE model cannot explain the strong evolution in space densities seen at low luminosities, whereas the PDE model fails to match the evolution at high luminosities.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** 2D contour plot showing the posterior distribution of $\rho_*(z)$ and $L_*(z)$ for the quiescent LERGs. For $0.5 < z \le 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 $\rho_*(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\sigma$ uncertainties on the evolved $\rho_*(z)$.</figcaption></figure>

a 2D contour plot of the posterior distribution of the LDE model fit in $L_*(z)$ and $\rho_*(z)$ parameter space for different redshift bins. For the $0.5 < z \le 1$ redshift bin, the contours shown are from the broken power law fit described above.

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}$ M$_\odot$ across $0.3 < z \le 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 \sim 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 to constrain the LDE model by scaling the $\rho_*(z = 0.75)$ from the broken power-law fit to predict $\rho_*(z)$ at higher redshifts, as shown by the dashed lines in the 2D contour plot in Fig. 12 for each

<figure><img src="image_15.png" /><figcaption>**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_{AGN}(M_*) = 0.018(\mathcal{M}_*/10^{11}\ \mathrm{M}_\odot)$; see the text) and integrating in stellar mass (down to $M_* = 10^9\ \mathrm{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 compiled by Best et al. (2014, taken from their fig. 6) and scaled to match the normalization in $f_{AGN}(M_*)$ 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.</figcaption></figure>

redshift bin, with the shaded regions representing the uncertainties in $\rho_*(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 $\rho_*(z)$ at higher redshifts, and fitted directly for $L_*(z)$ only for each of the three $1.0 < z \le 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\sigma$ hatched region at each redshift.

These results agree well with studies of the total LERG population out to $z \lesssim 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 \sim 2.5$. The LoTSS deep fields represents the deepest low-frequency radio continuum survey to date, covering 25 deg$^2$ 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 \le 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 \le 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 \sim 2.5$. We find relatively mild evolution within the redshift bins ($0.5 < z \le 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:

(i) The quiescent LERG LFs dominate the space densities of LERGs across essentially all luminosities at $z \le 0.75$.

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

(iii) 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$ Gyr.

(iv) The quiescent LERGs show a strong negative evolution in their LFs beyond $z \sim 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:

(i) The incidence of LERG activity in star-forming galaxies shows a much weaker stellar mass dependence ($\beta \approx 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.

(ii) 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 \le 1.5$, reaching a higher prevalence than that of the quiescent LERGs at the highest redshifts.

(iii) 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 \gtrsim 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 \sim 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/roh-itk-10/AGN_LF_Kondapally22. Other results presented in the paper are available upon reasonable request to the corresponding author.

## REFERENCES

Aihara H. et al., 2018, PASJ, 70, S8

Aird J., Coil A. L., Georgakakis A., 2019, MNRAS, 484, 4360

Allen S. W., Dunn R. J. H., Fabian A. C., Taylor G. B., Reynolds C. S., 2006, MNRAS, 372, 21

Ashby M. L. N. et al., 2009, ApJ, 701, 428

Astropy Collaboration, 2013, A&A, 558, A33

Astropy Collaboration, 2018, AJ, 156, 123

Baldry I. K. et al., 2012, MNRAS, 421, 621

Becker R. H., White R. L., Helfand D. J., 1995, ApJ, 450, 559

Best P. N., 2004, MNRAS, 351, 70

Best P. N., Heckman T. M., 2012, MNRAS, 421, 1569

Best P. N., Kauffmann G., Heckman T. M., Brinchmann J., Charlot S., Ivezić Ž., White S. D. M., 2005a, MNRAS, 362, 25

Best P. N., Kauffmann G., Heckman T. M., Ivezić Ž., 2005b, MNRAS, 362, 9

Best P. N., Kaiser C. R., Heckman T. M., Kauffmann G., 2006, MNRAS, 368, L67

Best P. N., Ker L. M., Simpson C., Rigby E. E., Sabater J., 2014, MNRAS, 445, 955

Bonato M. et al., 2021, A&A, 656, A48

Boquien M., Burgarella D., Roehlly Y., Buat V., Ciesla L., Corre D., Inoue A. K., Salas H., 2019, A&A, 622, A103

Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645

Brown M. J. I., Dey A., Jannuzi B. T., Brand K., Benson A. J., Brodwin M., Croton D. J., Eisenhardt P. R., 2007, ApJ, 654, 858

Brown M. J. I. et al., 2008, ApJ, 682, 937

Burgarella D., Buat V., Iglesias-Páramo J., 2005, MNRAS, 360, 1413

Butler A. et al., 2018, A&A, 620, A3

Butler A., Huynh M., Kapińska A., Delvecchio I., Smolčić V., Chiappetti L., Koulouridis E., Pierre M., 2019, A&A, 625, A111

Calistro Rivera G., Lusso E., Hennawi J. F., Hogg D. W., 2016, ApJ, 833, 98

Calistro Rivera G. et al., 2017, MNRAS, 469, 3468

Calistro Rivera G. et al., 2018, ApJ, 863, 56

Carnall A. C., McLure R. J., Dunlop J. S., Davé R., 2018, MNRAS, 480, 4379

Carnall A. C. et al., 2020, MNRAS, 496, 695

Carraro R. et al., 2020, A&A, 642, A65

Cattaneo A., Dekel A., Devriendt J., Guiderdoni B., Blaizot J., 2006, MNRAS, 370, 1651

Cattaneo A. et al., 2009, Nature, 460, 213

Ching J. H. Y. et al., 2017a, MNRAS, 464, 1306

Ching J. H. Y. et al., 2017b, MNRAS, 469, 4584

Cochrane R. K. et al., 2021, MNRAS, 503, 2622

Condon J. J., 1989, ApJ, 338, 13

Condon J. J., Cotton W. D., Greisen E. W., Yin Q. F., Perley R. A., Taylor G. B., Broderick J. J., 1998, AJ, 115, 1693

Croston J. H. et al., 2019, A&A, 622, A10

Croton D. J. et al., 2006, MNRAS, 365, 11

da Cunha E., Charlot S., Elbaz D., 2008, MNRAS, 388, 1595

Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827

de Ruiter H. R., Willis A. G., Arp H. C., 1977, A&AS, 28, 211

Delvecchio I. et al., 2017, A&A, 602, A3

Domínguez Sánchez H. et al., 2011, MNRAS, 417, 900

Duncan K. et al., 2014, MNRAS, 444, 2960

Duncan K. J. et al., 2018a, MNRAS, 473, 2655

Duncan K. J., Jarvis M. J., Brown M. J. I., Röttgering H. J. A., 2018b, MNRAS, 477, 5177

Duncan K. et al., 2019, ApJ, 876, 110

Duncan K. J. et al., 2021, A&A, 648, A4

Dunlop J. S., Peacock J. A., 1990, MNRAS, 247, 19

Eisenhardt P. R. et al., 2004, ApJS, 154, 48

Fabian A. C., 2012, ARA&A, 50, 455

Ferrarese L., Merritt D., 2000, ApJ, 539, L9

Fontana A. et al., 2009, A&A, 501, 15

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

Fritz J., Franceschini A., Hatziminaoglou E., 2006, MNRAS, 366, 767

Gao F. et al., 2020, A&A, 637, A94

Gaspari M., Ruszkowski M., Oh S. P., 2013, MNRAS, 432, 3401

Gaspari M., Brighenti F., Temi P., 2015, A&A, 579, A62

Gaspari M., Temi P., Brighenti F., 2017, MNRAS, 466, 677

Gehrels N., 1986, ApJ, 303, 336

Gendre M. A., Best P. N., Wall J. V., Ker L. M., 2013, MNRAS, 430, 3086

Gordon Y. A. et al., 2019, ApJ, 878, 88

Griffin M. J. et al., 2010, A&A, 518, L3

Gürkan G. et al., 2018, MNRAS, 475, 3010

Hardcastle M., 2018, Nat. Astron., 2, 273

Hardcastle M. J., Croston J. H., 2020, New Astron. Rev., 88, 101539

Hardcastle M. J., Evans D. A., Croston J. H., 2007, MNRAS, 376, 1849

Hardcastle M. J. et al., 2019, A&A, 622, A12

Heckman T. M., Best P. N., 2014, ARA&A, 52, 589

Hildebrandt H. et al., 2016, MNRAS, 463, 635

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Ilbert O. et al., 2013, A&A, 556, A55

Jannuzi B. T., Dey A., 1999, in Weymann R., Storrie-Lombardi L., Sawicki M., Brunner R., eds, ASP Conf. Ser. Vol. 191, Photometric Redshifts and the Detection of High Redshift Galaxies. Astron. Soc. Pac., San Francisco, p. 111

Janssen R. M. J., Röttgering H. J. A., Best P. N., Brinchmann J., 2012, A&A, 541, A62

Kaiser N. et al., 2010, in Stepp L. M., Gilmozzi R., Hall H. J., eds, Proc. SPIE Conf. Ser. Vol. 7733, Ground-based and Airborne Telescopes III. SPIE, Bellingham, p. 77330E

Kaviraj S. et al., 2017, MNRAS, 467, 4739

Kondapally R. et al., 2021, A&A, 648, A3

Kormendy J., Ho L. C., 2013, ARA&A, 51, 511

Lawrence A. et al., 2007, MNRAS, 379, 1599

Lonsdale C. J. et al., 2003, PASP, 115, 897

Maccagni F. M., Morganti R., Oosterloo T. A., Oonk J. B. R., Emonts B. H. C., 2018, A&A, 614, A42

Mahony E. K. et al., 2016, MNRAS, 463, 2997

Mandal S. et al., 2021, A&A, 648, A5

Mauch T., Sadler E. M., 2007, MNRAS, 375, 931

Mauduit J.-C. et al., 2012, PASP, 124, 714

McCheyne I. et al., 2022, A&A, in press

McLeod D. J., McLure R. J., Dunlop J. S., Cullen F., Carnall A. C., Duncan K., 2021, MNRAS, 503, 4413

McNamara B. R., Nulsen P. E. J., 2007, ARA&A, 45, 117
Merlin E. et al., 2018, MNRAS, 473, 2098

Merritt D., Ferrarese L., 2001, MNRAS, 320, L30

Mingo B., Hardcastle M. J., Croston J. H., Dicken D., Evans D. A., Morganti R., Tadhunter C., 2014, MNRAS, 440, 269

Mingo B. et al., 2022, MNRAS, 511, 3250

Mohan N., Rafferty D., 2015, PyBDSF: Python Blob Detection and Source Finder, Astrophysics Source Code Library. record ascl:1502.007

Morrissey P. et al., 2007, ApJS, 173, 682

Moustakas J. et al., 2013, ApJ, 767, 50

Murphy E. J., Momjian E., Condon J. J., Chary R.-R., Dickinson M., Inami H., Taylor A. R., Weiner B. J., 2017, ApJ, 839, 35

Muzzin A. et al., 2009, ApJ, 698, 1934

Muzzin A. et al., 2013, ApJ, 777, 18

Narayan R., Yi I., 1994, ApJ, 428, L13

Narayan R., Yi I., 1995, ApJ, 452, 710

Noll S., Burgarella D., Giovannoli E., Buat V., Marcillac D., Muñoz-Mateos J. C., 2009, A&A, 507, 1793

North E. V. et al., 2019, MNRAS, 490, 319

Ocaña Flaquer B., Leon S., Combes F., Lim J., 2010, A&A, 518, A9

Oke J. B., Gunn J. E., 1983, ApJ, 266, 713

Okuda T., Kohno K., Iguchi S., Nakanishi K., 2005, ApJ, 620, 673

Oliver S. J. et al., 2012, MNRAS, 424, 1614

Pacifici C. et al., 2016, ApJ, 832, 79

Padovani P., 2016, A&AR, 24, 13

Peterson J. R., Fabian A. C., 2006, Phys. Rep., 427, 1

Pilbratt G. L. et al., 2010, A&A, 518, L1

Poglitsch A. et al., 2010, A&A, 518, L2

Pracy M. B. et al., 2016, MNRAS, 460, 2

Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49

Retana-Montenegro E., Röttgering H. J. A., Shimwell T. W., van Weeren R. J., Prandoni I., Brunetti G., Best P. N., Brüggen M., 2018, A&A, 620, A74

Richards G. T. et al., 2006, AJ, 131, 2766

Rieke G. H. et al., 2004, ApJS, 154, 25

Ruffa I. et al., 2019a, MNRAS, 484, 4239

Ruffa I. et al., 2019b, MNRAS, 489, 3739

Sabater J., Best P. N., Argudo-Fernández M., 2013, MNRAS, 430, 638

Sabater J. et al., 2019, A&A, 622, A17

Sabater J. et al., 2021, A&A, 648, A2

Schawinski K. et al., 2014, MNRAS, 440, 889

Schmidt M., 1968, ApJ, 151, 393

Schreiber C. et al., 2015, A&A, 575, A74

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

Sharma P., McCourt M., Quataert E., Parrish I. J., 2012, MNRAS, 420, 3174

Shimwell T. W. et al., 2022, A&A, 659, A1

Shirley R. et al., 2021, MNRAS, 507, 129

Silva L., Maiolino R., Granato G. L., 2004, MNRAS, 355, 973

Simpson C., Westoby P., Arumugam V., Ivison R., Hartley W., Almaini O., 2013, MNRAS, 433, 2647

Smith D. J. B. et al., 2016, in Reylé C., Richard J., Cambrésy L., Deleuil M., Pécontal E., Tresse L., Vauglin I., eds, SF2A-2016: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics. French Soc. Astron. Astrophys., Lyon, p. 271

Smith D. J. B. et al., 2021, A&A, 648, A6

Smolčić V., 2009, ApJ, 699, L43

Smolčić V. et al., 2009, ApJ, 696, 24

Smolčić V. et al., 2017a, A&A, 602, A1

Smolčić V. et al., 2017b, A&A, 602, A6

Somerville R. S., Davé R., 2015, ARA&A, 53, 51

Speagle J. S., Steinhardt C. L., Capak P. L., Silverman J. D., 2014, ApJS, 214, 15

Stalevski M., Fritz J., Baes M., Nakos T., Popović L. Č., 2012, MNRAS, 420, 2756

Stalevski M., Ricci C., Ueda Y., Lira P., Fritz J., Baes M., 2016, MNRAS, 458, 2288

Sutherland W., Saunders W., 1992, MNRAS, 259, 413

Sweijen F. et al., 2022, Nat. Astron., 6, 350

Tadhunter C., 2016, A&AR, 24, 10

Tasca L. A. M. et al., 2015, A&A, 581, A54

Tasse C., Best P. N., Röttgering H., Le Borgne D., 2008, A&A, 490, 893

Tasse C. et al., 2021, A&A, 648, A1

Tremblay G. R. et al., 2016, Nature, 534, 218

van Haarlem M. P. et al., 2013, A&A, 556, A2

Vito F., Gilli R., Vignali C., Comastri A., Brusa M., Cappelluti N., Iwasawa K., 2014, MNRAS, 445, 3557

Werner M. W. et al., 2004, ApJS, 154, 1

Whitaker K. E., van Dokkum P. G., Brammer G., Franx M., 2012, ApJ, 754, L29

Whittam I. H., Prescott M., McAlpine K., Jarvis M. J., Heywood I., 2018, MNRAS, 480, 358

Williams W. L. et al., 2016, MNRAS, 460, 2385

Williams W. L. et al., 2018, MNRAS, 475, 3429

Williams W. L. et al., 2019, A&A, 622, A2

Wilman R. J. et al., 2008, MNRAS, 388, 1335

Wilson G. et al., 2009, ApJ, 698, 1943

Windhorst R., Mathis D., Neuschaefer L., 1990, in Kron R. G., ed., ASP Conf. Ser. Vol. 10, Evolution of the Universe of Galaxies. Astron. Soc. Pac., San Francisco, p. 389

Windhorst R. A., Fomalont E. B., Partridge R. B., Lowenthal J. D., 1993, ApJ, 405, 498

Yang G. et al., 2020, MNRAS, 491, 740

York D. G. et al., 2000, AJ, 120, 1579

Yuan F., Narayan R., 2014, ARA&A, 52, 529

## APPENDIX A: COMPARISON WITH WILLIAMS ET AL. (2018)

In Section 4, we compared the cosmic evolution of LERGs in LoTSS-Deep with other literature results and found that the Williams et al. LFs are systematically offset to lower space densities; here, we investigate the potential causes for this apparent discrepancy.

As detailed in Section 4, Williams et al. constructed their LFs using a sample of 243 LERGs which were selected using a radio-excess selection based on the FIRC and SED fitting from AGNFITTER. In principle, this selection method is similar to the one used by Best et al. (in preparation), however, in practice, the Williams et al. (2018) radio-excess criterion is more conservative, increasingly so with increasing redshift, and their criteria for separating HERGs and LERGs is also different (see Section 2.5). Therefore we reconstructed the LERG LFs with the LoTSS-Deep data using the classification scheme adopted by Williams et al. (2018). Specifically, we identified radio-excess sources making use of the FIRC of Calistro Rivera et al. (2017) using properties derived from AGNFITTER only instead of the consensus values, and we classified these sources as HERGs or LERGs using the Williams et al. (2018) $f_{AGN}$ parameter.

Fig. A1 shows comparison with Williams et al. (2018) when applying their source classification scheme to our LoTSS Deep data set (cyan squares; hereafter LoTSS-Deep W18). We find that the derived LFs now show better agreement at all redshift with the LFs derived in this work (within $2\sigma$) based on Best et al. (in preparation) classifications than with the Williams et al. (2018) results, especially at high radio luminosities and higher redshifts. This appears to be largely because we used an updated version of AGNFITTER compared to the one used in Williams et al. with improved models being used in the fitting process resulting in better fits in general. This has a significant effect on the classification of sources as HERGs or LERGs.

The reconstructed LoTSS-Deep W18 LFs do not agree well with ours at faint radio luminosities, at all redshifts. We find that this is due to a significant number of sources classified as radio-excess AGN
<figure><img src="image_16.png" /><figcaption>**Figure A1.** Comparison of the cosmic evolution of LERG LF with Williams et al. (2018) using different source classification schemes. We show the LoTSS-Deep LERG LF resulting from the Best et al. (in preparation) classification scheme (pink circles), and the resultant LoTSS-Deep LF when the source classification method of Williams et al. (2018) is applied to our sample (cyan squares; see also Section 4). The literature results from Williams et al. (2018) (yellow squares) and Best et al. (2014) (green triangles) are also shown.</figcaption></figure>

(and subsequently as LERGs) by our criteria, but as SFGs under the Williams et al. (2018) criteria. One reason for this discrepancy is that the two studies use different estimates of SFRs: consensus SFRs versus AGNFITTER based IR luminosities (and hence SFRs) used by Williams et al. (2018). Best et al. (in preparation) show that the SFRs reported by AGNFITTER are found to be systematically offset to higher values than those reported by the other three SED fitting routines used. This may be caused by a combination of a lack of energy balance in AGNFITTER unlike other routines, and differences between the various SED fitting routines used in this study. Another reason for the discrepancy is that the radio-excess selection used by Williams et al. (2018) is more conservative as their redshift-dependent FIRC correlation has a $\sigma = 0.527$, and so their $2\sigma$ cut corresponds to a larger radio-excess than that of Best et al. (in preparation). Furthermore, the use of an evolving FIRC by Williams et al. (2018) means that their selection becomes relatively more conservative with increasing redshift. The combination of these two factors can help explain the observed differences with Williams et al. (2018); a detailed SED fitting comparison to investigate the cause of this discrepancy is outside the scope of this paper.

## APPENDIX B: EFFECTS OF VARYING THE SOURCE CLASSIFICATION CRITERIA

In this section, we explore how deviations from our adopted criteria for source classification affect our main results on the evolution of the LERG LFs presented in Sections 4–6.

As discussed in Sections 3.3 and 3.4, the local and evolving LFs for the total radio-AGN population derived in LoTSS-Deep agree well with literature results (as do those of the SFGs by Cochrane et al., in preparation), giving us confidence in our separation of star-forming galaxies from radio-excess AGN across redshift. Furthermore, we have investigated increasing and decreasing the radio-excess criteria by 0.1 dex to find that both the star-forming and quiescent LERG LFs remain largely unchanged, with typical changes in space densities of <0.15 dex. In this section, we only investigate the modification of the selection of the so-called ‘SED AGN’, which will result in a change between sources classified as HERGs and LERGs. We note that it can be challenging to identify (optical) AGN of moderate to low luminosities using photometric data; this may be particularly important in star-forming galaxies as both AGN activity and star-formation contribute to the MIR emission, making it difficult to disentangle contributions associated with the different emission processes. If such sources also show radio-loud emission and the SED AGN component is missed, then based on our source classification criteria, these would be misclassified as LERGs (and likely as star-forming LERGs) instead of HERGs. Although this issue is inevitable in the absence of emission line diagnostics that allow a more robust selection of optical AGN, we consider the impact of this potential contamination on our results.

Fig. B1 shows the cosmic evolution of the quiescent LERGs, star-forming LERGs and HERGs between $0.5 < z \le 2.5$ based on our adopted source classification criteria from Best et al. (in preparation), plotted with triangle (red), square (blue), and cross (black) symbols, respectively. The shaded regions correspond to the $1\sigma$ uncertainties on the LFs. Also shown in Fig. B1 are the LFs for these three classes of AGN obtained using three variations of the ‘SED AGN’ selection criteria from Best et al. (in preparation): SC-1, SC-2, and SC-3. As outlined in Section 2.5, based on the output from each of AGNFITTER and CIGALE, Best et al. (in preparation) defined a value $f_{\text{AGN}}$ for each code, as the fraction of the mid-infrared luminosity associated with AGN emission compared to the host galaxy. This diagnostic (in particular the 16th percentile of $f_{\text{AGN}}$) forms a key part in the selection of ‘SED AGN’, and hence the separation of ‘radio-excess AGN’ into LERGs and HERGs. We arbitrarily reduced the thresholds for these parameters from both AGNFITTER and CIGALE by a factor of two (i.e. leading to an increase in the number of HERGs) and reconstructed the LFs, as shown by the SC-1 (dotted) lines.

We also implemented additional criteria to that adopted by Best et al. (in preparation) in the selection of ‘SED AGN’, using the best-fitting $f_{\text{AGN}}$ values from AGNFITTER and CIGALE in addition to the 16th percentile values. We used $f_{\text{AGN}}$ as we wanted to select objects with high $f_{\text{AGN}}$ values, but with large uncertainties (e.g. due to low SNR), that may be missed by the Best et al. (in preparation) criteria; this allowed us to test a more complete but possibly contaminated sample, compared to the clean but possibly incomplete ‘SED AGN’ sample of Best et al. (in preparation). In order to find an appropriate threshold in $f_{\text{AGN}}$, we considered first the distribution of sources in the $f_{\text{AGN}}(\text{AGNFITTER})$ versus $f_{\text{AGN}}(\text{CIGALE})$ 2D parameter space. Fig. B2 shows a scatter plot in this 2D space, with the red, blue, and black points corresponding to the quiescent LERGs, star-forming LERGs, and HERGs, respectively, as classified based on the Best et al. (in preparation) criteria. It is reassuring to see that although
<figure><img src="image_17.png" /></figure>

**Figure B1.** Evolution of the luminosity functions of star-forming and quiescent LERGs, and the HERGs in the LoTSS-Deep fields using different source classification criteria. The separation of quiescent and star-forming LERGs is performed based on the sSFR criteria defined in Section 5.1. The markers and associated shaded regions correspond to the LFs and $1\sigma$ uncertainties resulting from the adopted source classification criteria from Best et al. (in preparation), whereas the lines show the LFs obtained under three variations of the source classification criteria (see the text). In general, we find that the adopted quiescent and star-forming LERGs, which are the focus of this paper, show little change in space densities with changes in source classification selection, suggesting that the results in the main paper are robust against uncertainties in the source classification scheme adopted.

<figure><img src="image_18.png" /></figure>

**Figure B2.** 2D scatter plot showing the location of the quiescent LERGs, star-forming LERGs, and HERGs (in ELAIS-N1), in AGN fraction (CIGALE) – AGN fraction (AGNFITTER) parameter space over a broad redshift bin. The dashed lines along the x and y-axes show the $f_{\text{AGN}}(\text{AGNFITTER}) = 0.25$ or $f_{\text{AGN}}(\text{CIGALE}) = 0.15$ thresholds used to select additional ‘SED AGN’ for testing the source classification criteria.

the best-fitting $f_{\text{AGN}}$ parameters were not used in the selection of the ‘SED AGN’, the HERGs occupy a well-defined region of this 2D plane, corresponding to high $f_{\text{AGN}}$ values from both the SED fitting codes, albeit with some overlap with the star-forming LERG population. Moreover, the quiescent and star-forming LERGs also show very similar distributions in best-fitting $f_{\text{AGN}}$ values (from both AGNFITTER and CIGALE) – and distinct from the HERGs – suggesting that these are both drawn from the same parent population

and that the star-forming LERG population is not significantly affected by contamination from HERGs. Based on the location of the HERGs, and that of X-ray detected AGN (indicative of radiative-mode AGN) from the X-Boötes survey in the 2D scatter plot in Fig. B2, we defined an extra criteria of ‘SED AGN’ as sources with $f_{\text{AGN}}(\text{CIGALE}) > 0.15$ and $f_{\text{AGN}}(\text{AGNFITTER}) > 0.25$. We note that this selection criterion is applied in addition to the Best et al. criteria, such that there will always be an increase in the number of sources classified as HERGs; in the case adopted here, this leads to 258 more sources being classified as HERGs. The LFs resulting from this criteria (SC-2) are shown as dashed lines in Fig. B1.

Finally, we considered adopting a more extreme approach such that sources with either $f_{\text{AGN}}(\text{CIGALE}) > 0.15$ or $f_{\text{AGN}}(\text{AGNFITTER}) > 0.25$ are classed as ‘SED AGN’. As seen from Fig. B2, while this selection may result in a higher completeness of ‘SED AGN’, it will suffer from even more contamination; this is particularly the case at higher redshifts. This is due to a significant population of sources (mostly star-forming LERGs based on the Best et al. classification) with $f_{\text{AGN}}(\text{AGNFITTER}) > 0.25$ but significantly lower (up to a factor of 10 lower) values of $f_{\text{AGN}}(\text{CIGALE})$, which is a result of the significantly higher uncertainties on the AGNFITTER values. Nevertheless, the SC-3 criteria allows us to examine the extent to which even an extreme classification could potentially affect the LERG LFs (dash-dotted lines in Fig. B1).

As evident from Fig. B1, given both the small number of HERGs in the LoTSS Deep Fields (see also Fig. 1) and the small numbers of sources with a change in their classification based on SC-1 and SC-2 criteria, there is only a very small effect on the LERG LFs across all redshifts (though slightly more significant for the HERG LFs). The SC-1 and SC-2 lines agree with the originally determined
Best et al. HERG LFs within 1 − 2σ of our best estimate LFs at all but the faintest few luminosity bins; at the faint end, the change in the LFs is comparable to the field-to-field variation found in the LFs. For the quiescent and star-forming LERGs, both SC-1 and SC-2 variations result in very similar LFs across redshift. The quiescent LERG population in particular remains robust against changes to the source classification method adopted within the uncertainties across all of the redshift bins; this also holds even when they become the sub-dominant population of all LERGs at z > 1.5. The star-forming LERG LFs at z ≤ 1 based on SC-1 and SC-2 are lower by ~0.1 − 0.2 dex, however this has little to no effect on the total LERG LFs given that the quiescent LERGs dominate the space densities over most of the radio luminosities studied here; at higher redshifts, the differences with the SC-1 and SC-2 based LFs are consistent within 1 − 2σ of the uncertainties.

For SC-3, the space density of the HERGs (2734 sources) now increases substantially across all luminosities, and especially at the faint end at $z > 1.5$, as might be expected. For the LERG population, however, which is the focus of this paper, the reconstructed LFs at $z \le 1.5$ show good agreement with our best estimate of the LERG space densities; at higher redshifts, the space densities decrease, in particular at faint luminosities, reaching ~0.2 dex lower in the faintest bins.

We therefore conclude that the LERG LFs presented in the main sections of the paper based on the Best et al. (in preparation) source classification – and subsequent results on the LERG stellar mass fractions – are robust to uncertainties in the source classification criteria used.

This paper has been typeset from a TeX/LA TeX file prepared by the author.