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

R Kondapally, P.N Best, M Raouf, N.L Thomas, R Davé, S.S Shabala, H.J.A

Röttgering, M.J Hardcastle, M Bonato, R.K Cochrane, et al.

► To cite this version:

R Kondapally, P.N Best, M Raouf, N.L Thomas, R Davé, et al. Cosmic evolution of radio-AGN feedback: confronting models with data. Monthly Notices of the Royal Astronomical Society, 2023, 523 (4), pp.5292-5305. 10.1093/mnras/stad1813. hal-04138132

HAL Id: hal-04138132

https://hal.science/hal-04138132

Submitted on 21 Apr 2024

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
# Cosmic evolution of radio-AGN feedback: confronting models with data

R Kondapally, P.N Best, M Raouf, N.L Thomas, R Davé, S.S Shabala, H.J.A

Röttgering, M.J Hardcastle, M Bonato, R.K Cochrane, et al.

► To cite this version:

R Kondapally, P.N Best, M Raouf, N.L Thomas, R Davé, et al. Cosmic evolution of radio-AGN feedback: confronting models with data. Monthly Notices of the Royal Astronomical Society, 2023, 523 (4), pp.5292-5305. 10.1093/mnras/stad1813. hal-04138132

HAL Id: hal-04138132

https://hal.science/hal-04138132v1

Submitted on 21 Apr 2024

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
# Cosmic evolution of radio-AGN feedback: confronting models with data

Rohit Kondapally<sup>1,★</sup> Philip N. Best,<sup>1</sup> Mojtaba Raouf<sup>1,2</sup> Nicole L. Thomas<sup>1,3,4</sup> Romeel Davé<sup>1,1</sup>

Stanislav S. Shabala<sup>1,5</sup> Huub J. A. Röttgering,<sup>2</sup> Martin J. Hardcastle<sup>1,6</sup> Matteo Bonato,<sup>7</sup>

Rachel K. Cochrane<sup>1,8</sup> Katarzyna Małek,<sup>9,10</sup> Leah. K. Morabito<sup>1,3,4</sup> Isabella Prandoni<sup>1,7</sup> and
Daniel. J. B. Smith<sup>1,6</sup>

Accepted 2023 June 14. Received 2023 June 3; in original form 2023 April 12

## 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-active galactic nuclei (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_{\text{kin}}$ (i.e. the volume-averaged heating output). Compared to all radio-AGN, low-excitation radio galaxies dominate the feedback activity out to $z \sim 2.5$, with both these populations showing a constant heating output of $\Omega_{\text{kin}} \approx (4-5) \times 10^{32} \text{ 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_{\text{kin}}$, although their absolute normalization 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} \text{ 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.

**Key words:** galaxies: active – galaxies: evolution – galaxies: jets – radio continuum: galaxies.

## 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 disc (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 (LF; 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).

<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>Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the Netherlands</footnote>

<footnote n="3"><sup>3</sup>Centre for Extragalactic Astronomy, Department of Physics, Durham University, Durham DH1 3LE, UK</footnote>

<footnote n="4"><sup>4</sup>Institute for Computational Cosmology, Department of Physics, University of Durham, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="5"><sup>5</sup>School of Natural Sciences, Private Bag 37, University of Tasmania, Hobart, TAS 7001, Australia</footnote>

<footnote n="6"><sup>6</sup>Centre for Astrophysics Research, Department of Physics, Astronomy and Mathematics, University of Hertfordshire, College Lane, Hatfield AL10 9AB, UK</footnote>

<footnote n="7"><sup>7</sup>INAF-IRA, Via P. Gobetti 101, I-40129 Bologna, Italy</footnote>

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

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

<footnote n="10"><sup>10</sup>Aix Marseille Univ. CNRS, CNES, LAM, F-13284 Marseille, 13013, France</footnote>

<footnote>* E-mail: rohitkondapally@gmail.com</footnote>
# Cosmic evolution of radio-AGN feedback: confronting models with data

Rohit Kondapally<sup>1,★</sup> Philip N. Best,<sup>1</sup> Mojtaba Raouf<sup>1,2</sup> Nicole L. Thomas<sup>1,3,4</sup> Romeel Davé<sup>1,5</sup>

Stanislav S. Shabala<sup>1,5</sup> Huub J. A. Röttgering,<sup>2</sup> Martin J. Hardcastle<sup>1,6</sup> Matteo Bonato,<sup>7</sup>

Rachel K. Cochrane<sup>1,8</sup> Katarzyna Małek,<sup>9,10</sup> Leah. K. Morabito<sup>1,3,4</sup> Isabella Prandoni<sup>1,7</sup> and

Daniel. J. B. Smith<sup>1,6</sup>

Accepted 2023 June 14. Received 2023 June 3; in original form 2023 April 12

## 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-active galactic nuclei (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_{\text{kin}}$ (i.e. the volume-averaged heating output). Compared to all radio-AGN, low-excitation radio galaxies dominate the feedback activity out to $z \sim 2.5$, with both these populations showing a constant heating output of $\Omega_{\text{kin}} \approx (4-5) \times 10^{32} \text{ 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_{\text{kin}}$, although their absolute normalization 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} \text{ 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.

**Key words:** galaxies: active – galaxies: evolution – galaxies: jets – radio continuum: galaxies.

## 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 disc (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 (LF; 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).

<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>Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the Netherlands</footnote>

<footnote n="3"><sup>3</sup>Centre for Extragalactic Astronomy, Department of Physics, Durham University, Durham DH1 3LE, UK</footnote>

<footnote n="4"><sup>4</sup>Institute for Computational Cosmology, Department of Physics, University of Durham, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="5"><sup>5</sup>School of Natural Sciences, Private Bag 37, University of Tasmania, Hobart, TAS 7001, Australia</footnote>

<footnote n="6"><sup>6</sup>Centre for Astrophysics Research, Department of Physics, Astronomy and Mathematics, University of Hertfordshire, College Lane, Hatfield AL10 9AB, UK</footnote>

<footnote n="7"><sup>7</sup>INAF-IRA, Via P. Gobetti 101, I-40129 Bologna, Italy</footnote>

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

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

<footnote n="10"><sup>10</sup>Aix Marseille Univ. CNRS, CNES, LAM, F-13284 Marseille, 13013, France</footnote>

<footnote>* E-mail: rohitkondapally@gmail.com</footnote>
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 characterization 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, Evans & Croston 2007; Tasse et al. 2008; Smolčić 2009; Best & Heckman 2012; Janssen et al. 2012; Sabater, Best & Argudo-Fernández 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 disc 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, Ruszkowski & Oh 2013), and hence lack a radiatively-efficient accretion disc. However, recent deeper radio continuum surveys, probing over an order of magnitude fainter radio luminosities 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; Kondapally et al. 2022; Mingo et al. 2022; Whittam et al. 2022).

