# Cosmic evolution of radio-AGN feedback: confronting models with data

Rohit Kondapally,1 Philip N. Best,1 Mojtaba Raouf,2 Nicole L. Thomas,3,4
Romeel Davé,1 Stanislav S. Shabala,5 Huub J. A. Röttgering,2 Martin J.
Hardcastle,6 Matteo Bonato,7 Rachel K. Cochrane,8 Katarzyna Małek,9,10 Leah.
K. Morabito,3,4 Isabella Prandoni7 and Daniel. J. B. Smith6
1SUPA, Institute for Astronomy, Royal Observatory, Blackford Hill, Edinburgh,
EH9 3HJ, UK
2Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The
Netherlands
3Centre for Extragalactic Astronomy, Department of Physics, Durham University,
Durham DH1 3LE, UK
4Institute for Computational Cosmology, Department of Physics, University of
Durham, South Road, Durham DH1 3LE, UK
5School of Natural Sciences, Private Bag 37, University of Tasmania, Hobart,
TAS 7001, Australia
6Centre for Astrophysics Research, Department of Physics, Astronomy and
Mathematics, University of Hertfordshire, College Lane, Hatfield AL10 9AB UK
7INAF-IRA, Via P. Gobetti 101, 40129, Bologna, Italy
8Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue,
New York, NY 10010, USA
9National Centre for Nuclear Research, ul. Pasteura 7, 02-093 Warszawa, Poland
10Aix Marseille Univ. CNRS, CNES, LAM, Marseille, France
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Radio-mode feedback is a key ingredient in galaxy formation and evolution
models, required to reproduce the observed properties of massive galaxies in
the local Universe. We study the cosmic evolution of radio-AGN feedback out to
$z\sim 2.5$ using a sample of 9485 radio-excess AGN. We combine the evolving
radio luminosity functions with a radio luminosity scaling relationship to
estimate AGN jet kinetic powers and derive the cosmic evolution of the kinetic
luminosity density, $\Omega_{\rm{kin}}$ (i.e. the volume-averaged heating
output). Compared to all radio-AGN, low-excitation radio galaxies (LERGs)
dominate the feedback activity out to $z\sim 2.5$, with both these populations
showing a constant heating output of $\Omega_{\rm{kin}}\approx 4-5\times
10^{32}\,\rm{W\,Mpc^{-3}}$ across $0.5<z<2.5$. We compare our observations to
predictions from semi-analytical and hydrodynamical simulations, which broadly
match the observed evolution in $\Omega_{\rm{kin}}$, although their absolute
normalisation varies. Comparison to the Semi-Analytic Galaxy Evolution (sage)
model suggests that radio-AGN may provide sufficient heating to offset
radiative cooling losses, providing evidence for a self-regulated AGN feedback
cycle. We integrate the kinetic luminosity density across cosmic time to
obtain the kinetic energy density output from AGN jets throughout cosmic
history to be $\sim 10^{50}\,\rm{J\,Mpc^{-3}}$. Compared to AGN winds, the
kinetic energy density from AGN jets dominates the energy budget at $z\lesssim
2$; this suggests that AGN jets play an important role in AGN feedback across
most of cosmic history.

###### keywords:

galaxies: active – galaxies: evolution – radio continuum: galaxies – galaxies:
jets

††pubyear: 2022††pagerange: Cosmic evolution of radio-AGN feedback:
confronting models with data–B

##  1 Introduction

It is now widely believed that most galaxies in the Universe host a
supermassive black hole (SMBH) at their centres. Significant evidence over the
past two decades has shown the existence of a tight correlation between the
mass of the SMBH and that of its galaxy’s central disk (Gebhardt et al., 2000;
Merritt & Ferrarese, 2001), and between the growth of the SMBHs and that of
the galaxies (e.g. Kormendy & Ho, 2013). As these SMBHs accrete matter, during
which they are known as active galactic nuclei (AGN), they release vast
amounts of energy into their surroundings which can suppress star-formation
and thereby regulate subsequent growth of their host galaxies (e.g. Bower et
al., 2006; Croton et al., 2006). This so-called AGN feedback effect is thought
to result in the observed co-evolution between the SMBHs and their host
galaxies (e.g. Best et al., 2005, 2006; McNamara & Nulsen, 2007; Cattaneo et
al., 2009; Fabian, 2012; Kormendy & Ho, 2013; Heckman & Best, 2014; Hardcastle
& Croston, 2020).

AGN feedback was initially introduced within semi-analytical models of galaxy
formation and evolution to solve the cooling-flow problem and reproduce the
observed high-mass end of the galaxy luminosity function (e.g. Bower et al.,
2006; Croton et al., 2006). This feedback from AGN was related to heating of
the surrounding gas within the halo and preventing run-away star-formation,
thereby limiting the growth of the most massive galaxies. This mode of
feedback is commonly referred to as ‘radio-mode’ feedback. The implementation
of ‘radio-mode’ feedback within semi-analytical models was motivated by
observations of radio-jets from brightest cluster galaxies creating cavities
in hot gas within massive haloes (e.g. Boehringer et al., 1993; Fabian et al.,
2006).

Feedback from AGN now forms a key ingredient in semi-analytical and
hydrodynamical galaxy formation models, required to reproduce observed galaxy
properties (e.g. Bower et al., 2006; Croton et al., 2006; Schaye et al., 2015;
Croton et al., 2016; Kaviraj et al., 2017; Springel et al., 2018; Davé et al.,
2019). However, the implementation of the growth of black holes and AGN
feedback models varies across different simulations. Moreover, the physical
processes governing the growth of black holes and the formation and evolution
of AGN jets operate on much smaller scales than the resolution of
cosmological-scale hydrodynamical simulations, and are hence implemented using
sub-grid prescriptions. The growth of black holes coupled with the growth and
evolution of the gas and host galaxies can offer valuable insights, and when
compared to observations can help constrain our understanding of black hole
fuelling and feedback models.

Observationally, feedback from ‘radio-mode’ AGN stems from two types of radio-
loud AGN: high-excitation radio galaxies (HERGs) and low-excitation radio
galaxies (LERGs), classified based on the nature of their optical emission
line properties (e.g. Best & Heckman, 2012). Detailed characterisation of
these populations, primarily using wide-area radio surveys and optical
spectroscopy, has revealed differences in both the properties of the central
engine and the host galaxy properties. These observations show that LERGs have
more massive black holes and tend to be hosted in more massive, quiescent,
redder galaxies in richer environments than the HERGs (e.g. Best et al., 2006;
Hardcastle et al., 2007; Tasse et al., 2008; Smolčić, 2009; Best & Heckman,
2012; Janssen et al., 2012; Sabater et al., 2013; Heckman & Best, 2014; Mingo
et al., 2014; Tadhunter, 2016; Ching et al., 2017; Croston et al., 2019;
Hardcastle & Croston, 2020; Magliocchetti, 2022).

Studies of LERGs and HERGs in the nearby Universe have suggested that the
differences in the observed properties of LERGs and HERGs may arise from
differences in the Eddington-scaled accretion rates onto the SMBH (e.g. Best
et al., 2005; Hardcastle et al., 2007; Mingo et al., 2014; Delvecchio et al.,
2017; Hardcastle, 2018a; Hardcastle & Croston, 2020). HERGs are thought to be
undergoing radiatively-efficient accretion at high fractions of the Eddington-
scaled accretion rates, typically from cold gas, leading to the formation of
an optically-thick, geometrically thin accretion disk and torus structure
(e.g. Shakura & Sunyaev, 1973). LERGs on the other hand are fuelled at low
Eddington-scaled accretion rates by cooling hot gas within their haloes in an
advection dominated accretion flow (e.g. Narayan & Yi, 1995; Sharma et al.,
2012; Gaspari et al., 2013), and hence lack a radiatively-efficient accretion
disk. However, recent deeper radio continuum surveys, probing over an order of
magnitude fainter radio luminositites and higher redshifts, have found more
overlap in key host galaxy and AGN properties between the two modes,
suggesting that the picture may be more complicated (e.g. Whittam et al.,
2018; Mingo et al., 2022; Kondapally et al., 2022; Whittam et al., 2022).

The radio luminosity functions of the LERGs and HERGs evolve differently, with
the LERGs dominating the space densities at lower luminosities compared to the
HERGs at low redshifts (e.g. Best & Heckman, 2012; Gendre et al., 2013; Pracy
et al., 2016); this highlights the importance of studying both the LERG and
HERG populations separately as the precise origins of the differences in
fuelling and feedback between the HERGs and LERGs remain unclear. Determining
the cosmic evolution of these AGN, their host galaxies, and their feedback
effect are crucial in understanding their role in galaxy evolution across
cosmic time. Best et al. (2014) were the first to study the evolution of LERGs
and HERGs, separately, out to $z\sim 1$; they found that the HERGs showed a
strong evolution with redshift, whereas the LERGs showed an overall mild
evolution. Similar results were later found by Pracy et al. (2016) and Butler
et al. (2019). At low frequencies however, Williams et al. (2018) extended the
analysis of the luminosity functions of LERGs across $0.5<z\leq 2$ using 150
MHz LOw Frequency ARray (LOFAR; van Haarlem et al. 2013) observations of the
Boötes field, finding a strong decline in their space densities with
increasing redshift.