The radio LFs 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 LFs of LERGs across $0.5 < z \le 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 LFs out to $z \sim 2.5$ using the LOFAR Two-metre Sky Survey Deep Fields Data Release 1 (LoTSS-Deep DR1; Duncan et al. 2021; Kondapally et al. 2021; Sabater et al. 2021; Tasse 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 \text{ deg}^2$ (where there is overlap with high-quality multiwavelength 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 LF than deep narrow-area surveys (e.g. Smolčić et al. 2017). Kondapally et al. (2022) found that the LERG LFs showed mild evolution across $0.5 < z \le 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$ Gyr 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 LFs 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 multiwavelength data sets 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_m = 0.3$, $\Omega_\Lambda = 0.7$, and $H_0 = 70 \text{ km s}^{-1} \text{ 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 multiwavelength data

The radio catalogues are taken from the LoTSS Deep Fields DR1 (Sabater et al. 2021; Tasse et al. 2021). The radio observations consist
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 characterization 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, Evans & Croston 2007; Tasse et al. 2008; Smolčić 2009; Best & Heckman 2012; Janssen et al. 2012; Sabater, Best & Argudo-Fernández 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 disc 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, Ruszkowski & Oh 2013), and hence lack a radiatively-efficient accretion disc. However, recent deeper radio continuum surveys, probing over an order of magnitude fainter radio luminosities 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; Kondapally et al. 2022; Mingo et al. 2022; Whittam et al. 2022).

The radio LFs 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 LFs of LERGs across $0.5 < z \le 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 LFs out to $z \sim 2.5$ using the LOFAR Two-metre Sky Survey Deep Fields Data Release 1 (LoTSS-Deep DR1; Duncan et al. 2021; Kondapally et al. 2021; Sabater et al. 2021; Tasse 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 \text{ deg}^2$ (where there is overlap with high-quality multiwavelength 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 LF than deep narrow-area surveys (e.g. Smolčić et al. 2017). Kondapally et al. (2022) found that the LERG LFs showed mild evolution across $0.5 < z \le 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$ Gyr 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 LFs 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 multiwavelength data sets 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_m = 0.3$, $\Omega_\Lambda = 0.7$, and $H_0 = 70 \text{ km s}^{-1} \text{ 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 multiwavelength data

The radio catalogues are taken from the LoTSS Deep Fields DR1 (Sabater et al. 2021; Tasse 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 µJy beam⁻¹ at 150 MHz in the centres of each field, respectively. The LoTSS-Deep data set 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 Sabater et al. (2021); Tasse et al. (2021).

The three LoTSS-Deep fields were chosen due to the availability of deep, wide-field, multiwavelength imaging (see Kondapally et al. 2021, and references therein), making these fields ideal for characterizing 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, K, data in Boötes from Gonzalez et al. 2010), mid-infrared (Spitzer IRAC and MIPS surveys; Lonsdale et al. 2003; Eisenhardt et al. 2004; Ashby et al. 2009; Mauduit et al. 2012), 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, Willis & Arp 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 (~68 deg² in each field), the host galaxy identification process was limited to the areas in each field with the best available multiwavelength 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 ~25 deg² over the three fields (Kondapally et al. 2021). Photometric redshifts for the radio sources (and the underlying multiwavelength catalogues) were generated by Duncan et al. (2021) using a hybrid machine learning and template-based approach optimized for deep radio-continuum surveys (see Duncan et al. 2018a, b).

## 2.2 SED fitting and source classification

The multiwavelength 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, Buat & Iglesias-Páramo 2005; Noll et al. 2009; Boquien et al. 2019), and Multi-wavelength Analysis of Galaxy Physical Properties (MAGPHYS; da Cunha, Charlot & Elbaz 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 data set. In summary, Best et al. (2023)

used fits to the AGN accretion disc 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 dex (~3σ) in radio luminosity, over the level expected from star formation processes alone based on the consensus SFRs and SFR-radio luminosity relation for SFGs (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 LFs 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 that satisfy sSFR < 0.2/t<sub>H(z)</sub> are defined as quiescent galaxies, where t<sub>H(z)</sub> 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, 2020). In total, across the redshifts 0.5 < z ≤ 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 misclassified as LERGs), especially in the absence of X-ray data (which is not available in ELAIS-N1 and Lockman Hole). To test this, we utilized the X-Boötes Chandra survey (Kenter et al. 2005); this data set 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 ~5 per cent 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 LFs 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 LFs and results in this paper.
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 µJy beam⁻¹ at 150 MHz in the centres of each field, respectively. The LoTSS-Deep data set 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 Sabater et al. (2021); Tasse et al. (2021).

The three LoTSS-Deep fields were chosen due to the availability of deep, wide-field, multiwavelength imaging (see Kondapally et al. 2021, and references therein), making these fields ideal for characterizing 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, K, data in Boötes from Gonzalez et al. 2010), mid-infrared (Spitzer IRAC and MIPS surveys; Lonsdale et al. 2003; Eisenhardt et al. 2004; Ashby et al. 2009; Mauduit et al. 2012), 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, Willis & Arp 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 (~68 deg² in each field), the host galaxy identification process was limited to the areas in each field with the best available multiwavelength 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 ~25 deg² over the three fields (Kondapally et al. 2021). Photometric redshifts for the radio sources (and the underlying multiwavelength catalogues) were generated by Duncan et al. (2021) using a hybrid machine learning and template-based approach optimized for deep radio-continuum surveys (see Duncan et al. 2018a, b).

## 2.2 SED fitting and source classification

The multiwavelength 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, Buat & Iglesias-Páramo 2005; Noll et al. 2009; Boquien et al. 2019), and Multi-wavelength Analysis of Galaxy Physical Properties (MAGPHYS; da Cunha, Charlot & Elbaz 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 data set. In summary, Best et al. (2023)

used fits to the AGN accretion disc 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 dex (~3σ) in radio luminosity, over the level expected from star formation processes alone based on the consensus SFRs and SFR-radio luminosity relation for SFGs (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 LFs 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 that satisfy sSFR < 0.2/t<sub>H(z)</sub> are defined as quiescent galaxies, where t<sub>H(z)</sub> 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, 2020). In total, across the redshifts 0.5 < z ≤ 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 misclassified as LERGs), especially in the absence of X-ray data (which is not available in ELAIS-N1 and Lockman Hole). To test this, we utilized the X-Boötes Chandra survey (Kenter et al. 2005); this data set 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 ~5 per cent 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 LFs 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 LFs and results in this paper.
## 3 EVOLUTION OF RADIO-AGN FEEDBACK OVER COSMIC TIME

The observed 150-MHz radio LFs for the radio-AGN in LoTSS-Deep were presented by Kondapally et al. (2022). The LFs were calculated using the $1/V_{max}$ method and the evolution of the radio-excess AGN, LERGs, Q-LERGS (i.e. LERGs hosted in quiescent galaxies) was characterized. 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_W$, with typical values in the range of $f_W \sim 1-20$. Converted into 1.4-GHz luminosity (Heckman & Best 2014), this relation is given as

$$L_{\text{kin, sync}} = 4 \times 10^{35} (f_W)^{3/2} (L_{1.4\text{GHz}}/10^{25} \text{ W Hz}^{-1})^{0.86} \text{ W}, \quad (1)$$

where $L_{1.4\text{GHz}}$ is the 1.4-GHz radio luminosity, and $f_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_{cav} = f_{cav} pV$, where $f_{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-fitting relation, largely based on the results from Bîrzan et al. (2008) and Cavagnolo et al. (2010), given as

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

The normalization of this relation is found to be in good agreement with that of the Willott et al. (1999) relation when using $f_W = 15$ and $f_{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, Santoso & Godfrey 2012; Hardcastle & Krause 2013; Shabala & Godfrey 2013; Godfrey & Shabala 2016; Croston, Ineson & Hardcastle 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.

An additional issue is that the above scaling relation has been determined from characterization 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 preparation).

In the analysis that follows, we use the Heckman & Best (2014) relation (equation 2) with $f_{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. B2), 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 normalization by around 0.3 dex. This normalization uncertainty needs to be borne in mind when comparing the output of the simulations with the observational data.

### 3.2 Cosmic evolution of the kinetic luminosity density

The evolving radio LFs 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_{1.4\text{GHz}}, z) = L_{\text{kin}}(L_{1.4\text{GHz}}) \times \rho(L_{1.4\text{GHz}}, z) \text{ W Mpc}^{-3} \text{ logL}^{-1}. \quad (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
# 3 EVOLUTION OF RADIO-AGN FEEDBACK OVER COSMIC TIME

The observed 150-MHz radio LFs for the radio-AGN in LoTSS-Deep were presented by Kondapally et al. (2022). The LFs were calculated using the $1/V_{max}$ method and the evolution of the radio-excess AGN, LERGs, Q-LERGS (i.e. LERGs hosted in quiescent galaxies) was characterized. 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_W$, with typical values in the range of $f_W \sim 1-20$. Converted into 1.4-GHz luminosity (Heckman & Best 2014), this relation is given as

$$L_{\text{kin, sync}} = 4 \times 10^{35} (f_W)^{3/2} (L_{1.4\text{GHz}}/10^{25} \text{ W Hz}^{-1})^{0.86} \text{ W}, \quad (1)$$

where $L_{1.4\text{GHz}}$ is the 1.4-GHz radio luminosity, and $f_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_{cav} = f_{cav} pV$, where $f_{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-fitting relation, largely based on the results from Bîrzan et al. (2008) and Cavagnolo et al. (2010), given as

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

The normalization of this relation is found to be in good agreement with that of the Willott et al. (1999) relation when using $f_W = 15$ and $f_{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, Santoso & Godfrey 2012; Hardcastle & Krause 2013; Shabala & Godfrey 2013; Godfrey & Shabala 2016; Croston, Ineson & Hardcastle 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.

An additional issue is that the above scaling relation has been determined from characterization 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 preparation).

In the analysis that follows, we use the Heckman & Best (2014) relation (equation 2) with $f_{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. B2), 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 normalization by around 0.3 dex. This normalization uncertainty needs to be borne in mind when comparing the output of the simulations with the observational data.

## 3.2 Cosmic evolution of the kinetic luminosity density

The evolving radio LFs 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_{1.4\text{GHz}}, z) = L_{\text{kin}}(L_{1.4\text{GHz}}) \times \rho(L_{1.4\text{GHz}}, z) \text{ W Mpc}^{-3} \text{ logL}^{-1}. \quad (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
<figure><img src="image_2.png" /><figcaption>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-fitting 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 LF, and therefore most of the heating output comes from relatively high-luminosity sources.</figcaption></figure>

population, the best-fitting LF models were shifted to 1.4 GHz using the standard spectral index $\alpha = -0.7$ to obtain $\rho(L_{1.4\text{GHz}}, z)$ for each redshift bin; this was then convolved with $L_{\text{kin}}(L_{1.4\text{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. A1. 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 normalization, the heating output from Q-LERGs decline sharply, in line with the observed evolution of the LFs.

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

$$ \Omega_{\text{kin}}(z) = \int \Psi(L_{1.4\text{GHz}}, z) d \log L_{1.4\text{GHz}} \text{ W Mpc}^{-3}. \quad (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
<figure><img src="image_2.png" /><figcaption>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-fitting 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 LF, and therefore most of the heating output comes from relatively high-luminosity sources.</figcaption></figure>

population, the best-fitting LF models were shifted to 1.4 GHz using the standard spectral index $\alpha = -0.7$ to obtain $\rho(L_{1.4\text{GHz}}, z)$ for each redshift bin; this was then convolved with $L_{\text{kin}}(L_{1.4\text{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. A1. 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 normalization, the heating output from Q-LERGs decline sharply, in line with the observed evolution of the LFs.

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

$$ \Omega_{\text{kin}}(z) = \int \Psi(L_{1.4\text{GHz}}, z) \, d \log L_{1.4\text{GHz}} \, \text{W Mpc}^{-3}. \quad (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
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The evolution of the kinetic luminosity density, $\Omega_{\text{kin}}(z)$ compared to observations in literature (left-hand panel) and to simulations (right-hand panel). 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_{\text{cav}} = 4$ and the best-fitting 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_{\text{cav}} = 1$ and 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-hand 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, Fabian & Blundell (2013; their LERG-equivalent population) in the right-hand panel.</figcaption></figure>

on the extreme values of $f_{\text{cav}} = 1$ and 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$ dex, so smaller than these extreme values of $f_{\text{cav}}$.

The kinetic luminosity density for both the radio-excess AGN and the LERGs remains roughly constant across $0.5 < z \le 2.5$ with $\Omega_{\text{kin}} \approx (4-5) \times 10^{32} \text{ 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_{\text{kin}} \approx 2.7 \times 10^{32} \text{ W Mpc}^{-3}$ at $z = 0.75$ to nearly an order of magnitude lower by $z \sim 2.5$.

In Fig. 2 (left-hand panel), 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_{\text{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_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
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The evolution of the kinetic luminosity density, $\Omega_{\text{kin}}(z)$ compared to observations in literature (left-hand panel) and to simulations (right-hand panel). 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_{\text{cav}} = 4$ and the best-fitting 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_{\text{cav}} = 1$ and 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-hand 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, Fabian & Blundell (2013; their LERG-equivalent population) in the right-hand panel.</figcaption></figure>

on the extreme values of $f_{\text{cav}} = 1$ and 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$ dex, so smaller than these extreme values of $f_{\text{cav}}$.

The kinetic luminosity density for both the radio-excess AGN and the LERGs remains roughly constant across $0.5 < z \le 2.5$ with $\Omega_{\text{kin}} \approx (4-5) \times 10^{32} \text{ 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_{\text{kin}} \approx 2.7 \times 10^{32} \text{ W Mpc}^{-3}$ at $z = 0.75$ to nearly an order of magnitude lower by $z \sim 2.5$.

In Fig. 2 (left-hand panel), 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_{\text{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_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 LF at $z < 0.7$ and integrated this to find $\Omega_{\text{kin}} \sim 7 \times 10^{31}$ 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), and 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_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-hand panel), 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}_{BH,R}$) over cosmic time for the ‘radio-mode’ population from Croton et al. (2006, see their fig. 3). The $\dot{m}_{BH,R}(z)$ were then translated into a luminosity of the black hole assuming $L_{BH} = \eta \dot{m}_{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_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}_{BH,R}(z)$ and hence the black hole luminosities by $\kappa_R = 0.08$ (which we also apply here).

The resulting prediction from the SAGE model is shown in the right-hand 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 SFGs 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_{\text{Edd}} \equiv (M_{\text{BH}}/\dot{M}_{\text{BH,Edd}}) < \alpha_{\text{crit}}$), resulting in an advection-dominated accretion flow (Narayan & Yi 1994, 1995), and ‘cold-mode’ (or radiatively-efficient) accretion occurs at high rates ($f_{\text{Edd}} > \alpha_{\text{crit}}$) resulting in the formation of an optically-thin accretion disc (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_{\text{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_{\text{kin,hot}} = \eta \dot{M}_{\text{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 normalization that is lower by up to $\sim 0.5$ dex; this suggests that the SFGs 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 SFGs 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 idealized 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$ K, i.e. hot-mode) described by Bondi accretion (Bondi 1952), and the accretion from cold gas ($T < 10^5$ 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 Section 2.1).
inferred jet powers to calculate the jet kinetic LF at $z < 0.7$ and integrated this to find $\Omega_{\text{kin}} \sim 7 \times 10^{31} \text{ 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), and 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_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-hand panel), 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}_{BH,R}$) over cosmic time for the ‘radio-mode’ population from Croton et al. (2006, see their fig. 3). The $\dot{m}_{BH,R}(z)$ were then translated into a luminosity of the black hole assuming $L_{BH} = \eta \dot{m}_{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_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}_{BH,R}(z)$ and hence the black hole luminosities by $\kappa_R = 0.08$ (which we also apply here).

The resulting prediction from the SAGE model is shown in the right-hand 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 SFGs 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_{\text{Edd}} \equiv (M_{\text{BH}}/\dot{M}_{\text{BH,Edd}}) < \alpha_{\text{crit}}$), resulting in an advection-dominated accretion flow (Narayan & Yi 1994, 1995), and ‘cold-mode’ (or radiatively-efficient) accretion occurs at high rates ($f_{\text{Edd}} > \alpha_{\text{crit}}$) resulting in the formation of an optically-thin accretion disc (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_{\text{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_{\text{kin,hot}} = \eta \dot{M}_{\text{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 normalization that is lower by up to $\sim 0.5$ dex; this suggests that the SFGs 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 SFGs 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 idealized 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$ K, i.e. hot-mode) described by Bondi accretion (Bondi 1952), and the accretion from cold gas ($T < 10^5$ 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 Section 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_{\text{kin}} = 0.2 L_{\text{bol}}$, where $L_{\text{bol}} = 0.1 \dot{M}_{\text{BH}}c^2$ is the bolometric luminosity of the AGN, and $\dot{M}_{\text{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_{\text{kin}}$ that increases from $z=0$ and peaks at $\Omega_{\text{kin}} \approx 10^{33} \text{ 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_{\text{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 dashed line), which shows a similar evolution to the observed LERGs but with a systematically higher normalization by up to an order of magnitude. The results of Mocz et al. (2013) 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 normalization than 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.

## 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_{\text{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_{\text{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_{\text{kin}}(z)$ for this population out to $z_{\text{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

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Top panel: the cumulative kinetic energy density radiated from AGN in different forms across cosmic time (in units of J Mpc$^{-3}$). The observed kinetic luminosity densities for Fig. 2 are integrated from $z=5$ to 0, using extrapolation beyond the coverage of the data set (see Section 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 per cent. The green line shows the bolometric energy radiated from quasars calculated using the black hole mass densities (and accretion efficiency of 10 per cent) that was determined by Hopkins, Richards & Hernquist (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 the 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 panel: The fraction of the total (i.e. $z=0$) kinetic energy density emitted across cosmic time.</figcaption></figure>

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 0 in units of 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_{\text{kin}} \approx 10^{50} \text{ J Mpc}^{-3}$, with $\sim 50$ per cent of this emitted within the past $\sim 6$ Gyr; this agrees well with other measurements (e.g. Heckman & Best 2023). For Q-LERGs, the corresponding value is $E_{\text{kin}} \sim 7 \times 10^{49} \text{ J Mpc}^{-3}$, but they deposit most of their total kinetic energy at late times, $z \lesssim 0.5$ (Fig. 3; bottom panel).

For comparison, the green line shows the bolometric radiative energy density from AGN estimated by using an accretion efficiency of 10 per cent and the evolution of the black hole mass density
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_{\text{kin}} = 0.2 L_{\text{bol}}$, where $L_{\text{bol}} = 0.1 \dot{M}_{\text{BH}}c^2$ is the bolometric luminosity of the AGN, and $\dot{M}_{\text{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_{\text{kin}}$ that increases from $z=0$ and peaks at $\Omega_{\text{kin}} \approx 10^{33} \text{ 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_{\text{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 dashed line), which shows a similar evolution to the observed LERGs but with a systematically higher normalization by up to an order of magnitude. The results of Mocz et al. (2013) 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 normalization than 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.

## 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_{\text{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_{\text{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_{\text{kin}}(z)$ for this population out to $z_{\text{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

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Top panel: the cumulative kinetic energy density radiated from AGN in different forms across cosmic time (in units of J Mpc$^{-3}$). The observed kinetic luminosity densities for Fig. 2 are integrated from $z=5$ to 0, using extrapolation beyond the coverage of the data set (see Section 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 per cent. The green line shows the bolometric energy radiated from quasars calculated using the black hole mass densities (and accretion efficiency of 10 per cent) that was determined by Hopkins, Richards & Hernquist (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 the 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 panel: The fraction of the total (i.e. $z=0$) kinetic energy density emitted across cosmic time.</figcaption></figure>

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 0 in units of 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_{\text{kin}} \approx 10^{50} \text{ J Mpc}^{-3}$, with $\sim 50$ per cent of this emitted within the past $\sim 6$ Gyr; this agrees well with other measurements (e.g. Heckman & Best 2023). For Q-LERGs, the corresponding value is $E_{\text{kin}} \sim 7 \times 10^{49} \text{ J Mpc}^{-3}$, but they deposit most of their total kinetic energy at late times, $z \lesssim 0.5$ (Fig. 3; bottom panel).

For comparison, the green line shows the bolometric radiative energy density from AGN estimated by using an accretion efficiency of 10 per cent and the evolution of the black hole mass density
from Hopkins et al. (2007), who determined this by integrating their quasar LF 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 multiphase gas medium in samples of galaxies, with different studies finding values in the range 0.0001–1 per cent $L_{bol}$ (e.g. Fiore et al. 2017; Dall’Agnol de Oliveira et al. 2021; Kakkad et al. 2022, see Heckman & Best 2023 for a more extensive discussion). Given the large range of values observed from different studies for the different multiphase mediums, here, for illustrative purposes, we take values in the range 0.1–0.5 per cent $L_{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 per cent $L_{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 disc or torus structure (based on SED fitting of photometry). This data set has been used previously to study the evolution of the radio LFs of radio-excess AGN and LERGs by Kondapally et al. (2022).

We studied the implications of the observed LFs 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 LFs 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 \le 2.5$, with both populations at $\Omega_{\text{kin}} \approx (4-5) \times 10^{32} \text{ W Mpc}^{-3}$. At $z \le 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, SFGs 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 SFGs, 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_{\text{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 SFGs. SIMBA predicts $\Omega_{\text{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 overpredict the heating at $z \le 1.5$, and underpredict it at higher redshifts compared to the observations. The Q-LERGs in SIMBA also overpredict the heating output at $z \le 1.5$ but match well with observations at higher redshifts. Although these differences are at the 0.3–0.4 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_{\text{kin}}$. Radio-excess AGN (and LERGs) output a total of $E_{\text{kin}} \approx 10^{50} \text{ 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 \le 0.5$. We compared this kinetic energy density from AGN jets to that output by AGN winds, estimated using measurements of the bolometric quasar LF, to find that AGN jets dominate the total energy budget over AGN winds at $z \le 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 LFs 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 benefited from the following recent major funding sources: CNRS-INSU, Observatoire de Paris and Université d'Orléans, France; BMBF, MIWF-NRW, MPG, Germany; Science
from Hopkins et al. (2007), who determined this by integrating their quasar LF 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 multiphase gas medium in samples of galaxies, with different studies finding values in the range 0.0001–1 per cent $L_{bol}$ (e.g. Fiore et al. 2017; Dall’Agnol de Oliveira et al. 2021; Kakkad et al. 2022, see Heckman & Best 2023 for a more extensive discussion). Given the large range of values observed from different studies for the different multiphase mediums, here, for illustrative purposes, we take values in the range 0.1–0.5 per cent $L_{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 per cent $L_{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 disc or torus structure (based on SED fitting of photometry). This data set has been used previously to study the evolution of the radio LFs of radio-excess AGN and LERGs by Kondapally et al. (2022).

We studied the implications of the observed LFs 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 LFs 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 \le 2.5$, with both populations at $\Omega_{kin} \approx (4-5) \times 10^{32} \text{ W Mpc}^{-3}$. At $z \le 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, SFGs 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 SFGs, 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_{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 SFGs. SIMBA predicts $\Omega_{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 overpredict the heating at $z \le 1.5$, and underpredict it at higher redshifts compared to the observations. The Q-LERGs in SIMBA also overpredict the heating output at $z \le 1.5$ but match well with observations at higher redshifts. Although these differences are at the 0.3–0.4 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_{kin}$. Radio-excess AGN (and LERGs) output a total of $E_{kin} \approx 10^{50} \text{ 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 \le 0.5$. We compared this kinetic energy density from AGN jets to that output by AGN winds, estimated using measurements of the bolometric quasar LF, to find that AGN jets dominate the total energy budget over AGN winds at $z \le 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 LFs 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 benefited 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 that have helped improve and clarify this 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'). KM 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 LF 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/roh-itk-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 H. et al., 2019, PASJ, 71, 114

Anglés-Alcázar D., Davé R., Faucher-Giguère C.-A., Özel F., Hopkins P. F., 2017, MNRAS, 464, 2840

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

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

Astropy Collaboration, 2018, AJ, 156, 123

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., 2005, MNRAS, 362, 25

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

Best P. N. et al., 2023, MNRAS, 523, 1729

Bírzan L., Rafferty D. A., McNamara B. R., Wise M. W., Nulsen P. E. J., 2004, ApJ, 607, 800

Bírzan L., McNamara B. R., Nulsen P. E. J., Carilli C. L., Wise M. W., 2008, ApJ, 686, 859

Boehringer H., Voges W., Fabian A. C., Edge A. C., Neumann D. M., 1993, MNRAS, 264, L25

Bondi H., 1952, MNRAS, 112, 195

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

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

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

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

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

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

Cavagnolo K. W., McNamara B. R., Nulsen P. E. J., Carilli C. L., Jones C., Bírzan L., 2010, ApJ, 720, 1066

Chambers K. C. et al., 2016, preprint (arXiv:1612.05560)

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

Croston J. H., Ineson J., Hardcastle M. J., 2018, MNRAS, 476, 1614

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

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

Croton D. J. et al., 2016, ApJS, 222, 22

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

Dall'Agnol de Oliveira B. et al., 2021, MNRAS, 504, 3890

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

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. J. et al., 2021, A&A, 648, A4

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

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

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 F. et al., 2017, A&A, 601, A143

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

Gebhardt K. et al., 2000, ApJ, 539, L13

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

Godfrey L. E. H., Shabala S. S., 2016, MNRAS, 456, 1172

Gonzalez A. H. et al., 2010, American Astronomical Society Meeting Abstracts #216, p. 415.13

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

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

Hardcastle M. J., 2018b, MNRAS, 475, 2768

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

Hardcastle M. J., Krause M. G. H., 2013, MNRAS, 430, 174

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

Heckman T. M., Best P. N., 2023, Galaxies, 11, 21

Hopkins P. F., Quataert E., 2011, MNRAS, 415, 1027

Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ, 654, 731

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

Ineson J., Croston J. H., Hardcastle M. J., Mingo B., 2017, MNRAS, 467, 1586

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

Kakkad D. et al., 2022, MNRAS, 511, 2105

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

Kenter A. et al., 2005, ApJS, 161, 9

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

Kondapally R. et al., 2022, MNRAS, 513, 3742

Körding E. G., Fender R. P., Migliari S., 2006, MNRAS, 369, 1451

Körding E. G., Jester S., Fender R., 2008, MNRAS, 383, 277

Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
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 that have helped improve and clarify this 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'). KM 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 LF 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/roh-itk-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 H. et al., 2019, PASJ, 71, 114

* Anglés-Alcázar D., Davé R., Faucher-Giguère C.-A., Özel F., Hopkins P. F., 2017, MNRAS, 464, 2840

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

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

* Astropy Collaboration, 2018, AJ, 156, 123

* 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., 2005, MNRAS, 362, 25

* 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

* Best P. N. et al., 2023, MNRAS, 523, 1729

* Bîrzan L., Rafferty D. A., McNamara B. R., Wise M. W., Nulsen P. E. J., 2004, ApJ, 607, 800

* Bîrzan L., McNamara B. R., Nulsen P. E. J., Carilli C. L., Wise M. W., 2008, ApJ, 686, 859

* Boehringer H., Voges W., Fabian A. C., Edge A. C., Neumann D. M., 1993, MNRAS, 264, L25

* Bondi H., 1952, MNRAS, 112, 195

* 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

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

* 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

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

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

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

* Cavagnolo K. W., McNamara B. R., Nulsen P. E. J., Carilli C. L., Jones C., Bîrzan L., 2010, ApJ, 720, 1066

* Chambers K. C. et al., 2016, preprint (arXiv:1612.05560)

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

* Croston J. H., Ineson J., Hardcastle M. J., 2018, MNRAS, 476, 1614

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

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

* Croton D. J. et al., 2016, ApJS, 222, 22

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

* Dall'Agnol de Oliveira B. et al., 2021, MNRAS, 504, 3890

* 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

* 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. J. et al., 2021, A&A, 648, A4

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

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

* 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 F. et al., 2017, A&A, 601, A143

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

* Gebhardt K. et al., 2000, ApJ, 539, L13

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

* Godfrey L. E. H., Shabala S. S., 2016, MNRAS, 456, 1172

* Gonzalez A. H. et al., 2010, American Astronomical Society Meeting Abstracts #216, p. 415.13

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

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

* Hardcastle M. J., 2018b, MNRAS, 475, 2768

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

* Hardcastle M. J., Krause M. G. H., 2013, MNRAS, 430, 174

* 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

* Heckman T. M., Best P. N., 2023, Galaxies, 11, 21

* Hopkins P. F., Quataert E., 2011, MNRAS, 415, 1027

* Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ, 654, 731

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

* Ineson J., Croston J. H., Hardcastle M. J., Mingo B., 2017, MNRAS, 467, 1586

* 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

* Kakkad D. et al., 2022, MNRAS, 511, 2105

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

* Kenter A. et al., 2005, ApJS, 161, 9

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

* Kondapally R. et al., 2022, MNRAS, 513, 3742

* Körding E. G., Fender R. P., Migliari S., 2006, MNRAS, 369, 1451

* Körding E. G., Jester S., Fender R., 2008, MNRAS, 383, 277

* 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

McNamara B. R., Nulsen P. E. J., 2007, ARA&A, 45, 117

Magliocchetti M., 2022, A&AR, 30, 6

Marconi A., Risaliti G., Gilli R., Hunt L. K., Maiolino R., Salvati M., 2004, MNRAS, 351, 169

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

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

Merloni A., Heinz S., 2008, MNRAS, 388, 1011

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

Mocz P., Fabian A. C., Blundell K. M., 2013, MNRAS, 432, 3381

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

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

O'Sullivan E., Giacintucci S., David L. P., Gitti M., Vrtilek J. M., Raychaudhury S., Ponman T. J., 2011, ApJ, 735, 11

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

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

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

Raouf M., Shabala S. S., Croton D. J., Khosroshahi H. G., Bernyk M., 2017, MNRAS, 471, 658

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

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

Schaye J. et al., 2015, MNRAS, 446, 521

Shabala S., Alexander P., 2009, ApJ, 699, 525

Shabala S. S., Godfrey L. E. H., 2013, ApJ, 769, 129

Shabala S. S., Santoso J. S., Godfrey L. E. H., 2012, ApJ, 756, 161

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., 2019, A&A, 622, A1

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

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

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

Soltan A., 1982, MNRAS, 200, 115

Springel V. et al., 2018, MNRAS, 475, 676

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

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

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

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

Thomas N., Davé R., Jarvis M. J., Anglés-Alcázar D., 2021, MNRAS, 503, 3492

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 R. J., Shabala S. S., 2015, ApJ, 806, 59

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

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

Whittam I. H. et al., 2022, MNRAS, 516, 245

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

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

Willott C. J., Rawlings S., Blundell K. M., Lacy M., 1999, MNRAS, 309, 1017

## APPENDIX A: RADIO-AGN LUMINOSITY FUNCTIONS USED FOR ESTIMATING KINETIC HEATING RATES

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_{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 \le 1$, $1 < z \le 1.5$, $1.5 < z \le 2$, and $2 < z \le 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. A1. 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. A1, 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 \le 2.5$; however, when split by the 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^*(z)}{(L^*(z)/L)^\beta + (L^*(z)/L)^\gamma}, \quad (A1) $$

where $L^*(z)$ is the characteristic luminosity at redshift $z$, $\rho^*$ 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^*(z)$ and $L^*(z)$ for LERGs (taken from Kondapally et al. 2022) and radio-excess AGN are shown in Table A1. For the Q-LERGs, we directly fitted a broken power law to the $0.5 < z \le 1.0$ LF to derive empirical power-law slopes. Using these slopes, we fitted the evolution at higher redshifts by fixing $\rho^*(z)$ at each redshift bin in accordance with the evolution of the available quiescent host-galaxies and only fitting for $L^*(z)$ (see Kondapally et al. 2022 for details); the parameters for the best-fitting Q-LERG LFs are also shown in Table A1 for completeness.
Lawrence A. et al., 2007, MNRAS, 379, 1599

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

McNamara B. R., Nulsen P. E. J., 2007, ARA&A, 45, 117

Magliocchetti M., 2022, A&AR, 30, 6

Marconi A., Risaliti G., Gilli R., Hunt L. K., Maiolino R., Salvati M., 2004, MNRAS, 351, 169

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

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

Merloni A., Heinz S., 2008, MNRAS, 388, 1011

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

Mocz P., Fabian A. C., Blundell K. M., 2013, MNRAS, 432, 3381

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

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

O'Sullivan E., Giacintucci S., David L. P., Gitti M., Vrtilek J. M., Raychaudhury S., Ponman T. J., 2011, ApJ, 735, 11

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

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

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

Raouf M., Shabala S. S., Croton D. J., Khosroshahi H. G., Bernyk M., 2017, MNRAS, 471, 658

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

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

Schaye J. et al., 2015, MNRAS, 446, 521

Shabala S., Alexander P., 2009, ApJ, 699, 525

Shabala S. S., Godfrey L. E. H., 2013, ApJ, 769, 129

Shabala S. S., Santoso J. S., Godfrey L. E. H., 2012, ApJ, 756, 161

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., 2019, A&A, 622, A1

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

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

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

Soltan A., 1982, MNRAS, 200, 115

Springel V. et al., 2018, MNRAS, 475, 676

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

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

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

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

Thomas N., Davé R., Jarvis M. J., Anglés-Alcázar D., 2021, MNRAS, 503, 3492

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 R. J., Shabala S. S., 2015, ApJ, 806, 59

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

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

Whittam I. H. et al., 2022, MNRAS, 516, 245

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

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

Willott C. J., Rawlings S., Blundell K. M., Lacy M., 1999, MNRAS, 309, 1017

## APPENDIX A: RADIO-AGN LUMINOSITY FUNCTIONS USED FOR ESTIMATING KINETIC HEATING RATES

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_{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 \le 1$, $1 < z \le 1.5$, $1.5 < z \le 2$, and $2 < z \le 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. A1. 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. A1, 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 \le 2.5$; however, when split by the 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^*(z)}{(L^*(z)/L)^\beta + (L^*(z)/L)^\gamma} \quad (A1) $$

where $L^*(z)$ is the characteristic luminosity at redshift $z$, $\rho^*$ 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^*(z)$ and $L^*(z)$ for LERGs (taken from Kondapally et al. 2022) and radio-excess AGN are shown in Table A1. For the Q-LERGs, we directly fitted a broken power law to the $0.5 < z \le 1.0$ LF to derive empirical power-law slopes. Using these slopes, we fitted the evolution at higher redshifts by fixing $\rho^*(z)$ at each redshift bin in accordance with the evolution of the available quiescent host-galaxies and only fitting for $L^*(z)$ (see Kondapally et al. 2022 for details); the parameters for the best-fitting Q-LERG LFs are also shown in Table A1 for completeness.
<figure><img src="image_5.png" /><figcaption>**Figure A1.** Cosmic evolution of AGN LFs from LoTSS-Deep across $0.5 < z \le 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 the text) is evolved in both density and luminosity (the shaded regions represent the uncertainty on this fit; see the text).</figcaption></figure>

<table><caption>**Table A1.** Results from modelling the evolution of three populations of AGN LFs, plotted in Fig. A1.</caption><thead><tr><th rowspan="2">z</th><th colspan="3">Radio-excess AGN</th><th colspan="3">LERGs</th><th colspan="3">Quiescent LERGs</th></tr><tr><th>log<sub>10</sub>ρ<sup>*</sup>(z)</th><th>log<sub>10</sub>L<sub>*</sub>(z)</th><th>χ<sup>2</sup><sub>v</sub></th><th>log<sub>10</sub>ρ<sup>*</sup>(z)</th><th>log<sub>10</sub>L<sub>*</sub>(z)</th><th>χ<sup>2</sup><sub>v</sub></th><th>log<sub>10</sub>ρ<sup>*</sup>(z)</th><th>log<sub>10</sub>L<sub>*</sub>(z)</th><th>χ<sup>2</sup><sub>v</sub></th></tr></thead><tbody><tr><td>0.5 &lt; z &le; 1.0</td><td>-5.53<sup>+0.07</sup><sub>-0.08</sub></td><td>26.31<sup>+0.14</sup><sub>-0.12</sub></td><td>7.22</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>5.59</td><td>-6.37<sup>+0.17</sup><sub>-0.24</sub></td><td>27.03<sup>+0.40</sup><sub>-0.22</sub></td><td>3.57</td></tr><tr><td>1.0 &lt; z &le; 1.5</td><td>-5.89<sup>+0.08</sup><sub>-0.09</sub></td><td>26.77<sup>+0.16</sup><sub>-0.14</sub></td><td>5.83</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>5.56</td><td>-6.92</td><td>27.33<sup>+0.04</sup><sub>-0.04</sub></td><td>3.39</td></tr><tr><td>1.5 &lt; z &le; 2.0</td><td>-5.71<sup>+0.07</sup><sub>-0.08</sub></td><td>26.76<sup>+0.13</sup><sub>-0.12</sub></td><td>7.98</td><td>-5.52<sup>+0.08</sup><sub>-0.08</sub></td><td>26.35<sup>+0.12</sup><sub>-0.12</sub></td><td>6.28</td><td>-7.28</td><td>27.91<sup>+0.06</sup><sub>-0.06</sub></td><td>2.28</td></tr><tr><td>2.0 &lt; z &le; 2.5</td><td>-5.41<sup>+0.09</sup><sub>-0.09</sub></td><td>26.09<sup>+0.13</sup><sub>-0.13</sub></td><td>4.98</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>3.45</td><td>-7.61</td><td>27.93<sup>+0.12</sup><sub>-0.12</sub></td><td>5.89</td></tr></tbody></table>

*Notes.* 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^*(z)$ at each redshift is fixed based on the available quiescent host galaxies at that redshift, and the $L_*(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 \le 1$ (with slopes of $\beta = -2.88$ and $\gamma = -0.55$; see Kondapally et al. 2022).

## APPENDIX B: COMPARISON OF JET POWER INFERENCE METHODS

There exist many scaling relations for translating monochromatic radio luminosities ($L_{1.4\,GHz}$) into jet kinetic powers ($L_{\text{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
<figure><img src="image_5.png" /><figcaption>**Figure A1.** Cosmic evolution of AGN LFs from LoTSS-Deep across $0.5 < z \le 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 the text) is evolved in both density and luminosity (the shaded regions represent the uncertainty on this fit; see the text).</figcaption></figure>

<table><caption>**Table A1.** Results from modelling the evolution of three populations of AGN LFs, plotted in Fig. A1.</caption><thead><tr><th rowspan="2">z</th><th colspan="3">Radio-excess AGN</th><th colspan="3">LERGs</th><th colspan="3">Quiescent LERGs</th></tr><tr><th>log<sub>10</sub>ρ<sup>*</sup>(z)</th><th>log<sub>10</sub>L<sub>*</sub>(z)</th><th>χ<sup>2</sup><sub>v</sub></th><th>log<sub>10</sub>ρ<sup>*</sup>(z)</th><th>log<sub>10</sub>L<sub>*</sub>(z)</th><th>χ<sup>2</sup><sub>v</sub></th><th>log<sub>10</sub>ρ<sup>*</sup>(z)</th><th>log<sub>10</sub>L<sub>*</sub>(z)</th><th>χ<sup>2</sup><sub>v</sub></th></tr></thead><tbody><tr><td>0.5 &lt; z &le; 1.0</td><td>-5.53<sup>+0.07</sup><sub>-0.08</sub></td><td>26.31<sup>+0.14</sup><sub>-0.12</sub></td><td>7.22</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>5.59</td><td>-6.37<sup>+0.17</sup><sub>-0.24</sub></td><td>27.03<sup>+0.40</sup><sub>-0.22</sub></td><td>3.57</td></tr><tr><td>1.0 &lt; z &le; 1.5</td><td>-5.89<sup>+0.08</sup><sub>-0.09</sub></td><td>26.77<sup>+0.16</sup><sub>-0.14</sub></td><td>5.83</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>5.56</td><td>-6.92</td><td>27.33<sup>+0.04</sup><sub>-0.04</sub></td><td>3.39</td></tr><tr><td>1.5 &lt; z &le; 2.0</td><td>-5.71<sup>+0.07</sup><sub>-0.08</sub></td><td>26.76<sup>+0.13</sup><sub>-0.12</sub></td><td>7.98</td><td>-5.52<sup>+0.08</sup><sub>-0.08</sub></td><td>26.35<sup>+0.12</sup><sub>-0.12</sub></td><td>6.28</td><td>-7.28</td><td>27.91<sup>+0.06</sup><sub>-0.06</sub></td><td>2.28</td></tr><tr><td>2.0 &lt; z &le; 2.5</td><td>-5.41<sup>+0.09</sup><sub>-0.09</sub></td><td>26.09<sup>+0.13</sup><sub>-0.13</sub></td><td>4.98</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>3.45</td><td>-7.61</td><td>27.93<sup>+0.12</sup><sub>-0.12</sub></td><td>5.89</td></tr></tbody></table>

*Notes.* 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^*(z)$ at each redshift is fixed based on the available quiescent host galaxies at that redshift, and the $L_*(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 \le 1$ (with slopes of $\beta = -2.88$ and $\gamma = -0.55$; see Kondapally et al. 2022).

## APPENDIX B: COMPARISON OF JET POWER INFERENCE METHODS

There exist many scaling relations for translating monochromatic radio luminosities ($L_{1.4\,GHz}$) into jet kinetic powers ($L_{\text{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
<figure><img src="image_6.png" /><figcaption>**Figure B1.** Comparison of different scaling relations between monochromatic radio luminosity (1.4 GHz) and the jet powers in the literature, as highlighted in Table B1 (Willott et al. 1999; Körding, Fender & Migliari 2006; Körding, Jester & Fender 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_{\text{cav}} = 4$, which shows good agreement in its normalization with that of Willott et al. (1999) when using $f_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$ 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}$ W Hz$^{-1}$, which shows a higher normalization than the other relations considered.</figcaption></figure>

impacts on the cosmic kinetic luminosity density results presented in Section 3.

Fig. B1 shows the different scaling relations between 1.4-GHz radio luminosity and jet kinetic powers. These scaling relations are also listed in Table B1. 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$ 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_W$ which is expected to lie in the range $f_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_W =$

15 (e.g. Heckman & Best 2014; Smolčić et al. 2017), this relation produces a similar normalization to that of the cavity based estimates at typical luminosities of $L_{1.4\,\text{GHz}} \sim 10^{25}$ 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_{1.4\,\text{GHz}} \gtrsim 10^{24}$ W Hz$^{-1}$), finding a similar slope to the Willott et al. (1999) relation with a normalization that is consistent with a lower $f_W$ value of $f_W \sim 4$ (see Table B1).

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. B1, 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_W \sim 4$.

Also shown in Fig. B1 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 normalization. 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. B2). 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 normalization (by $\sim 0.3$ dex) such that it predicts higher kinetic powers than the Heckman & Best (2014) scaling relation used in this study over all radio luminosities.

In Fig. B2, we show the cosmic evolution of the kinetic luminosity density of AGN within LoTSS-Deep (see Section 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_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_{1.4\,\text{GHz}} \lesssim 10^{24}$ 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_{150\,\text{MHz}} \sim 10^{26}-10^{27}$ 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_{\text{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 normalization of the kinetic luminosity densities.
<figure><img src="image_6.png" /><figcaption>**Figure B1.** Comparison of different scaling relations between monochromatic radio luminosity (1.4 GHz) and the jet powers in the literature, as highlighted in Table B1 (Willott et al. 1999; Körding, Fender & Migliari 2006; Körding, Jester & Fender 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_{\text{cav}} = 4$, which shows good agreement in its normalization with that of Willott et al. (1999) when using $f_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$ 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}$ W Hz$^{-1}$, which shows a higher normalization than the other relations considered.</figcaption></figure>

impacts on the cosmic kinetic luminosity density results presented in Section 3.

Fig. B1 shows the different scaling relations between 1.4-GHz radio luminosity and jet kinetic powers. These scaling relations are also listed in Table B1. 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$ 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_W$ which is expected to lie in the range $f_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_W =$

15 (e.g. Heckman & Best 2014; Smolčić et al. 2017), this relation produces a similar normalization to that of the cavity based estimates at typical luminosities of $L_{1.4\,\text{GHz}} \sim 10^{25}$ 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_{1.4\,\text{GHz}} \gtrsim 10^{24}$ W Hz$^{-1}$), finding a similar slope to the Willott et al. (1999) relation with a normalization that is consistent with a lower $f_W$ value of $f_W \sim 4$ (see Table B1).

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. B1, 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_W \sim 4$.

Also shown in Fig. B1 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 normalization. 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. B2). 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 normalization (by $\sim 0.3$ dex) such that it predicts higher kinetic powers than the Heckman & Best (2014) scaling relation used in this study over all radio luminosities.

In Fig. B2, we show the cosmic evolution of the kinetic luminosity density of AGN within LoTSS-Deep (see Section 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_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_{1.4\,\text{GHz}} \lesssim 10^{24}$ 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_{150\,\text{MHz}} \sim 10^{26}-10^{27}$ 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_{\text{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 normalization of the kinetic luminosity densities.
<table><caption>**Table B1.** Different literature scaling relations for converting monochromatic radio luminosities ($L_{1.4\,GHz}$) to jet kinetic powers ($L_{\text{kin}}$) that are compared.</caption><thead><tr><th>Reference</th><th>Scaling relation</th></tr></thead><tbody><tr><td>Willott et al. (1999)</td><td>log<sub>10</sub>(L<sub>kin</sub>/W) = 0.86 log<sub>10</sub>(L<sub>1.4GHz</sub>/10<sup>25</sup> W Hz<sup>-1</sup>) + 37.37 + log<sub>10</sub>(f<sub>W</sub>/15)</td></tr><tr><td>Körding et al. (2008)</td><td>log<sub>10</sub>(L<sub>kin</sub>/W) = 0.71 log<sub>10</sub>(L<sub>1.4GHz</sub>/10<sup>25</sup> W Hz<sup>-1</sup>) + 37.72</td></tr><tr><td>Cavagnolo et al. (2010)</td><td>log<sub>10</sub>(L<sub>kin</sub>/W) = 0.75 log<sub>10</sub>(L<sub>1.4GHz</sub>/10<sup>25</sup> W Hz<sup>-1</sup>) + 38.06</td></tr><tr><td>O'Sullivan et al. (2011)</td><td>log<sub>10</sub>(L<sub>kin</sub>/W) = 0.63 log<sub>10</sub>(L<sub>1.4GHz</sub>/10<sup>25</sup> W Hz<sup>-1</sup>) + 37.76</td></tr><tr><td>Shabala & Godfrey (2013)</td><td>log<sub>10</sub>(L<sub>kin</sub>/W) = 0.8 log<sub>10</sub>(L<sub>1.4GHz</sub>/10<sup>25</sup> W Hz<sup>-1</sup>) + 35.12 + log<sub>10</sub>(1+z) + 0.58 log<sub>10</sub>(D/kpc)</td></tr><tr><td>Heckman & Best (2014)</td><td>log<sub>10</sub>(L<sub>kin</sub>/W) = 0.68 log<sub>10</sub>(L<sub>1.4GHz</sub>/10<sup>25</sup> W Hz<sup>-1</sup>) + 37.45 + log<sub>10</sub>(f<sub>cav</sub>/4)</td></tr><tr><td>Ineson et al. (2017)</td><td>log<sub>10</sub>(L<sub>kin</sub>/W) = 0.89 log<sub>10</sub>(L<sub>1.4GHz</sub>/10<sup>25</sup> W Hz<sup>-1</sup>) + 36.05</td></tr></tbody></table>

*Notes.* For Willott et al. (1999) and Heckman & Best (2014), the relations listed have explicitly assumed $f_W = 15$ and $f_{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.

<figure><img src="image_7.png" /><figcaption>**Figure B2.** 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. B1. 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.</figcaption></figure>

This paper has been typeset from a *TEX/LATEX* file prepared by the author.
<table><caption>**Table B1.** Different literature scaling relations for converting monochromatic radio luminosities ($L_{1.4GHz}$) to jet kinetic powers ($L_{kin}$) that are compared.</caption><thead><tr><th>Reference</th><th>Scaling relation</th></tr></thead><tbody><tr><td>Willott et al. (1999)</td><td>log<sub>10</sub>(L<sub>kin</sub>/W) = 0.86 log<sub>10</sub>(L<sub>1.4GHz</sub>/10<sup>25</sup> W Hz<sup>-1</sup>) + 37.37 + log<sub>10</sub>(f<sub>W</sub>/15)</td></tr><tr><td>Körding et al. (2008)</td><td>log<sub>10</sub>(L<sub>kin</sub>/W) = 0.71 log<sub>10</sub>(L<sub>1.4GHz</sub>/10<sup>25</sup> W Hz<sup>-1</sup>) + 37.72</td></tr><tr><td>Cavagnolo et al. (2010)</td><td>log<sub>10</sub>(L<sub>kin</sub>/W) = 0.75 log<sub>10</sub>(L<sub>1.4GHz</sub>/10<sup>25</sup> W Hz<sup>-1</sup>) + 38.06</td></tr><tr><td>O'Sullivan et al. (2011)</td><td>log<sub>10</sub>(L<sub>kin</sub>/W) = 0.63 log<sub>10</sub>(L<sub>1.4GHz</sub>/10<sup>25</sup> W Hz<sup>-1</sup>) + 37.76</td></tr><tr><td>Shabala & Godfrey (2013)</td><td>log<sub>10</sub>(L<sub>kin</sub>/W) = 0.8 log<sub>10</sub>(L<sub>1.4GHz</sub>/10<sup>25</sup> W Hz<sup>-1</sup>) + 35.12 + log<sub>10</sub>(1+z) + 0.58 log<sub>10</sub>(D/kpc)</td></tr><tr><td>Heckman & Best (2014)</td><td>log<sub>10</sub>(L<sub>kin</sub>/W) = 0.68 log<sub>10</sub>(L<sub>1.4GHz</sub>/10<sup>25</sup> W Hz<sup>-1</sup>) + 37.45 + log<sub>10</sub>(f<sub>cav</sub>/4)</td></tr><tr><td>Ineson et al. (2017)</td><td>log<sub>10</sub>(L<sub>kin</sub>/W) = 0.89 log<sub>10</sub>(L<sub>1.4GHz</sub>/10<sup>25</sup> W Hz<sup>-1</sup>) + 36.05</td></tr></tbody></table>

*Notes.* For Willott et al. (1999) and Heckman & Best (2014), the relations listed have explicitly assumed $f_W = 15$ and $f_{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.

<figure><img src="image_7.png" /><figcaption>**Figure B2.** 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. B1. 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.</figcaption></figure>

This paper has been typeset from a TEX/LATEX file prepared by the author.