Recently, Kondapally et al. (2022) studied the evolution of the total radio-
AGN and LERG luminosity functions out to $z\sim 2.5$ using the LOFAR Two-metre
Sky Survey Deep Fields Data Release 1 (LoTSS-Deep DR1; Tasse et al. 2021;
Sabater et al. 2021; Kondapally et al. 2021; Duncan et al. 2021; Best et al.
2023. LoTSS-Deep forms one of the deepest wide-field radio continuum surveys
and covers a sky area of $\sim$ 25 $\rm{deg^{2}}$ (where there is overlap with
high-quality multi-wavelength data), and detects $>11\,000$ radio-AGN
(including $>10\,000$ LERGs); this allowed them to probe much fainter radio
luminosities than many previous studies (e.g. Best et al., 2014; Pracy et al.,
2016; Williams et al., 2018; Butler et al., 2019) and study the low-luminosity
AGN in detail while also tracing a wide range of galaxy environments and
better constraining the bright-end of the luminosity function than deep
narrow-area surveys (e.g. Smolčić et al., 2017). Kondapally et al. (2022)
found that the LERG luminosity functions showed mild evolution across
$0.5<z\leq 2.5$ (the differences with the Williams et al. (2018) results were
found to be largely due to differences in the source classification schemes
employed; see Kondapally et al. 2022). When split by host galaxy type
(quiescent versus star-forming), they found that the quiescent LERGs showed a
strong decline in their space densities with increasing redshift; LERGs hosted
by star-forming galaxies (SFGs) become more prominent across all luminosities
at $z\gtrsim 1$. Moreover, they found that quiescent LERGs showed a roughly
constant duty-cycle over the past $\sim$ 10 Gyrs and that the strong negative
evolution of the quiescent LERGs was in accordance with the space density of
massive quiescent galaxies as their hosts. The radio-mode feedback models
within simulations are often tuned to balance radiative cooling losses in
massive quiescent hosts, in order to prevent gas cooling (and hence star
formation) in these systems; therefore the different evolution seen from the
different modes of radio-AGN has interesting implications for these feedback
models.

In this paper, we study the cosmic evolution of radio-mode feedback by using
the evolving radio luminosity functions to calculate the kinetic luminosity
densities of radio-AGN and LERGs from LoTSS-Deep; this traces the volume
integrated total heating output in the form of mechanical (kinetic) power of
the radio jets. The observational measurements are then compared with
predictions from hydrodynamical and semi-analytical galaxy formation models.
The paper is structured as follows. In Section 2 we describe the radio and
multi-wavelength datasets used and the selection of radio-AGN and LERGs.
Section 3 presents a comparison of the cosmic kinetic luminosity density to
other observations and simulations. In Section 4, we calculate the total
kinetic energy density output across cosmic history in the form of AGN jets to
study the importance of jet AGN feedback. We present our conclusions in
Section 5. Throughout this work, we use a flat $\Lambda$CDM cosmology with
$\Omega_{\rm{m}}=0.3$, $\Omega_{\Lambda}=0.7$, and
$H_{\rm{0}}=70~{}\rm{km\,s^{-1}\,Mpc^{-1}}$, and for calculating radio
luminosities, a radio spectral index $\alpha=-0.7$ (where
$S_{\nu}\propto\nu^{\alpha}$).

##  2 Data

###  2.1 Radio and multi-wavelength data

The radio catalogues are taken from the LoTSS Deep Fields DR1 (Tasse et al.,
2021; Sabater et al., 2021). The radio observations consist of deep repeated
LOFAR High Band Antenna (HBA) observations (with an angular resolution of 6
arcsec) of the ELAIS-N1, Lockman Hole, and Boötes fields that reach an rms
sensitivity of 20, 22, and 32 $\rm{\mu Jy\,beam^{-1}}$ at 150 MHz in the
centres of each field, respectively. The LoTSS-Deep dataset forms the deepest
wide-field radio-continuum survey at low-frequencies to date; this makes it
ideal for studying the cosmic evolution of the faint radio-AGN population. The
details of the radio calibration, imaging, and comparison to other radio-
continuum surveys are presented by Tasse et al. (2021); Sabater et al. (2021).

The three LoTSS-Deep fields were chosen due to the availability of deep, wide-
field, multi-wavelength imaging (see Kondapally et al., 2021, and references
therein), making these fields ideal for characterising the physical properties
of the radio sources. In summary, the three LoTSS-Deep fields have coverage
from the ultra-violet (GALEX; Morrissey et al., 2007), optical (PanSTARRS
Medium Deep Survey; Chambers et al. 2016, Hyper Suprime-Cam Subaru Strategic
Program; Aihara et al. 2019, and the NOAO Deep Wide-Field Survey; Jannuzi &
Dey 1999), near-infrared (UK Infrared Telescope Deep Sky Survey in ELAIS-N1
and Lockman Hole; Lawrence et al. 2007, and J, H, Ks data in Boötes from
Gonzalez et al. 2010), mid-infrared (Spitzer IRAC and MIPS surveys; Lonsdale
et al. 2003; Mauduit et al. 2012; Eisenhardt et al. 2004; Ashby et al. 2009,
and far-infrared (Herschel; Oliver et al. 2012) wavelengths.

The identification of the host galaxies of the radio-detected sources, and the
association of radio components (where the radio source finder had not
correctly grouped physical sources in the catalogue) was performed by
Kondapally et al. (2021). The host galaxies were identified using a
combination of the statistical likelihood ratio method (de Ruiter et al.,
1977; Sutherland & Saunders, 1992) and a visual classification scheme (see
also Williams et al., 2019), whereas the source association was performed
using visual classification only. While the radio data covers a much larger
area ($\sim 68\,{\rm{deg^{2}}}$ in each field), the host galaxy identification
process was limited to the areas in each field with the best available multi-
wavelength data; this process resulted in host galaxies for $>$ 97 per cent of
the radio-detected sources, with the final radio catalogue consisting of 81
951 radio sources, across $\sim 25\,{\rm{deg^{2}}}$ over the three fields
(Kondapally et al., 2021). Photometric redshifts for the radio sources (and
the underlying multi-wavelength catalogues) were generated by Duncan et al.
(2021) using a hybrid machine learning and template-based approach optimised
for deep radio-continuum surveys (see Duncan et al., 2018a; Duncan et al.,
2018b).

###  2.2 SED fitting and source classification

The multi-wavelength photometry and photometric redshifts were used to perform
spectral energy distribution (SED) fitting for all of the radio-detected
sources(Best et al., 2023). For each source, the SED fitting process was
performed using four different routines: 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 et al., 2005; Noll et al., 2009; Boquien et al.,
2019), and Multi-wavelength Analysis of Galaxy Physical Properties (magphys;
da Cunha et al., 2008).

The output of this SED fitting process was used to perform source
classifications and derive stellar masses and star-formation rates (SFRs) for
the full radio dataset. In summary, Best et al. (2023) used fits to the AGN
accretion disk and torus models present in both AGNFitter and cigale to define
a parameter, ‘AGN fraction’, which represents the fraction of the mid-infrared
emission arising from AGN as compared to the host galaxy components. This ‘AGN
fraction’, along with a comparison of the goodness of fit from the SED fitting
codes with and without AGN components, were used to identify the so-called
‘optical’ (SED) AGN (also known as radiative-mode AGN). For such sources, the
cigale stellar masses and SFRs were adopted, as these were found to have lower
scatter than those from AGNfitter, while the masses estimated from magphys or
bagpipes for these sources may be inaccurate due to the lack of an AGN
component in these codes. For sources that did not show signs of a radiative-
mode AGN, the average of the bagpipes and magphys results were generally used
to derive ‘consensus’ stellar masses and SFRs as these codes include better
sampling of the range of potential stellar populations. Then, to identify the
radio-AGN, Best et al. (2023) selected sources whose radio emission
predominantly arose from the AGN by identifying sources that showed an excess
of $>0.7\,\rm{dex}$ ($\sim 3\sigma$) in radio luminosity, over the level
expected from star-formation processes alone based on the consensus SFRs and
SFR-radio luminosity relation for star-forming galaxies (e.g. Gürkan et al.,
2018; Smith et al., 2021).

HERGs are the radio-loud subset of the radiative-mode AGN population (i.e.
also display signs of an AGN in their SED); these are identified as sources
showing radio-excess and an SED AGN. The LERGs are identified by the presence
of radio-jets only; these were hence identified as sources that display a
radio-excess AGN but were not selected as AGN based on their SEDs. The radio
luminosity functions for these LERG and HERG samples were constructed and
discussed by Kondapally et al. (2022).

In this study, we follow Kondapally et al. (2022) in selecting quiescent
galaxies using the specific star-formation rates (sSFRs) of galaxies, such
that sources which satisfy $\rm{sSFR}<0.2/t_{\rm{H}}(z)$ are defined as
quiescent galaxies, where $t_{\rm{H}(z)}$ is the age of the Universe at
redshift $z$; this criterion was found to be broadly consistent with quiescent
galaxies selected using rest-frame UV colour-colour diagrams (e.g. Pacifici et
al., 2016; Carnall et al., 2018; Carnall et al., 2020). In total, across the
redshifts $0.5<z\leq 2.5$ analysed in this study, our sample consists of 9485
radio-excess AGN, of which 8409 are LERGs (with 2974 of these LERGs being
hosted in quiescent galaxies; hereafter quiescent LERGs or Q-LERGs).

We note that the identification of radiative-mode AGN by Best et al. (2023)
may be incomplete (leading to HERGs being mis-classified as LERGs), especially
in the absence of X-ray data (which is not available in ELAIS-N1 and Lockman
Hole). To test this, we utilised the X-Boötes Chandra survey (Kenter et al.,
2005); this dataset was also used to identify X-ray detected AGN during the
source-classification process in the Boötes field (see Duncan et al., 2021;
Best et al., 2023). We find that only $\sim 5\%$ of the radiative-mode AGN in
Boötes were identified as AGN using X-ray data alone. Therefore, while the
lack of wide-field X-ray observations in ELAIS-N1 and Lockman Hole may result
in a small fraction of mis-classifications, Best et al. (2023) found that the
fraction of sources classified as LERGs and HERGs across the three fields were
consistent with each other. Furthermore, Kondapally et al. (2022) also
constructed the LERG luminosity functions in each of the three LoTSS-Deep
fields separately, finding good agreement across the fields. These results
indicate that the lack of X-ray data in ELAIS-N1 and Lockman Hole does not
have a significant effect on the derived luminosity functions and results in
this paper.

##  3 Evolution of radio-AGN feedback over cosmic time

The observed 150 MHz radio luminosity functions (LFs) for the radio-AGN in
LoTSS-Deep were presented by Kondapally et al. (2022). The LFs were calculated
using the $1/V_{\rm{max}}$ method and the evolution of the radio-excess AGN,
LERGs, Q-LERGS (i.e. LERGs hosted in quiescent galaxies) was characterised.
The observed LFs and the best-fitting evolution models derived by Kondapally
et al. (2022) are shown in Appendix A, with the radio-excess AGN LFs
recomputed over different redshift bins for the analysis in this paper. In
this section, we consider the implications of the evolution of the observed
radio-AGN LFs from Kondapally et al. (2022) on the evolution of the amount of
energy deposited into their host galaxies by the radio-AGN across cosmic time.
We also perform a comparison with simulations in physical space by using
scaling relations to convert observed radio luminosities into jet kinetic
powers.

###  3.1 Kinetic powers of radio-AGN

The kinetic energy carried in the jets can provide a significant energetic
output in the form of work done on the surrounding environment; this so-called
mechanical (kinetic) power from the jets can be orders of magnitude larger
than the monochromatic radio luminosities. The jet kinetic power is often
estimated from observed radio luminosities by using a scaling relation. One
such relation often used in literature (e.g. Smolčić et al., 2017; Butler et
al., 2019) was derived by Willott et al. (1999) who determined the kinetic
powers by using minimum energy arguments to estimate the energy stored within
lobes, and combined this with estimates of radio source lifetimes and energy
losses from inflating the radio source. Uncertainties in our knowledge of the
physics of radio jets and their composition, along with departures from the
minimum energy condition result in significant uncertainties in the
calibration of this relation; Willott et al. (1999) combined all these
uncertainties into a single parameter, $f_{\rm{W}}$, with typical values in
the range of $f_{\rm{W}}\sim 1-20$. Converted to 1.4 GHz luminosity (Heckman &
Best, 2014), this relation is given as:

$L_{\rm{kin,\,sync}}=4\times
10^{35}\left(f_{\rm{W}}\right)^{3/2}\left(L_{\rm{1.4\,GHz}}/10^{25}\,\rm{W\,Hz^{-1}}\right)^{0.86}\,\rm{W},$
(1)

where $L_{\rm{1.4\,GHz}}$ is the 1.4 GHz radio luminosity, and $f_{\rm{W}}$ is
the uncertainty parameter on the calibration of the scaling relation.

Another method of calculating jet kinetic powers is based on studying the
cavities in the hot gas created as the radio jets plough through the
surrounding material, which can be observed in the X-rays (Boehringer et al.,
1993). The jet kinetic power can be estimated from the radio luminosity by
considering the energy required to expand the lobes and inflate the X-ray
cavities of pressure $p$ with a volume $V$, $E_{\rm{cav}}=f_{\rm{cav}}\,pV$,
where $f_{\rm{cav}}=4$ is the commonly adopted value (e.g. Cavagnolo et al.,
2010; Best et al., 2014; Pracy et al., 2016; Butler et al., 2019),
corresponding to the $pV$ of work done, and $3pV$ of energy stored in the
relativistic particles; a good correlation is seen between the cavity powers
and 1.4 GHz radio luminosities (e.g. Bîrzan et al., 2004; Cavagnolo et al.,
2010; Timmerman et al., 2022). Heckman & Best (2014) derived a best-fit
relation, largely based on the results from Bîrzan et al. (2008); Cavagnolo et
al. (2010), given as

$L_{\rm{kin,\,cav}}=7\times
10^{36}f_{\rm{cav}}\left(L_{\rm{1.4\,GHz}}/10^{25}\,\rm{W\,Hz^{-1}}\right)^{0.68}\rm{W}.$
(2)

The normalisation of this relation is found to be in good agreement with that
of the Willott et al. (1999) relation when using $f_{\rm{W}}=15$ and
$f_{\rm{cav}}=4$ (e.g. Heckman & Best, 2014; Smolčić et al., 2017), however
the Heckman & Best (2014) relation has a shallower slope which will result in
higher jet powers at low luminosities.

It is important to note that the above relations have a large scatter that is
dominated by the systematic effects and assumptions about the unknown physics
of the radio sources. Therefore a simple conversion between radio luminosity
and kinetic power is likely not appropriate; even for jets of a consistent
kinetic power, the radio luminosity varies over the lifetime of a radio
source, and depends on the energy density and magnetic field of the radio
lobes and hence the environment into which the radio lobes are expanding, and
is also influenced by the assumed spectral index (see Shabala et al., 2012;
Shabala & Godfrey, 2013; Hardcastle & Krause, 2013; Godfrey & Shabala, 2016;
Croston et al., 2018; Hardcastle, 2018b; Hardcastle et al., 2019). However,
while a relation based on radio luminosity alone may not be accurate for
inferring jet kinetic powers for individual sources, we are interested in the
heating output from the AGN at a population level; the use of either of the
above scaling relationships for typical values of the uncertainty parameters
(see below) should therefore provide a reasonable mean value.

Figure 1: Cosmic evolution of the heating rate as a function of radio
luminosity in four redshift bins for LERGs (black), Q-LERGs (red), and radio-
excess AGN (blue). These are calculated by convolving the best-fit LFs with
the jet kinetic powers (estimated using equation 2), with the shaded regions
showing uncertainties on the best-fitting LF models. The heating rate function
peaks at the break of the luminosity function, and therefore most of the
heating output comes from relatively high luminosity sources.

An additional issue is that the above scaling relation has been determined
from characterisation of sources at low-redshifts and it is possible that the
radio luminosity of the sources of a given jet power (and hence the above
relation) will evolve with redshift (e.g. due to cosmic evolution of magnetic
field strengths, or increasing inverse-Compton losses); a well-constrained
relation out to high redshifts is, however, still lacking. Recently,
Hardcastle et al. (2019) used the projected linear source sizes, redshifts,
and radio luminosities to infer the jet kinetic powers of resolved AGN out to
$z\sim 0.7$ from the first data release of the wide-area LoTSS survey
(Shimwell et al., 2019) using the dynamical model from Hardcastle (2018b).
High-resolution imaging using the LOFAR international baselines for LoTSS-Deep
(e.g. Sweijen et al., 2022) will provide more robust sizes and help extend the
jet power inference models to higher redshifts (Hardcastle et al. in prep.).

In the analysis that follows, we use the Heckman & Best (2014) relation
(Equation 2) with $f_{\rm{cav}}=4$ and assume that this local relation is
applicable over the full redshift range studied, but note that this might lead
to systematic errors. In Appendix B, we study the impact of using different
scaling relations for estimating jet powers on our results. We find that the
volume-integrated heating rate across cosmic time predicted from these two
relations (Equations 1 and 2) is in good agreement (see Fig. 6), while the use
of some other relations in the literature would also produce almost identical
trends in the cosmic evolution of the volume-integrated heating rate, but with
typical changes to the overall normalisation by around 0.3 dex. This
normalisation uncertainty needs to be borne in mind when comparing the output
of the simulations with the observational data.

Figure 2: The evolution of the kinetic luminosity density,
$\Omega_{\rm{kin}}(z)$ compared to observations in literature (left) and to
simulations (right). Measurements for all LERGs (black points), Q-LERGs (red
triangles), and radio-excess AGN (blue open circles) in the LoTSS Deep Fields
are shown on both panels. The shaded regions in each case represent the
uncertainties on the kinetic luminosity density based on uncertainties in
fitting the evolution of the LFs. The kinetic luminosity density was
calculated using equation 2 with $f_{\rm{cav}}=4$ and the best fit LDE models
for the three AGN populations (see Appendix A). The black lines show the
systematic shift that would be obtained for the total LERG heating rate using
the values $f_{\rm{cav}}=1$ and $f_{\rm{cav}}=10$ that account for the
systematic uncertainties in the estimate of kinetic powers. Observational
results for the LERG population from Best et al. (2014) and Butler et al.
(2019), along with radio-excess AGN from Smolčić et al. (2017), and radio-AGN
from Hardcastle et al. (2019) are shown in the left panel. Our measurements
are compared to predictions from sage (Croton et al., 2016), radio-sage (Raouf
et al., 2017), simba (Davé et al., 2019; Thomas et al., 2021), and Mocz et al.
(2013; their LERG-equivalent population) in the right panel.

###  3.2 Cosmic evolution of the kinetic luminosity density

The evolving radio luminosity functions of AGN from Kondapally et al. (2022;
also shown in Appendix A) can be combined with equation 2 to estimate the
kinetic heating rate as a function of radio luminosity (also known as the
specific heating rate function), given as

$\Psi(L_{\rm{1.4\,GHz}},z)=L_{\rm{kin}}(L_{\rm{1.4\,GHz}})\times\rho\left(L_{\rm{1.4\,GHz}},z\right)\,\rm{W\,Mpc^{-3}\,logL^{-1}}.$
(3)

Using the above best-fitting models for the evolution of the LFs, we
calculated the specific kinetic heating rate for the radio-excess AGN, the
total LERG, and the Q-LERG populations separately. For each population, the
best-fitting LF models were shifted to 1.4 GHz using the standard spectral
index $\alpha=-0.7$ to obtain $\rho(L_{\rm{1.4\,GHz}},z)$ for each redshift
bin; this was then convolved with $L_{\rm{kin}}(L_{\rm{1.4\,GHz}})$
(calculated using equation 2) to determine the kinetic heating rate as a
function of radio luminosity. The resulting curves for the radio-excess AGN,
LERGs, and Q-LERGs are shown in Fig. 1 with the shaded regions corresponding
to uncertainties in the modelling of the evolution of the LFs. For all three
populations, the heating rate function peaks at high radio luminosities, near
the break in the LFs; the location of this peak occurs at slightly higher
luminosities at higher redshifts. At $z<1$, the majority of the heating rate
output from LERGs comes from the Q-LERGs across all radio luminosities, as
would be expected given the LFs observed in Fig. 4. The heating rate from
Q-LERGs peaks at higher radio luminosities than that of the other populations.
This may in part be due to the uncertainty in modelling the break luminosity
in the LFs for the LERGs and radio-excess AGN; for Q-LERGs, the characteristic
space density at each redshift is fixed based on the available host galaxies,
resulting in one less free parameter at each redshift interval. At higher
redshifts, while the heating rates from LERGs show little evolution in shape
and normalisation, the heating output from Q-LERGs decline sharply, in line
with the observed evolution of the LFs.

The specific heating rate function from Eq. 3 can be integrated to estimate
the total kinetic luminosity density (also known as the kinetic heating rate),
$\Omega_{\rm{kin}}(z)$. This is given, in units of $\rm{W\,Mpc^{-3}}$, as

$\Omega_{\rm{kin}}(z)=\int\Psi(L_{\rm{1.4\,GHz}},z)\,\rm{d}\log
L_{\rm{1.4\,GHz}}\,\rm{W\,Mpc^{-3}}.$ (4)

For each redshift bin in Fig 1, the specific heating rate function is
integrated as a function of the radio luminosity down to
$10^{-4}\,L_{\star}(z)$ to obtain the redshift evolution of the integrated
kinetic luminosity density for each of the total LERGs (black circles), the
Q-LERGs (red triangles), and the radio-excess AGN (blue open circles), as
shown in Fig. 2. In each case, the shaded region corresponds to the
uncertainties in the kinetic luminosity density based on uncertainties in
fitting the evolution of the LFs. The narrow black lines show the systematic
shift that would be produced for the total LERGs based on the extreme values
of $f_{\rm{cav}}=1$ and $f_{\rm{cav}}=10$; this illustrates the effect of
uncertainties in the calibration of the radio-luminosity to jet kinetic power
relationship on the absolute values obtained; changes due to different
calibration relations discussed in Appendix B are typically $\lesssim
0.3\,\rm{dex}$, so smaller than these extreme values of $f_{\rm{cav}}$.

The kinetic luminosity density for both the radio-excess AGN and the LERGs
remains roughly constant between $0.5<z\leq 2.5$ with
$\Omega_{\rm{kin}}\approx 4-5\times 10^{32}\,\rm{W\,Mpc^{-3}}$ for both
populations at $z=0.75$; this indicates that the dominant source of heating
from the radio-loud population comes from the LERGs (rather than the HERGs)
across redshift, and that at least out to $z\sim 2.5$, the importance of
feedback from LERGs is broadly uniform. In contrast, for the Q-LERG
population, the kinetic luminosity density decreases steadily with increasing
redshift, from $\Omega_{\rm{kin}}\approx 2.7\times 10^{32}\,\rm{W\,Mpc^{-3}}$
at $z=0.75$ to nearly an order of magnitude lower by $z\sim 2.5$.

In Fig. 2 (left), we compare our observations with other measurements of the
kinetic luminosity densities from the literature. The estimates obtained for
the ‘jet-mode AGN’ population by Best et al. (2014) for $z<1$, also derived
using the Heckman & Best (2014) jet power scaling relation with
$f_{\rm{cav}}=4$ are shown by green symbols. The results from Best et al.
(2014) show an increase out to $z\sim 0.6$ and then show a decrease in the
kinetic luminosity density, resulting in better agreement with our Q-LERG
population, within errors, than the total LERG population. One potential
reason for the better agreement with the quiescent-LERGs could be due to the
method by which the ‘jet-mode AGN’ were selected by Best et al. (2014); these
were identified using emission line ratio diagnostics to select sources with
relatively low emission line fluxes or equivalent widths from [O ii] or [O
iii] lines. Jet-mode AGN hosts with considerable star-forming activity will
result in higher equivalent widths of these spectral lines and therefore the
spectroscopic classification used by Best et al. (2014) may result in a sample
similar to the quiescent-LERG population identified in this study.

The grey dash-dotted line illustrates the kinetic luminosity density for the
LERGs obtained using a pure density evolution model for the evolution of the
LFs as determined by Butler et al. (2019) out to $z\sim 1.3$, using the
Cavagnolo et al. (2010) relation. The results from Butler et al. (2019) show a
steady increase with increasing redshift and agree well with our estimates for
the total LERG population, even when extrapolating their results to higher
redshifts. We note that even if the Best et al. (2014) selection is similar to
our Q-LERGs, the good agreement between the Best et al. (2014) and Butler et
al. (2019) results is expected at $z<0.6$, where the total LERG population is
dominated by Q-LERGs. The dashed purple line shows the results for all radio-
AGN from Smolčić et al. (2017) determined using the Willott et al. (1999)
relation with $f_{\rm{W}}=15$; their results agree well with our measurements
for the evolution of the radio-excess AGN kinetic luminosity density across
redshift.

We also compare to results from Hardcastle et al. (2019), who used data from
LoTSS DR1 to calculate the jet kinetic powers for radio-AGN by incorporating
the projected source sizes using the analytic model from Hardcastle (2018b).
Hardcastle et al. (2019) used their inferred jet powers to calculate the jet
kinetic luminosity function at $z<0.7$ and integrated this to find
$\Omega_{\rm{kin}}\sim 7\times 10^{31}\,\rm{W\,Mpc^{-3}}$ (shown as a blue
star in Fig. 2). This result is systematically lower (by a factor of $\sim
2-3$) than other observations shown in Fig. 2 at this redshift, however this
offset is found to be due to the different methods of estimating jet kinetic
powers. Similarly, Turner & Shabala (2015) found a lower kinetic luminosity
density using their dynamical model than the results presented in this work.
Unlike Hardcastle et al. (2019), our study and those by Best et al. (2014);
Smolčić et al. (2017); Butler et al. (2019) use radio luminosity to jet-power
scaling relationships which generally predict higher jet powers at low
luminosities (see fig. A4 of Hardcastle et al., 2019, and Appendix B). We note
that if we use the Willott et al. (1999) scaling relationship with
$f_{\rm{W}}=4$ instead, which will result in lower jet powers at low
luminosities compared to the Heckman & Best (2014) relation adopted in this
study, our measurement for the kinetic luminosity density of radio-excess AGN
would match well with the results from Hardcastle et al. (2019).

###  3.3 Comparison of kinetic luminosity density with simulations

In Fig. 2 (right), we compare the observed kinetic luminosity densities with
predictions from recent simulations. We first compare with predictions from
the Semi-Analytic Galaxy Evolution (sage) model (Croton et al., 2016), which
provides an enhancement to the previous model by Croton et al. (2006). Croton
et al. (2006) considered the prescription of SMBH growth and AGN feedback in
two classes: “quasar mode” and “radio mode” in their terminology. We used the
predictions of the black hole accretion rate density ($\dot{m}_{\rm{BH,R}}$)
over cosmic time for the “radio-mode” population from Croton et al. (2006, see
their fig. 3). The $\dot{m}_{\rm{BH,R}}(z)$ were then translated into a
luminosity of the black hole assuming
$L_{\rm{BH}}=\eta\dot{m}_{\rm{BH,R}}c^{2}$, where $\eta=0.1$ is the standard
black-hole accretion efficiency and $c$ is the speed of light. This luminosity
generated from the accretion process is considered equivalent to the kinetic
luminosity density output from the LERGs derived in this study. As part of an
update to this model, Croton et al. (2016) introduced a “radio mode
efficiency” parameter, $\kappa_{\rm{R}}$ within sage to provide a more
realistic treatment of AGN feedback cycle by attempting to couple the heating
provided by the AGN with the cooling. As a result, in addition to $\eta=0.1$,
Croton et al. (2016) scale their $\dot{m}_{\rm{BH,R}}(z)$ and hence the black
hole luminosities by $\kappa_{\rm{R}}=0.08$ (which we also apply here).

The resulting prediction from the sage model is shown in the right panel of
Fig. 2 (dot-dot dashed blue line). The radio-excess AGN measurements show a
similar shape to the sage model prediction, but offset to higher values by
nearly an order of magnitude. This could be due to calibration errors in the
scaling relation, or due to some fraction of the heating from radio-AGN (in
particular from very extended sources) being deposited on scales larger than
that useful to offset cooling losses. Moreover, radio-mode feedback in the
sage model is scaled to a level required to provide AGN heating that can
offset radiative cooling, however, at higher redshifts in particular, the LERG
activity and heating output occurs predominantly in star-forming LERGs; the
feedback cycle from LERGs within star-forming galaxies remains unclear. The
kinetic luminosity density of the Q-LERG population is also systematically
higher than the sage model at all but the highest redshift bin; these results
therefore suggest that the energy output by the Q-LERG population alone is
sufficient to offset the radiative cooling losses within sage at least out to
$z\sim 2.5$, thereby providing evidence for a self-regulating AGN feedback
cycle in quiescent galaxies hosting a LERG.

We then compare our observations to predictions from another semi-analytical
model, radio-sage (Raouf et al., 2017), that builds upon the sage model to
focus on intermittent AGN jet activity, providing an improved modelling of
radio-mode feedback (Raouf et al., 2017). The accretion rate of matter onto
the black hole is approximated by the Bondi-Hoyle relation (Bondi, 1952), with
two accretion states implemented in radio-sage. The ‘hot-mode’ accretion
occurs at low Eddington-scaled accretion rates ($f_{\rm
Edd}\equiv\frac{\dot{M}_{\rm BH}}{\dot{M}_{\rm BH,Edd}}<\alpha_{\rm{crit}}$),
resulting in an advection dominated accretion flow (Narayan & Yi, 1994, 1995),
and ‘cold-mode’ (or radiatively-efficient) accretion occurs at high rates
($f_{\rm{Edd}}>\alpha_{\rm{crit}}$) resulting in the formation of an
optically-thin accretion disk (Shakura & Sunyaev, 1973). The ‘hot-mode’ (i.e.
with low Eddington-scaled accretion rates) sources are taken to be analogous
to the LERGs following observed differences in the accretion rate properties
of LERGs and HERGs (e.g. Best & Heckman, 2012); here, we used
$\alpha_{\rm{crit}}=0.03$ to select LERGs from radio-sage (Shabala &
Alexander, 2009). The quiescent galaxies in radio-sage are selected using the
same sSFR criterion as the observations, with the SFR calculated based on the
average over the past $\sim 280$ Myr (the time-step duration between each
redshift snapshot). The output jet power as a result of this accretion is then
estimated following $L_{\rm{kin,hot}}=\eta\dot{M}_{\rm{BH}}c^{2}$, where c is
the speed of light and $\eta$ is the accretion efficiency that is constrained
in radio-sage to be 0.35 using observations as described by Raouf et al.
(2017).

The radio-sage predictions for the LERG and Q-LERG populations (dashed black
and red lines, respectively) are also shown in Fig. 2. The integrated kinetic
luminosity densities for LERGs from radio-sage are in excellent agreement with
our observations, showing a gradual increase with redshift. The predictions
for Q-LERGs decline with redshift, showing a similar slope to the observations
but with a normalisation that is lower by up to $\sim 0.5\,\rm{dex}$; this
suggests that the star-forming galaxies in radio-sage play an important role
in radio-AGN feedback at all times. Given the good agreement with the total
LERG population, this also implies that radio-sage may not match the observed
space density of the Q-LERGs, with a significantly higher fraction of the
heating output from LERGs being performed in star-forming galaxies within
radio-sage; we will investigate this directly by comparing the radio LFs and
other host galaxy properties from simulations in future work. In addition,
radio-sage is able to model large sources which heat gas out to well beyond
the cooling radius; therefore more jet power is required for the same amount
of feedback within radio-sage compared to the more idealised feedback
efficiency model in sage. The observed size distribution of sources, using
high-resolution LOFAR data, can be compared to predictions from radio-sage in
a future paper.

We also compare our results to predictions from the simba cosmological
hydrodynamical simulation, which incorporates a unique sub-grid prescription
for black hole growth and feedback (Davé et al., 2019). simba employs a two-
mode accretion model, with the accretion from hot gas ($T>10^{5}\,\rm{K}$;
i.e. hot-mode) described by Bondi accretion (Bondi, 1952), and the accretion
from cold gas ($T<10^{5}\,\rm{K}$; i.e. cold-mode) being described by a
torque-limited accretion model (Hopkins & Quataert, 2011; Anglés-Alcázar et
al., 2017). Following Thomas et al. (2021), who studied the properties of
radio-AGN within simba, LERGs are identified as sources where the Bondi
accretion (or accretion via hot gas) dominates the black hole accretion rate
(based on the dichotomy in accretion rates observed for the LERGs and HERGs;
Best & Heckman 2012). This classification is performed based on the average
accretion rate (from both hot and cold gas) over the past 50 Myr period.
Subsequently, quiescent galaxies (and hence Q-LERGs) in simba are identified
using the same sSFR criterion as in the observations (see Sect. 2.1).
Following the prescription of Körding et al. (2008) used by Thomas et al.
(2021) to estimate radio-source observables, the jet kinetic powers for the
LERGs and Q-LERGs within simba are calculated as
$L_{\rm{kin}}=0.2\,L_{\rm{bol}}$, where
$L_{\rm{bol}}=0.1\dot{M}_{\rm{BH}}c^{2}$ is the bolometric luminosity of the
AGN, and $\dot{M}_{\rm{BH}}$ is the black hole accretion rate density.

The predictions from simba for both the LERGs and Q-LERGs (dotted black and
red lines, respectively) are also shown in Fig. 2. For LERGs, simba predicts
$\Omega_{\rm{kin}}$ that increases from $z=0$ and peaks at
$\Omega_{\rm{kin}}\approx 10^{33}\,\rm{W\,Mpc^{-3}}$ by $z\sim 1$ before
declining at higher redshifts. Although these predictions show a different
redshift evolution than our observations, displaying higher
$\Omega_{\rm{kin}}$ values at lower redshifts and lower values at higher
redshifts, the simba predictions for LERGs are typically within a factor of
two of our observations. The Q-LERGs within simba show a similar slope to the
observations, matching them at high redshifts but producing more heating
output by a factor of $\sim 3$ at $z\lesssim 1.5$ than observed. These results
suggest that in order to reproduce the massive galaxy population at $z=0$,
simba requires more heating from AGN jets than observed at $z\lesssim 1.5$ (as
determined for the observations using the Heckman & Best (2014) jet power
scaling relationship).

Finally, we also show the prediction for the evolution of radio-mode feedback
from Mocz et al. (2013) who develop a model for the evolution and energetic
output of supermassive black holes based on accretion only. Using the derived
Eddington luminosity distributions from this model, they predict the kinetic
luminosity density for different modes of AGN, with the LERG equivalent
population being referred to as the ‘low kinetic’ (LK) mode AGN. Their LK mode
AGN results are displayed in Fig. 2 (black dash-dotted line), which shows a
similar evolution to the observed LERGs but with a systematically higher
normalisation by up to an order of magnitude. The Mocz et al. (2013) results
are also higher than other observations and similar studies that directly use
radio LFs and a scaling relationship for kinetic powers to compute the kinetic
luminosities (e.g. Körding et al., 2008; Merloni & Heinz, 2008). Mocz et al.
(2013) propose that this systematic offset is due to faint sources missing
from the observed radio LFs and due to a fraction of bright, older sources
that may fall below the flux limit due to the decrease in radio luminosity as
a source ages. The latter effect, however, is also incorporated in the radio-
sage simulation, where the dynamics and radio luminosity evolution of sources
across their lifetimes are modelled, which shows a much lower normalisation
that is in agreement with our observations. Furthermore, as indicated in Fig.
1, the contribution of low luminosity sources to the integrated kinetic
luminosity density is not significant.

Figure 3: Top: The cumulative kinetic energy density radiated from AGN in
different forms across cosmic time (in units of $\rm{J\,Mpc^{-3}}$). The
observed kinetic luminosity densities for Fig. 2 are integrated from $z=5$ to
$z=0$, using extrapolation beyond the coverage of the dataset (see Sect. 4) to
determine the kinetic energy output from AGN jets for the three groups of
radio-AGN. The circle shows the total energy radiated by BHs using the local
BH mass density from Marconi et al. (2004) and an accretion efficiency of 10%.
The green line shows the bolometric energy radiated from quasars calculated
using the black hole mass densities (and accretion efficiency of 10%) that was
determined by Hopkins et al. (2007). The green shaded region shows the range
of possible values for kinetic energies from AGN driven winds determined by
scaling the Hopkins et al. (2007) line (see text). Our results indicate that
while AGN winds may be more important at early times, AGN jets play a more
important role in AGN feedback at $z\lesssim 2$. Bottom: The fraction of the
total (i.e. $z=0$) kinetic energy density emitted across cosmic time.

##  4 Kinetic energy density from AGN jets across cosmic history

We can study the total kinetic energy output from the radio-AGN in the form of
jets by integrating the kinetic luminosity density across cosmic time from
Fig. 2. To do this, we extrapolated the observed $\Omega_{\rm{kin}}$ to higher
redshifts (i.e. $z>2.5$) by assuming that our results follow the evolution
observed by Smolčić et al. (2017) out to $z_{\rm{max}}\sim 5$ (the approximate
redshift limit of measurements made by Smolčić et al. 2017) for the radio-
excess AGN and LERGs; for the Q-LERGs, due to the strong evolution with
redshift, we linearly extrapolated our measurements of $\Omega_{\rm{kin}}(z)$
for this population out to $z_{\rm{max}}=5$. We extrapolated our results to
lower redshifts ($z<0.5$) for all three groups of AGN by again using the
values from Smolčić et al. (2017); at these low redshifts, the LERGs will
dominate the total radio-excess AGN population (as shown by the agreement
between the Smolčić et al. (2017) and the Best et al. (2014) LERG data at
$z<0.7$), with the Q-LERGs dominating the total LERG population. The kinetic
energy density output from the AGN jets was then obtained by integrating these
extrapolated kinetic luminosity density curves with respect to time.

In Fig. 3, we show the cumulative integral (with respect to time) of these
extrapolated kinetic luminosity density curves from $z=5$ to $z=0$ in units of
$\rm{J\,Mpc^{-3}}$. The results for radio-excess AGN, LERGs, and Q-LERGs are
shown by the same colours as in Fig. 2. The total time-integrated kinetic
energy density in the form of jets from radio-AGN (and for the LERG subset
given their similarities in the kinetic heating rates) across cosmic history
is $E_{\rm{kin}}\approx 10^{50}\,\rm{J\,Mpc^{-3}}$, with $\sim$ 50 per cent of
this emitted within the past $\sim 6\,\rm{Gyr}$; this agrees well with other
measurements (e.g. Heckman & Best, 2023). For Q-LERGs, the corresponding value
is $E_{\rm{kin}}\sim 7\times 10^{49}\,\rm{J\,Mpc^{-3}}$, but they deposit most
of their total kinetic energy at late times, $z\lesssim 0.5$ (Fig. 3; bottom).

For comparison, the green line shows the bolometric radiative energy density
from AGN estimated by using an accretion efficiency of 10% and the evolution
of the black hole mass density from Hopkins et al. (2007) who determined this
by integrating their quasar luminosity function using the Soltan argument
(Soltan, 1982); the present day value of this agrees well with integrated
measurements of the local black hole mass density from Marconi et al. (2004,
black circle). As discussed by Heckman & Best (2023), only a portion of this
bolometric luminosity is in the form of kinetic energy from AGN-driven winds.
This fraction can be estimated by using the typical kinetic energy outflow
rates from the multi-phase gas medium in samples of galaxies, with different
studies finding values that range from $0.0001\%-1\%\,L_{\rm{bol}}$ (e.g.
Fiore et al., 2017; Kakkad et al., 2022; Dall’Agnol de Oliveira et al., 2021,
see Heckman & Best 2023 for a more extensive discussion). Given the large
range of values observed from different studies for the different multi-phase
mediums, here, for illustrative purposes, we take values in the range
$0.1-0.5\%\,L_{\rm{bol}}$ and multiply this by the bolometric radiative energy
density determined from Hopkins et al. (2007), giving the range indicated by
the green shaded region. We note that the upper end of this range is
consistent with the median value of $0.5\%\,L_{\rm{bol}}$ used by Heckman &
Best (2023). These results indicate that AGN winds may have a comparable or
even more important role in kinetic feedback at early cosmic times, but that
the AGN jets may have been more important for feedback since $z\sim 2$.

##  5 Conclusions

In this paper, we have presented the cosmic evolution of radio-AGN feedback
out to $z\sim 2.5$, and compared these to predictions from recent
cosmological-scale simulations. The observational data come from a sample of
radio-AGN from the LoTSS Deep Fields, which consists of 9485 radio-excess AGN,
of which 8409 are LERGs. The radio-excess AGN were identified based on sources
with a $>3\sigma$ excess in their radio emission over that expected from star-
formation, and the LERGs identified as the subset of radio-excess AGN that do
not show signs of having an accretion disk or torus structure (based on SED
fitting of photometry). This dataset has been used previously to study the
evolution of the radio luminosity functions of radio-excess AGN and LERGs by
Kondapally et al. (2022).

We studied the implications of the observed luminosity functions of radio-
excess AGN, LERGs, and LERGs hosted by quiescent galaxies (Q-LERGs) on cosmic
radio-AGN feedback. To do this, we translated the observed radio luminosities
into jet kinetic powers using a scaling relationship. This was combined with
the best-fitting models for the evolution of luminosity functions to compute
the cosmic evolution of the kinetic luminosity density of radio-AGN (a proxy
for the amount of heating provided by the radio-AGN jets) over cosmic time. We
find that feedback from LERGs dominates the radio-mode feedback with the
heating output for both the radio-excess AGN and LERGs being roughly constant
across $0.5<z\leq 2.5$, with both populations at $\Omega_{\rm{kin}}\approx
4-5\times 10^{32}\,\rm{W\,Mpc^{-3}}$. At $z\lesssim 1$, most of the kinetic
luminosity density of the total LERG population is output by Q-LERGs, beyond
this redshift, the contribution from Q-LERGs decreases by an order of
magnitude out to $z\sim 2.5$; this indicates that at higher redshifts star-
forming galaxies contribute significantly to radio-AGN feedback. Our
observations provide new measurements for the LERGs and Q-LERGs at $z\gtrsim
1$, while the measurements for the radio-AGN show good agreement with other
studies in the literature. We also find that the relatively high-luminosity
AGN dominate the overall heating output.

We then compared the observed evolution of the kinetic luminosity densities
with predictions from various recent simulations. A comparison with the sage
model finds that the Q-LERGs alone may deposit sufficient energy into their
surrounding environment to balance the radiative cooling losses, providing
evidence for a self-regulating AGN feedback cycle in these systems out to
$z\sim 2.5$. In contrast, the mechanisms driving feedback in LERGs hosted by
star-forming galaxies, which dominate at early times, remain unclear at
present. The predictions from radio-sage for the LERGs show a good match to
the observations; while the evolution of $\Omega_{\rm{kin}}(z)$ for the
Q-LERGs from radio-sage shows a similar slope to the observations, it is
offset to systematically lower values, suggesting a much higher contribution
of heating output from LERGs in radio-sage coming from those hosted in star-
forming galaxies. simba predicts $\Omega_{\rm{kin}}(z)$ for LERGs and Q-LERGs
that agree with local observations at $z=0$, but rise out to $z\sim 1$ before
declining at higher redshifts. The LERGs in simba over-predict the heating at
$z\lesssim 1.5$, and under-predict it at higher redshifts compared to the
observations. The Q-LERGs in simba also over-predict the heating output at
$z\lesssim 1.5$ but match well with observations at higher redshifts. Although
these differences are at the $0.3-0.4\,\rm{dex}$ level, these results suggest
that in order to reproduce the observed massive galaxy population, the AGN jet
feedback model employed in simba appears to require more heating output from
the AGN across cosmic history than that observed.

We integrated the kinetic luminosity density across cosmic time to obtain the
kinetic energy per unit volume output from AGN jets over cosmic history,
$E_{\rm{kin}}$. Radio-excess AGN (and LERGs) output a total of
$E_{\rm{kin}}\approx 10^{50}\,\rm{J\,Mpc^{-3}}$. The Q-LERGs are found to be
the dominant source of this kinetic energy output in the form of AGN jets,
with the Q-LERGs outputting most of their kinetic energy since $z\lesssim
0.5$. We compared this kinetic energy density from AGN jets to that output by
AGN winds, estimated using measurements of the bolometric quasar luminosity
function, to find that AGN jets dominate the total energy budget over AGN
winds at $z\lesssim 2$; this indicates that AGN jets may play a more important
role in kinetic feedback across much of cosmic history.

Our results indicate that while radio-AGN feedback plays an important role in
galaxy evolution since at least $z\sim 2.5$, and simulations are broadly able
to reflect this, black hole growth and AGN feedback models within simulations
require modifications to reproduce the observed evolution of radio-AGN
feedback in detail. To gain more insights into this, it is vital to not only
study the total amount of feedback within simulations, but to also compare the
feedback from different subsets of the radio-galaxy population, and the radio
luminosity functions and host galaxy properties of different modes of radio-
AGN; we will focus on studying this in a subsequent paper. In future,
improvements to the observations will also be offered by deep high-resolution
LOFAR imaging, enabling robust measures of source sizes out to high redshifts;
this will enable improved dynamical modelling of the radio sources and a
better-constrained jet power scaling relationship as a function of redshift.

## Acknowledgements

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. We thank the anonymous referee for their
comments which have helped improve and clarify the paper.

RK and PNB acknowledge support from the UK Science and Technology Facilities
Council (STFC) via grant ST/V000594/1. MJH and DJBS acknowledge support from
the UK STFC [ST/V000624/1]. MB and IP acknowledge support from INAF under the
SKA/CTA PRIN “FORECaST” and PRIN MAIN STREAM “SAuROS” projects, as well as
through the Large Grant 2022 funding scheme (project “MeerKAT and LOFAR Team
up: a Unique Radio Window on Galaxy/AGN co-Evolution”). K.M. has been
supported by the National Science Centre (UMO-2018/30/E/ST9/00082). LKM and
NLT acknowledge support from the Medical Research Council [MR/T042842/1]. This
research made use of Astropy, a community-developed core Python package for
astronomy (Astropy Collaboration et al., 2013, 2018) hosted at
http://www.astropy.org/, and of Matplotlib (Hunter, 2007).

## Data Availability

The luminosity function observational data used in this study come from the
analysis by Kondapally et al. (2022), which can be accessed in a machine-
readable format at the following repository:
https://github.com/rohitk-10/AGN_LF_Kondapally22. The simulations data are
based on analysis performed with the simba simulations (Davé et al., 2019;
Thomas et al., 2021), and with radio-sage (Raouf et al., 2017). The simba
snapshots can be found at http://simba.roe.ac.uk. Other data presented in the
paper are available upon request to the corresponding author.

## References

  * Aihara et al. (2019) Aihara H., et al., 2019, arXiv e-prints, p. arXiv:1905.12221
  * Anglés-Alcázar et al. (2017) Anglés-Alcázar D., Davé R., Faucher-Giguère C.-A., Özel F., Hopkins P. F., 2017, MNRAS, 464, 2840
  * Ashby et al. (2009) Ashby M. L. N., et al., 2009, ApJ, 701, 428
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Best & Heckman (2012) Best P. N., Heckman T. M., 2012, MNRAS, 421, 1569
  * Best et al. (2005) Best P. N., Kauffmann G., Heckman T. M., Brinchmann J., Charlot S., Ivezić Ž., White S. D. M., 2005, MNRAS, 362, 25
  * Best et al. (2006) Best P. N., Kaiser C. R., Heckman T. M., Kauffmann G., 2006, MNRAS, 368, L67
  * Best et al. (2014) Best P. N., Ker L. M., Simpson C., Rigby E. E., Sabater J., 2014, MNRAS, 445, 955
  * Best et al. (2023) Best P. N., et al., 2023, MNRAS, 523, 1729
  * Bîrzan et al. (2004) Bîrzan L., Rafferty D. A., McNamara B. R., Wise M. W., Nulsen P. E. J., 2004, ApJ, 607, 800
  * Bîrzan et al. (2008) Bîrzan L., McNamara B. R., Nulsen P. E. J., Carilli C. L., Wise M. W., 2008, ApJ, 686, 859
  * Boehringer et al. (1993) Boehringer H., Voges W., Fabian A. C., Edge A. C., Neumann D. M., 1993, MNRAS, 264, L25
  * Bondi (1952) Bondi H., 1952, MNRAS, 112, 195
  * Boquien et al. (2019) Boquien M., Burgarella D., Roehlly Y., Buat V., Ciesla L., Corre D., Inoue A. K., Salas H., 2019, A&A, 622, A103
  * Bower et al. (2006) 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
  * Burgarella et al. (2005) Burgarella D., Buat V., Iglesias-Páramo J., 2005, MNRAS, 360, 1413
  * Butler et al. (2019) 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 et al. (2016) Calistro Rivera G., Lusso E., Hennawi J. F., Hogg D. W., 2016, ApJ, 833, 98
  * Calistro Rivera et al. (2017) Calistro Rivera G., et al., 2017, MNRAS, 469, 3468
  * Carnall et al. (2018) Carnall A. C., McLure R. J., Dunlop J. S., Davé R., 2018, MNRAS, 480, 4379
  * Carnall et al. (2020) Carnall A. C., et al., 2020, MNRAS, 496, 695
  * Cattaneo et al. (2009) Cattaneo A., et al., 2009, Nature, 460, 213
  * Cavagnolo et al. (2010) Cavagnolo K. W., McNamara B. R., Nulsen P. E. J., Carilli C. L., Jones C., Bîrzan L., 2010, ApJ, 720, 1066
  * Chambers et al. (2016) Chambers K. C., et al., 2016, arXiv e-prints, p. arXiv:1612.05560
  * Ching et al. (2017) Ching J. H. Y., et al., 2017, MNRAS, 469, 4584
  * Croston et al. (2018) Croston J. H., Ineson J., Hardcastle M. J., 2018, MNRAS, 476, 1614
  * Croston et al. (2019) Croston J. H., et al., 2019, A&A, 622, A10
  * Croton et al. (2006) Croton D. J., et al., 2006, MNRAS, 365, 11
  * Croton et al. (2016) Croton D. J., et al., 2016, ApJS, 222, 22
  * Dall’Agnol de Oliveira et al. (2021) Dall’Agnol de Oliveira B., et al., 2021, MNRAS, 504, 3890
  * Davé et al. (2019) Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827
  * Delvecchio et al. (2017) Delvecchio I., et al., 2017, A&A, 602, A3
  * Duncan et al. (2018a) Duncan K. J., et al., 2018a, MNRAS, 473, 2655
  * Duncan et al. (2018b) Duncan K. J., Jarvis M. J., Brown M. J. I., Röttgering H. J. A., 2018b, MNRAS, 477, 5177
  * Duncan et al. (2021) Duncan K. J., et al., 2021, A&A, 648, A4
  * Eisenhardt et al. (2004) Eisenhardt P. R., et al., 2004, ApJS, 154, 48
  * Fabian (2012) Fabian A. C., 2012, ARA&A, 50, 455
  * Fabian et al. (2006) Fabian A. C., Sanders J. S., Taylor G. B., Allen S. W., Crawford C. S., Johnstone R. M., Iwasawa K., 2006, MNRAS, 366, 417
  * Fiore et al. (2017) Fiore F., et al., 2017, A&A, 601, A143
  * Gaspari et al. (2013) Gaspari M., Ruszkowski M., Oh S. P., 2013, MNRAS, 432, 3401
  * Gebhardt et al. (2000) Gebhardt K., et al., 2000, ApJ, 539, L13
  * Gendre et al. (2013) Gendre M. A., Best P. N., Wall J. V., Ker L. M., 2013, MNRAS, 430, 3086
  * Godfrey & Shabala (2016) Godfrey L. E. H., Shabala S. S., 2016, MNRAS, 456, 1172
  * Gonzalez et al. (2010) Gonzalez A. H., et al., 2010, in American Astronomical Society Meeting Abstracts #216. p. 415.13
  * Gürkan et al. (2018) Gürkan G., et al., 2018, MNRAS, 475, 3010
  * Hardcastle (2018a) Hardcastle M., 2018a, Nature Astronomy, 2, 273
  * Hardcastle (2018b) Hardcastle M. J., 2018b, MNRAS, 475, 2768
  * Hardcastle & Croston (2020) Hardcastle M. J., Croston J. H., 2020, New Astron. Rev., 88, 101539
  * Hardcastle & Krause (2013) Hardcastle M. J., Krause M. G. H., 2013, MNRAS, 430, 174
  * Hardcastle et al. (2007) Hardcastle M. J., Evans D. A., Croston J. H., 2007, MNRAS, 376, 1849
  * Hardcastle et al. (2019) Hardcastle M. J., et al., 2019, A&A, 622, A12
  * Heckman & Best (2014) Heckman T. M., Best P. N., 2014, ARA&A, 52, 589
  * Heckman & Best (2023) Heckman T. M., Best P. N., 2023, Galaxies, 11, 21
  * Hopkins & Quataert (2011) Hopkins P. F., Quataert E., 2011, MNRAS, 415, 1027
  * Hopkins et al. (2007) Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ, 654, 731
  * Hunter (2007) Hunter J. D., 2007, Computing in Science and Engineering, 9, 90
  * Ineson et al. (2017) Ineson J., Croston J. H., Hardcastle M. J., Mingo B., 2017, MNRAS, 467, 1586
  * Jannuzi & Dey (1999) Jannuzi B. T., Dey A., 1999, in Weymann R., Storrie-Lombardi L., Sawicki M., Brunner R., eds, Astronomical Society of the Pacific Conference Series Vol. 191, Photometric Redshifts and the Detection of High Redshift Galaxies. p. 111
  * Janssen et al. (2012) Janssen R. M. J., Röttgering H. J. A., Best P. N., Brinchmann J., 2012, A&A, 541, A62
  * Kakkad et al. (2022) Kakkad D., et al., 2022, MNRAS, 511, 2105
  * Kaviraj et al. (2017) Kaviraj S., et al., 2017, MNRAS, 467, 4739
  * Kenter et al. (2005) Kenter A., et al., 2005, ApJS, 161, 9
  * Kondapally et al. (2021) Kondapally R., et al., 2021, A&A, 648, A3
  * Kondapally et al. (2022) Kondapally R., et al., 2022, MNRAS, 513, 3742
  * Körding et al. (2006) Körding E. G., Fender R. P., Migliari S., 2006, MNRAS, 369, 1451
  * Körding et al. (2008) Körding E. G., Jester S., Fender R., 2008, MNRAS, 383, 277
  * Kormendy & Ho (2013) Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
  * Lawrence et al. (2007) Lawrence A., et al., 2007, MNRAS, 379, 1599
  * Lonsdale et al. (2003) Lonsdale C. J., et al., 2003, PASP, 115, 897
  * Magliocchetti (2022) Magliocchetti M., 2022, A&ARv, 30, 6
  * Marconi et al. (2004) Marconi A., Risaliti G., Gilli R., Hunt L. K., Maiolino R., Salvati M., 2004, MNRAS, 351, 169
  * Mauch & Sadler (2007) Mauch T., Sadler E. M., 2007, MNRAS, 375, 931
  * Mauduit et al. (2012) Mauduit J.-C., et al., 2012, PASP, 124, 714
  * McNamara & Nulsen (2007) McNamara B. R., Nulsen P. E. J., 2007, ARA&A, 45, 117
  * Merloni & Heinz (2008) Merloni A., Heinz S., 2008, MNRAS, 388, 1011
  * Merritt & Ferrarese (2001) Merritt D., Ferrarese L., 2001, MNRAS, 320, L30
  * Mingo et al. (2014) Mingo B., Hardcastle M. J., Croston J. H., Dicken D., Evans D. A., Morganti R., Tadhunter C., 2014, MNRAS, 440, 269
  * Mingo et al. (2022) Mingo B., et al., 2022, arXiv e-prints, p. arXiv:2201.04433
  * Mocz et al. (2013) Mocz P., Fabian A. C., Blundell K. M., 2013, MNRAS, 432, 3381
  * Morrissey et al. (2007) Morrissey P., et al., 2007, ApJS, 173, 682
  * Murphy et al. (2017) 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
  * Narayan & Yi (1994) Narayan R., Yi I., 1994, ApJ, 428, L13
  * Narayan & Yi (1995) Narayan R., Yi I., 1995, ApJ, 452, 710
  * Noll et al. (2009) Noll S., Burgarella D., Giovannoli E., Buat V., Marcillac D., Muñoz-Mateos J. C., 2009, A&A, 507, 1793
  * O’Sullivan et al. (2011) O’Sullivan E., Giacintucci S., David L. P., Gitti M., Vrtilek J. M., Raychaudhury S., Ponman T. J., 2011, ApJ, 735, 11
  * Oliver et al. (2012) Oliver S. J., et al., 2012, MNRAS, 424, 1614
  * Pacifici et al. (2016) Pacifici C., et al., 2016, ApJ, 832, 79
  * Pracy et al. (2016) Pracy M. B., et al., 2016, MNRAS, 460, 2
  * Raouf et al. (2017) Raouf M., Shabala S. S., Croton D. J., Khosroshahi H. G., Bernyk M., 2017, MNRAS, 471, 658
  * Sabater et al. (2013) Sabater J., Best P. N., Argudo-Fernández M., 2013, MNRAS, 430, 638
  * Sabater et al. (2021) Sabater J., et al., 2021, A&A, 648, A2
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Shabala & Alexander (2009) Shabala S., Alexander P., 2009, ApJ, 699, 525
  * Shabala & Godfrey (2013) Shabala S. S., Godfrey L. E. H., 2013, ApJ, 769, 129
  * Shabala et al. (2012) Shabala S. S., Santoso J. S., Godfrey L. E. H., 2012, ApJ, 756, 161
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Sharma et al. (2012) Sharma P., McCourt M., Quataert E., Parrish I. J., 2012, MNRAS, 420, 3174
  * Shimwell et al. (2019) Shimwell T. W., et al., 2019, A&A, 622, A1
  * Smith et al. (2021) Smith D. J. B., et al., 2021, A&A, 648, A6
  * Smolčić (2009) Smolčić V., 2009, ApJ, 699, L43
  * Smolčić et al. (2017) Smolčić V., et al., 2017, A&A, 602, A6
  * Soltan (1982) Soltan A., 1982, MNRAS, 200, 115
  * Springel et al. (2018) Springel V., et al., 2018, MNRAS, 475, 676
  * Sutherland & Saunders (1992) Sutherland W., Saunders W., 1992, MNRAS, 259, 413
  * Sweijen et al. (2022) Sweijen F., et al., 2022, Nature Astronomy, 6, 350
  * Tadhunter (2016) Tadhunter C., 2016, A&ARv, 24, 10
  * Tasse et al. (2008) Tasse C., Best P. N., Röttgering H., Le Borgne D., 2008, A&A, 490, 893
  * Tasse et al. (2021) Tasse C., et al., 2021, A&A, 648, A1
  * Thomas et al. (2021) Thomas N., Davé R., Jarvis M. J., Anglés-Alcázar D., 2021, MNRAS, 503, 3492
  * Timmerman et al. (2022) Timmerman R., van Weeren R. J., Botteon A., Röttgering H. J. A., McNamara B. R., Sweijen F., Bîrzan L., Morabito L. K., 2022, A&A, 668, A65
  * Turner & Shabala (2015) Turner R. J., Shabala S. S., 2015, ApJ, 806, 59
  * Whittam et al. (2018) Whittam I. H., Prescott M., McAlpine K., Jarvis M. J., Heywood I., 2018, MNRAS, 480, 358
  * Whittam et al. (2022) Whittam I. H., et al., 2022, MNRAS, 516, 245
  * Williams et al. (2018) Williams W. L., et al., 2018, MNRAS, 475, 3429
  * Williams et al. (2019) Williams W. L., et al., 2019, A&A, 622, A2
  * Willott et al. (1999) Willott C. J., Rawlings S., Blundell K. M., Lacy M., 1999, MNRAS, 309, 1017
  * da Cunha et al. (2008) da Cunha E., Charlot S., Elbaz D., 2008, MNRAS, 388, 1595
  * de Ruiter et al. (1977) de Ruiter H. R., Willis A. G., Arp H. C., 1977, Astronomy and Astrophysics Supplement Series, 28, 211
  * van Haarlem et al. (2013) van Haarlem M. P., et al., 2013, A&A, 556

##  Appendix A Radio-AGN luminosity functions used for estimating kinetic
heating rates

Table 1: Results from modelling the evolution of three populations of AGN LFs, plotted in Fig. 4. The radio-excess AGN and LERG LFs are modelled by fitting a luminosity and density evolution (LDE) model at each redshift (independently) based on the local radio-AGN LF of Mauch & Sadler (2007). For the Q-LERGs, the $\rho^{\star}(z)$ at each redshift is fixed based on the available quiescent host galaxies at that redshift, and the $L_{\star}(z)$ is then allowed to evolve freely at each redshift based on a broken power law fit to the observed LF at $0.5<z\leq 1$ (with slopes of $\beta=-2.88$ and $\gamma=-0.55$; see Kondapally et al. 2022.) $z$ | Radio-excess AGN | LERGs | Quiescent LERGs
---|---|---|---
| $\log_{10}\rho^{\star}(z)$ | $\log_{10}L_{\star}(z)$ | $\chi^{2}_{\nu}$ | $\log_{10}\rho^{\star}(z)$ | $\log_{10}L_{\star}(z)$ | $\chi^{2}_{\nu}$ | $\log_{10}\rho^{\star}(z)$ | $\log_{10}L_{\star}(z)$ |  $\chi^{2}_{\nu}$
${0.5}<z\leq{1.0}$ | $-5.53_{-0.08}^{+0.07}$ | $26.31_{-0.12}^{+0.14}$ | 7.22 | $-5.39_{-0.07}^{+0.07}$ | $25.98_{-0.11}^{+0.12}$ | 5.59 | $-6.37_{-0.24}^{+0.17}$ | $27.03_{-0.22}^{+0.40}$ | 3.57
${1.0}<z\leq{1.5}$ | $-5.89_{-0.09}^{+0.08}$ | $26.77_{-0.14}^{+0.16}$ | 5.83 | $-5.63_{-0.09}^{+0.08}$ | $26.24_{-0.13}^{+0.14}$ | 5.56 | -6.92 | $27.33_{-0.04}^{+0.04}$ | 3.39
${1.5}<z\leq{2.0}$ | $-5.71_{-0.08}^{+0.07}$ | $26.76_{-0.12}^{+0.13}$ | 7.98 | $-5.52_{-0.08}^{+0.08}$ | $26.35_{-0.11}^{+0.12}$ | 6.28 | -7.28 | $27.91_{-0.06}^{+0.06}$ | 2.28
${2.0}<z\leq{2.5}$ | $-5.41_{-0.09}^{+0.09}$ | $26.09_{-0.13}^{+0.13}$ | 4.98 | $-5.13_{-0.11}^{+0.12}$ | $25.61_{-0.14}^{+0.14}$ | 3.45 | -7.61 | $27.93_{-0.12}^{+0.12}$ | 5.89
Figure 4: Cosmic evolution of AGN luminosity functions from LoTSS-Deep across
$0.5<z\leq 2.5$. The LFs for radio-excess AGN (blue open circles), LERGs
(black circles), and quiescent LERGs (Q-LERGs; red triangles) are shown. The
corresponding solid lines at each redshift show the best-fitting LFs for the
three populations, in which the lower-$z$ LF shape (see text) is evolved in
both density and luminosity (the shaded regions represent the uncertainty on
this fit; see text).

The observed radio LFs used in this study, for the radio-excess AGN, LERGs,
and Q-LERGs are taken from the analysis by Kondapally et al. (2022). In
summary, the LFs were calculated using the $1/V_{\rm{max}}$ method, which
accounted for the radio flux-density incompleteness and noise variation across
the field. The radio flux-density completeness corrections were estimated on a
field-by-field basis by performing a large suite of simulations involving the
injection and recovery of mock sources onto the image plane; we refer to
Kondapally et al. (2022) for the details of this process. The 150 MHz radio
luminosities were computed assuming a spectral index $\alpha=-0.7$ (e.g.
Calistro Rivera et al., 2017; Murphy et al., 2017).

Kondapally et al. (2022) computed the LERG and Q-LERG LFs over the redshift
range: $0.5<z\leq 1$, $1<z\leq 1.5$, $1.5<z\leq 2$, and $2<z\leq 2.5$. They
computed the radio-excess AGN LFs over different redshift bins to readily
compare with available observational measurements for that population. Here,
we recompute the radio-excess AGN LFs over the same redshift bins as that for
the LERGs, which are shown in Fig. 4. Both the radio-excess AGN LFs and LERG
LFs show excellent agreement with the observed LFs from other studies (Best et
al., 2014; Smolčić et al., 2017; Butler et al., 2019) over the full redshift
range studied as shown by Kondapally et al. (2022). From Fig. 4, it is also
evident that across all radio luminosities and redshifts studied, the LoTSS-
Deep radio-excess AGN sample is dominated by LERGs; the HERGs form a minority
of the radio-excess AGN population. As noted by Kondapally et al. (2022), we
observe little-to-no evolution in the LERG LFs across $0.5<z\leq 2.5$, however
when split by host galaxy type, we see a strong decline in the space densities
of the quiescent-LERGs (Q-LERGs, red triangles) with increasing redshift.
Kondapally et al. (2022) showed that this decline is in line with the decrease
in the availability of quiescent hosts.

The evolution of the radio-AGN LFs were modelled at each redshift as a broken
power law of the form

$\rho(L,z)=\frac{\rho^{\star}(z)}{\left(L^{\star}(z)/L\right)^{\beta}+\left(L^{\star}(z)/L\right)^{\gamma}},$
(5)

where $L^{\star}(z)$ is the characteristic luminosity at redshift $z$,
$\rho^{\star}$ is the characteristic space density at redshift $z$, and
$\beta$ and $\gamma$ are the bright and faint end slopes, respectively. For
both the radio-excess AGN and the LERGs, a broken power-law fit at each
redshift was performed by fixing the bright and faint end slopes,
$\beta=-1.27$ and $\gamma=-0.49$, respectively, following the broken power-law
fit to the local radio-AGN LF by Mauch & Sadler (2007); we then fitted for a
combined luminosity and density evolution (LDE) model at each redshift, as
detailed by Kondapally et al. (2022). The resulting fitted $\rho^{\star}(z)$
and $L^{\star}(z)$ for LERGs (taken from Kondapally et al. 2022) and radio-
excess AGN are shown in Table 1. For the Q-LERGs, we directly fitted a broken
power law to the $0.5<z\leq 1.0$ LF to derive empirical power law slopes.
Using these slopes, we fitted the evolution at higher redshifts by fixing
$\rho^{\star}(z)$ at each redshift bin in accordance with the evolution of the
available quiescent host-galaxies and only fitting for $L^{\star}(z)$ (see
Kondapally et al. 2022 for details); the parameters for the best-fitting
Q-LERG LFs are also shown in Table 1 for completeness.

##  Appendix B Comparison of jet power inference methods

Table 2: Different literature scaling relations for converting monochromatic radio luminosities ($L_{1.4\,\rm{GHz}}$) to jet kinetic powers ($L_{\rm{kin}}$) that are compared. For the Willott et al. (1999) and Heckman & Best (2014), the relations listed have explicitly assumed $f_{\rm{W}}=15$ and $f_{\rm{cav}}=4$, as used in this study; relations for other factors can be calculated using the additional term in each case. The Shabala & Godfrey (2013) relation includes a dependence on redshift and source size. Reference | Scaling Relation
---|---
Willott et al. (1999) |  $\log_{10}(L_{\rm{kin}}/\rm{W})=0.86\log_{10}\left(\frac{L_{\rm{1.4\,GHz}}}{10^{25}\,\rm{W\,Hz^{-1}}}\right)+37.37+\log_{10}\left(\frac{\mathit{f}_{\rm{W}}}{15}\right)$
Körding et al. (2008) | $\log_{10}(L_{\rm{kin}}/\rm{W})=0.71\log_{10}\left(\frac{L_{\rm{1.4\,GHz}}}{10^{25}\,\rm{W\,Hz^{-1}}}\right)+37.72$
Cavagnolo et al. (2010) | $\log_{10}(L_{\rm{kin}}/\rm{W})=0.75\log_{10}\left(\frac{L_{\rm{1.4\,GHz}}}{10^{25}\,\rm{W\,Hz^{-1}}}\right)+38.06$
O’Sullivan et al. (2011) | $\log_{10}(L_{\rm{kin}}/\rm{W})=0.63\log_{10}\left(\frac{L_{\rm{1.4\,GHz}}}{10^{25}\,\rm{W\,Hz^{-1}}}\right)+37.76$
Shabala & Godfrey (2013) | $\log_{10}(L_{\rm{kin}}/\rm{W})=0.8\log_{10}\left(\frac{L_{\rm{1.4\,GHz}}}{10^{25}\,\rm{W\,Hz^{-1}}}\right)+35.12+\log_{10}(1+\textit{z})+0.58\log_{10}\left(\frac{\textit{D}}{\rm{kpc}}\right)$
Heckman & Best (2014) | $\log_{10}(L_{\rm{kin}}/\rm{W})=0.68\log_{10}\left(\frac{L_{\rm{1.4\,GHz}}}{10^{25}\,\rm{W\,Hz^{-1}}}\right)+37.45+\log_{10}\left(\frac{\mathit{f}_{\rm{cav}}}{4}\right)$
Ineson et al. (2017) | $\log_{10}(L_{\rm{kin}}/\rm{W})=0.89\log_{10}\left(\frac{L_{\rm{1.4\,GHz}}}{10^{25}\,\rm{W\,Hz^{-1}}}\right)+36.05$
Figure 5: Comparison of different scaling relations between monochromatic
radio luminosity (1.4 GHz) and the jet powers in the literature, as
highlighted in Table 2 (Willott et al., 1999; Körding et al., 2006, 2008;
Cavagnolo et al., 2010; O’Sullivan et al., 2011; Heckman & Best, 2014; Godfrey
& Shabala, 2016; Ineson et al., 2017; Raouf et al., 2017). The 150 MHz radio
luminosities (assuming $\alpha=-0.7$) are also shown. The Heckman & Best
(2014) relation is shown for $f_{\rm{cav}}=4$, which shows good agreement in
its normalisation with that of Willott et al. (1999) when using
$f_{\rm{W}}=15$, and with the O’Sullivan et al. (2011) relation. For
illustrative purposes, the Shabala & Godfrey (2013) relation is shown for a
source of size $D=100\,\rm{kpc}$ and $z=2$. The median relation for the jet
model from radio-sage is taken from Raouf et al. (2017) and extrapolated
linearly beyond $\sim 10^{26}\,\rm{W\,Hz^{-1}}$, which shows a higher
normalisation than the other relations considered.

There exist many scaling relations for translating monochromatic radio
luminosities ($L_{\rm{1.4\,GHz}}$) into jet kinetic powers ($L_{\rm{kin}}$)
derived using different methods (e.g. Willott et al., 1999; Bîrzan et al.,
2004; Cavagnolo et al., 2010; O’Sullivan et al., 2011; Shabala & Godfrey,
2013; Godfrey & Shabala, 2016), each with their own set of uncertainties on
the calibration (see Best et al., 2014; Godfrey & Shabala, 2016; Hardcastle &
Croston, 2020). In this section, we compare the commonly used scaling
relations and assess their impacts on the cosmic kinetic luminosity density
results presented in Section 3.

Fig. 5 shows the different scaling relations between 1.4 GHz radio luminosity
and jet kinetic powers. These scaling relations are also listed in Table 2.
Throughout the main sections of this paper, we have used the scaling relation
from Heckman & Best (2014); this method is based on determining the work done
by the radio jets in inflating cavities using X-ray observations (e.g. Bîrzan
et al., 2004, 2008; Cavagnolo et al., 2010). O’Sullivan et al. (2011) and
Heckman & Best (2014) extended the analysis to fainter radio luminosities,
albeit for a small number of systems in total at low redshift ($z<0.04$), and
found slopes $\sim 0.6-0.8$ with a typical scatter of $\sigma\approx
0.7\,\rm{dex}$.

Another commonly used relation in the literature is that based on theoretical
arguments; Willott et al. (1999) used the minimum energy argument to estimate
the energy stored within radio lobes using estimates of radio source lifetimes
and the efficiency with which the jet power is converted to synchrotron
emission as the radio lobes expand, as described in Section 3. They combined
all the uncertainties (departure from minimum energy, jet plasma composition,
and energy stored in radiating particles) in their calibration into a single
parameter, $f_{\rm{W}}$ which is expected to lie in the range
$f_{\rm{W}}=1-20$. Willott et al. (1999) find a steeper slope ($\sim$ 0.9)
than the cavity power based methods, but for the typical value of
$f_{\rm{W}}=15$ (e.g. Heckman & Best, 2014; Smolčić et al., 2017), this
relation produces a similar normalisation to that of the cavity based
estimates at typical luminosities of $L_{\rm{1.4\,GHz}}\sim
10^{25}\,\rm{W\,Hz^{-1}}$ (Heckman & Best, 2014). More recently, Ineson et al.
(2017) estimated the jet powers using internal energies of the lobes for a
sample of bright Fanaroff-Riley II (FR-II) sources
($L_{\rm{1.4\,\rm{GHz}}}\gtrsim 10^{24}\,\rm{W\,Hz^{-1}}$), finding a similar
slope to the Willott et al. (1999) relation with a normalisation that is
consistent with a lower $f_{\rm{W}}$ value of $f_{\rm{W}}\sim 4$ (see Table
2).

Shabala & Godfrey (2013) used a sample of FR-II sources to derive a jet-power
scaling relationship based on both the radio luminosity and the source-size
(which is taken to be a proxy for the age of the radio source). In Fig. 5, we
show the Shabala & Godfrey (2013) relation assuming a source-size of 100 kpc
at $z=2$; this relation predicts lower jet powers than the other relations,
but is found to be consistent with the Willott et al. (1999) relation when
assuming $f_{\rm{W}}\sim 4$.

Also shown in Fig. 5 is the median relationship found for the jet model
employed within radio-sage (Raouf et al., 2017); this has a notably steeper
slope than the Heckman & Best (2014) and O’Sullivan et al. (2011) relation,
similar to that of the Willott et al. (1999) relation, however with a higher
normalisation. Therefore, radio-sage AGN on average would be expected to
produce higher kinetic powers at high luminosities and hence result in higher
kinetic luminosity density across redshift (see Fig. 6). The Körding et al.
(2008) relation for jet powers based on observations of X-ray binaries is also
shown; this relation has a similar slope to that of Heckman & Best (2014),
however with a higher normalisation (by $\sim 0.3\,\rm{dex}$) such that it
predicts higher kinetic powers than the Heckman & Best (2014) scaling relation
used in this study over all radio luminosities.

Figure 6: The cosmic evolution of the kinetic luminosity density for radio-
excess AGN (blue) and Q-LERGs (red) assuming different scaling relations for
converting monochromatic radio luminosities into jet powers in Fig. 5. The
darker shades of colours correspond to measurements based on cavity power
method and the lighted shades correspond to methods based on the synchrotron
emission.

In Fig. 6, we show the cosmic evolution of the kinetic luminosity density of
AGN within LoTSS-Deep (see Sect. 3.2) when assuming different scaling
relations for converting radio luminosities to kinetic jet powers. The solid
lines with symbols (the same as in Fig. 2) show the measurements based on the
Heckman & Best (2014) relation that is used in this study (see Section 3) for
radio-excess AGN (blue) and Q-LERGs (red); the ‘All-LERGs’ are not shown on
this plot for clarity. Other line styles show how the kinetic luminosity
density changes for these two populations when different scaling relations are
used for the two AGN populations.

Measurements based on cavity power scaling relations (i.e. from Cavagnolo et
al., 2010; O’Sullivan et al., 2011; Heckman & Best, 2014) are shown by darker
colour shades compared to estimates based on radio source modelling by Willott
et al. (1999); Körding et al. (2008); Ineson et al. (2017); Raouf et al.
(2017). The cavity power based estimates from O’Sullivan et al. (2011) and
Heckman & Best (2014) agree well with each other and with the results from
Willott et al. (1999, for the value of $f_{\rm{W}}=15$). The results from
Cavagnolo et al. (2010) predict systematically higher kinetic luminosity
densities than the other cavity based methods (by a factor of $\sim$ 2) due to
the steeper slope; the Cavagnolo et al. (2010) sample of 21 systems covered
relatively low radio luminosities ($L_{\rm{1.4\,GHz}}\lesssim
10^{24}\,\rm{W\,Hz^{-1}}$) compared to both O’Sullivan et al. (2011) and
Heckman & Best (2014) resulting in poorer constraints and hence extrapolation
at high luminosities. This can have an important effect as the heating output
peaks near the break in the LFs, i.e. at $L_{\rm{150\,MHz}}\sim
10^{26}-10^{27}\,\rm{W\,Hz^{-1}}$ (see Fig. 1). Both the Körding et al. (2008)
and Raouf et al. (2017; radio-sage) models predict systematically higher
$\Omega_{\rm{kin}}(z)$ than the Heckman & Best (2014) relation used in this
study. The shape of the cosmic evolution seen for the various relations used
is similar, with the different relations largely introducing a difference in
overall normalisation of the kinetic luminosity densities.